15MB - City Research Online - City, University of London [PDF]

results), it has the disadvantage that solution for the friction factor requires iteration. This makes it less attractiv

0 downloads 7 Views 15MB Size

Recommend Stories


City Research Online - City, University of London
You have to expect things of yourself before you can do them. Michael Jordan

City Research Online - City, University of London
Make yourself a priority once in a while. It's not selfish. It's necessary. Anonymous

City Research Online - City, University of London
Almost everything will work again if you unplug it for a few minutes, including you. Anne Lamott

City Research Online - City, University of London
Forget safety. Live where you fear to live. Destroy your reputation. Be notorious. Rumi

City Research Online - City, University of London
If you want to become full, let yourself be empty. Lao Tzu

City Research Online - City, University of London
You can never cross the ocean unless you have the courage to lose sight of the shore. Andrè Gide

City Research Online - City, University of London
At the end of your life, you will never regret not having passed one more test, not winning one more

City Research Online - City, University of London
Be grateful for whoever comes, because each has been sent as a guide from beyond. Rumi

City Research Online - City, University of London
You often feel tired, not because you've done too much, but because you've done too little of what sparks

City Research Online - City, University of London
The best time to plant a tree was 20 years ago. The second best time is now. Chinese Proverb

Idea Transcript


              

City Research Online

City, University of London Institutional Repository Citation: Kimambo, C. (1996). Modelling of linebreak in high-pressure gas pipes. (Unpublished Doctoral thesis, City University London)

This is the accepted version of the paper. This version of the publication may differ from the final published version. Permanent repository link:

http://openaccess.city.ac.uk/7934/

Link to published version: Copyright and reuse: City Research Online aims to make research outputs of City, University of London available to a wider audience. Copyright and Moral Rights remain with the author(s) and/or copyright holders. URLs from City Research Online may be freely distributed and linked to. City Research Online:

http://openaccess.city.ac.uk/

[email protected]

MODELLING OF LINEBREAK IN HIGH-PRESSURE GAS PIPES

by Cuthbert ZM Kimambo

A thesis submitted to City University in partial fulfilment of the requirement for the degree of Doctor of Philosophy

July, 1996

Department of Mechanical Engineering. & Aeronautics

City University London

TABLE OF CONTENTS Page 2

TABLE OF CONTENTS LIST OF TABLES

6

LIST OF FIGURES

7

ACKNOWLEDGMENTS

10

DECLARATION

11

ABSTRACT

12

NOMENCLATURE

13

CHAPTER

16

1: INTRODUCTION

1.1

TRANSPORTATION OF NATURAL GAS IN PIPELINES

16

1.2

RUPTURE AND BLOWDOWN OF GAS PIPELINES

17

1.3

COMPUTER CODES FOR ANALYSIS OF FLUID FLOW IN PIPES

19

1.4

PREVIOUS WORK AT CITY UNIVERSITY

22

CHAPTER 2: THEORETICAL DEVELOPMENT OF THE BASIC EQUATIONS FOR UNSTEADY FLOW

27

2.1

INTRODUCTION

27

2.2

MAJOR ASSUMPTIONS AND SIMPLIFICATIONS

32

2.2.1 Dimensionof Flow

32

2.2.2 Flow Phase

34

2.2.3 Fluid Structure Interaction

36

2.2.4 Minor Lossesand Changesin Cross-section

38

DERIVATION OF THE BASIC EQUATIONS FOR UNSTEADY FLOW OF A COMPRESSIBLE FLUID IN A PIPE

39

THERMODYNAMIC

48

2.3 2.4

AND TRANSPORT PROPERTIES

2.4.1 Introduction 2.4.2 Thermal Equationsof Statefor Real Gases 2.4.3 TheoreticalApproach to the Equation of State 2.4.3.1 Introduction to the Basic Principlesand Equations

48 50 53 53

2.4.3.2 The Principle of Corresponding States 2.4.3.3 Extension of the Principle of Corresponding States

56

2.4.3.4 CommercialComputer Software 2.4.3.4.1 IUPAC 2.4.3.4.2 QUANT

58 58 58

2.4.3.4.3 PREPROP

59

2.4.4 Other Approachesfor Real-GasMixtures 2.4.5 Caloric Equationsof State 2

57

60 60

2.5

2.6

2.7

FRICTIONAL FORCE

62

2.5.1 Introduction 2.5.2 Steady Flow Friction Factor 2.5.3 Flow-Dependent Friction Factor 2.5.4 Frequency-Dependent Friction Factor 2.5.5 Friction Factor for Two-Phase Homogeneous Flow 2.5.6 Friction Factor With Respect to Fluid Structure Interaction 2.5.7 Approximation of Friction Factor When Solving the Basic Equations

62 63 64 68 69 70 70

HEAT TRANSFER

71

2.6.1 Introduction 2.6.2 Heat Transfer Process 2.6.3 Calculation of Heat Transfer

71 71 74

OTHER APPROACHESTO THE BASIC EQUATIONS

79

CHAPTER 3: METHODS FOR SOLUTION OF THE BASIC EQUATIONS 82 3.1

INTRODUCTION

82,

3.2

NUMERICAL METHODS OF SOLUTION 3.2.1 Finite-DifferenceMethods 3.2.1.1 GeneralDescription 3.2.1.2 Explicit Finite-DifferenceMethods 3.2.1.3 Method of Characteristics 3.2.1.4 Lax-Wendroff Second-OrderTwo-Step Method 3.2.1.5 Implicit Finite-DifferenceMethods 3.2.2 Finite-ElementMethods 3.2.3 Flux-DifferenceSplitting Schemes

83

3.2.4 Method of Lines

94

3.2.5 Wave-PlanMethod

96

A REVIEW OF SOME NUMERICAL STUDIES

97

3.3.1 GeneralStudieson TransientFlow

97

3.3

3.4

4.2

83 84 86 88 89 90 92

3.3.2 Studies on Pipeline Rupture and Blowdown

109

DISCUSSION OF NUMERICAL METHODS

118

CHAPTER 4: THE COMPUTER

4.1

83

MODEL

FOR LINEBREAK

ANALYSIS

122

SIMPLIFICATION OF THE BASIC EQUATIONS OF FLOW

122

4.1.1 Shouldthe SmallTermsbe Neglected? 4.1.2 Solution Procedure 4.1.3 Representationof Non-Linear Terms

122 124 126

INITIAL AND BOUNDARY CONDITIONS

126

4.2.1 Introduction 4.2.2 SteadyStateAnalysis

126 128

3

4.2.2.1 Basic Equations for Steady State Flow 4.2.2.2 Incompressible Flow

128 130

4.2.2.3 IsothermalCompressibleFlow 4.2.2.4 Adiabatic CompressibleFlow 4.2.2.5 Non-isothermalNon-adiabaticCompressibleFlow

130 132 139

4.2.3 Break Boundary Conditions

140

NUMERICAL SOLUTION OF THE BASIC EQUATIONS FOR UNSTEADY FLOW

146

4.3.1 Introduction

146

4.3.2 Method of Characteristics

148

4.3.3 MacCormack Second-orderTwo-step Method 4.3.4 Warming-Kutler-Lomax Third-order Method

179 182

4.4

VARIABLE GRID SIZE

184

4.5

APPLICATION OF THE QUANT SOFTWARE FOR

4.3

THERMODYNAMIC

AND TRANSPORT PROPERTIES OF FLUIDS

187

4.6

COMPUTER CODES

192

4.7

PREPARATION OF GENERAL GAS AND SYSTEM DATA

205

CHAPTER 5: A REVIEW OF SOME EXPERIMENTAL AND NUMERICAL DATA

212

5.1

INTRODUCTION

212

5.2

LABORATORY EXPERIMENTS

213

5.2.1 Description of the Shock Tube Test 5.2.2 Review of SomeLaboratory Experiment

213 214

5.3

REVIEW OF SOME FULL-SCALE PIPELINE EXPERIMENTS

217

5.4

REVIEW OF SOME NUMERICAL DATA

220

5.5

SELECTION OF TEST DATA

222

CHAPTER

6: VALIDATION

OF THE COMPUTER

MODEL

6.1

VALIDATION

6.2

COMPARISON OF COMPUTER MODEL PREDICTIONS

6.3

226 226

PROCEDURE

WITH EXPERIMENTAL RESULTS

227

6.2.1 Foothills Test Data

227

6.2.2 British Gas Test Data

231

6.2.3 SNGSO Test Data 6.2.4 API Test Data

237 244

DISCUSSION OF VALIDATION RESULTS

249

4

CHAPTER

7: CASE STUDY: THE SONGO SONGO-DAR NATURAL GAS PIPELINE

ES SALAAM 272

7.1

A BRIEF OVERVIEW OF TANZANIA'S

7.2

DESCRIPTION OF THE SONGO SONGO GAS DEVELOPMENT PROJECT

274

7.2.1 Background Information 7.2.2 Gas Properties

274 275

7.2.3 Preliminary Pipeline Design

277

COMPUTER SIMULATION OF A LINEBREAK IN THE SONGO SONGO-DAR ES SALAAM GAS PIPELINE

281

DISCUSSION OF SIMULATION RESULTS

286.

7.3 7.4

CHAPTER 8: RECOMMENDATIONS CHAPTER 9:

ENERGY SECTOR

FOR FURTHER WORK

CONCLUSIONS

272

288 291

REFERENCES AND BIBLIOGRAPHY

297

APPENDICES

309

A

GeneralEquations

309

B

Expressionsfor Equation of State

310

C

Expressionsfor Friction Factor

314

D

List of Programmesand Sub-routines

318

E

Programme Listings (in a floppy disk)

5

LIST OF TABLES Table 2.1

McAdams Constantsfor Natural Convectionto Air

Table 4.1

PossibleBoundary Conditions

Table 4.2

General Gas and System Data Layout

Table 6.1

Foothills Tests GeneralData

Table 6.2

AverageGasCompositionfor Foothills Tests

Table 7.1

Average(Mol. %) Composition of Songo SongoGas

Table 7.2

Songo SongoGasProperties

Table 7.3

ProjectedGasDemandfor the Three Scenariosfor Years 2011 and 2016

Table 7.4

EstimatedGasVelocities and Pipe Sizes

Table 7.5

Pipe Parameters

Table 7.6

Block Valve Parameters

6

LIST OF FIGURES Figure 2.1

Flow Pattern for Air-water Mixture Flowing in a Pipe

Figure 2.2

SketchShowing the Control Volume at the Beginningof a Time Step

Figure 2.3

SketchIllustrating Equation 2.2

Figure 2.4

Conical Control Volume Illustrating Body Forces

Figure 2.5

Control Volume Illustrating Heat and Work Transfer

Figure 2.6

Distributionof Pressure,TemperatureandVelocity Along the Pipe Diameter

Figure 2.7

RadialTemperatureDistributionAcrossa PipeSectionand its Surroundings

Figure 3.1

Finite-differenceGrid Illustrating the Two-step Lax-Wendroff Method

Figure 4.1

Grid Illustrating SteadyState Analysis

Figure 4.2

Pressureand Velocity Variations Illustrating Flatt's Decompression Model

Figure 4.3

Decompression Curve at the Break Boundary

Figure 4.4

Grid IllustratingHybrid Method of Characteristics Solutionat Interior Points

Figure 4.5

Solution Domain Using Interior Points Only

Figure 4.6

Grid Illustrating Hybrid Method of Characteristics Solution at Boundary Points

Figure 4.7

Finite-differenceGrid Illustrating the MacCormackMethod Solution

Figure 4.8

Finite-differenceGrid Illustrating the Warming-Kutler-Lomax Method Solution

Figure 4.9

Grid Illustrating Hybrid Method of CharacteristicsSolution at Interior Boundary Points Without Interpolation

Figure 4.10

Grid Illustrating Hybrid Method of CharacteristicsSolution at Interior Boundary Points With Interpolation

Figure 4.11

Flow Chart for Tiley's ComputerProgramme

Figure 4.12

Flow Chart for the Main Computer Code for the New Model

Figure 4.13

Flow Chart for Sub-programmeSTEAD

Figure 4.14

Flow Chart for Sub-programmesTRANS and BREAK for Transient Analysis

Figure 4.15

Flow Chart for a Typical Method of CharacteristicsProgramme

Figure 4.16

Flow Chart for Sub-routinesS..MOC.. 1 and S..MOC..2

Figure 4.17

Dimensionsof the Test Sectionof the Pipe

Figure 4.18

Diametersof the Pipe in the Test Section

Figure 4.19

Test Sectionof the Pipe After the Break Illustrating a Variable Grid

Figure 6.1

Schematicof Foothills Test NABTFI 7

Figure 6.2

Schematicof Foothills Test NABTF7

Figure 6.3

p-t Curves for Foothills Test NABTFI East Comparison with First-order Method of Characteristics Predictions

Figure 6.4

p-t Curves for Foothills Test NABTFI East Comparison with MacCormack Method Predictions

Figure 6.5

PressureWave PropagationSpeedfor Foothills Test NABTF 1 East

Figure 6.6

p-t Curvesfor Foothills Test NABTF7 West Comparisonwith First-order Method of CharacteristicsPredictions

Figure 6.7

Pressure Wave Propagation Speed for Foothills Test NABTF7 West

Figure 6.8

Schematic of British Gas Tests BGT1, BGT2 and BGT3

Figure 6.9

p-t Curvesfor British GasTest BGTI Comparisonwith First-order Method of CharacteristicsPredictions

Figure 6.10

p-t Curvesfor British GasTest BGT1 Comparisonwith Second-order Method of CharacteristicsPredictions

Figure 6.11

p-t Curvesfor British GasTest BGT1 Comparisonwith MacCormack Method (Alternative 1) Predictions

Figure 6.12

p-t Curvesfor British GasTest BGT2 Comparisonwith MacCormack Method (Alternative 1) Predictions

Figure 6.13

p-t Curvesfor British GasTest BGT2 Comparisonwith MacCormack Method (Alternative 2) Predictions

Figure 6.14

p-t Curvesfor British GasTest BGT2 Comparisonwith MacCormack Method (Alternative 3) Predictions

Figure 6.15

p-t Curves for British Gas Test BGT3 Comparison with MacCormack Method (Alternative 3) Predictions

Figure 6.16

p-t Curvesfor British GasTest BGT3 Comparisonwith WarmingKutler-Lomax Method Predictions

Figure 6.17

p-t Curvesfor British GasTestBGT3 Comparisonwith First-order Method of CharacteristicsPredictions(Ax=0.01 m)

Figure 6.18

p-t Curves for British Gas Test BGT3 Comparison with Second-order Method of Characteristics Predictions (Ax=O. OIm)

Figure 6.19

p-t Curvesfor British GasTest BGT3 Comparisonwith Secondorder Method of CharacteristicsPredictions(Ax=0.1m)

Figure 6.20

PressureWave PropagationSpeedCurvesfor BMI Tests

Figure 6.21

PressureWave PropagationSpeedCurvesfor University of Calgary Tests

Figure 6.22

Schematicof SNGSOTest

Figure 6.23

u-t Curvesfor SNGSO Test Comparisonwith First-order Method of CharacteristicsPredictions

8

Figure 6.24

Method First-order for SNGSO Comparison Curves Test of with p-t Characteristics Predictions

Figure 6.25 Figure 6.26

Velocity Head Curvesfor SNGSOTest Comparisonwith First-order Method of CharacteristicsPredictions Schematicof Alberta PetroleumIndustry Tests APIT1 and APIT2

Figure 6.27

Schematicof Alberta PetroleumIndustry Tests APIT3

Figure 6.28

Mass Flow Rate for Alberta PetroleumIndustry Test APITI Comparison with Second-orderMethod of CharacteristicsPredictions

Figure 6.29

Pressure at Intact End for Alberta Petroleum Industry Test APITI Comparison with Second-order Method of Characteristics Predictions

Figure 6.30

Mass Flow Rate for Alberta Petroleum Industry Test APIT2 Comparison

with First-order Method of CharacteristicsPredictions Figure 6.31

Pressureat Intact End for Alberta PetroleumIndustry Test APIT2 Comparisonwith First-order Method of CharacteristicsPredictions

Figure 6.32

Mass Flow Rate for Alberta PetroleumIndustry Test APIT3 Comparison with First-order Method of CharacteristicsPredictions

Figure 6.33

Pressureat Intact End for Alberta PetroleumIndustry Test APIT3 Comparisonwith First-order Method of CharacteristicsPredictions

Figure 7.1

Map of Africa Showing the Location of Tanzania

Figure 7.2

Map of TanzaniaIndicating the Location of the SongoSongo Gas DevelopmentProject

Figure 7.3

Kilwa Kivinje-Dar es SalaamGasPipelineLocation Map

Figure 7.4

PipelineRoute from Songo SongoIsland to Dar es Salaam

Figure 7.5

Line PressureProfile for Low Scenario

Figure 7.6

Line PressureProfile for Medium Scenario

Figure 7.7

Block Valves Location

Figure 7.8

Line PressureProfile for High Scenario Comparison with Model Predictions

Figure 7.9

Line Pressure Profile for Simulated Break in the Songo Songo-Dar es Salaam Pipeline

Figure 7.10

Mass Flow Rate Profile for Simulated Break in the Songo Songo-Dar es Salaam Pipeline

9

ACKNOWLEDGMENTS individuals to and I wish to expressmy sinceregratitude and appreciation all organizationswho contributedto the successfulcompletionof this study. I wish to expressspecialthanksand appreciationto my supervisor,Prof. A. R. D. Thorley,first for providingthe researchtopic; for his continuousencouragement;guidance; life. in beyond academic support, even matters and criticisms; I alsowish to extendparticularthanksto the GermanGovernmentwho, through the GermanAcademicExchangeService(DAAD) provided generousfinancial support for the study programmeand also for the purchaseof the QUANT software for thermodynamic and transport propertiesof fluids. Particular thanks also go to the University of Dar es Salaamfor granting study leave, support and encouragement during the study period.

It is not possibleto mention every individual and organizationwho assistedbut I amvery muchindebtedto Mr. ReneFlatt of the SwissFederalInstitute of Lausanne,whose untiring contributionin variousacademicaspectsof the project hasbeeninvaluable.In this regard I also wish to mentionProf. StephenRichardsonof Imperial College, London. I wish to thank all organizationswhich generouslyprovided various literature and data used in this study. In particular I would like to mention British Gas Plc, the Energy ResourcesConservation Board, and the NOVA Gas Transmission both of Alberta, Canada for providing experimentaldata; and the Ministry of Water, Energy and Mineral, and Hardy BBT Ltd. in Dar es Salaamfor providing information about the proposed Songo Songo to Dar es Salaam natural gas pipeline.

10

DECLARATION I grant powers of discretion to the University Librarian to allow the thesisto be in in copied whole or part without further referenceto the author. The permissioncovers only single copies made for study purposes, subject to normal conditions of acknowledgment.

11

ABSTRACT Although there are many computercodesavailablefor analysisof fluid transients, is limited. linebreak known be few their to to scope situationsand applicable only. a are Thereis, therefore,still a big potentialfor developmentwork in the subject. Discrepancies between different models which have been developed have mainly centred on the assumptions used in developing the basic partial differential equations of flow, and simplifications;the thermophysicalmodel used;representationof various terms subsequent in the equationssuchasthe friction term; andthe numericalmethod of solution of the basic partial differential equations. A previous model developedby Tiley (1989), overestimated the actual wave speeds and had problems of instability of the solution. A new approach, in which the three basic partial differential equation of flow are derived, based on the assumption of an unsteady quasi-one-dimensional flow of a real gas through a rigid constant cross-section area pipe, and using the Gamma Delta method is used. No further simplification is made on the basic equations. Significant improvements have been made on the type of equation of state, thermodynamic model, heat transfer approximation and friction factor representation. The QUANT software for thermodynamic and transport properties of real gases is used. A flow dependentexplicit equation of Chen (1979) is used to calculate the frictional force and heat transfer is calculated using the concept of recovery factor and adiabatic wall temperature. Numerical solution of the basic equations is performed using the third-order WarmingKutler-Lomax method, the second-order MacCormack method and the method of characteristics. A pc based computer coding with the C language is used.

The QUANT softwarehassuccessfully beenincorporatedwith the programme. The full benefitsof the softwarecould not be realisedwith linebreakproblemsdue to limitation of the range within which it gives output at present, but satisfactory results have neverthelessbeenattained. An improved and more accurateway of calculatingthe break boundarycondition has beenused. A non-uniform grid spacinghasbeenused,which allow fine grid spacing in the vicinity of the break in order to enableaccuratemodelling of the rapid transients occurringin that part. Two differentmodelsfor calculatingthe heat transfer i.e. one for the caseof pipesexposedto the atmosphereand buried pipeshavebeenincorporatedwith the model. Experimental data from full-scale pipeline tests is used to validate the computer models. Results from the computer model simulations show good agreement with the experimental data. The MacCormack method has been found to be unsuitable for modelling transient flow following linebreak in high-pressure gas pipelines. The method of characteristics has proved to be the method of solution for such applications. A better understanding of the flow following a break in high-pressure gas pipes is achieved, especially the decompression behaviour at the break boundary. Data gathered from feasibility studies conducted in the late 1980's for a pipeline in Tanzania is used to validate the steady state analysis model and to simulate a linebreak in the pipeline. Results of the computer simulation are discussedand recommendationsmade on the suitability the pipeline design.

Additional work is recommendedon refining and further testing of the computer programmesandusingthe Gamdelepsmethod which covers all the three phasesregion i.e. gas, liquid and gas/liquid.

12

NOMENCLATURE Description

Units

A

Cross-sectionareaof pipe

m-

a

Wave speed

m/s

c

Pipe wall thickness

in

C

Courant number

Cp

Specificheat of gas at constantpressure

J/kgK

Cv

Specificheat of gasat constantvolume

J/kgK

d

Pipe diameter

m

D

Depth of pipe

m

E

Energy of the system

i

e

Specificinternal energyof gas

J/kg

F

Force

N

f

Friction factor

g

Gravitationalacceleration

ßs2

H

Vertical height/piezometricheight

m

h

Specificenthalpyof gas

J/kg

hBL

Heattransfercoefficientthroughtheboundarylayer

WATR

hcA

Convectiveheat transferto the atmosphere

W

I

Inventory/mass of gas in the pipeline

kg

k

Thermal conductivity

W/UK

L

Length of pipeline/wettedperimeter

m

M

Molecular weight of gas

l

Ma

Mach number

m

Mass flow rate of gas

kg/s

NL

Avogadro's number

-

Symbol

p

Static pressure of gas

Pa

Pr

Prandtl number

Q

Heat transferrate per unit volume

J/m3s

q

Pseudo-viscosity variable

N/m2

r

Radial distanceof the pipe

m

R

Specificgas constant

J/kg

13

Rc

Recovery factor

-

Re

Reynolds number

-

S

Specificentropy of gas

J/kgK

St

Stanton number

-

T

Temperature of gas

K

t

Time

S

U

Overall heat transfercoefficient

Walk

u

Flow velocity of gas

mIs

V

Volume flow rate of gas

m3/s

W

Work done by the system

i

x

Horizontal distancealong the pipe

m

Z

Compressibilityfactor of gas

-

Greek Symbols

a

Polytropic alphacoefficient

-

ß

Polytropic beta coefficient

-

A

Polytropic gammacoefficient Small changein the quantity

DE

Changein energyof the system

J

S

Polytropic delta coefficient

-

e

Polytropic epsaloncoefficient

-

e

Roughness of pipe

m

6

Angle of inclination of pipe to horizontal

MdiWs

K

Ratio of specificheats Coefficient of dynamicviscosity

kg/ms

v

Kinematic viscosity

m2/s

p

Density of gas

kg/m3

v

Specificvolume of gas

m3/kg

(D

Viscous dissipationfunction

-

x

Thermodynamic quantity/dryness fraction

-

ýr

Conical angleof the pipe

I

cp 0

Free parameter

=

Heat flow into the pipe per unit length of pipe per unit time Frictional force per unit length of pipe

Jim s

w

14

N/m

Subscripts

A

Atmosphere

aw

Adiabatic-wall

BL

Boundary layer

c

Critical

CA

Convectionto atmosphere

D

Darcy's

e

Equalization

eff

Effective

F

Fanning's

m

surroundingmedia

o

Stagnation/initial

p

At constant pressure/pipe

pw

Pipe wall

R

Reduced(ratio to critical value)

r

radial

S

Isentropic

wo

Outer pipe wall

wi

Inner pipe wall

Superscripts

-

mean

15

CHAPTER

I

INTRODUCTION 1.1 TRANSPORTATION OF NATURAL GAS IN PIPELINES Transportation of natural gas in pipelines takes place at very high pressures (typically 35-150bar, though higher pressures are also used). In main transmission lines, the gas mixture normally exists in one phase (liquid phase), thus avoiding problems such as multi phaseflow and reduction in pipe capacity. These high pressure pipelines, normally cover is USA in largest long distances. the The the one, world natural gas pipeline network very is former USSR kilometres. The the length 1.5 than network total million of more with a second largest in the world.

In the United Kingdom, the national transmission system

kilometres, kilometres 12000 5000 with pressures of up plus another covers approximately to 7bar in the regional distribution systems. In the Netherlands, the high-pressure system kilometres. length 6000 lines total transmission over of with a consists of a series of Following recent discoveries of natural gas in the South Eastern Coast of Africa, plans are highest to to transport the to market with areas. gas now well underway construct pipelines potential. In Tanzania, a 200km pipeline is under construction for transportation of natural in 900km from Salaam Mozambique, Songo Songo Island Dar to pipeline will a es and gas be constructed from Pande to the Republic of South Africa.

The amountsof natural gas stored in the transmissionlines are enormous. In the Netherlandssystemdescribedabove,the yearly flow amountsto over 7x 1010kilogramme of naturalgas. To give a feelingof the amountsand their qualities, a 145 kilometre typical high-pressure pipelinewould containabout 7000 tonnesof natural gas,which is equivalent to 100MW of heat power for 40 days. With such large amountsof naturalgas contained in pipelinesystems,potentialhazardssuchas pollution, fire and economiclossesdue to the spillingof the gas in caseof a rupture are of great concernand needto be thoroughly and accurately assessed. This will in turn be the basis for designing and implementing preventivemeasuresfor suchhazards.One exampleof such caseswas in Alberta, Canada, where it was necessaryto accuratelymodel sour gas (natural gas containing hydrogen sulphide)pipeline rupture, in order to set guidelinesfor land developmentpolicy. At present, there are over a dozen computer codes for analysingtransient flow behaviourof fluids in pipelines.Somework which hasbeendone in the past three decades, on pipe rupture and blowdown in gas transmissionlines, is summarisedin Section 1.2. 16

1.2 RUPTURE AND BLOWDOWN

OF GAS PIPELINES

It hasbeenstatedin section 1.1 that the amountsof natural gas stored in the transmission lines are enormousand the potential hazards,such as pollution, fire and economiclosses due to the spillingof the gas in caseof a rupture, are of great concern. Assessmentof the flow situationfollowing a linebreakin suchpipelinesis crucial, in order to be able to design and implementpreventivemeasuresand to plan corrective action for such hazards. Most pipeline ruptures occur accidentallydue to causessuch as material failure; faultyoperation;damageto pipelinesdue to excavationwork; poor or lack of maintenance, anddamageof pipelinesby falling objectsfrom platforms, shipsand icebergs. Within the last decade,the world has witnessedseveraldisastersresulting from accidentalpipeline ruptures. Theseincludethe Piper Alpha disasteron the night of 61hJuly 1986 and more recently(April andMay 1995)in Russiaand South Korea respectively. The South Korean rupture was causedby pipeline damage,due to excavationwork and it resultedin many deaths.The Russianpipelinerupture is thought to havebeencausedby pipe failure due to corrosion. After the Russianpipelinerupture, flamescould be seenup to 25000ft high and severalsquaremiles of forest were set on fire. Pipe blowdown is a controlled venting to is the atmosphere and usuallydoneprior to shutdownor repair. In this report, both rupture and blowdown are referredto as linebreak. Analysisof the flow following a linebreakin a gas pipelineis done using computer modelswhich havebeenvalidatedwith the best availableexperimentaldata. A review of experimentaldata which is availableon linebreakhas shown that it consistsof data from shock tube tests; laboratory experimentson short sectionsof pipes;linebreaksimulation experiments on full-scale pipelines such as those reported by Bisgaard, Sorensenand Spangenberg (1987)andVan Deeen& Reintsema(1983); rupture and blowdown tests on sectionsof full-scale pipeline such as the Foothills Pipe Lines (Yukon) Ltd. (1981); and measurementson full-scale pipeline systemduring accidentalrupture such as the Piper Alpha disasterduring which somedata was recorded. One essentialrequirementbefore performing computer modelling of a given pipeline

ruptureis that the test data, specificgas data and pipeline systemdata must be preparedin the form requiredby the computerprogramme. Often this involves making a number of assumptions andsimplifications,suchasassumingsomeparametersto be constantfor some interval. Computer software such`as QUANT, PPDS-IUPAC, ASPEN PLUS and

17

PREPROP for calculation of thermophysical properties of fluids are available and make the processmuch easier and accurate, while maintaining a sound computation speed. System data such as pipeline dimensions,is usually included in experimental data. However, some variables such as grid size, friction factor, Stanton number etc. need to be determined. A suitable grid size needs to be chosen in order to produce stable results and to adequately model the physical behaviour such as the shock wave fronts. Transient analysis in the vicinity of the break is very crucial. However, it is often not easy to obtain the required experimental data and with the required accuracy in this area. Consequently, the only available option is to use computer models. In the more recent studies, such as those by Picard and Bishnoi (1989), Tiley (1989) and Chen, Richardson and Saville (1992) variable grid sizeswere used in the vicinity of the break, in order to enable close monitoring of the expansion waves.

The flow of gas following a rupture in a high-pressuregas pipeline is one caseof fluids in pipelines.Suchflows are governedby the well rapidtransientflow of compressible known setof threepartialdifferentialequationsderived from the principlesof conservation of mass,momentumandenergy. The gas propertiesare describedby an equationof state. These,togetherwith appropriateauxiliary conditions, determinethe mathematicalstate of the gas. Many assumptionsand simplificationsare involved in the processof formulating and manipulatingtheseequations. All theseaspectsare presentedin this thesis. Whena ruptureoccurs,the operatingparameterssuchas pressure(p), flow velocity (u), density(p) andtemperature(T) at different parts of the systemvary considerably. For example, T could fall to such an extent that it could render part of the pipeline unusable after the rupture, due to changesin it's material properties. It is thereforevery important to be able to determinethesevariations in order to be able to decide which parts of the pipeline should not be usedand thus avoid further disaster. In modelling a linebreakin a gaspipeline system,questionssuchaswhat is the magnitudeof the flow rate at the broken pipe end, how quickly will it diminish with time, what is the pressure,what is the temperatureandhow are thesegoing to affect the rest of the systemhaveto be answered. Anotherimportantparameteris the speedof propagation of pressurewavesin the system, which also hasa bearingon the fracture behaviourof the pipe material.

18

Resultsof transientanalysisin rupturedpipelines,whether producedby experiments in models, are or computer normally presented a set of graphs. The more commonlyused arethe pressure,temperature,flow velocity,massflow rate, pressurewave speedand mass of gas in the pipelineas functionsof position and time.

1.3 COMPUTER CODES FOR ANALYSIS OF FLUID FLOW IN PIPES Nearlythirty computercodes,for the analysisof fluid transientsituationsin gaseousand/or multi-phase flows, are known to have been used; and some of these are available commercially.TheseincludePIPENET, PIPESIM, PIPESIM PLUS, SURGE, EXPRES, PIPETE, IMF, IPSA, TOFFEA, SIMPLE-2P, PISO-2P, PIPETRAN, TRANSFLOW, TRAC, FLASH, RODFLOW, RELAP, FEAT, OLGA, PLAC, FLUENT, FLOW-3D, PHOENICS,PRO-II, HYSIM, FLOWMASTER, TGNET, PIPE 1, PIPE2, BLOWDOWN andPIPEBREAK. In addition,someorganisations,suchas British Gas [Mallinson (1996)] andGASUNZEof theNetherlands[De Bakker (1993)], have developedtheir own models suited to their own pipelineoperations. FLUENT, FLOW-3D and PHOENICS are general purpose codes for multidimensionaltwo-phaseflow, but limitedto dispersedflow or particle tracking. Algorithms for separatedand stratified flow analysisinclude IMF, IPSA, TOFFEA, SIMPLE-2P and PISO-2P. OLGA and PLAC are one-dimensionalflow codes which do not assume homogeneous two-phaseflow and are therefore suitedto low-speedtwo-phasefluid flow. FLASH, RELAP, RODFLOW andTRAC arenormallyusedfor Loss-Of-Coolant-Analysis (LOCA) in nuclearpower plants. They are therefore suited for situationswhere the twophase fluid is primarily a steam-water mixture. FLOWMASTER, of the British Hydromechanics ResearchGroup (BHRG), can handlefluid systemswhere both rapid and slow changes occur in the boundary condition. FEAT is based on a finite-element numericalmethod of solution. TGNET, PIPE1, PIPE2 and BLOWDOWN are among the group of codes for analysisof transientsfollowing high-pressure gaspipelineruptures. PIPE1 and PIPE2 were developedby Flatt (1985-1989) at the SwissFederalInstitute of Technology, Lausanne. PIPET is for perfect gaseswhereasPIPE2 is for real gases. It is known that at least one companyhasoffered to acquirethe codes. BLOWDOWN is a computer code which can be used to simulate the emergency blowdown of vessels and pipelines containing hydrocarbons. The programme can predict pressure, temperature and multi-phase composition within a vessel or a pipeline; temperatureof the wall; and efflux: all as 19

functions of time. Development of the code [Richardson (1993-96)] started in 1984 at Imperial College, London; and the first version of the code was produced in 1988. It has have it is testing, then not yet and claimed which since undergone constant modifications involved any correction, but have extended considerably the range of problems which it can handle. An independentstudy commissioned to compare BLOWDOWN with PRO-II and HYSIM concluded that it is BLOWDOWN alone which is capable of dealing with systems is BLOWDOWN liquid' than nowone-vessel. containing condensate or comprising more used by many gas and oil companies for the simulation of depressurisation systems. Applications have included a large number of individual vessels on offshore platforms and onshore installations, a number of multiple vessels connected to a common blowdown header,complete platforms and sub-seapipelines. BLOWDOWN is being used throughout the SHELL BLOWDOWN

Group of Companies.

British Gas has decided that it will only use

for depressurisation calculations.

BLOWDOWN

can also be used to

simulate the repressurisation of a system (which is effectively just depressurisation in reverse) and the whole programme is common to repressurisation and depressurisation.

The main physicaldifferencebetweenthe two is that the fluid tends to get warmer (and muchwarmerif the processis fast enough)in the former, while in the latter it tends to get colder (and much colder if the depressurisationis fast enough). PIPEBREAKis a computerprogrammewhich was developedand is being usedby British GasPlc for linebreakanalysis.The GASUNIE pc model for hazardassessment can model gas outflow for complex pipeline networks with different elementsand different outflow scenarios,all in one model. It can model linebreak,venting and leakage. It can handleelementsdefining different boundaryconditions,which can represent,for example, the simplifiedbehaviourof a compressor. The emphasisin developingthis model was put on user friendliness,robustnessand the ability to model complex networks. The basic relations are solvedusing implicit finite-differenceschemesand a graphicaluser interface makesinputting of the network very easy. An estimateof the accuracyis claimedto be 5 to 20 per cent,andin the GASUNIE context, this is consideredto be sufficient for hazard analysispurposes. There are numerousstudiesin which attemptswere made to develop computer codesfor analysisof transientsin rupturedhigh-pressuregas pipelines. Perhapsthe earliest publishedmodelfor the analysisof gaslinebreakis that by Sens,Jouveand Pelletier (1970). Therehavesincethenbeenseveralattemptsto developcomputer modelsfor gas linebreak analysis,including that by Tiley (1989) which was the starting point of this study. Sens, Jouve and Pelletier (1970) claimedthat results from their model were identical with their 20

experimental results.

Arrison, Hancox, Sulatisky and Banerjee (1977) compared

blowdown data for from RODFLOW of a their code with experimental predictions heat heated two loop two sections and containing two pumps, recirculating water in figure-of-eight geometry. a exchangers arranged

Groves, Bishnoi and Wallbridge

(1978) developed a computer model to calculate decompression wave velocities in natural (1976). data Groves the The of experimental model was validated with gas pipelines. Cheng and Bowyer (1978) used their quasi-one-dimensional unsteady compressible fluid flow code to simulate two cases of transient flow. Transients caused by a sudden pipe Lyczkowski, duct left hand three the system were predicted. steam side of a rupture at Grimesey and Solbrig (1978) presented comparative results of their alternating gradient method with analytical results and also their experimental results.

A study by Alberta PetroleumIndustry, GovernmentEnvironmentalCommittee (1978) reported on two existing hydrogen sulphide isopleth (constant concentration of hydrogen sulphide lines) prediction models, one blowdown model and a simplified blowdown model developedduring the study. Predictions from the blowdown models were Industry, in [Alberta Petroleum a subsequent study validated with experimental results Government Environmental Committee (1979)]. Knox, Atwell, Angle, Willoughby and Dielwart (1980), presentedresults of their theoretical model and compared them with their Svrcek (1980) data. Cronje, Bishnoi presented their adiabatic model and and experimental include data (1976). Other Groves its the studies of experimental compared results with those by Fannelop and Ryhming (1982), Van Deen and Reintsema (1983), Bisgaard, Sorensenand Spangenberg (1987), Lang and Fannelop (1987), Picard and Bishnoi (1988 Sj, (1991), (1989), (1989), Kunsch, Botros, Jungowski Weiss Lang 1989), They en and and Imide (1993). (1991 Olorunmaiye Fannelop 1995) and and and and Whereas the physical phenomena of a pipeline system response following a rapid transient event is well explained for the other cases such as pump trip, rapid valve closure etc., the contrary applies for the caseof a linebreak, despite all the studies mentioned in this section. It is stated that when a break occurs in a high-pressure pipeline, the pressure drops virtually instantaneously at the break and rarefaction waves are transmitted up and down the pipeline and are rapidly dissipated when the fluid in the pipe is a gas. Expansion waves are created in the section of the pipe upstream of the break and flow reversal occurs in the section of the pipe downstream of the break. This study investigates further the linebreak

phenomenon,especially at the break boundary with the aim of achieving a better understandingof the phenomenonand a more accurateand stablemodel.

21

The processof developinga computer model for analysisof transientflow of gas in pipelinesentailsthree main steps,namelyformulation of the basicgoverning equations, numericalsolution of the equationsand finally validation of the model with experimental data. Discrepanciesbetweendifferent modelswhich havebeendevelopedfor analysisof transientsfollowing a linebreakin gas pipelineshave mainly centred on the assumptions in developing the basic partial differential equations of flow, and subsequent made simplifications; the thermophysicalmodel used; representationof various terms of the equationssuchasthe friction term etc.; and the numericalmethodsusedfor solution of the basicpartial differential equations. Following a critical review of the previous study by Tiley (1989) a new approach, basedon the GammaDelta method developedby Flatt (1989), is used. The three basic partial differential equationsof flow are derived for unsteadyquasi-one-dimensionalflow of real gasesthrough a non-rigid variable cross-sectionarea pipe. The only subsequent simplification are assumptionsof constant cross-sectionarea and rigid pipe walls. The QUANT software for thermodynamicand transport propertiesof real gasesis used. The softwareis basedon the virial equationof stateand also containsthe isentropicgammaand delta coefficients(ys and SS)required for the GammaDelta method. A flow dependent explicit equation is usedto calculatethe friction losses. Numerical solution of the basic equationsis doneusing the methodof characteristics;the second-orderMacCormack and third-order Warming-Kutler-Lomax explicit finite-differencemethods.

1.4 PREVIOUS WORK AT CITY UNIVERSITY This work [Tiley (1989)] was an attempt to developa computer model which can simulate a linebreak occurring in a gas pipeline. The procedure involved derivation of the basic simultaneous partialdifferentialequationsof motion which mathematicallymodel unsteady fluid flow in a pipeline,by assuminga one-dimensionalhomogeneousflow, constantvalue steadyflow friction factor,constantvalue Stanton numberand neglectingminor lossesand changesin cross-sectionareaof the pipe. The equationswere solved numericallyusing the methodof characteristics.A computercodewaswritten in FORTRAN 77 and a Gould PN 9005 mainframecomputerwas usedto solve the mathematicalmodel. Finally, the model was tested using experimentaldata. It is reported that although the theoretical model successfullysimulatedthe rapidexpansionof a gas following a linebreakin a high-pressure pipeline,comparisonof the model resultswith experimentaldata which was available,did

22

highlight someareasof concern. Theseinclude overestimatingactual wave speedsand problemsof instability of the solution. Further work was recommendedin the two areas mentionedaboveand also further refinementand testing of the model. In this study,the whole approachusedby Tiley (1989) hasbeencritically reviewed. The subjectswhich havebeenexplored include the following: (a)

One-dimensionalflow assumption.

(b)

Homogeneousflow assumption.

(c)

Neglecting minor lossesand changesin cross-sectionareaof the pipe.

(d)

Classificationof the partial differential equations.

(e)

Use of a constantvalue steadyflow friction factor to calculatefriction losses.

(f)

Use of a constantvalue Stanton number to account for heat transfer (simplified model of heat transfer).

(g)

Use of the method of characteristicsfor numerical solution of the basic partial differential equations.

(h)

Use of FORTRAN 77 and Gould N 9005 mainframecomputer.

(i)

Estimationof the gasconstantsspecificheat at constantpressure(Cr), specificgas constant(R), critical temperature(T,,) and critical pressure(Ps).

(j)

Physicalmodellingof a linebreaksituation.

(k)

Assumptionthat pipe walls are inelastic.

(1)

Determinationof compressibilityfactor. (4n) Limitation in the amountand scopeof experimentaldata. (n)

Simplificationof the basicequations.

The one-dimensional flow approximationwas usedon the basisthat for high Reynoldsnumberflows, asin gas transmissionlines, the approximationhasbeenshown to be very good for steadyand slowly varying flows. This was despitethe fact that in the region of influence on flow of bendsand fittings and particularly where there are large, rapid changesin conditions,as in the caseof a linebreak,larger discrepanciesare expected from the one-dimensional approximation. An alternativeto this approximationis to apply the method used in turbulent boundary-layertheory which is lengthy and complicated. Someworkers haveusedthis method. A review of their work hasbeendone in order to find out the limitationsof usingthis method and the significanceof errors introduced in the model by assuminga one-dimensionalflow.

23

Homogeneousone-dimensionalflow meansthat property distribution is uniform be is to includes This the profile, which assumptionabout velocity cross-section. over flow. The homogeneous by Tiley The in detail assumeda study presently. examined more factor friction flow investigated. Darcy's was used this steady assumptionwas validity of flows. It for defined factors friction been basis had transient was gas there the that no on little flow for factor friction flow transient error very causes stated that the use of steady it but frequency low flow was also and amplitude, variations are of relatively when the be disturbances large, may error are occurring, a significant rapid admitted that when incurred. The latter fact was consideredwhen selectingthe steadyflow friction factor and flow (Reynolds in in From the this of review study, made, also subsequentcalculations. friction, it is strongly felt that the use of constantvalue number)andfrequency-dependent in introduce flow factor for the friction flow transient error considerable quite could steady inaccuracies in fact be the the observedwith causes of oneof major modelresultsandcould the model. A considerableamount of work on frequency-dependentfriction in transient fluid flow simulationhasbeenidentified. This hasbeenreviewedbut not usedin this study becausethe effectsof unsteadyfriction, in addition to steadyfriction, are very small. The factor for friction friction factor issues by Tiley twothe such as concerning raised other phaseflow and approximationof the friction term havealso beencritically reviewed. The threepartialdifferentialequationsdescribingthe transientflow were classified in fact differential "first-order hyperbolic "non-linear they equations", are partial while as in hyperbolic differential However, this partial equations". variation quasi-linear impact had the classification of equations no on the solution and results obtained. An investigationhasbeenmadeon the relative magnitudesof the minor lossesand changesin in in cross-section area a practical situation order to assessthe significanceof the error introducedby neglectingthem. The decisionon whether or not to neglectthem hasbeen reviewed basedon the findings of the above investigationand other factors affecting the developmentand executionof the model. A constantvalueStantonnumberwasusedto accountfor heat transferon the basis that variationsin Stantonnumberwith flow rate are not sufficient to warrant the additional computationinvolved,andthat the heattransferterm in the basicequationsis comparatively small. Also in deriving the equation for the heat transfer term, it was assumedthat in a rapid transient, the pressurechangesoccur instantaneously,allowing no time for heat transferto take placebetweenthe pipe and the surroundings(adiabaticflow). It hasbeen argued by others that no models are known for the general case of heat transfer in non-steady flows, where there is no thermal equilibrium between the pipe and the 24

literature has The revealedthat there hasbeenprevious works surroundings. searchmade done in A has been Section 2.6 and, as a result, the heat this these review on subject. of transfermodel hasbeenimproved. The Berthelot equationof statewas usedto determinethe compressibilityfactor, since it was said to producecomparativelyaccurateresults for gasesand vapours at low temperatures. It was statedpreviouslyby Thorley and Tiley (1987), that sincethe terms involvingcompressibilityfactor in the basicequationsareusually relatively small, a complex equationof statewould be uneconomicalin terms of computer time. Therefore, only the relatively simple equationswere examined. Also equations of state were preferred to compressibilitychartsbecausethe former can easilybe programmedinto a computer and canbe solvedfor derivatives. The equationof state for a real gas and the thermodynamic identity were used before applying the Berthelot equation of state. In this study, other approaches andequationsof statewhich could be usedto produce better results havebeen considered. The isobaricheat constant(Cr), isocharicheat constant(C'), specific gas constant (R), critical temperature (T) and critical pressure (pJ of the gases used in various experimentsemployedto validate the model were estimatedby assumingconstantvalues for pressureandtemperaturerangesencounteredfor CP;dividing the universalgas constant by the mean molecular weight of gas for R; Grieves and Thodes method for T,; and Rausnitzand Gunn methodfor p, A recommendationwas madefor this to be one of the areasto be experimentedin, in an attempt to improve the model. In this work, different methodsof estimatingtheseparametershavebeenreviewed and the QUANT software for thermodynamicandtransport propertiesof fluid hasbeenchosen. Also test data in which as manyas possiblegas constantsare specified,hasbeenselected. The assumptionthat the pipe walls are inelasticwas madewithout any discussion to justify its validity. It is known that the transientbehaviourof a pipeline systemdepends upon the elastic propertiesof the conduit (including conduit size,wall thicknessand wall material), external constraints(including type of supports) and the freedom of conduit movement in the longitudinal direction. Thesehavebeendiscussedin Section 2.2.3. In deriving the basic equationsthe terms a(pu2A)/ax and a(puA)/ät were neglectedon the basisthat they are very small. Studiesare known which have arguedfor the retention of each and both of theseterms. These have been reviewed in this study and it has been decidedthat they shouldboth be retainedbecausethe error intoducedby neglectingthem is very significant.

25

Tiley reviewed some of the popular numericalmethodsused for modelling fluid transientsandselectedthe methodof characteristics. The main merits of this method were saidto be the ability to handlesmalltime stepsrequiredfor rapid transients;relatively more accuracy,andability to posethe boundaryconditions properly. It was also statedthat the method of characteristicsis more suited to single pipelines. Since the conclusionof the studyby Tiley, there hasbeenmore studiesemployingand comparingdifferent numerical methods of solution available. It was therefore found necessaryto update the review carried out by Tiley, in order to decideon the most suitablenumericalmethod. A

novel feature of Tiley's model, compared to those which were previously availableis that it allowsfor the possibilityof flow reversaldownstreamof the break as well ashandlinggrid size reduction in the vicinity of the break. Both thesefeatureshavebeen adopted in the model developed in this study. Tiley's main code consists of two programmes,one performinga transientanalysison a given shock tube or single pipe, and the other converting the required section of the numericaloutput of the first programme into graphicalform.

26

CHAPTER 2

THEORETICAL DEVELOPMENT OF THE BASIC EQUATIONS FOR UNSTEADY FLOW 2.1 INTRODUCTION Unsteady flow of compressiblefluids in pipelines is describedby a set of three partial differential equations,derived from the principles of conservationof mass (continuity (equation law Newton's of motion second of or equation), conservationof momentum First Law of or energy and conservation of motion or momentum equation) Thermodynamics(energyequation). The gas propertiesare describedby an equationof determine These the mathematical together auxiliary conditions, state. with appropriate in involved Many the processof the assumptionsand simplifications are state of gas. formulating and manipulatingtheseequations. The form of the equationvaries with the assumptions made regarding the condition of operation of the gas pipeline and any simplifications made. Based on the above the equations may, in the more general firstbe hyperbolic; linear, and or classification, quasi-linearor non-linear; parabolic or second-order. It is generally preferred to keep the equations as simple as possible, without significantlyreducingthe accuracyof results in a particular model, in order to economise on computationallabourandtime andalsoto minimisethe computer memory requirement. This rule is very effectivefor modelling of slow transients,but there seemsto be not much room for manoeuvrein the caseof rapid transients,and especiallylinebreakproblems. In the latter case,the equationshaveto be usedalmost without any simplification in order to achievethe required accuracycriteria. In most of the studiesreviewed,it is generallyfelt that researchershavemadevarious assumptionsand simplifications;usedone form of the basicequationsor the other; representedthe various terms in particular ways; in order to suit a particular methodof solution or application and vice versa. Suwan and Anderson (1992)arguedthat alternativeformulations,interpolations,friction force representation,or time integration,which maybe appropriatefor parabolicproblems,will all violate the basic information carrying physics of a hyperbolic problem. Most cases of unsteady one-dimensionalflow, where disturbancepropagationvelocities do not vary significantly, by quasi-linearhyperbolicpartial differential equationsfor continuity and arecharacterised momentum. On the other hand, complex phenomenasuch as stratified and intermittent stratified-bubble(slug) flows require a two-dimensionaltransientanalysis.

27

Various analyticaltechniqueshavebeenused to reducethe numberof equations before employing the relevant numerical procedure. Van Deen and Reintsema (1983), for example, introduced a technique which reduces the energy equation to a single parameter-in-massequation without the assumption of isothermal or isentropic flow. It was stated in section 1.3, that the processof developing a computer model for transient analysis of gas pipelines entails three main steps, namely formulation of the basic governing equations, solution of the equations and finally validation of the model with experimental

data. This chapterdealsentirelywith the first step. Essentiallythe sameapproachas used by They (1989) is used but an improvement is made on the assumptionsused and representationof the various terms, with an overall objective of achievingmore accurate and numericallystableresults. Tiley (1989)assumeda one-dimensional homogeneousflow in an inelasticpipe and neglectedminor lossesand changesof cross-sectionareaof the pipe. Three simultaneous non-linear(Though the equationshavebeenclassified by Tiley (1989) as non-linear,they are in actual fact quasi-linear)partial differential equationswhich mathematicallymodel pressuretransientsin a non-perfectgaswere derived. A constantvalue steadyflow Darcy's friction factor was used to calculate frictional losses. Heat transfer into the pipe was approximatedby a constantvalue Stanton number. The Berthelot equationof state was consideredto be the most suitablefor this application. A formulation using the general equation of state for a perfect gas and the thermodynamicidentity from the Joule-Kelvin effect was used. Since the developmentof the Tiley model and even before, there has been a considerable amountof researchactivity on theoretical modelling of fluid transients,and in particular those aspectsrelatedto linebreakin gas transmissionsystems. A great number of these have been thoroughly and critically reviewed and consideredin this study. Followingthis reviewa new approachbasedon the gammadeltamethoddevelopedby Flatt (1989), is used. The three basic partial differential equationsof flow are derived for flow of real gasesthrough a rigid constantcross-section unsteadyquasi-one-dimensional area pipe. No further simplification is made on the basic equationsand the QUANT software for thermodynamicand transport propertiesof real gasesis used. The software is basedon the virial equationof state and also containsthe coefficientsrequired for the gammadelta method. A flow dependentexplicit equationis usedto calculatethe friction losses.Numericalsolutionof the basicequationsis effectedusing explicit finite-difference methods.

28

There are three basic approaches by which fluid flow problems are formulated; namely the control volume, integral or large-scale analysis; infinitesimal system, differential or small-scale analysis; and experimental or dimensional analysis. In all these cases, the three basic conservation laws of mechanics i. e. conservation of mass, momentum and

energyplusthe thermodynamic equationof state and associatedboundaryconditions must be satisfied. A control volume is a finite region chosencarefully with open boundaries which mass,momentumand energyare allowed to cross. A balanceis madebetweenthe incoming and outgoing fluid and the resultant changeswithin the control volume, and details of the flow are normally ignored. In the infinitesimalapproach,the conservation laws are written for an infinitesimal system of fluid in motion and become the basic differential equationsof the fluid flow. To apply them to a specific problem, one must integrate these equations mathematically over the volume, subject to the boundary conditionsof the particularproblem. Exact analyticalsolutionsare often possibleonly for very simple geometriesand boundary conditions, otherwise numerical integration on a digital computeris used. Experimentalstudy may be full-scale or laboratory based. The differential approachcan be usedfor any problem but in practice the lack of mathematicaltools and the inability of digital computersto model small-scaleprocesses makes the application of the differential approach rather limited. Even computer programmessometimesfail to provideaccuratesimulationbecauseof either the inadequate storageor inability to model the finely detailedflow structure characteristicsof irregular geometriesor turbulentflow patterns. The dimensionalanalysisapproachcan similarly be applied to any problem but lack of time and money and generality often makes the experimentalapproachlimited. The control-volume analysisis consideredthe most useful of all the three approachesas far as practical engineeringapplicationsare concerned. It givesusefulresultsin a reasonableamountof time, though sometimesgross and crude. In applicationssuchasmodellingof a linebreakin high-pressuregas pipelines,the differential infinitesimal or approachis normally usedand later validatedwith experimentalresults. In fluid flow problemsthe dependentvariables,i.e. p, p, u etc., are functions of the independentvariablesLe. t and x, in the caseof one-dimensionalflow. In any given flow situation,the determinationby experimentor theory of the fluid propertiesas a function of position and time is consideredto be the solution to the problem. There are two distinct fundamentalmethodsof specifyingthe flow field, namelythe Eulerian description and the . Lagrangian description. In the Eulerian form, the independentspacevariablesrefer to a coordinatesystemassumedto be fixed in spaceor translatingand through which the fluid is moving. The flow is characterised by a time-dependent velocity field which is to be found 29

by solving the initial value problem. In the Lagrangian form, the independentspace in fluid fixed the the motion and all and to undergoing system variablesrefer a coordinate distortion of the fluid, so that the particlesof the fluid are permanentlyidentified by their Lagrangian variables, while their actual positions in space are among the dependent variablesto be solvedfor. Although the Eulerian and Lagrangian forms are essentially equivalent, the Lagrangian form gives more information i.e. it tells where each bit of fluid came from initially and has the virtue that conservation of mass is automatic. This results in for (1957) in Ritchmyer the that stated considerablygreater accuracy some problems. in for is form Lagrangian the generallypreferred someproblems a one space abovereasons, Lagrangian in For time, the method two and or more spacevariables problems variable. decreases difficulties. In the seriouslyas usually accuracy particular, encounters serious time goeson, dueto distortions,unlessa new Lagrangianpoint-net is definedfrom time to time, which requiresdifficult and usually rather inaccurateinterpolations. From this point is form is form Eulerian However, Eulerian the the almost useless moreattractive. of view, if thereareinterfacesbetweenfluids having different thermodynamicpropertiesbecauseit found is be kind fluid for to telling at a given of provides no simple mechanism which instantat a givennet-point(for problemsin one spacevariable,sucha mechanismcan easily be provided). Many schemeshad been tried, some combining the features of the LagrangianandEulerianformsbut therewas no satisfactoryuniversalmethod that had been found for generalmulti-dimensionalproblems. Fashbaughand Widawsky (1972) stated that the Eulerian formulation is usually used in steady-statefluid flow problems and the Lagrangian formulation is used more extensively in unsteadyflow and is more desired for solution of shock propagation problems. Also the fact that the Lagrangian formulation yields more information, such as where each bit of fluid originates initially, facilitates locating temperature contact surfaces propagated in the pipe. The study compared two methods namely the pseudo-viscosity and the Lax-Wendroff numerical method in modelling the effects of duct area change on shock strength in onedimensional viscous air flow.

It was concluded that a one-dimensional variable area

Lagrangeanalysisis adequatefor predicting shock flow through a duct areaincreaseof up to at leastarearatiosof 10 to 1. In addition,the Lagrangianformulation using the pseudoviscosity numericalmethod of solution was recommendedand preferred to the Eulerian formulationwith the Lax-Wendroff numericalmethod of solution for accuratepredictions in very long ducts, where an appropriatefunction for the friction factor shouldbe used.

jo

Ames (1977) stated that until the early 1960's, the only known stable difference for less for form Eulerian than those the Lagrangian the accurate were approximations form. Consequently, the Lagrangian forms were preferred. Both systems have disadvantage Eulerian The the the system arises of major same approximately complexity. form different density. Lagrangian The fluids interfaces (shocks) of occur separating when does not have the spatial coordinate mesh fixed in advance and may require refinement of the mesh as computation advances. This possibility of regridding arises since the Lagrange form is constructed so that mass between two successive mesh points is approximately fluid bounded (1988) White that certain of sharply stated numerical analyses conserved. flow, such as the motion of isolated fluid droplets are very conveniently computed in Lagrangian coordinates. Other more recent writers, including Batchelor (1992) argued that

the Lagrangiantype of specificationis usefulin certainspecialcontextsbut it leadsto rather cumbersomeanalysisand in generalis at a disadvantagein not giving directly the spatial gradientsof velocity in the fluid. In fluid dynamicsmeasurements,the Eulerian method is the most suitable. To the probewould haveto movedownstreamat the fluid simulatea Lagrangianmeasurement, particlespeeds.However, this is sometimesdone in oceanographicmeasurements,where flow metersdrift along with the prevailing current. The Eulerian formulation has almost exclusivelybeenusedin all recentstudiesandliterature,evenin such caseswhere according to the abovediscussion,the Lagrangiandescriptionwould seemmore appropriate. For the sake of convenience,the continuity equation has been derived based on the Eulerian approach, while derivation of the momentum and energy equations is based on the Lagrangianapproach. The equationsof flow of fluid in a pipe can be quite complex, especiallyfor cases such as high-pressuregas pipeline ruptures. This is more so if all aspectsof the flow are to be taken into consideration. However, the equationscan be handledconvenientlyby makinga setof assumptions andsimplificationsregardingthe condition of the flow. In this particular casefor example,one hasgot to decidewhether to use a one-dimensionalor flow approximation,singlephaseor multi-phaseflow, include structural multi-dimensional effects of the system,consider minor lossesand changesin cross-sectionarea etc. The inclusion of viscous and thermal effects has beenwell discussedby among others Tiley (1989) and Issa (1970), and therefore need not be repeatedhere. It is therefore being assumedthat in the equationsbeing formulated,viscous and thermal effects are included i. e. non-isothermalviscous flow. Further simplificationscould be made,for exampleby neglecting body forces etc. Various workers have made different assumptionswhen 31

developingthe basic flow equations,dependingon their particular flow situations. They (1989)assumeda one-dimensional homogeneous flow in an inelasticpipe, neglectingminor lossesandchangesin cross-sectionarea,when developinga model for pressuretransients in a rupturedhigh-pressuregas pipeline. The sameassumptionswere madeby Goldwater and Fincham(1981) and Osiadacz(1987). They all arrived at the sameset of equations. White (1988) gave a detaileddescriptionof the derivation of similar equationsfor threedimensionalflow andMoe andBendiksen(1993)usedequationsfor three-dimensionaltwophaseflow to develop their model. Chen, Richardsonand Saville (1993) presenteda set of equationsfor a transienthomogeneoustwo-phaseflow model.

2.2 MAJOR ASSUMPTIONS AND SIMPLIFICATIONS 2.2.1 DIMENSION

OF FLOW

In one-dimensional flow, the components of the fluid velocity in the circumferential and radial directions are ignored. A flow can be assumed to be one-dimensional if the rate of changeof gas properties normal to the streamline direction is negligible compared with the rate of changealong the streamline. This means that over any cross-section of the pipe all the gas properties can be assumedto be uniform. The assumption of one-dimensional flow gives satisfactory solutions to many problems where either the cross-section area changes slowly along the path of the stream of gas, the radius of curvature of the pipe is large compared with its diameter or the shape of the velocity and temperature profiles are approximately constant along the pipe. For one-dimensional flow of a fluid; p, p, u etc. are only functions of t and x. The one-dimensional model enables derivation of the basic equations of pipeline flow in a simple fashion. In real fluid flow situations, especially in the case of high-pressure natural gas flow in pipelines, the flow cannot be truly one-dimensional because viscous effects will produce a velocity profile across the pipe with the local velocity zero at the pipe wall and reaching a maximum in the centre. Moreover, the flow is turbulent so that there are random motions

superimposedupon the meanflow. Ansari (1972) investigatedthe influenceof including radialflow on the solution of unsteadypipe flow equations. It was shown that neglecting the radialvelocitycan lead to substantialerror in the determinationof the axial velocity of flow. The effectof the assumptionof negligibleradialflow was also investigatedin. relation to anotherassumptionwhich is commonlymade,i.e. the gradient of the axial velocity in the

32

in in flow, is direction, the to gradient radial negligibly small compared viscous axial direction. It was found that for slow transientsthe latter assumptionis valid while the former is not and vice versafor the caseof rapid transients. The departurefrom the onedimensionalflow assumptionwill be even more pronouncedwhere there are bendsand fittings in the pipeline. Nevertheless, despite the foregoing, the one-dimensional for be has been in to steadyand transmission good very shown systems approximation gas in highlarge in flow flows. When the are and rapid, such as variations slowly varying from larger discrepancies the one-dimensional are expected ruptures, pressuregaspipeline approximation.A morerigorousapproachis to apply the integratedequationmethod used in the turbulent boundarylayer theory. Therearemanystudiesandcomputermodelsdevelopedfor multi-dimensionalflow flow FLOWfew based A as such models computercodes on multi-dimensional analysis. 3D arecommerciallyavailable.However,thesemodelscan only be appliedto somespecial caseswith confidence.Shin(1978)discussedthe extensionof the one-dimensionalmethod of characteristicsto two space dimensions. The extended method uses the same simplificationsand retainsthe similar simplicity and efficiency as in the conventionalone. The two-dimensionalmethodis applicableto both the cartesianand axisymmetricsystems includes the conventionalmethod as a specialcase. Wylie (1982) presenteda similar and flow fluid for low-velocity transient two-dimensional numerical method analysis of problems.

The method contains similarities to the one-dimensional method of

but does follow the traditionalcharacteristicstheory for two-dimensional characteristics not problems. The most up to date and comprehensivestudy is that by Moe and Bendiksen (1993), which presentedthe physicalbasisof a new type of multi-dimensionaltwo-fluid model, particularly suited for transient flow problems. A basic differencebetween this modeland the other modelsis in the solution procedure,aiming in particular at improved predictions of transientproblems. The numericalschemeis basedon an extensionof an earlier one-dimensionalmodel and employsan implicit finite-differencescheme. It follows from the abovediscussionthat the three-dimensional is approximation the natural method for the linebreak problem being modelled. However, for the sake of simplicity, the one-dimensionalmodel is used with the aim of extending it to threedimensionsin the future. As shown in the abovementionedstudies,the extensionof onedimensionalmodelsto two- and three-dimensionalmodelscould, in principle, be achieved especiallyif caution is madein representationof the other terms of the equations. 33

2.2.2 FLOW PHASE When a high-pressurenatural gas pipelines is ruptured, the temperature falls hence condensationoccurs. A multi-phaseflow situation is therefore likely to be encountered. Multi-phaseflow in generalis a very complicatedphenomenon.It is convenientlyclassified into flow regimeswhich may be further grouped into two main types, namelydispersed (particles,bubblesand droplets) and separated(stratified, annularand elongatedbubbles) flows. More complexflow regimesoften occur as combinationsof these,suchas stratified and annular flows with entrainmentand slug flow. Flow regimes that are commonly encountered are stratified (wavy or smooth), annular dispersed,intermittent (slug) and dispersedbubbleflows. Theseare shown in Figure 2.1. In relatively long pipelines,with possiblylargepressurelosses,severalof theseflow regimesmay exist simultaneouslyas a result of changingin situ flow rates and physical properties of the fluids. A gas-liquid mixture may be treated as a pseudo-fluid,if the mixture and its motion may be treated as homogeneous. Dispersedflow regimeshavebeenquite extensivelystudied recentlyusing two- and three-dimensionalmodels. For two-phaseflow, however, at leastthree generalpurpose codes are available namely FLUENT, FLOW-3D and PHOENICS. For separatedor stratified flow, there has been severalstudies,including those by Oliemans(1987); Wu, Pots, Holdenberg & Meerhoff (1987); and Moe & Bendiksen (1993). In the multidimensionalMoe-Bendiksenmodel,the generaltwo-fluid equationshavebeenappliedwith the assumptionof a singlepressurefield. The modellingof constitutive laws at the interface is not generalbut focusedon separatedor stratified flows. A volume equationis applied for the pressure, enabling a direct two-step solution procedure. Oliemans (1987) investigated the accuracyof two-phaseflow trunk line predictionsby a one-dimensional steady model for stratified wavy flow in horizontal and inclined pipes. Wu, Pots, Holdenberg and Meerhoff (1987) investigatedthe exact location of the transition from stratifiedto non-stratified flow for a high-pressuregas-condensatesystemin a large pipe both experimentallyand using the Taitel-Dukler mechanistic model and the Wallis onedimensionalwave theory. According to the one-dimensionallinear wave stability criterion for high-pressure flow in an 8 inch pipe,unstable(non-stratified) flow occurs at much higher gas-condensate liquid fractionsthanaccordingto the Taitel-Dukler flow regime map. This is confirmed by

34

Superficial Water Velocity, Vgl NÖÖOO

b

,

ft /sec

N 2r Zi



00

:1 z': 0 z p

On 0G0-. o

OC

--b-

0o

0 000 o

s mG y

1 r: 49, ....

0

ý,.

w

co

00 o0.

0o

rlt

vcD

G7 0 N

o

.,

(CD O C)

tj

rt. .,

n ---P

ýý Gý

--º

oC'

-'

-4' o N CD 0

o co

O

C

ýt

""s Tt

Rf

ýý ,'fi t

l4

'ýjý

ý1

C C

Source: Govier and Aziz (1972) Fig. 2.1 Flow Pattern for Air-water Mixture Flowing in a Pipe

35

measurements in a pipe using natural gas and live condensate at a pressure of 75 bar. Experimentally, at the theoretical transition to unstable waves, stratified flow shifts to a flow pattern featuring pipe wall wetting and liquid entrainment.

The transition to

intermittent flow was observed at still larger liquid loadings. It was concluded that a proper theoretical description of this flow regime transition most likely requires

a more

sophisticated modelling with due account being taken of two-dimensional effects, liquid entrainment and non-linearity. For the intermittent flow boundary, some promising results have been obtained by determining the unstable solution of the one-dimensional transient stratified model.

In analysingpressuretransientsin bubbly air-water mixtures; Padmanabhan, Ames andMartin (1978)useda homogeneous modelwhich consistsof one-dimensionalequations of conservationof massfor each of the phasesand conservationof momentumfor the mixture. BhallamudiandChaudhry(1990) developeda third-order accurateexplicit finitedifferenceschemefor transientflows in one-dimensionalhomogeneousgas-liquid mixtures in pipes. The gas-liquid mixture was treated as a pseudo-liquid. Chen, Richardsonand Saville (1993) used a two-phase model called the two-fluid model, in which the local instantaneous conservationequationsareformulated for eachphase. In the model by Tiley (1989), the fluid was assumedto be homogeneousand no provision was made in the representationof the various terms such as friction, for the possibility of a two-phase mixturebeingpresent..For the sakeof simplicity, the sameapproachis usedin this study. However, significant errors could be introduced in the solution because of this simplification.

2.2.3 FLUID STRUCTURE

INTERACTION

The classicalwaterhammertheory only predicts the extremeloading on a system long as as it is rigidly anchored. When a piping systemhas certain degreesof freedom severe deviationfrom the classicaltheorymay occur due to motion of the system.Pressurewaves exertforceswhich causea compliant systemto move. As a result of the motion, pressure waves are formed. This phenomenonis known as Fluid Structure Interaction. Fluid StructureInteractionis essentiallya dynamicphenomenon,the interactionbeing causedby dynamicforceswhich act converselyon fluid and pipe. The forces in sucha piping system are classifiedin two groups namelydistributed forces and local forces. Typical examples 36

of distributed forces are fluid pressureand friction. In the caseof fluid pressure,rapid pressurefluctuationscausea pipeto expandor contract therebycreatingaxial stresswaves in the pipewall. The stresswavesin return generatepressurefluctuations in the enclosed fluid, resulting in a coupling known as Poisson Coupling. Similarly, friction force is responsiblefor friction coupling. In most practical systems,frictional coupling is weak. Local forces act on specific points in a system,such as elbows, teesor valves and cause axial and lateral motions which generate pressure waves in the fluid resulting in an interactioncalledjunction coupling. Junctioncouplingis generallydominantcomparedwith Poissonandfriction couplings. Lawooij and Tijsselling (1990) conducteda study of Fluid StructureInteractionin compliantpipingsystems.They concludedthat the most significant couplingmechanismsare Poissonandjunction coupling. A simpleguidelinewhich states interaction is important is formulated in terms of three time scalesi.e., the time in when which the pressureis built up, the eigenperiodsof the structure and the time scaleof the waterhammerwaves. Theyalsoconcludedthat Poissoncoupling is important for the fluid when significantmodesof interaction are dominatedby stiffness. To consider the motion of the pipe, three main displacementsare distinguished namely axial displacement,lateral displacementand rotation. Therefore in general, the dynamics of a pipe systemare influencedby four wave families, i. e. axial, flexural and torsional waves in the pipe wall and pressurewaves in the fluid. Theseco-exist during transience and have different degreesof influence on the transient behaviour. Fluid Structure Interaction in compliant piping systemsis modelledby extended waterhammer theory for the fluid and by beam theory for the pipes. In buried pipelines, the lateral restraint is usually sufficient to ensurethat the overall behaviour is dominatedby axial effects. Hence most analysesof Fluid Structure Interaction have focused on the propagationof axialwavesi. e. pressurewavespropagatingalong the walls of the pipe. In suspendedpipelineshowever, few lateral restraintsexist and account should be taken of flexural and torsional waves. In a studyon Fluid StructureInteractionin flexiblecurvedpipes, Stittgen (1990)concludedthat the influenceof structural motion

and Zielke

on the internal pressureis smaller

than the time dependentinternal pressurebut can still be of influencein particular cases. They also concludedthat pressurewaves in flexible pipes are quite dependenton viscoelastic wall properties. Rachid and Stuchenbruck (1990) stated that the mechanical propertiesof the material,the temperatureand the degreeof stiffnessof the supports,have 37

a significant influence on the response of the system. Stiffer pipe materials i. e. those with mechanical damping will experience higher pressures than more compliant pipe systems. They also found that Poisson effects can induce significant piping motion especially in long pipe reachesand cause high frequency peaks for pressure and stress in visco-elastic pipes, e.g. up to 25% above the rigid pipe model. However, mechanical damping of the pipe material tends to degrade the solid wave front quite quickly.

In most cases, high

frequenciesgeneratedby the Fluid Structure Interaction mechanism are virtually eliminated before the first fluid cycle.

In practicethe valueof the sonicvelocityin the fluid in a pipeis influencedby the elasticityof the confiningwalls and the compressibilityof the fluid. As the elasticity of the wall materials increases,the effective value of sonic velocity decreases. This effect is commonly neglectedfor typical high-pressuregas pipelines. However, some workers including Wood, Dorsch and Lightner (1966), Zielke (1968), Hirose (1971) and Beauchemin and Marche (1992) have all taken into consideration the effect of pipe elasticity. This has resulted in an additional term in the continuity equation. In the Beauchemin-Marchemodel, no simplification was made on the basic equationsand in addition the effect of variable cross-sectionarea was included. Despite this, the model could still be implementedwith acceptablecomputationalcosts. In the study by They (1989) the effects of Fluid Structure Interaction were not considered. The sameapproachis followed here in order to maintain simplicity of the modelandalsobecausethe review madedoesnot concludeany significantjustification for incorporating Fluid Structure Interaction in the model. A consideration of the wall elasticityis takeninto accountwhendevelopingthe basicequationsof fluid flow. It should be noted however, that in the study by Tiley (1989) the pipe was assumedto be inelastic i.e. elasticityof the pipe walls is negligiblecomparedwith the compressibilityof the fluid.

2.2.4 MINOR LOSSES AND CHANGES IN CROSS-SECTION Anothereffectworth consideringin the formulationof thebasicequationof the flow, is that arisingfrom minor lossesandchangesin cross-sectionarea. Swaffield (1967 and 1968-69) investigated the influenceof bendson fluid transientswhich were propagatedby a rapid valve closure in a long pipeline. He used a short dimensionalanalysisand a seriesof experimentaltestswhich were designedto record the pipe diameterto thicknessratio, pipe 38

elasticity to fluid bulk modulus ratio, restraint applied to the pipe, bend radius of curvature to pipe bore ratio, included angle of the bend and initial flow velocity. He concluded that the influence of a bend is solely dependent on its geometry. Empirical formulae were developedto expressthe reflection and transmission of a transient at a bend in terms of the bend radius of curvature to pipe bore ratio and included angle. Otwell, Wiggert and Hatfield (1985) investigated the effect of elbow restraint on pressure transients. They observedthat if an elbow was fully restrained,there was no alteration of a pressure transient travelling through the elbow. However, if an elbow was not fully restrained, there was a significant alteration of the liquid transient. The alteration was related to the direction and amplitude of the motion of the elbow and it was therefore thought to be dependent on the mechanicalcharacteristicsof the piping and pipe support structure. A numerical model was developed, which accurately predicted pressure and velocity responses.

Sincein this casewe are dealingwith relatively long and straight lines (rather than networks), the assumptionthat minor lossesare negligible is reasonable. Tiley (1989) assumedthat the abovelosseswere smallcomparedwith the distributedfrictional lossesand therefore neglectedthem. The sameassumptionis used in this study. Generally,most workershaveassumedthat the cross-sectionareaof the pipe is not a function of the axial distanceor varies slowly. In somespecialcasessuch as in cooling ducts etc. where the cross-sectionareavaries appreciablywith distance,the effect hasbeenincorporatedinto the basic equations. However, in somestudieson compressibleflow, including thoseby Zielke (1968) and Flatt (1989), the equationsare basedon a variable changein crosssectionarea. Tiley (1989) neglectedchangesin cross-sectionareaof the pipe. Also in this study, it is assumedthat the cross-sectionarea of the pipe is not a function of the axial distance. However, the basic equationsof flow have been derived with provision for varying cross-sectionareawith axial distance.

2.3

DERIVATION OF THE BASIC EQUATIONS FOR UNSTEADY FLOW OF A COMPRESSIBLE FLUID IN A PIPE

Derivation of the basicequationsdescribinga homogeneouscompressiblefluid flow in a pipefrom first principlesis very well known. For one-dimensionalflow in an inelasticpipe with minor lossesandchangesin crosssectionalareaneglected,They (1989) and Osiadacz (1987) derived the three conservationequations,in the conventionalform using different 39

derived (1977) Zucrow both Nevertheless, they similar arrived at equations. approaches. the basic governing equations for unsteady one-dimensional and quasi-one-dimensional isentropic flow and also for generalised flow. White (1988) described the derivation of similar equations but for the case of a three-dimensional flow.

Literatureon derivation from first principlesof more complexphenomenasuchas the one being proposed in this study is rather scarce,but the equations concerned have been stated by various workers who used them. In this study, the basic equations for unsteady quasi-one-dimensionalflow of real gases through non-rigid non-constant area pipes, using the gamma delta method, are used. The gravity term, is also included. These equations are valid for three fluid forms namely, real (or perfect) gas, homogeneous liquid/vapour mixture and liquid.

This form of the equations simplifies considerably computer codes and is

conceived to deal with problems where several of the three fluid forms appear simultaneously. The approach is rigorous from the thermodynamic point of view and takes into account the compressibility of the fluid and considers the vapour as a real gas. Formulation of the energy equation for unsteady flow of fluid in pipes has commonly contained either specific internal energy or specific enthalpy. Each of the above mentioned dependent variables is related to the other dependent variables p, p and T by a caloric equation of state which is often a complicated non-linear empirical correlation in integral form.

This procedure sometimes involves as many as 20 or more gas dependent

coefficients. With the help of the two non-dimensional coefficients y and S, it is possible to eliminate the specific internal energy and specific enthalpy from the energy equation, resulting in considerable computing economy.

The threebasicequationsof conservationfor unsteadyflow are derived from first principles,assumingthat the cross-sectionareaof the pipe varies with t as well as with x. The approach used in deriving the basic equationsis very similar to that used by Flatt (1993a). The major differencebetweenthis derivation and that of Flatt is that the gravity term has beenincluded,whereasFlatt neglectedthe gravity term. The centre line of the pipe is assumedto be stationary. Dependingon which one is more convenient,either the Euler concept which is denotedby (c3.. /at) and (a../ax) or the Lagrangeconcept of full derivativeswhich is denotedby (d../dt) areused. It is possibleto transformfrom one of the two approachesto the other, by using the following mathematicalequivalence:

-"-"u-

dt

(2.1)

at

ax 40

In this derivation, the control volume is regardedto be fixed in space(according to the Euler concept) and attachedto moving particles of the fluid (according to the Lagrange concept). Continuity

equation:

The continuity equationis derived using the Eulerian approach. An infinitesimal control is is in depicted for Fig. 2.2 used. volume a pipe with varying cross-sectionarea,which

LZ.,

v

s

®v

ýj svv

DUCT

It+

,..

WALL

AT

dt UCT

WALL

AT

t

LA --u

------------------

Coritrol

Volume

\

Fig.2.2 Sketch Showingthe Control Volume at the Beginning of a Time Step Massbalanceis performedacrossthe control volume. The total massof the control volume with length dx after time dt is the sum of massflux into the control volume, massleaving the control volume and mass flux through the fictitious stationary pipe wall. In equation form, this is represented as follows: A+A l

2x/

1e at

[APu -A pu

a(axu)A '

where

L= wettedperimeter Ur= radial velocity of the pipe

41

"(

)Pur (2.2)

x

at

BA

t+clt

at

at tt

x

dt

at

Fig. 2.3 Sketch Illustrating Equation (2.2)

following Thetwo parameters through the arerelated equation: aA

(2.3)

Lur .

at

Neglecting the small terms i.e. those containing (dx)2, substituting equation (2.3) and dividing throughout by Adx, equation(2.2) reducesto the following equation: ap

aApu

at

A ax

ap

ax

au

aAp

ax

A at

Rearrangingand using equation(2.1) we get the following equation: 8p.

dp

äp

at

ax

.

au_äA _p

dt

8x

äx

pu_öAp A at A

ý2 4ý

Let the last two terms of equation (2.4) be represented by the notation 1;. Therefore

--P

aAl

LA'

at AA

ax

(2.5)

For rigid and cylindrical pipes ý=0. - The continuity equationis obtainedby substituting equation(2.5) into equation(2.4), resulting into the following equation:

ä?

.u

ap .

(2.6)

PL

42

Momentum equation:

The momentumequationis derivedaccording to the Lagrangeconcept. The forces acting on the control volume are shown in Fig. 2.4.

al)

dx

PA -

Fig. 2.4 Conical Control Volume Illustrating Body Forces Applying Newton's second law of motion to the control volume, in x direction, the following equationis obtained: (A pdx)

[PA

du

. pA dt

.

a(ý'A)dx ax

dx 9 dx odl pAgsin .L ax

Substitutingthe value of dl = dx/Cosl, and simplifying we get the following equation: Ap dx

du

ap --A

dt

aA

dx

-p

ax

ax

dx .pA

8x

dx --

dx cos 4r

0 dx pAgsin -

Dividing throughout by Adx we get the following equation: du r dt

8p -_

er

-

w A cos liº

9 Pgsin -

Using equation(2.1) for du/dt anddividingthroughoutby p, we get the following equation: all at

ap '1

p ax

. ti

au ax

.-"-

pAcos 4r

grin e

For a uniform diameterpipe, Cosgr= 1. 43

(2.7)

Energy equation: The Lagrangian approach is used in deriving the energy equation, and heat transfer into the control volume in x direction is neglected. This assumption is very common for gaseous flows.

il

dx

dt

x.

e

Fig. 2.5 Control Volume Illustrating Heat and Work Transfer Applying the first law of thermodynamics to the control volume, the following equation results: 0. pAdxde pAudt

Au.

a(Mu)

dr-

ax

pdxdt

pLdxu,.

dt

-pAgudxsin0(2.8)

The specificinternalenergy at stagnationconditions,ea,is given by the following equation: eo =e+

'/2u2

(2.9)

Dividing throughoutequation(2.8) by Adx and using equation(2.3), we get the following equation: de dt

°u

aP

pu _ 8x A

äA 8x

P

au cox

P aA A at

n_pg A

rt me

Rearrangingequation(2.6) for au/ax we get the following equation: 44

(2.10 )

au

lp1dp1

8A

ax

dt

A at

u öA A öx

(2.11)

Similarly, rearranging equation (2.7) for ap/caxwe get the following equation: du

ap ax

"_" ar A cos 4r

.-p

(2.12)

pgsine

Substituting equation(2.9), (2.11) and (2.12) into equation (2.10) we get the following equation: de p.

(

du pu-

dt f

du pu-

-l

dt

(pdt pdp

dt

pöA + Aat

""

(a u

Acos'4r

puaA . A ax

pugsin0

paA _ Aat

.0_ A

pu äA --

A ax pugsin0

Simplifying and rearrangingresultsin the following equation: de dt

cri Acos 4r

pdp p dt

A

(2.13)

A

The specific internal energy is related to specific enthalpy through the following thermodynamicrelationship:

e-h-p

(2.14 ) P

Differentiating equation(2.14) with respectto t we get the following equation: de dt

A1 di

dp

(2.15)

.pdp dt p dt p2

Both the continuity and momentum equations i. e. equations (2.6) and (2.7) contain

only p, u and p as the dependentvariables and neither of these equationscontain the variable h. It would therefore be convenient to eliminate h from equation (2.15) by is it in This the the replacing with variablesused continuity and momentumequations. achievedby using the gammadelta method, which was developedby Flatt (1993b) and which is -derivedin the following part of this section. For small changesin the domain as is the casein this model, y and 8 are assumedto be constant. 45

In deriving the basic relationship for the gamma delta method, the polytropic coefficients of Dzung (1944) which are described in Section 2.4.5, for the particular case of isentropic process are used. The coefficients are as follows: as "

(2.16) v(aT)

Rs-

(2.17) ar

p

s

ap

Ys'

-v

av)

p

(2.18) s

The three coefficients are related to each other by the following relationship: Ps = -as Ys

(2.19)

Flatt (1989) derived a relationshipfor the coefficient 8s which is expressedas follows: 8s.

1-

1 -1-vIaTJ Tl öv as

(2.20) s

Another expression relating the specific enthalpy to the isentropic alpha and beta coefficientswas derived by Flatt (1985b) and is as follows: dh = (1- as) v dp + ßs p dv

(2.21)

From equation(2.20) 1 -as.

as

I "1as-

1

as-

(2.22)

1

and from equations(2.19) and (2.20) YS

PS -

(2.23)

as -1

Substitutingequations(2.22) and(2.23) into equation(2.21) we get the following equation: as dh .

as-

YS 1

vdp .

8S-

46

1

p dv

Rewriting using p we get the following relationship: as dh -as-1

YS

1

dp 6sp

p

(2.24)

dp

1 p2

is following (2.24) to t, the Differentiating equation equation obtained: with respect A

as

dt

as -1p

Ys

1 dp

as -1

dt

pdp dt p2

(2.25)

Substitutingequation(2.25) into equation(2.15) andsubsequentlysubstitutingthe resulting following into (2.13) the equation: we get equation equation as

dp

Ys

Ss -1

dt

Ss -1p

p

dp

dp

dt

dt

p dp p dt

p

dp

p

dt

wu Acos 4r

Q A

Simplifying we get the following equation: 6s - (6s - 1) dp 6s- 1 dt .

(ys

1

p

dp

wry Acos ilr

dt

6s -1p

p

(2.26)

A

The speedof soundis relatedto ys by the followingequation: a2.

(2.27)

YSP P

Substituting equation (2.27) into equation (2.26), using the equivalencepresentedin following: (2.1) the to transforms the and rearranging; energyequation equation ap

at

f

uap

ax

_Q2raP

.

lJlJ at

uaPl

(as_

1)

Ira

A

ax

.ul

(2.28)

cos4'

The continuityequation(2.6) containsthe partialderivativeof p with respectto time, while the momentumequation(2.7) containsthe partial derivative of u with respectto time. For convenience of numerical solution of the three equation of conservation, the energy equationis rewritten suchthat the term containingpartial derivative with respectto time, would be that of p only. This condition could be achievedby substitutingthe continuity equation(2.6) into equation(2.28). The resulting equation,which is usedin the numerical solution is as follows:

47

ap

at

It ap

ax

au

2

das 1) äx . -

+aP

1(

l\ n+

(aIt _)i J.

Cos

(2.29)

a

Equation (2.29) is simpler than equation (2.28). However, in deriving the characteristic and compatibility equationsfor the numerical method of characteristics, equation (2.28) has to be used instead of equation (2.29) because the latter fails to produce a unique solution. This is explained further in Section 4.3.2.

2.4

THERMODYNAMIC

AND TRANSPORT PROPERTIES

2.4.1 INTRODUCTION The threebasicequationof conservation[Equations (2.6),(2.7) & (2.29)] must be written with dependentvariablesin sucha way that the solution of pressure,velocity, densityand temperaturecould be obtained.This is in general,commonlydone by using the equations of stateandsomeotherthermodynamicrelationships,which are definedby the model used to representthe transientevent. An equationof state is an expressionwhich interrelates propertiesof substances, andwe distinguishbetweentwo typesof equationsof statenamely the thermal and the caloric equationsof state. A gas which obeysthe ideal gas thermal equationof state(EquationB-Iin the Appendix)is saidto be thermally perfect; and if it has constantspecificheats,it is saidto be caloricallyperfect(the terms ideal gas and perfect gas are sometimesusedas synonyms). Sections2.4.2,2.4.3 and 2.4.4 deal with the thermal equationsof statewhile section2.4.5 dealswith the caloric equationsof state. Therearemanydifferentwaysto approachthis problem,for exampleVan Deen and Reintsema (1983) and Tiley (1989) used the general equation of state for a real gas, P-

ZpRT

(2.30)

andthe thermodynamicidentity basedon the Joule-Kelvin effect [Zemanskyand Dittman (1983)]; dh

.Cp

dT

T81d e.

1p

(2.31)

Pp

to obtainthreeequationscontainingthe compressibilityfactor, Z, and its partial derivatives. The latter were then substitutedin terms of pressure,temperatureand velocity; using the Berthelot thermal equationof state. 48

Generallythe thermodynamicprocessduring a transientevent can be represented by two major categorisations,namelythe perfect or real gas model and the isentropic or non-isentropic(polytropic) decompressionmodel. Thereis a wide selectionof equations in in The to of statevarying accuracyandcomplexity. choiceof which equation use a given situation is determinedby the type of gas that is being modelled, the level of accuracy required, the temperatureand pressurerangesthat are likely to be encounteredand the amountof availablecomputerspacethat caneconomicallybe used. This subjectis covered in the following sectionsof this chapter,and someof the availableequationsof state are presentedin Appendix B. The perfect gas equationof state holds at low pressuresin the vapour and gas regions. The perfect-gasmodel, whether isentropic or non-isentropicis perhapsthe simplestand thus it hasbeencommonpractice,when predicting transientsin a ruptured gas pipeline to assumethat the fluid will behaveas a perfect gas. The basic premisein makingthis approximationis that any non-idealeffect that may occur will only affect the initial stagesof decompression,and for the majority of the releaseprocess,the fluid will behaveas a perfect gas. Significantdeparturescan, however occur betweenthe perfect gas and the real fluid decompressionbehaviourdue to condensationeffects. Resultsfrom studiesby PicardandBishnoi (1988 and 1989) and Flatt (1985-1989) show that the perfect-gasmodel can underestimatethe transientpressurein gas pipeline the massflow rate of gas escapingfrom the rupturesby morethan20%, andunderestimate rupturedpipe by up to 50%. One of the major weaknessesof the perfect-gasmodel is its inability to predict the sudden drop of sound velocity that occurs at the onset of condensation. However, the perfect gas assumption has been used by several workers including Lyczkowski, Grimesey & Solbrig (1978); Lang & Fannelop (1987); Kunsch, Sj oen & Fannelr p (1991); Lang (1991); Chen, Richardson & Saville (1992); and Olorunmaiye & Imide (1993). Lang (1991) admitted that the assumption of a perfect gas leads to errors in the flow rate of about 25% but made the assumption in order to demonstratethe possible use of the spectral method. Kunsch, Sjuen and Fannelop (1991) made the assumptionof an ideal gasjust for the sake of simplicity. They expressed the nonideality of the gas by a compressibility factor, which was taken to be 0.82 for the pressure levels encountered in the investigation. It was claimed that the important results such as mass flow rates were qualitatively well described by the ideal gas model. Another aspect to be considered in using: an equation of state for natural gas is the composition of the different substancesthat make up the gas. For example, Flatt (1985), instead of regarding

49

differing locally C. H. diverse composition, methane of gases natural gas as consisting of (CH4) was used. Approximately 90 to 95% of natural gas consists of this substance whose for is known. The perfect-gas model not considered thermodynamic characteristics are into decompression takes instead but model, which also a real-fluid application in this study, is the mixed nature of natural gas, used. consideration

2.4.2 THERMAL

EQUATIONS

OF STATE FOR REAL GASES

The term equationof staterefersto the equilibrium relation, in the absenceof specialforce fields,betweenpressure,volume, temperatureand composition of a substance,whether it be quite pure or in a uniform mixture. In terms of a functional relation it is expressedas, (2.32)

(p, v, T, x) -0 ,

An equation of state may be applied to gases, liquids and solids. An equation of state for high important is ionisation dissociation for very where and which a real gas accounts temperaturesare encountered. Over five dozen equations of state are known to exit, which Van Der Waals include liquid-vapour These liquid, the the regions. vapour and represent (1873), Clausius (1880) Dietrici (1899), Berthelot(1899-1903),

Wohl (1914), Keyes

(1917), NBS-National Bureau of Standards (1923) Beattie-Bridgman (1928), Keyes, Smith-Gerry (1936), Lennard-Jones-Devonshire (1937), Benedict-Webb-Rubin (1940), Hlrschfelder-Bird-Spotz (1949), Redlich-Kwong (1949), Bloomer-Rao (1952), Martin-Hou (1955), Hirschfelder-Buehler-McGee-Sutton (1958), Pings-Sage(1959), Storbridge (1962), Costolnick-Thodos (1963), Flory-Orwall-Vrij

(1964), McCarty-Stewart (1965), Goodwin

(1967), Martin (1967), Soave-Redlich-Kwong (SRK) (1972), Peng-Robinson (PR) (1976), Van Reet-Skogman (1987) and the Kamerlingh-Onnes equations of state. There are two general approaches to the development of an equation of state, namely the theoretical approach and the empirical or semi-theoretical approach. The theoretical approach is basedon either the kinetic theory or statistical mechanics, involving intermolecular forces. The empirical approach has been covered in detail by, among others, Martin (1967). The empirical approach is more of a convenient method of interpolation if it is applied to a pure substance. It can become dangerously inaccurate when used in unfamiliar situations. Empirical equationshave been used more widely than theoretical ones becauseof lack of data to enablethe application of the latter equations to various practical

50

situations and fluids [Münster (1969)].

However, the theoretical approach gives high

accuracy and a basis for a wide application. In the searchfor an equation of state, one must first make a decision on the amount and kind of data that will be required to obtain the parameters in the equation; range of density to be covered; and the precision with which the pressure, volume and temperature (p,v,T) data are to be represented. For some applications it is desirable to have an equation which can be obtained from a minimum of data i. e. the critical temperature, pressure and volume; or the boiling point; or some molecular parameter. In other applications, an equation is sought which will represent a large amount of experimental p-v-T data within the precision of the experiment. In all pressure-explicit equations, density is paramount. If the same order of precision is to be maintained, a relatively simple short equation will suffice for low densities,whereas a long complicated equation will be required if coverage is to include both high and low densities. High precision equations have many arbitrary constantswhose values depend primarily upon the density range and in a minor way upon the temperature range. For example, an equation representing precisely up to a fiftieth of the critical density may need only two constants, while four or five constants will be necessary when in the region of half the critical density. To continue up to the critical density at least a half dozen or more constants will be required and twice that many will be needed if the goal is one and a half, or two times, the critical density. The Keyes equation of state produces data within the experimental precision for densities of up to a little over half the critical density, while the BWR equation of state, used by Kiuchi (1993), gives good precision for a density range well past the critical density. The Bloomer-Rao equation of state is essentially the same as the BWR equation of state in which the temperature coefficients have been modified to fit nitrogen data. The Keyes-Smith-Gerry equation of state can achieve extremely high precision in calculating thermodynamic properties of steam, but is valid only to a third of the critical density. The Hirschfelder-Buehler-McGee-Sutton

and Costolnick-Thodos equations of state cover a

wide range of densities and simultaneously maintain a precision which is acceptablefor practical applications. In these equations the p-v-T plateau of data is divided into regions and separate equations are written for the different regions. The Flory-Orwall-Vrij

and

Lennard-Jones-Devonshireequations of state possesshigh precision but the range of their applicability is limited to the liquid region where the densities are always significantly greater than critical densities.The Van Der Waals, Clausius, Dietrici, Berthelot, Wohl and

51

infinitely from density the Redlich-Kwong equations of state cover whole range of from deviations large fairly but liquid, experimental with attenuated gas to compressed inaccurate is it but low is pressure, results. The Van Der Waals equation quite accurate at for is the many Dietrici The point critical equation reliable near near the critical point. isotherm far from in incurred the fluids. However, critical other regions errors are organic Berthelot The largely for be it temperatures. equation hence varying used cannot and low for temperatures. at vapours gases and results accurate produces comparatively

and

Thorley and Tiley (1987) and Tiley (1989) considered the Van Der Waals, Dietrici Berthelot equations for application in the analysis of transients in a ruptured high-

Picard be found Berthelot the and to the suitable. most equation pressuregas pipeline and Bishnoi (1987) used the Peng-Robinson (PR) and the Soave-Redlich-Kwong(SRK) fluids in thermodynamic to the single-phase speed of sound calculate equations of state in the two-phase region of a multiand consisting of pure components and mixtures for their SRK PR Both ability the noted of state were equations and component mixture. to provide quite accurate vapour densities and generated reliable equilibrium ratios. However, neither was able to provide liquid densities with the same level of acceptability SRK better the PR than equation of state. the equation of state was somehow although Both the SRK and PR equations of state were weak in predicting sound velocities in liquids, i. for have to vapours e. application shown some promise although eachwas considered to by Kiuchi less 10%. The Van Reet-Skogman than state, used of equation with errors (1993), is considered adequate for reasonable ranges of pressure and temperature.

It is generallyrecommended therefore, that for caseswhere the fluid is transported Also, the simple are required. sophisticated equations at conditions near critical point, more fit data far fail by Redlich-Kwong to the the state equations of state such as equation of 5% density, Even the the the to order of of errors experimental precision. critical within up data fit intolerable. is These can equations of more can occur which considered completely quite precisely over limited rangesbut they definitely are not the tools to correlate good pv-T data within the experimental precision. They cannot therefore be seriously regarded as more than a qualitative tool for predicting p-v-T behaviour. The Beattie-Bridgman equation of state, which contains five adjustable constants, represents with some accuracy the whole range above the triple point. The Martin equation of state, used by Flatt (1985), is the sameas the Martin-Hou equation of state with the addition of the exponential terms, which permit the former equation to go to about 2.3 times the critical density as opposed

52

to 1.5 times the critical density reached by the latter. The equation of Martin has the valuable property that it is also solvable explicitly

in respect of temperature. For the

practical application, the extra expenditure in computation time with the real gas variant was 13% as compared with the corresponding ideal gas computation.

In deciding on the use of one of the equationsof state, or the other, it has been common practice to choose the simplest possible equation. Various workers have justified their choice in different ways. For example Flatt (1985) argued that, in view of the relative simplicity

of the one-dimensional flow model used, it was not necessary to demand too

high an accuracy from the equation of state and therefore used the Martin equation of state. However, he later admits [Flatt (1993-96)] that the use of the Martin equation of state and an empirical C(p) correlation is inconsistent from the thermodynamic point of view, and therefore recommends the use of a theoretical equation of state of the type h= h(p, S). Thorley and Tiley (1987) argued that since the two terms in which the compressibility factor and its derivatives appeared i. e. T

az

Z arp

and

p

az

z

ap

are usually relatively small, a complex equationof state would be uneconomicalin terms of computertime and thereforeused the Berthelot equation of state. In this study the theoretical approachto equationsof stateis used.,Though more complex, it is as accurateas experimentaldata. As a result of this high accuracy, the theoretical approachto equationsof state is well suited for analysisof high-pressuregas transients and it has been used under similar circumstancesby Richardsonand Saville (1991),with satisfactoryresults. As far as possible,real gas propertiesand data are used in the model.

2.4.3 THEORETICAL 2.4.3.1 Introduction

APPROACH

TO EQUATION

OF STATE

to the Basic Principles and Equations

In this approach,the macroscopicpropertiesof the fluid such as pressureand temperature areessentiallydeterminedby the interaction of very many particles(atoms, molecules),in contrast to properties definable for single atoms or molecules, even though the latter properties are usually measuredin systemswith very large number of particles. The 53

help be idea is that the aggregate of which a probability can constructed, underlying with certainprobabilisticstatementscan be derived about the systemunder consideration. The theoreticalapproachto the developmentof the equationsof state hasdevelopedover many years,utilisingboth the principlesof classicaland quantummechanics.From the viewpoint of statistical mechanics,a system is consideredto consist of an enormous number of molecules,eachof which is capableof existing in a set of statesat different energylevels. The moleculesareassumedto interactwith one anotherby meansof collisionsor by forces causedby fields. A systemof moleculesmay be perceivedto be isolatedor, in somecases, consideredto be embeddedin a set of similar systemsor ensembleof systems. Concepts of probabilityareapplied,andthe equilibriumstateof the systemis assumedto be the state of highestprobability. The fundamentalproblemis to find the numberof moleculesin each of the molecularenergystates(population of the states)when equilibrium is reached. It is not the intention of this study to go into the detailsof theseprinciplesbut rather to find a useableform of an equationof statebasedon this approach,that will representthe situation beingmodelledasaccuratelyaspossible. For further detailson theseprinciples,the reader is referred to Münster (1969). On the basisof this theoretical approach,the microscopicdescriptionof a system involvesspecifyingmanyquantities,which are not suggestedby our senseperceptionsand cannotbe measurede.g. configurational,residual,partition,moleculardistribution functions etc. The relationbetweenthe macroscopicand microscopicpoints of view lies in the fact that the few directlymeasurable properties,whosespecificationconstitutesthe macroscopic description are really averagesover a period of time of a large numberof microscopic characteristics. For example,the macroscopicquantity, pressure,is the averagerate of changeof momentumdueto all the molecularcollisions madeon a unit area. In statistical mechanics,an ideal gas is definedas a systemof non-localisedparticleswhoseinteraction energycanbe neglectedrelativeto the total energyof the system. The laws of ideal gas are limiting laws at zero pressure. Their practical significanceis due to the fact that they frequently give sufficiently good approximationsfor simple gasesat ordinary pressures. Deviations from the ideal gas laws, which are due to intermolecularforces, ariseat finite pressures thus the necessityfor real gas considerations. For real gases at ordinary pressures, at leasttwo effectsappear,which can only be includedby the addition of special constants. Theseconstantsare known as the volume correction, which takes account of the volume of the moleculesthemselvesand the pressurecorrection; which takes account 54

of the intermolecular forces of attraction. The Van der Waals equation of state was the first to take into account these considerations. The quantitative results of the Van der Waals equation are generally not satisfactory, but nevertheless represent the behaviour of real gases in a fairly satisfactory qualitative manner. The Beattie-Bridgman equation of state proved itself the best from the viewpoint of convenience of computation.

The most generalequation of state is the Kammerlingh-Onnesequation of state (EquationB-17). In this formulation, the equationis representedas a power seriesin v'. The equationis commonlycalled the virial equationof stateand the coefficientsB, C. D. E, etc. are known as the second,third, fourth, fifth, etc. virial coefficients respectively. The coefficientsvary with temperature,anddependon the nature of the gas and are related to molecularpropertiesof the gas. In the pressurerangefrom zero to about forty standard atmospheres,only the first two terms in the expansionare significant. In general, the greater the pressurerange, the larger the number of terms in the virial expansion. The equationis usedto accountfor the behaviourof realgasesespeciallyat high pressures.The is equation tedious to use but one of its merits is that once the virial coefficients are determined,the behaviour of any real gas can be predicted. An averagereduced virial equationof statecanbe obtainedby introducing the reducedstatevariablesp1,T, and v, in placeof p, T and v. The reducedvirial equationof stateis: 2468

D, Ce

pryr

(1ýBrKe+C/e yr vý

Te

yr

F ýE/Ce v'

e) 8 Vr

(2.33)

where the reducedstatevariableare P,

"pT, P,

(2.34)

"TV-v Tc

vc

and the coefficientsBr, C1,Dr etc. are of the form b2 B,

bl "

-

b3

b4

bs (2.35)

Tr

Tr

Tr

7"r

and the critical coefficient Ke

NLRTe -

(2.36) pcvc

55

The importance of the reduced form of the equation of state

(NL is Avogadro's number).

(Principle of Corresponding States-Refer to Section 2.4.3.2) is due to the fact that many substances have almost the same thermal equation of state in the reduced state variables. The value of K,, evaluatedby Münster (1969) from equation (2.36) above, when pr Tr yr 1 is K. 3.43. The remaining 25 numerical coefficients in equation (2.33) are available from Landolt-Bornstein Tables. Also availablefrom this reference are the coefficients of the Van der Waals equation for various substances.

2.4.3.2 The Principle of Corresponding States

The principle of correspondingstatesis basedon the fact that any equationof state f(p, v, T, A, B, C) =0

(2.37)

which does not contain more than three independentconstantscan by insertion of the critical data be put into the form 60. equation film Re>104,0.7x x PO

ut

x Ax

x

x Ax

x

PQ "

PxeX . 1- (xe

x

x) x

(4.135)

ux

x)

1- (xA

Pxex'

(4.134)

Px

(xe

uxex'

x Ax

x)

(4.136)

px x

Fluid propertiesat position R: If XRa u Positive 3. uf

Upsti-earzi

Pipe

tFi

brealc m

jiH

a V w A

a Direotiora

of

1DCt8 tn[o]

flow

b

(b)

Section

of

Pipe

Dowrtstream

tine

BraI<

Q

Time

3 0 Q Di.

tancs

Fig. 4.19 Variable Grid Mesh Layout be takento ensurethat the respectivedata files for the last time step of the previousbatch are kept intact before starting the next batch analysis. Physical Properties of the Pipe Material: ID[21, ID[10J AND ID[11] ID[2] is the thermal conductivity of the pipe materialand it is usedin the non-isothermal non-adiabaticsteadystate analysisprogrammesand in the transientanalysisprogrammes to calculatethe heattransferthroughthe pipewall. The pipe wall thickness(ID[ 11]) is also required and used in exactly the same way as ID[2].

The inner pipe wall roughness

(ID[10]) is usedby all the steadystateand transientanalysisprogrammesto calculatethe frictional force. Physical Properties of the Surrounding and ID[141

Environment:

ID[1], ID[31,1D[121, ID[13]

The outer pipewall temperatureand temperatureof the surroundingatmosphere(ID[ 13]) are used by all the non-isothermal non-adiabatic steady state and transient analysis programmesto calculate the heat transferthrough the pipe wall and from,the pipe to the surrounding environment. The pipe wall temperatureneedsto be specified only at the 209

upstream end. In situations where this data is not available, an estimated value is used. The estimatedvalue should be based on the temperature of the surroundings or of the gas, and allow for a reasonable temperature drop. If the value of ID[1] estimated is not accurate enough and also for points other than the upstream end, the programmes will calculate the precise value through iterations.

The thermal conductivity of the surroundingmedia(ID[13]) and the depth of the pipe below the surface(ID[12]) apply only in caseswhere the pipe is buried and for nonisothermalnon-adiabaticprogrammes.They areusedto calculatethe heat transfer through the medium, from the pipe outer wall to the surroundingatmosphere. In most casesthe pipes are buried either under ground or water. The surroundingmediain thesecasesare the mediain which the pipesareburied, and ID[ 12] is the depth of the pipe in the medium, below the atmosphere. Pressure of the surroundings(ID[14]) is the pressureof the place at which the broken pipe discharges. If the broken pipe is open to the atmosphere,the surrounding pressurewill be the local atmosphericpressure. In somecases,for examplepipesburied under water, the broken pipe would dischargethe gas under water and the surrounding pressure will be higher than the atmosphericpressure. In such case, the surrounding is pressure given by the following equation:

ID[14] = PA+ (P,,.g.ID[12])

(4.298)

wherepAisatmosphericpressureand pwis densityof the water and ID[12] is the depth of the pipe in the water. Physical Properties of the Gas: ID[15], ID[16], ID[17], ID[20], ID[21], ID[22], ID[23], ID[25], ID[26], ID[271, ID[281, ID[30] and ID[31] Among the above data the primary data which must be specified in all the steady state and transient analysis programmes is the initial temperature, pressure and density of the gas at the upstream end i. e. ID[15], ID[16] AND ID[21] respectively and mass flow rate of gas before the break (ID[17]). With programmeswhich use the QUANT software i. the e. non-

isothermal non-adiabatic steady state and all transient analysis programmes it is not to specifyboth T and p. The most practical way is to specifyT and calculate p necessary using the QUANT software, althoughthe opposite could also be used. This appliesonly when the domain range of the parameterslies within the range in which the QUANT software delivers output. Otherwiseboth T and p haveto be specified.

210

The fluid properties u (ID[20]), n (ID[22]), Z (ID[23]) all at the initial state of the gas and M (ID[28]) are not required for the non-isothermal non-adiabatic steady state and the transient analysis programmes. They are used by the other steady state analysis programmes,and the are not all required by each of the programmes. For the programmes using the QUANT software, all the fluid properties mentioned in this paragraph can be calculated from it and need not necessarily be specified in the data file SYSTEM. DAT.

The fluid propertiesat the final conditions, are the propertiesof the gas at the state in which the gas is expected to be after the transient event has ended. The properties which need to be specified at the final condition are p (ID[25]), (ID[30])

and K (ID[31]).

T (ID[26]),

a (ID[27]),

Z

These together with the pressure and temperature of the

surroundings (ID[14] and ID[15] respectively) are used by the programme for transient analysis after the break as the properties of the gas at the final state, and to which the programme must finally converge.

211

CHAPTER 5

A REVIEW OF SOME EXPERIMENTAL AND NUMERICAL DATA 5.1 INTRODUCTION Experiments on high-pressuregas pipeline rupture are commonly divided in two categories, namely laboratory experiments and full-scale pipeline experiments. Laboratory experiments have traditionally been performed using shock tubes although other variations are also used. Full-scale pipeline experiments representmore realistically the situation being modelled, but are very expensive, hazardous and entail significant practical problems. Data of this type is therefore very scarce compared with laboratory data, and is in most cases restricted. Laboratory experiments, on the other hand, are more practicable but do not accurately representthe full-scale situation. However, results from shock tubes tests, especially those include (to tubes the effects of heat transfer, change in cross-sectional using modified shock area, etc.) are good enough and have been used in the absence of full-scale experimental data.

Laboratory experiments have traditionally been performed using shock tubes although other variations are also used. An alternativeto performing full-scale pipeline experimentswould be to wait until a pipeline rupturesaccidentallyand collect the required data. But sinceaccidentsoccurunexpectedly,areinherentlyunsafe,and often when people ill are preparedto collectthe requireddata,it is not possibleto collect all the requireddata. The other alternative data available, is results from other similar computer models. However, this type of data should be used with great caution. Experimentaldata on linebreak problems and especiallyfull-scale-experimentaldata is becoming increasingly available. Very often, experimentaldata is availablein a form which can not be used directly without modifications,assumptions,extrapolationsand interpolations. Resultof transientanalysisin ruptured pipelines,whether producedby experiments or computermodels,are normally presentedin a set of graphs. The more commonlyused arethe p-x, T-x, u-x, p-t, T-t, u-t, I-t, m-t and p-a curves. A review of experimentaldata which is availableon pipeline rupture has shown that it consistsof data from shock tube tests; laboratory rupture experimentson short sections of pipes; linebreak simulation experiments on full-scale pipelines such as those reported by Bisgaard, Sorensenand Spangenberg (1987)andVan Deeen& Reintsema(1983); rupture and blowdown testson 212

sectionsof full-scale pipeline such as the Foothills Pipelines(Yukon) Ltd. (1981); and measurements on full-scale pipeline system during accidental rupture such as the Piper

Alpha disaster. Oneessentialrequirementbeforeperformingcomputermodellingof a given pipeline ruptureis that the test data, specificgas data and pipeline systemdata must be preparedin the form requiredby the computer programme. Often this involves making a numberof assumptions andsimplifications,suchasassumingsomeparametersto be constantfor some interval. Computersoftware for calculationof thermophysicalpropertiesof fluids such as PPDS-IUPAC,QUANT, ASPEN PLUS andPREPROPareavailableand makethe process mucheasierandaccurate,while maintaininga soundcomputation speed. Systemdata such aspipelinedimensions,is usually includedin experimentaldata. However, somevariables such as grid size, friction factor, Stanton numberetc. needto be determined. A suitable grid size needsto be chosenin order to meet the required quality of results. Transient analysisin the vicinity of the break is very crucial. However, it is often not easyto obtain the requiredexperimentaldataandwith the requiredaccuracyin this area. In that case,the only availableoption is to use computer models. Consequently,a review of the available computer modelsis made,from which suitableresults are to be selectedfor comparison with the model being developed.

5.2 LABORATORY

EXPERIMENTS

5.2.1 DESCRIPTION OF THE SHOCK TUBE TEST A simple shock tube consists of a tube of constant cross-section, in which a diaphragm initially separatestwo bodies of gas at different pressures. Rapid removal of the diaphragm generatesa flow of short duration containing waves of finite amplitudes separated by quasisteady regions. Initially, after diaphragm removal, a shock wave travels into the low

pressuregaswhile an expansionor rarefactionwavetravels into the high pressuregas. The quasi-steadyflow regions inducedbehind thesewavesare separatedby a contact surface across which pressureand velocity are equal,but density and temperatureare in general different. The shockheatingof the low pressuregasand the expansioncooling the highof pressuregas,permita very wide rangeof flow temperatures to be achieved. There is a very wide rangeof researchapplicationsof shock tubes,including the study of the propagation of pressurewavesin unsteadyfluid flow situations. A large numberof resultsfrom shock 213

tube experimentsexist and some of these are specifically on the modelling of unsteady fluid flow, similar to that occurring in ruptured high-pressure pipelines.

Stronger shock waves can be producedwith various modificationsto the simple shocktube. Thesemodificationsinclude driver gas heatingand areareduction from driver to drivensections. The theory and performanceof simpleshock tubeshasbeendescribed by Glass(1958); and Hall (19580) covers modified shock tubes.

5.2.2 REVIEW OF SOME LABORATORY EXPERIMENTS There is a wide range of experimental results from shock tubes available. These have been presented and discussed by Glass (1958), Issa (1970), Thorley & Tiley (1987) and Tiley (1989).

There seems to be not much published recently, on shock tube modelling of

ruptures in high-pressure gas pipelines. Most of the recent work is based on full-scale pipeline experiments, which are becoming increasingly available.

Basedon a three point criteria, describedin Section 5.5, Tiley (1989) selectedthe Groves-Bishnoi-Wallbridge (1978) [referredto asthe Grovesdata] shock tube data and the British Gasshocktubedata[JonesandGough (1981)]. The Groves data was also usedby Cronje,BishnoiandSvrcek(1980)to simulategaspipeline rupture, with successfulresults. In using the Groves data, Tiley (1989) had to makevarious assumptionsregardingshock tube material (in order to estimatethe friction factor and Stanton number);the effective rupture time of the diaphragm and the accuracy and sensitivity of the measuringand recording devices. Therefore, no assessmentcould be made of the experimentalerrors incurred. No suchproblemswere experiencedwith the British Gas data. Lyszkowski,GrimeseyandSolbrig(1978) usedexperimentalresultsfrom ideal gas shocktubetests and blowdown of an ideal gas from a 13ft (approximately4m) long pipe. In the caseof the shocktube tests,the length of the pipe was 20ft (approximately6m) and the initial pressureratio at the middle was 7.9 to 1 [114.7Ib/in' abs.to 14.7Ibrn'-abs.(7.9 to 1 bar)]. The velocity was uniformly zero and the internal energy was uniformly 80.62Btu/Ib(187.5kJ/kg)initially. A diatomicgas having a specific heat capacityratio (K) of 1.4 was used. Resultswere presentedin the form of graphsof pressure,velocity and internal energy; all as functions of the pipe length. The same gas was used in the blowdownexperiment.The ratio of the initial pressurein the pipe to ambientpressurewas

214

2.83 and the initial pressure was 1,000 Ib/in2 abs (69 bar). Results were presented in the form of closed end pressure and break velocity as functions of time after the break.

Ilic (1987)determinedexperimentally the rate of transientmassdischargeof a twophasemixturethrougha shortpipe. An initially stagnantand saturatedcolumn of Freon-12 liquid stored in a 0.25 litre glass vesselwas used. Blowdown was initiated through a 3.2mm bore brass tube, 50mm long by opening a quick acting ball valve located downstream from the pipe. Strain gaugepressuretransducerswith natural frequencyof about 20kHz were usedto measurevapour pressurein the spaceabovethe liquid in the glassvesselandthe pressureof the liquid at the exit of the dischargepipe. A Bourdon tube compoundpressuregaugewasusedto indicatethe initial and final receiverpressures.The fluid temperaturewas measuredat the exit plane of the pipe with a 0.2mm diameter chromel/alumelthermocouplewhosebare hot junction was machinedflush with the pipe bore andthermallyinsulatedfrom the pipe wall by a resin plug (1.6mm diameter)in which the thermocouplewireswereembedded.A similarthermocouplemountedon the discharge pipe wall showedthat the temperaturedrop was insignificantthroughout the blowdown period, exceptduring the vapour stage. The pressuretransducerand temperaturesignals were recordedon a high speedchart recorder. Resultshave beenpresentedin terms of photographstaken at 0.33 secondsintervals,dischargingvesselmassinventory variation with time and dischargetube exit characteristicsand driving pressurehistory. An initial delayperiod was observedbefore the start of ebullition by the formation of vapour at the bottom of the vessel.Discretebubbleseventuallycombinedto form a slug of vapour which displacedthe descendingliquid surfaceand causedconsiderablemixing in its wake. At the endof blowdown,only vapourissuedfrom the vessel. The criterion of choking was taken to be the independenceof the pipe exit pressurefrom the receiverpressurechanges. This condition was confirmedby the simultaneousmeasurementof the fluid exit temperature, which showedsimilar trends to the exit pressure. Haque, Richardson, Saville and Chamberlain(1990) reported on experiments conducted on rapid depressurisationof large pressurevesselsat the Imperial College, London. The aim of the experiments wasto enableunderstandingof the physicalprocesses involvedduringblöwdown. The experiments involved the depressurisationof three vessels ranging in diameter from 5 to 110cm, with length to diameter ratios from 10 to 3 respectively. Most of the measurements were madeon vertical vesselsblowndown from the top andfrom the bottom. However,measurements werealso madeon horizontal vessel 215

blowdown through an exit on the axis of the end closure of the vessel. Measurements recordedincludedpressure,temperatureand composition;all as a function of time during the blowdown process. The vesselswere instrumentedwith pressuretransducersto measurethe internal pressureand thermocouplesto measuretemperaturewithin the gas spacein the vesselandthe surfacetemperaturesof the vesselwall, both inside and outside. Pressurisationwas with nitrogen at 150 bar and depressurisationwas through an orifice down to atmosphericpressure. A data logger was used to record all the pressuresand temperaturesonce every 3 seconds. Haque, Richardson, Saville, Chamberlain and Shirvill (1992) presented and discussedexperimentalmeasurements which were taken to validate their computer model, BLOWDOWN. The experimentswere conductedusing different sizedvesselsoriented verticallyandhorizontally;containinga rangeof different fluids; and with blowdown from the top, bottom andsidethroughchokesof varioussizes.The experimentswere conducted usingthreevesselsof differentsizes,in orderto checkthe effect of vesselsize on the model predictions. The model showed that the predictions are scale-independent.The three vesselsusedwerenamedasVessel1 to 3; andhad lengthsof 3.240m, 1.524mand 0.671m; insidediametersof 1.130m,0.273m and 0.040m; and thicknessof 59mm,25mm and 5mm respectively.Transducerswereusedto measurethe pressurein the vesselsto an estimated accuracyof 10.2 bar. Pressuregaugeswere also usedto give direct measurements.Barewire thermocouples wereusedto measurethe temperatureof the fluid within the vesseland alsoof the insideandoutsidewalls to an estimatedaccuracyof f0.5K. It was also possible to withdraw samplesof the fluid from the top and bottom of the vesselat arbitrary stages during blowdown, and to measure the composition of the samples using a mass spectrometeranda gaschromatograph.In addition, a windowed port was attachedto the upperpart of the vesselanda mirror set at an angleof 45° to the vertical within the vessel. A video camerawasthenusedto view the fluid within the vesselboth horizontally, across the upper part of the vesseloccupiedby gas; and vertically downwards, from the part occupied by gas to the part occupiedby liquid. In this way condensationin the gas and evaporation in the liquid could be monitored directly and continuously during the blowdown. All measurementswere transmitted to data logging systemsand thence to micro-computersfor subsequentdata reduction. The experimentswere carried out at the Imperial College(for Vessel2) and the British Gas test site at Spadeadam(for Vessel 1).

216

ExperimentsusingVessel1 at Spadeadam wereconductedon mixtures of methane, ethane and propane;together with some on nitrogen (for comparisonwith the Imperial Collegeexperiments).Someeighteenand fifteen setsof experiments,which were referred to as SI to S18 and 11 to 115 were conducted at Spadeadamand Imperial College respectively.The experimentscovereda wide rangeof different fluid compositions,vessel blowdown directions College Imperial The and choke sizes. at orientation, experiments were,for safetyreasons,confined to non-inflammablebut representativefluids. Nitrogen was usedasa representativegas phasesinceits critical propertiesare in the samerangeas thoseof methane;andcarbondioxidewas used as a representativecondensablephasesuch aspropane.The initial pressurewas,except in two cases120 bar for the Spadeadamtests; and 50 bar for the Imperial College tests. Blowdown times were of the order of 1,500 seconds for the Spadeadamexperiments and 100 seconds for the Imperial College experiments.Threeexperiments,namely11,S9 and S12; all of which were conductedwith a vertical orientation of the vesselsand blowdown at the top position; were selectedas representativecomparisonsfor the BLOWDOWN model. However, in the caseof the pipelinerupturemodel,experimentscarriedout at a horizontalorientation of the vesselsare morerepresentative.Richardson(1993-96)reportson further experimentalblowdown tests on a 1.6m long, 0.3m diameter vesselat Imperial College; a 3.Om high, 1.3m diameter vessel conductedby SHELL at Spadeadam;and a 0.3m high, 0.15m diametervesselat Imperial College for a numberof flashingflow tests.

5.3 REVIEW OF SOME FULL-SCALE PIPELINE EXPERIMENTS Perhaps the earliestfull-scale pipelinerupture tests are those reported by Sens,Jouve & Pelletier(1970). The test was conductedon an 11,800mlong 0.1065minternal diameter naturalgaspipelineat a pressureof 31.4bar. The resultsin form of p-t and u-t curves were compared with calculatedresults and said to be identical. However, this is one of the simplestcasesof naturalgasflow, sincethe pressurewas relatively low and measurements were taken away from the vicinity of the break. In most cases,the pressurewould be so high that the modelling would involve handling of liquid and gas phases. In a test programme conducted by Alberta Petroleum Industry Government Environmental Committee (1979), an existing 168.3mmoutside diameterpipeline which was typical of sour gas lines in the province was used. The test section,approximately4.0km in length, 217

was burst at the midpoint.

In addition, one rupture was performed on a pipeline of

323.9mm diameter and approximately 7.1 km long. Of the results presented, only the p-t and m-t curves are relevant for this case.

Tiley (1989) used resultsfrom somefull-scale testson high-pressuregas pipelines including those reported by Jones& Gough (1981) and Foothills Pipelines(Yukon) Ltd. (1981). As for the resultspresentedby Jonesand Gough (1981), the tests were conducted diameter 1.22m Tests, and sectionswere with naturalgasandusingshortsectionsof pipe. both The to long 50m pressurised pipe was at ends. sections reservoir with approximately Pressure initiated test 90bar the the section. of at centre and a crack was approximately historieswere recordedat pointseitherside of the break, and resultswere presentedas p-a for Gas; by behalf British for BMI, tests Similar tests and of on carriedout results; curves. initial In for SHELL; the Gas these British cases, that were also presented. conducted by Foothills bar The 120 140 respectively. and experimentsreported pressureswere Pipeline (Yukon) Ltd. (1981) were carried out between 1979 and 1981 at the Northern AlbertaTest Facility. Short sectionsof about 120mlong and 1.4mand 1.2mdiameterpipe between 74bar known to and composition and pressurised were chargedwith naturalgasof 87 bar. Fracturewas initiated at the centre of the test sectionand resultswere presented in the form of p-t and p-a curves. Van Deen& Reintsema(1983)usedexperimentaldata from GASUNIE to validate their theoreticalmodel. In the experiment,a linebreakwas simulatedby rapidly openinga Another lower the test to set of pipe a parallel pipe at a pressure. connected valve which datafrom full-scalepipeline simulation of a rupture is that reported and usedby Bisgaard, Sorensenand Spangenberg(1987), and which was carried out in 1979on a 77.33km gas pipeline from Neustadt through Sörzen to Unterföhring in Germany. Gradle (1984) presentedblowdown curves for natural gas i. e. graphsof blowdown time versuspipeline volume (t-V curves), for different types of pipeline and at line pressuresranging from 13.8bar to 103.4bar. Knox, Atwell, Angle, Willoughby and Dielwart (1980) presented results of sour gas pipeline rupture experiments. Two pipelines of a 3.8km, 168mm diameteranda 7.1km,323mmdiameterwere rupturedat the midpoints of the test sections. m-t curves were presentedfor the two pipelinesat an initial pressureof 69bar and gas temperatureof 10°C. The total massof gas releasedwas also measured.Wilson (1981) alsoreportedon full-scalerupture tests performedin 1978by Alberta Petroleum Industry, GovernmentEnvironmentalCommittee(1979). Botros, Jungowskiand Weiss(1989) used 218

the data obtained during the blowdown of a straight gas pipeline section and a three-unit compressor station for comparison with their models. Results were presented in the form of p-t curves.

Richardsonand Saville (1991) used experimentaldata from one of three ruptures lines of connectedto the Piper Alpha platform, which was monitored sufficiently well to bore length information; data line 100m from and of and also a give useful validatory 150mm, containing LPG. In the latter case,initial pressuresof up to 21bar, and initial temperatures of about293K were used. The blowdown was done through orifices ranging in equivalentdiameterfrom 0.0im to 0.15m(full-bore). Resultswere presentedin the form of p-t and p-T curves,at the intact end of the line; and m-t, M-p curves. Someof these results were also used by Chen, Richardsonand Saville (1992). Chen, Richardsonand Saville (1993) reported on a 100mlong pipe with an internal diameterof 150mm,which was fully rupturedat oneend. Pressureat the intactend and total massin the pipelinewere recorded as a function of time. The pipe was filled with 95% propane and 5% butane mixture at a pressureof l lbar and temperatureof 20°C. Much more data resultingfrom laboratory experimentsand also full-scale tests,but which hasnot beenpublished,exist. However, the data was not availablefor this study. Mallinson(1996)reportsof recentexperimentaldata, which havebeenacquiredby British Gas in addition to the onesreported by Jonesand Gough (1981). The new data were obtained through full-scale pipeline tests, in collaboration with. several other companies.Mallinson(1996) furtherstatesthat the cost of the experimentswas very high, and as such the results are both valuableand also commerciallysensitive. This data was therefore not availablefor this study. Morrow (1996)reportson other recentfull-scalepipelineventing tests, which were performedby a gas transportation company,in order to evaluateleak detectionsystems. The pipelinesystemconsistedof parallelpipelineswith interconnectionsor cross-oversthat could be openedor closed. The interconnectionswere spacedat regular intervalsbetween compressorstations. A partial linebreakwas simulatedby venting gas through a smaller diameter branch line terminatedby a remotely actuatedquick acting relief valve. A full linebreakwas not simulated. P and T data were recorded. This information was received too late to be able to follow up on the possibility of the data being madeavailable.

219

5.4 REVIEW

OF SOME NUMERICAL

DATA

Sens, Jouve and Pelletier (1970) claimed that results from their model was identical with their experimental results described in Section 5.3.

Arrison, Hancox, Sulatisky and

Banerjee (1977) compared predictions by their RODFLOW code with experimental data for blowdown of a recirculating water loop containing two pumps, two heated sections and two heat exchangersarranged in a figure-of-eight geometry. Results were presented as p-t and T-t curves at different sections of the loop and different break sizes.

Groves, Bishnoi and Wallbridge (1978) developeda computermodel to calculate decompression wavevelocitiesin natural gas pipelines. The model was validatedwith the experimentaldata of Groves (1976). Resultswere presentedas p-a curves. Chengand Bowyer (1978)usedtheir quasi-one-dimensional unsteadycompressiblefluid flow code to simulatetwo casesof transientflow. Transientscausedby a suddenpipe rupture at the left hand side of a three duct steam systemwere predicted. Results were presentedas p-t curves. Lyczkowski, Grimeseyand Solbrig (1978) presentedcomparativeresultsof their alternating gradient method with analytical results and also their experimental data describedin Section3.2. The resultswere presentedas p-x, u-x and e-x in an ideal shock tube 5msafter the break; p-t and u-t at the closedend of the pipe for an ideal gas; and p-t for a blowdown of steam-watermixture pipe. A study by Alberta Petroleum Industry, GovernmentEnvironmentalCommittee (1978) reported on existing two isopleth prediction models, one blowdown model and a simplified blowdown model developed during the study. Predictions from the blowdown models were presented as m-t curves and validated with experimental results in a subsequent study [Alberta Petroleum Industry Government Environmental Committee (1979)].

Knox, Atwell, Angle, Willoughby and Dielwart (1980), presented M-t curves

resulting from their theoretical model and compared them with their experimental data reported in Section 5.3. Cronje, Bishnoi and Svrcek (1980) presented graphs of results from their adiabatic model and compared them with the experimental data of Groves (1976). Results were presentedas p-a curves for pure methane and argon. Wilson (1981) presented an M-t and exit pressure ration against time for the full-scale pipeline rupture experiments described by Knox, Atwell, Angle, Willoughby and Dielwart (1980).

Fannelopand Ryhming (1982) presentedseveralgraphswhich model the pressure profilesandreleaseratesfollowing the ruptureof a gas pipeline. No experimentaldata was

220

break high Van Deen in the took the the and pressure end. simulation place at used,and Reintsema(1983) comparedresults from their numerical model with full-scale pipeline experiments whichweredesignedto simulateda linebreak. Graphsof pressureagainsttime were presentedfor the two simulation experiments. Flatt (1985 and 1986), presented data from his However, was no experimental predictionsof model. variousgraphsresulting include S-x, his The to a-x and u-x, p-x, used validate modelpredictions. graphspresented m-x; and m-t curves. Bisgaard,Sorensenand Spangenberg(1987) comparedresultsfrom the linebreak in simulation experimentsreported Section 5.3 with results of their numerical model. Resultswerepresentedasp-t curvesfor the abovecaseand also for simulationof a rupture in a straightpipelineusingtheir numericalmodel. Lang and Fannelop(1987) presentedm-t curves,resulting from simulationswith three different methodsand comparedthem with thosecalculatedby FannelopandRyhming(1982). They also presentedp-x and u-x curves for various times after the rupture. Picard and Bishnoi (1988) compared results from their computer models with experimentaldatafrom the Northern Alberta Burst Tests [Foothills Pipeline(Yukon) Ltd. (1981)]. Graph of pressureratio p/pL (where PLis the pressureat the initial condition) versusexpansionwavevelocitywerepresented.Picardand Bishnoi (1989) used their three modelsto demonstratethe importanceof real-fluid behaviourin modelling of high pressure gaspipelineruptures. A casewasconsideredwherebyafter a suddenrupture of a 168.3mm diametersour-gaspipeline,an upstreamemergencyshutdown valve closes,restricting the blowdown length to 1000m. Pressureof the gas was I lobar. Resultswere presentedas m-t, p-t, u-t a-t, and T-t curves. They (1989) used the shock tube data used by Groves, Bishnoi & Wallbridge (1978)i.e. that of Groves(1976)andthat of Jones& Gough (1981); and also the full-scale pipeline experimentalresults of Foothills Pipeline (Yukon) Ltd. (1981) to validate her computer model. Resultswere presentedas p-a and p-t curves. Botros, Jungowski and Weiss (1989) presenteda seriesof p-t profiles resulting from their models predictions. Comparisonswere made with their field measurementsreported in Section 5.3. Lang (1991) presentedp-x, m-x and m-t curves,resulting from his numericalmodel. Kunsch, Sjeen and Fannelop(1991) presentedm-t, I-t, u-t, p-t and p-x curves produced by their computer model for rupture of sub-seapipeline connectedto processequipmenton the platform through a vertical segment. 221

Richardson and Saville (1991) presented a series of graphs comparing predictions from their model BLOWDOWN with experimental results. The results consist of p-t, T-t and I-t curves at the intact end of the pipe; m-t and mass efflux curves against time at the ruptured end; and massflux through the orifice with upstream pressure. Chen, Richardson and Saville (1992) presented m-t, p-t curves at both the open and intact ends of the pipe. The graphs consist of results obtained from the BLOWDOWN

model and from

experiments. Haque, Richardson, Saville, Chamberlain and Shirvill (1992) presented p-t, T-t and m-t following blowdown of a gas vessel. BLOWDOWN

The data calculated using the

model was compared with experimental results. Chen, Richardson and

Saville (1993) presentedp-t curves at the intact end and m-t curves of the line predicted by their homogeneous two-phase model and also resulting from experiments.

OlorunmaiyeandImide (1993) presentedm-t curves, producedby their model and comparedthemwith resultscalculatedby Flatt(1986) and also Lang and Fannelop(1987).

5.5 SELECTION

OF TEST DATA

In her selection of experimentaldata for comparisonwith the theoretical model, Tiley (1989) usedthe following criteria: (i) The variablesrequiredby the programmemust either be given in the experimental data or be calculatedfrom it. (ii) The experimentalresults must be of a form that can be directly comparedwith the theoreticalcomputeroutput. (iii) Details of the apparatusand procedureare necessaryin order to be able to assess the experimentalerror and evaluatethe resultsobtained. The samecriteria will be used in this study. All the recent data described in Section 5.3 was not available for this study, despite the big effort made to secure it. The discussion of the data which could be obtained and its screening for suitability in this study is covered in this section.

Whena breakoccursin a high-pressuregas pipeline,the pressureof the gas drops virtually instantaneouslyat the break and rarefactionwavesare transmittedup and down the pipeline and rapidly dissipatedwhen the fluid in the pipe is a gas. A flow reversal occurs in the sectionof the pipe down streamthe break. In order to model thesewaves properly,a reducedgrid sizeis normallyrequired in the vicinity of the break. Mass release rates from pipeline ruptureswas studiedby Wilson (1981). The behaviourof gas during 222

its initial release,expansionand plume rise was studied by carrying out an exact momentum analysisof the expandingjet outside the pipe rupture, and developing a more sophisticated model for both the mass release rate and the plume rise of the momentum jet. The most hazardous portion of the gas released during a pipeline rupture is emitted during the first 10 to 20 seconds, when the high initial mass flow causes the largest downwind concentrations. There have been other but simpler models for predicting the initial release rates from high-pressuregas pipeline rupture, including those by Flatt (1986) and Kunsch, Sjoen and Fannelop (1991). Another factor which affects the release rates is the rupture mode. This was studied by Knox, Atwell, Angle, Willoughby and Dielwart (1980). In most cases a full bore rupture is assumed. Wilson (1981) also studied the heat transfer to the moving gas through the pipe walls and observed from experiments almost an isothermal condition throughout the length of the pipe, except for about the last 200 diameters during which the rapidly accelerating flow near the pipe exit was moving too quickly to gain heat from the pipe walls.

To be ableto producecomparable results,a modeldevelopedfor suchas above conditions and assumptions,should be validated with data obtained from experiments performed under similar conditions. Criteria such as the one usedby They (1989), and listed in this section, is very useful in selectingexperimentaldata for comparisonwith theoreticalmodels. The reviewof experimentaldata availablehasrevealedthat the criteria are very seldommet, and consequentlya numberof test setshave to be used in order to validate a particular model. Also someassumptionshaveto be madefor somedata, thus reducingthe accuracyof the data. For example,in using the Grovesdata, Tiley (1989) had to makevariousassumptions regardingshocktube materialin order to estimatethe friction factor and Stanton number;the effective rupture time of the diaphragmand the accuracy be andsensitivityof the measuringand recording devices. Therefore, no assessment could madeof the experimentalerrors incurred. Four differentcategoriesof data for validating computer modelshavebeenreview in this chapter. The main disadvantageof laboratory experimentaldata, including shock tube data, is that of its much smaller scale comparedto full-scale pipelines. The small diameter effects observed by Groves, Bishnoi and Wallbridge (1978), thermal effects includedby amongothers Tiley (1989) and frictional effects discussedby Flatt (1986) are someof the importantfactorsagainstlaboratoryexperimentsdata. Under such a situation, manyof the assumptions suchasfrictionless,adiabaticflow etc. which would not normally apply to full scalepipelines,hold true. This meansthat a model developedunder such 223

assumptionsand validated with such data would produce comparable results, but it would full-scale to applied a pipeline. However, one the positive not produce correct results when is be high laboratory that used, thus compression ratios could aspects of experiments, producing strong shock waves. On the other side, simulated rupture experiments represent but in full the effects encountered a pipeline rupture, a more realistic presentation of all by in Van Deen For the example experiments used operate under very small pressureratios. and Reintsema(1983) and Bisgaard Sorensenand Spangenberg (1987) the pressures varied by approximately 0.2bar in 400 secondand 2.3bar in 2 hours. In a typical high-pressure gas pipeline rupture, the pressureat the break would be reduced to near ambient pressure within full-scale fraction Simulated pipeline experiments therefore represent rather a of a second. slower transients compared to actual rupture experiments.

Full-scale pipeline rupture tests provide the most suitable data for validating numericalmodels. But evenwith full-scalepipelineexperimentaldata, often measurements do not existfor the sectionof the pipelinein the vicinity of the break and at the time interval requiredto showall the importantvariationsof the properties.Somemodelsfor simulating pipelineruptureswere claimedto havepassedsimplybecausethe datausedto validate them based was on laboratory experimentswhich could not representall the effects existing in a full-scale pipeline. The claim by Haque, Richardson,Saville, Chamberlainand Shirvill (1992) that their model predictionswere scale-independent, should be treated with great caution, since the study was on vesselsand not pipelines. In such casesthe effects mentionedaboveare not significant. Numerical resultsfrom other modelsare only useful asa qualitativecomparisonwhendevelopinga computermodel and preliminaryevaluation and comparisonof results. They should not be regardedas final validatory evidence. Although as seenin this review, there is plenty of experimentaldata for validation of computer models, it is not easily available. When available,it is often in a form of printed graphswhich requiresone to translateit into numericaldata. In somecases,the specificationof the gas and test condition is not fully given and evenwhen given, it is not always easyfor workers to obtain the propertiesof suchmixtures. All theseare someof the factors that significantly affect the accuracyand validation of computer models for rupture in natural gas pipelines. Even with a working model, there is a needto regularly update the gas composition sinceas observedby Foothills Pipeline(Yukon) Ltd. (1981), it changeswith time.

224

Computermodelsfor simulationof transientflow of gasfollowing a rupture in highpressurepipelinesmustbe validatedwith suitableexperimentaldata. Although all types of datadiscussedin this paperare useful at somestageof developmentof a computer model, full-scale pipeline rupture tests provide the most suitableand reliabledata for validating computermodelsfor suchevents. Wheneverpossiblethis categoryof data shouldbe used as the final and most reliable validatory test for computer models. The review of experimentaldata availablehas reveal that a numberof test setshaveto be usedin order to provide all the necessarydata required to fully validate a particular model. Also some assumptionshave to be made for some data, thus reducing the accuracy of the data. Although as seen in this review, there is plenty of experimentaldata for validation of computer models,it is not easilyavailable. Basedon this review, four setsof experimentaldata are selectedfor validation of this computer model. The data is that reported by Jones& Gough (1981), Foothills Pipeline (Yukon) Ltd. (1981), Alberta Petroleum Industry GovernmentEnvironmental Committee(1979)and Sens,JouveandPelletier(1970). The test data is referredto in this study as the British Gas, Foothills, API and SNGSO test data respectively. More details about the tests and the data is are presentedin Chapter6.

225

CHAPTER 6

VALIDATION 6.1 VALIDATION

OF THE COMPUTER MODEL

PROCEDURE

It was decidedin Chapter 5 that four setsof full scalepipeline experimentaldata on pipe sectionsof varyinglengths,diametersand operatingconditions will be usedto validate the computermodelpredictions.The tests,which arereferredto as British Gas, Foothills, APT data and SNGSO;arebriefly describedin this chapter. The necessary which gas and system was providedis presentedin tabular form and the experimentaland predictedresults from the computermodelsarepresentedin graphicalform. The resultsare discussedin Section 6.3. The procedureusedin validating the computer model is as follows: (i) All the data provided concerningthe gas and the test systemis stored in relevant files in the QUANT softwareand transientanalysisprogramme. Data which is not is provided, estimatedbasedon someknown parameterof the systemor of similar systems. (ii) A datafile of thermodynamicand transport propertiesof the fluid which covers all the range of dependablevariables(p, p and T) encounteredduring the transient event is producedusing the QUANT software. (iii)

Computersimulationis performedfor the test sectionsdescribedin section4.6 and numericalresults are stored in a data file at intervals for the dependantvariablest and x sufficient to produce good graphicaloutput.

(iv) The computerprogrammeDSORT is usedto processthe data generatedin (iii) in to a form which can be used by standardgraphical packagesto plot the require graphs. (v)

A commercial graphical software EXCEL is used to plot the required graphs, combiningboth the experimentaldatageneratedin (v) and the data predictedby the computer model, which is generatedin (iv). Both setsof data are plotted in the samegraph in order to provide good comparison.

(vi) Experimental data which are available in graphical form are converted into numericaldata and stored in ASCII files which can be usedby standardgraphical packagesto plot graphs.If the experimentaldata is available in numericalform, in ASCII datafiles,thenthis step is not required. However, if the data file is too big,

226

data. In the be this to the all (iv) study, also experimental applied then step may data used to validate the computer models are availablein the form of printed graphs. Four setsof experimentaldata havebeenselectedto validate the computer been has behaviour is believed data used. to transient Only representrapid which model. For all the four setsof data, the pipes are effectively horizontal, and are exposedto the the based to heat exposed transfer Therefore pipes the on model only atmosphere. is atmosphere used.

6.2 COMPARISON OF COMPUTER MODEL PREDICTIONS WITH EXPERIMENTAL RESULTS 6.2.1 FOOTHILLS

TEST DATA

The testswhich were reportedby Foothills Pipeline(Yukon)Ltd.(1981), were carried out between December 1979 and April 1981 at the Northern Alberta Burst Test Facility (NABTF). A total of six tests, which are denoted as NABTFI, NABTF3, NABTF4, NABTF5, NABTF6 andNABTF7, werecarriedout andwere reported. The main purpose behaviour fracture the the of the to the test on effect of gas composition of was examine 1.4m diameters 1.2 243m lengths Shorts total were and and of approximately of of a pipe. bar. between 87 74 known to and composition and pressurised chargedwith naturalgasof Fracture was initiated at the centre of the test sectionby detonatingan explosivecutter. Resultswerepresentedin the form of pressurehistoriesand timing wire data showing the is broken data for The the tip pipe which are presented each section of crack position. denoted as West and East. In relation to the computer model they are referred to as downstreamand upstreamsectionsrespectively. Out of the results presented,of relevanceto this study are the p-t and a-p curves. Although the data can be usedto someextent to validate computer modelsfor linebreak analysis, it was not intendedfor that purpose. The major reasonwhich makesthis data unsuitablefor validating line break models is the fact that the fracture was designedto lengths. direction This the the propagatealong of pipe covering someconsiderable axial break it boundary, difficult break the this to the using model where makes especially model boundaryis assumedto be fixed in the x-t plane.The shorter the axial distancecoveredby the fracturepropagationthe more suitablethe resultsarefor validating the computer model. 227

Two major weaknesseshave been observed in the FOOTHILLS data. The first is

that of p-t curvescrossingeachother. The is secondthat of non-identicalcurves for the two identical halvesof the broken pipe section. During the depressurisationof a broken broken higher intact that the than the the end and the at pipe, end remains pressureat pressuregradient follows the same direction throughout the depressurisationprocess. Therefore, it is not possiblefor the pressureat a point in the pipe which is further away from the break to be lower than that of a point which is nearerto the broken end, during the first transient of the pressurewave following the break. The fact that this rule is not followed by the FOOTHILLS tests results,indicatesthat the data presentedin this report must be usedwith extremecaution. The stateof the gas inside the pipe was stationarybefore the break. Assuminga in length is throughout the temperature the there uniform of a symmetry all the pipe, propertiesof the fluid and the geometryof the pipe about the crosssection at the middle of the pipe. The symmetry should be maintainedeven after the break, if the fracture

is both the sameon sectionsof the pipe,unlessif therearestrongexternal propagation factorssuchaswind whichwould affectthe out flow of gasfrom the brokenends. In the FOOTHILLS results,this condition is also not well followed.

Outer diameter(mm) Wall thickness

1422 13.71

1219 15.24

1219 13.71

P[kPa]]

7446

8687

8687

23-26 Gr.483

(-3)-(-4) Gr.483

18.5 Gr.483

86.59 6.8 4.03 0.262 0.421

85.36 8.22 4.34 0.182 0.278

i-C5H12

85.36 7.68 4.46 0.238 0.331

84.7 8.21 4.38 0.201 0.235

85.19 8.07 4.4 0.203 0.3

85.71 7.94 4.27 0.214 0.311

0.057

0.029

0.032

0.029

0.029

0.033

n-C5H12

0.034

0.028

0.032

0.03

n-C, H14

0.029

0.03

0.008

0.013

0.011

0.008

0.01

0.009

CO2

0.076

1.56

0.049

2.212

1.177

1.1409

1.71 0.016 0016

1.56

1.804 0.013 0011

2.212 2.212 2212 1

1.177 1.177 1-1771

1.1409 1.1409 1-1409

T [C] Pipe Material Gascomposition [%] CH4 C2H6 C,H= i-C4H, 0 n-C4H,o

N2 Ar lo,

1.56 _

Source:Foothills Pipeline(Yukon) Ltd (1981) Table 6.1: Foothills Test GeneralData

228

1422 13.71 74446

1219 15.24

1219 13.71

8697

8143

1819 (-4)-(-5) (-5)-(-6) Gr.483 Gr.483 Gr.483

Two test results, namely NABTFI

EAST and NABTF7 WEST, are selected for

validation of the computer model. The result NABFTI EAST are used to validate the model for flow reversalin the downstreamsection of the broken pipe and the NABFT7 WEST results are used for the upstreamsection.Both the NABTFI and NABTF7 tests were performed with relatively short axial fracture propagation(approximately5m and 18m respectively). The length of the axial fracture propagationis more than three times higher in NABFT7 test than in the NABFTI test. The test results produced for test NABTF7 WEST seemto be more promising for validating of the model. Two main reasonsfor this arethe lower timing wire velocityused for NABFT7 comparedto that used for NABFT 1; and the long distanceaway from the broken end (18.28m comparedwith 4.2m for NABTFI)used asthe shortestdistanceto presentthe results. The shematicof the Foothills test NABTFI and NABTF7 are presentedin Figs. 6.1 and 6.2 respectively.

I

SECTION OF PIPE

GASFILL

UPSTREM OF THE BREAK

UNE

+

BROKEN END

UPSTREAMEND

BROKEN END

DQ HSTRUW ENO

SECTION OF PIPE

A

DOWNSTREM OF THE BREAK

PRESSURE TRANSDUCER POSITIONS DISTANCEor PRESSURETRANSDUCER FROM UPSTREAMEND

Fig. 6.1 Schematicof Foothills Test NABTFI 229

121.5m

la 0 1219mm

t

0

0

08

P76W

PT5W2

PT5WI

Q

UU P74W

PT3W

L

P

PTRW

GASFlu UNE

SECTION OF PIPE

UPSTREM OF THE BREAK

BROKENEND

UPSTREAMEND

DOWNSTREAM END

BROKEN END

SECTIONOF PIPE DOWNSTREAIOF THE BREAK

A

PRESSURE TRANSDUCER POSITIONS DISTANCE OF PRESSURE TRANSDUCER FROM UPSTREAM END

Fig. 6.2 Schematicof Foothills Test NABTF7

Sufficient data was provided concerning the test gas, to enable determination of the rest of the required fluid properties using the QUANT software. The former are summarisedin Table 6.1. The average composition of the gas used in the seven test which is given in Table 6.2 is used as the input data to the QUANT software. A data file FLUID. DAT

is generated using the routine DFGEN, to cover the whole range of the

dependent variables(p, T and p) which will be covered by the transient event.

Results produced by the computer model are summarisedin Figs.6.7 to 6.11. togetherwith the experimentalresultspresentedby FoothillsPipelines(Yukon) Ltd. (1981). A grid spacingof Ox = 0.1m and At=0.0001s, at the broken end is used in the numerical model and a variable grid spacingis used. 230

0

TSM 0 0

H FS 0 0

C3Hg i-C,H, o n-C4HIO

0 1 2

0 0 0

4.40 0.20 0.30

i-C5H12

1

0

0.03

n-CSH12 n-C6H14 CO2 N2 Ar

2 1 0 0 0

0 0 0 0 0

0.03 0.01 0.07 1.70 0.08

Chemical Formula CH4 CA

85.10 8.00

Table 6.2 Average GasCompositionfor Foothills Test

6.2.2 BRITISH GAS TEST DATA The tests referred to as British Gas Tests are those which were reported by Jonesand Gough (1981). They include a seriesof experimentswhich were carried out using short disc length bursting (36.6m) diameter, 120ft 4in (0.1016m) at and a sections of with pipe by Gas. developed The British to the tests theoretical aim of was validatea oneend. model The computer model, DECAY, could be usedto predict the decompressionbehaviourof any gas mixture from any starting conditionsof pressureand temperature. Threegasmixtures,namelymethanelethane, and a typical rich gas methane/propane mixture, were used. The decompressionwas performedfrom pressuresrangingfrom 70 to 125 bar (7. OMPato 12.5MPa) and temperaturesranging from 13 to 30°C. Gas of

knowncomposition waspumpedintothepipeuntil the requiredconditions(p andT) were attained. The bursting disc was explosively ruptured and the resulting decompression behaviourwas recordedby pressuretransducersat a numberof locations along the pipe. p and T curves together with the appropriate details on the gas composition, initial conditionsof p and T were presented.The curvesinclude both the experimentalresults and predictionsby the British GasDECAY model. Also presentedin the report, are four setsof experimentalp-a curveswhich were obtained from BMI using a 20 ft (6.1m) long 4in (0.1016m) diameter pipe; full-scale fracture tests carried out by British Gas for Shell and North America sponsors;and fullscalefracture tests carried out by BMI for Artic GasPipelinecompanies..For thesetests, only the p-a datais presentedtogetherwith predictedresultsusing the British Gas DECAY 231

(ii) Position 14.17m From the Break

(i) P osition 4.21m From the Break

8000--

Experiment MOC-1

---

Experiment M6

---

7000 6000to

8000 7000 6000-

ä H

it

5000 4000

5000

d

ooo 3000 2000

3000

a

2000 1000

1000

0

0 500

0

..

0

1000

1000

500

Time [ms]

Time [ms]

(iii) Position 25.02m From the

(iv) Position 36.06m From the

Break

Break --

8000

Y

Experiment - MOC-1

7000

6000

,...

6000

5000--

ic 4000--

4000-ä

---

8000

7000

N

3000

3000

2000

2000

1000

1000

0

0 0

500

0

1000

(v) Position 46.57m From the Break

(vi) Position 101.44m From the Break

Experiment

8000

---

7000

MOC-1

ä

7000

6000

6000

5000

5000

4000

4000

3000

3000

2000

2000

1000

1000

0

0 500

Experiment

8000 v

0

1000

500

Time [ms]

Time [ms]

ä

Experiment MOC-1

A

0

1000

---

500

MOC-1

1000 Time [ms]

Time [ms]

Fig. 6.3 p-t Curves for Foothills Test NABTFI East Comparison With First-order Method of Characteristics Predictions

232

(i) Position 4.21m From the Break aua.8000 .

Experiment MCC-1

---

6000 I

y 'A

ý

ä

4000 3000

---

5000-

ý 3000

ýý..

2000

Experiment MCC-1

8000 arc 7000 , 6000

7000

a

(ii) Position 14.17m From the Break

\

2000 1000

1000

0-

0

1

0

1000

500

0

1000

500

Time [ms]

Time [ms]

(iii) Position 25.02m From

(iv) Position 36.06m From the

the Break

Break

8000 .

0- 7000

.

Experiment MCC-1

ý 5 y

U) y

--

1

6000

-MCC-1

5000

5000--

a

Experiment

^m J 7000

4000-3000

3000

2000

2000

1000

11000--

0

0 0

1000

500

0

500

Time [ms]

1000 Time [ms]

(v) Position 46.57m From the

(vi) Position 101.44m From the

Break

Break

8000 I

--

7000

8000

Experiment -MCC-1



6000 N

Experiment MCC-1

6000

5000

y

5000--

0. 4000

4000

3000

3000

2000

2000

1000

1000

0-

0

i

0

500

1000 Time [ms]

0

500

1000 Time imsi

Fig. 6.4 p-t Curves for Foothills Test NABTFI East ComparisonWith MacCormackMethod Predictions 233

-4

- Experiment

(dX=0.1m) -f-MCC-1 -, º-MOC-1 (dX=0.5m)

N

500 a N

200 100

0

2000

1000

0

3000

4000

5000

7000

6000

8000

Pressure JkPal

Fig. 6.5 Pressure Wave Propagatio n Speed for Foothills Test NABTFI East (ii) Position PT3W

(i) Position PT2W 7 Y

9000

... d

8000

-

Experiment MOC-1

---

9000--

d

7000 1

CL

---

8000 1

7000

N

I

ä

6000 -

Experiment MOG1

0000 5000

5000 4000

4000

\

3000

3000

:: ::

2000

2000

1000

1000

0

00

200

400

600

0

200

400

Time [ms]

Time (ms]

(iv) Position PT5W1

(iii) Position PT4W ä

Experimen MOG1

---

8000

Y 8000

7000

7000

6000

6000

5000

5000

4000

4000

3000

3000

2000

2000

1000

1000

0

0 0

200

400

i 600

---

0

600

Eýe riment MOG1

200

400

600

Time [ms]

Time [ms]

Fig. 6.6 p-t Curves for Foothills Test NABTF7 West Comparison With First-order Method of Characteristics Predictions

234

(vi) Position

(v) Position PT5W2 9000-8000

9000--

7000

7000 2 a

6000--

a

Experi Cl; -M -ent

ä--0 " 8000

Experiment MOC-1

---

PT6W

6000--

5000

5000

4000

4000

3000

3000

2000

2000

1000

1000

0

0 0

200

400

600

0

200

400

Time [ms]

Time [ms]

9000ýö 0. 8000 7000

ä 5000

---

ent MOC-1

4000 3000 2000 1000

0 0

100

200

300

400

500

Timeims]

(vii) Position PTRW

Fig. 6.6 p-t Curves for Foothills Test NABTF7 West Comparison

With First-order Method of Characteri stics Predictions

4

--f- Experiment --E-MOC-1 (dX=0.5m)

350 300

250 ö, 200 N

150

3

100 50

0 0

2000

4000

600

6000

8000

10000

Pressure JkPaj

Fig. 6.7 Pressure Wave Propagation Speed for Foothills Test NABTF7 West

235

developed by The the two other models models. were computer existing model and other Exxon Production Research Company and the University of Calgary, and were based on in developing Gas British the the to theoretical model. one used approach a similar

No furtherdetailswere givenon the full-scaletests,apartfrom the gas compositions be full-scale initial T. The therefore the used to results cannot conditions of p and and by Jones The the this presented results experimental pattern of computermodel. validate by Foothills Pipelines is in Gough (1981) than that consistent presented general more and (Yukon) Ltd (1981). The problem of p-t curves for different positions on the pipe in data in is Section 6.2.1, the so present which was reported not crossing each other, describedin this section. The threetestsperformedby British Gasusing a multi-constituent (typical)naturalgasmixture areusedto validatethis computermodel. The initial pressures for the testwere 70bar,100barand 125bar,and the initial temperatureswere 29°C 25°C , and 25°C respectively. The three tests are referred to as British Gas Test 1,2 and 3 respectively.Theseresultsare usedto validate the computer model only in relation to the p-t curves. As for the p-a curves,only the BMI (PRUDHOE BAY 1), BMI (PRUDHOE BAY 2) andthe Universityof Calgary testsare used.These tests are referred to as BMI 1, BM12 andUCT respectively. There are two main reasonsfor this selectionof data, apart from that of availability of all the necessarydata to enablecomputer modelling. The two fact firstly, the that the multi-constituentgasmixture representsmore closely reasonsare: a typical natural gas mixture. The secondreasonis one which is of convenience. The procedureusedto preparethe fluid property data file using the QUANT software, (refer to Sections4.5 and 6.1.1) hasthe major deficiencyin that manualoperation is necessary for each set of input data. Consequently,this procedure is very tedious and time consuming.Test dataon gaswith compositionswhich are closest to that of the Foothills tests(referto Section6.2.1), is selectedso that the samefluid property data file could be used. The compositiondatausedfor the gasmixtures in all the tests selectedis presented in Table6.2 andresultsproducedby the computermodel together with their corresponding experimentaldata, are presentedin Figs. 6.9 to 6.21. The samegrid spacingat the break is used as for- the Foothills tests, including a variable grid spacing. A grid spacing of Ax=0.1m and it ß. 0001s is used for all the tests,at the broken end. The Schematicof the British Gas tests is presentedin Fig. 6.8.

236

36.6m

_a

EE 9E8F

101.6mm

mmm

DISC BURSTING

CAS FILL UNE

A

PRESSURE TRANSDUCER POSITIONS DISTANCE OF PRESSURETRANSDUCER FROM UPSTREAM END

Fig. 6.8 Schematicof British GasTests BGT I, BGT2 and BGT3

6.2.3 SNGSO TEST DATA

The experimentaldata used to validate this model so far, were carried out using rather data full-scale The was significanceof using pipeline experimental short pipe sections. discussedin Chapter5. In this section,and also the one which follows full-scale pipeline in experimentsareused order to obtain a better validatory evidencefor the model. In this sectionthe test resultswhich were reported by Sens,Jouve, and Pelletier (1970) are used The test wascarriedout on a 11,800mlong 0.1065m internal diameternatural gas pipeline,at a pressureof 3.14MPa, on the S.N. G.S.O gas system. A sectionof a pipe was duplicatedand takenout of service temporarily. One end of the pipe could be connected to the systemor isolated,while the other end was equippedwith a venting assemblywith a cylindrical passageof 0.1m diameter.The assemblypermitted venting of the pipe to the atmosphereby eitherburstinga discor by rapid operation of a stop-cock. It was observed that the two methodsof ventingthe pipewere indistinguishablefrom one another when the decompressionbehaviour at a point 6.159m from the broken end were considered. However,substantialdifferenceswere observedat the broken end during the first few milli secondsinstants after the break occurred. No details were provided regardingthe gas composition, but it is assumedthat the gas used was a typical natural gas mixture of compositioncloseto that of the Foothillstests(referto Table 6.2). Consequently,the same data file for the fluid propertiesasusedfor the Foothills test is used. Comparativemodelpredictionandexperimentalresultsare presentedin Figs. 6.23, 6.24 and 6.25. A much bigger grid spacingat the break (ten times that usedfor the 237

Experiment PT 1 80

---

Experiment PT 2

""---"

Experiment PT 3

MOC-1(dx=0.1m) PT 1 - MOC-1(dx=0.lm) PT 2 MOC-1(dx=0.1m) PT 3

----

-(u 1

'-`

""`

`" "`

60 a

50

.1".

40 30 20

10 0iiii 05

10

---I 15

20

IIi 25

30

35

Time [ms]

Fig. 6.9 p4 Curves for British Gas Test BGTI Comparison with First-order Method of Characteristics Predictions

ExperimentPT 1 ExperimentPT 2 ExperimentPT 3 MOC-2(dx=0.1m)PT 1

--""-"""

80

-"-" 70

ti

MOC-2 (dx=0.1 m) P12 MOC-2 (dx=0.1 m) PT 3

----

60

'ý CL



" ".

5040 30 20

10 0ItIIIIF, 05

10

15

20

25

30

35

Time [ms]

Fig. 6.10 p-t Curves for British Gas Test BGTI Comparison with Second-order Method of Characteristics Predictions

238

Experiment PT 1

80

Experiment PT 2

---

" Experiment PT 3

I-

(V

,

_-........

ý.

a+

----

mýý.

----

MCC-1 (dx=0.1 m) PT 2 MCC-1 (dx=0. lm)PT3

60-50

" I fax=v.

nný ri

". ,\\

-

ti

4^

1 ý-

30 20 10 0iIIIII

10

05

15

20

25

30

35

Time (ms]

Fig. 6.11 p-t Curves for British Gas Test BGTI Comparison with MacCormack Method (Alternative 1) Predictions

ExperimentPT 1 ExperimentPT 2 --ExperimentPT 3 "--"-" (dx=0.lm)PT 1 -"-"MCC"1 (dx=0.1m)PT2 ----MCC-1 -- - MCC-1(dx=0.1m) PT 3

120

T A ý.

100

"-"--ý

ý1

a

1'

60 '

4 20

0i 0

F1 10

20

ii 30

40

Time (ms]

Fig. 6.12 p-t Curves for British Gas Test BGT2 Comparison with MacCormack Method (Alternative 1) Predictions

239

Experiment PT 1 ------"

120

Experiment PT 2 Experiment PT 3

lm)PT1 ----MCC-2(dx=0. MCC-2 (dx=0.1 m) PT 2 ----

MCC-2(dx=0.1m) PT 3

100N

°:

'""

80.1 ý

60 40

201

oa

2iº

20

10

0

40

30 Time [ms)

Fig. 6.13 p-t Curves for British Gas Test BGT2 Comparison with MacCormack Method (Alternative 2) Predictions

I ä

120

ExperimentPT I ExperimentPT 2 --ExperimentPT 3 -----" 1m)PT1 ---"MCC-3(dx=O.

loo

(dx=O.tm)PT2 -"--MCC-3 MCC-3(dx=0.t m) PT 3

`

80

_. ...

,

1 "1

60 It 40

20

0i 0

ri 10

+ 30

20

40 Time (ms)

Fig. 6.14 p-t Curves for British Gas Test BGT2 Comparison with MacCormack Method (Alternative 3) Predictions

240

Experiment PT 1 ---

Experiment PT 2

140

---"-"

Experiment PT 3

120

MCC-3(dX=O.lm) PT I -------MCC-3(dX=0.1m)PT2

------ý'.

MCC-3 (dX=0.1 m) PT 3 looP.

801 60-- 1 40 20

0

10

05

15

20

25

30

35

Time (ms]

Fig. 6.15 p-t Curves for British Gas Test BGT3 Comparison with MacCormack Method (Alternative 3) Predictions

--"""""

ExperimentPT 1 ExperimentPT 2 ExperimentPT 3

tm)PT1 -"-"WKL(dX=O. WKL (dX=O.tm) PT 2 -""3 PT WKL (dX=0.1 m) -

140

120

"ý 11 t 11

`

ýs

20 0i

1i 05

10

II1I1 15

20

25

30

35

Time Imsj

Fig. 6.16 p-t Curves for British Gas Test BGT3 Comparison with Warming-Kutler-Lomax Method Predictions

241

Experiment PT 1

140

Experiment PT 2

---

n rimentPT3

----"F

ä1""""---"

MOC-1(dX=0.01m) PT I "'ý.

100--l. 80

----MOC-1 -

(dX=0.01m)PT2

MOC-1 (dX=0.01 m) PT 3



t\

---'

40

"_ ý-

--

-

._....

"-

.... -

20

0

15

10

05

20

25

30 35 Time [ms]

Fig. 6.17 pA Curves for British Gas Test BGT3 Comparison With First-order Method of Characteristics Predictions (dX=O.Olm)

ExperimentPT I ExperimentPT 2

---"""" 140

----MOC-2 ""

120

''

"'

a

100



80

(dX=0.01m) PT I

MOC-2 (dX=0.01 m) PT 2

.ý-"--

11; Ä

ExperimentPT 3

-..,

--

--

MOC-2 (dX=0.01 m) PT 3

'.

\

Ni

v

ýý_"tiý 60

20 0i1Ii11I

05

10

15

20

25

35 30 Time Imsj

Fig. 6.18 p-t Curves for British Gas Test BGT3 Comparison with Second-order Method of Characteristics Prediction (dX=0.01m)

242

Experiment PT 1

140-ý...

120-

..

Experiment PT 2

---

_

----"

H 100-

'ý.

ExperimentPT 3

(dX=0.lm)PT I ----MOC-2 MOC-2(dX=0.1m) PT 2 --"-

:1

MOC-2 (dX=0.1 m) PT 3

20 0 10

05

15

20

25

30

35

Time [ms]

Fig. 6.19 p-t Curves for British Gas Test BGT3 Comparison with Second-order Method of Characteristics Predictions (dX=0. lm)

I -+ -Experiment BMI-1 1-f-Experiment BMI-2

500

--A-MOC-2

400 N

(dX=0.05m)

300

200 100

0 0

20

40

60

80

100

120

Pressure [bar]

Fig. 6.20 Pressure Wave Propagation Speed Curves for BMI Tests ExperimentUCT -"--f- MOC-2(dX=0.1m) -A-MOC-2 (dX=0.01m)

400 350 Ný

250 200 150

goo so 0

0

0.2

0.4

0.6

0.8 Pressure Ratio

Fig. 6.21 Pressure Wave Propagation Speed Curves for University of Calgary Test

243

Foothills and British Gas Test) is used. The reasonfor using such a big spacing is to minimise CPU time. The analysisin this caseis performedfor a much longer time (25s) in British Foothills hundreds tens the the gas and of milliseconds used and comparedwith testsrespectively.No problemswere experiencedwith the numerical algorithm due to the A broken large (Ax=1m, At=0.0001s). the this variable selectionof end grid spacingat grid spacingis used and the largest ix=512m.

IA

t

100mm

BURSTINGDISC

GAS FlU. UNE

PRESSURETRANSDUCER POSITIONS DISTANCE OF PRESSURE TRANSDUCER FROM THE BURSTING DISC

Fig. 6.22 Figure 6.22 Schematicof SNGSO Test

6.2.4 API TEST DATA The second set of full-scale experimental data is that which was reported by the Alberta Petroleum Industry, Government Environmental Committee (1979).

The tests were

for the phase a field verification programme of isopleth prediction performed as second techniques. It involved a full-scale field programme to investigate, among other things, the behaviour of gas emissions resulting from a pipeline rupture. The rupture mode was the most critical variable examined. Two tests setups were used at the test site, which is located in the Western gas field in Southwestern Alberta. The first and second tests are referred to here as APIT1 and APIT2.

The tests were carried out using an existing

168.3mm outside diameter pipeline which was typical of sour gas pipelines in the province at pressuresof 6.9 and 3.45MPa respectively. The tests section was approximately 4.0km long. It was burst at the mid point. The third test is referred to in the report as APIT3. It was performed on a 323.9mm outside diameter and approximately 7.1km long pipeline. Also in this test the pipeline was pressurisedto 6.9MPa pressureand ruptured at the middle. The specification of the pipes used in the two tests was provided. As with SNGSO test, no

244

(ii) Position 2250m From the Break

(i) Position 1500m From the Break Experiment 16

16

12

12

12

10

10

10--

v ;

, v 3

6

k

4

2

2

2--

0ý 0

0

0-

20

30

16 14

---

4

U.

0

10

20

30

10

0

16

Experiment MOC-1

14 '

---

16

Experiment MOC-1

14

12

Experiment MOC-1

---

12--

0

3 üý.

10--

ö

8

6

d 3

6..

6

4

4

4

2

2

2

01 0

0

0

30

0

10

Time (s]

20

8

30



0

20

Time (s]

(vii) Position 6000m From the Break

30

(vi) Position 5250m From the Break

8

20

20

Time (s]

(v) Position 4500m From the Break

12

10

.

Time (s]

(iv) Position 3750m From the Break

LL

6

4

10

MOC-1

---

14

Time (s]

d

16

MOC-1

---

14

6--

ö

Experiment

14

v

° 2 -

Experiment

MOC-1

---

(iii) Position 3000m From the Break

40 Time (s]

(viii) Position 6750m From the Break

(ix) Position 7500m From the Break N

16 P

14

---

16

Experiment MOC-1

14

---

ExperimentMOC-1

16-8

14

12

' 12

"Z 12

3

10

10

10

M

8

3

0 W

6

---

U.

8 6

6

4

4

4--

2

2

2--

0

0

0

0

20

40

Experiment MOC-1

0

20

40

Time (s]

Time (s]

t--; 0

20

40

Time (s]

Fig. 6.23 u-t Curves for SNGSO Test Comparison With First-order Method of Characteristics Predictions

245

32

ä

Experimen MOC-1

---

Experiment

32

31

31

30

30

a

MOC - 1

---

-

ä

28

28

27

27

27

26

26

26

25

25

25 0

40

20

40

0

(vi) Position 5250m From the Break

32 Ä a

30

32

31

a

31

30

E

30

d

29

a Experiment MOC-1

28 ---

d

29 28

Experiment MOC-1

---

a

Experiment 28

27

27

27

26

26

26

25

25

25

0

20

40

0

20

40

0

---

MOC-1

10

20

Time [s]

Time [s]

(ix) Position 7500m From the Break

32

32

31

31

31

30

30--

30--

v

N 04 N

28

---

Experimen

Experiment

MOC-1

MOC-1

28

---

a

27

26

26

26

25

25

25

40

0

20

Time [s]

40 Time [s]

---

28

27

20

Experiment

R 29

27

0

30 Time [s]

(viii) Position 6750m From the Break

(vii) Position 6000m From the Break 32

30

20

Time [s]

(v) Position 4500m From the Break

32 31

10

Time [s]

(iv) Position 3750m From the Break

Vgl

, 31

28

20

MOC-1

---

30

Time [s]

a

ý

29

0

Experiment

32

29

Q. 29

V ä d

(iii) Position 3000m From the Break

(ii) Position 2250m From the Break

(i) Position 1500m From the Break

0

10

MOC-1

20

30

Time [s]

Fig. 6.24 p-t Curves for SNGSO Test Comparison With First-order Method of Characteristics Predictions

246

(ii) Position 2250m From the Break

(i) Position 1500m From the Break

Experiment MOC-1

Experiment ---

MOC-1

--16 14 12 w 10 8 6 tLr__4 d >4 2 0 0 10 20

16 c 14 90 12 d 10 x8

>2

6 4 0

30

0

10

(iii) Position 3000m From the Break Experiment ---ö 16 14 v 12 x 10 8 g v 4 2 0 0

20

30 Time is]

Time [s]

MOC-1

+-+--a 10 20 30 Time

(iv) Position 3750m

(v) Position 4500m From

(vi) Position 5250m

From the Break

the Break

From the Break

16 g

14

V

12 --

41 X

---

Experiment MOC-1

ä

6

MOGt

---

L

>

14

ä

---

12--

10

10

MOC-1

d x 8

6

6 >

L6

4

Experiment

12

8

8

16

Experiment

14

x

10

16

4 L4 2

0 0

-ýý 10 20 30

0

10

20

(vii) Position 6000mFrom the Break

Time is]

(i x) Position 7500m From the Break

(viii) Position 6750m From the Break

16 14

Expenme E "

12 X

0

-1 10 20 30

Time (s]

Time is]

d C, =

0

1 30

---

MOGt

10

"W ä 16 14

Experiment

=

MOG1

.

8

>

12

---

14

>

4 2 0

10

,

MOC-1

---

10

8 >

6

6

4

4

2

2

0

0 0

0

Expeimerýt

'too 12 v

8 6-

16

ä

10 20 30 Time (s]

10

20

30

Time (s]

0

10 20

30

Time [s]

Fig. 6.25 Velocity Head for SNGSO Test ComparisonWith First-order Method of Characteristics Predictions 247

data was provided for the gas composition. The same data which was used in the previous sections is used for the fluid properties. Due to the long run times required, a coarse grid spacing is used in order to reduce the CPU time. A grid spacing of Ox=5m and At=0.005s is used for APIT1 and APIT2 tests, and a grid spacing Ax=1Om and Ot=0.01s is used for APIT3 test. A variable grid spacing is used in all the computer predictions.

It was observed [Alberta Petroleum Industry, Government Environmental Committee(1979)] that the rupturemodewasunpredictableand difficult to control. It was found that under similar test conditions any of three different rupture modescould occur. The threerupturemodesare bell opening,ring-off and a combinationof the two in which fracture is bell bell In is propagateda short the one end and other ring-off. rupture, the distancealongthe longitudinalaxis with the pipe remainingintact. In ring-off rupture, the fracture is propagated a short distance along the longitudinal axis and then turned circumferentialaroundthe pipe.The rupture sectionwould either be blown completelyout of the line or remainattachedto the pipe by a small tab at either end. A pressuresensorwas placedwithin 1mof the rupturepoint in order to ensurethat *obtained were in the critical flow zone. It was observed that the the measurements instrumentswere subjectedto extremeforces by the exiting gas jet and often they were blown off the line if, the rupture mode was a ring-off. When the rupture mode was a bell opening,the instrumentswere unaffected. A total of elevenexperimentswere conducted for measuringflow rate data, six of which could not provide completedata becausethe instrumentto measurethe critical flow were lost duringthe blowdown. The remainingfive testsproducedsufficientdatato validate the computer models. Data from three tests (two bottom ring-off andonesidebell) were combinedto form a set of experimentaldata for the 168.3mmdiameter pipe at an initial pressureof 6.9MPa (APIT1). Another test, in which the rupture was top ring-off was used for the samepipeline but pressurisedto 3.45MPa (APIT2). The lastof the five good tests was carried on with the 323.9mmdiameterpipe. The line was pressurised to 6.9MPa and the rupture mode was bottom bell. This provided the datafor test APIT3. The datafor tests APIT 1, APIT2 and APIT3, which is described in this sectionandthe correspondingpredictionfrom the computer model developedin this study, are presentedin Figs 6.28 to 6.33.

248

2000m

- ----------

--

------

----------

- --- - ---------------

- --- - ---------------

-----

- --

----

$ 168.3mm BREAK

us FlLL UNE

A

PRESSURE TRANSDUCERPOSITIONS

Fig. 6.26

Schematicof Alberta PetroleumIndustry Tests APITI and APIT2

3550m

S 323.9mm

4

BREAK

CASnu urge

POSITIONS PRESSURETRANSDUCER

Fig. 6.27

Schematicof Alberta PetroleumIndustry Tests APIT3

6.3 DISCUSSION OF VALIDATION GENERAL

RESULTS

DISCUSSION

FannelopandRyming(1982) definedthe different time regimesfollowing a linebreak,each in Knox, described Also Section different 3.3.2. These are requiringa methodof solution. Atwell, Willoughby and Dielwart (1980) discussedthe different rupture modes, which both late different In the time regimes, this early and study, approaches. require modelling data defined by (1982), Ryming Fannelop the various are modelled using and which were both for in The theoretical this numerical presented and modelsare used chapter. same time regimes.Unlessthe pipe diameterand/or cross-sectionareais specifieddifferently at 249

200-180-API Model

160-140

---

Exprimental

......

MOC-2 (dX=5m)

120

LL 100 80 60--

40-20--

00 50

200

150

100

Time (s]

Fig. 6.28 Mass Flow Rate for Alberta Petroleum Industry Test APIT1 Comparison with Second-order Method of Characteristics Prediction

7000 ý'" 1,, `" \

6000 a, v

d

I---

"'".

......

Experimental API Model MOC-2 (dX=5M)

4000--

3000-\\ 2000-1000 1_1 0 0

50

100

150

200

250

Time [s]

Fig. 6.29 Pressure at the Intact End for Alberta Petroleum Industry Test APITI Comparison with Second-order Method of Characteristics Prediction t

250

100

90 o, 80-

--0

API Model

--W

Exprimental MOC-1 (dX=10m)

A

D

70

ý-

60 50. 401. 30 20 10

50

0

100

250

200

150

Time (s]

Fig. 6.30 Mass Flow Rate for Alberta Petroleum Industry Test APIT2 Comparison with First-order Method of Characteristics Prediction

3500 Experimental

-0 -f-API

3000 --

Model

--A-MOC-1

(dX=10m)

2500 N N

a 2000 d rn

1500 1000 500

0

-------ý 0

50

100

150

250

200

Time (s]

Fig. 6.31 Pressure at the Intact End for Alberta Petroleum Industry Test APIT2 Comparison with First-order Method of Characteristics Predictions

251

700- -

0

API Model

--01-- Exprimental --Ar-MOC-1 (dX=10m)

600 w

500 U.

U, 400 G

300

200

100

0 0

50

100

150

200

_I

----

250

300

350

400

450

500

Time (s]

Fig. 6.32 Mass Flow Rate for Alberta Petroleum Industry Test APIT3 Comparison with First-order Method of Characteristics Prediction

Experimental -0 Model -ý-API --Ar-MOC-1 (dX=10m)

m^ 8000 --19 7000 6000

-5000 -a. 4000 -3000 -2000 -1000 -0 0

50

100

150

200

250

300

350

400

Fig. 33 Pressure at the Intact End for Alberta Petroleum Industry ielsl Is] APIT3 Comparison with First-order Method of Characteristics Prediction

252

450

brokenend,full bore rupturemodeis assumed.The modelis limited to analysisof the flow following a linebreakin a straight section of pipeline. Various aspectsof the model and basedon resultsobtained from the computer modelling the predictedresultsarediscussed, which is describedin sections6.2.1 to 6.2.4 and previoustrial runs. The symbolsused in the graphsare definedas follows: MOC-1 MOC-2 -

First-order method of characteristics

MCC-1 MCC-2 -

MacCormackmethod using alternative 1

MCC-3 WKL -

MacCormackmethod using alternative3

Second-ordermethod of characteristics

MacCormackmethod using alternative2

Warming-Kutler-Lomax method

Two setsof experimentaldatawhich were usedby Tiley (1989), havebeenused in this study. These are the British Gas test data and [Jones and Gough (1981)] and the Foothillstest data[FoothillsPipeline(Yukon) Ltd. (1981)]. A better comparisonhasbeen in achieved this study with the two setsof data than in the study by Tiley (1989). In the study by They (1989), the length of the time step, grid size and break conditions were for in simulation, varied each order to optimize convergencetowards a stablesolution. In this study, the ratio betweenthe time step and the distancegrid is selectedsuchthat the Courant-Friedrichs-Levy stabilitycriterionis satisfied.In this uniform grid spacing,the time step varies proportionally as the distancegrid varies. For the test data used in this study, the maximum value of speedof sound and hencethe maximumflow velocity was above 400m/s,but neverit neverreached500m/s. Therefore, in order to ensurethat the stability criterionis satisfied,a constantratio of ix/At=1000m/s was used. This ratio could handle situationsin which the flow velocity and wave speedare each500m/s,without failing the stability criterion. In every calculation step the programme checks to ensurethat the Courant-Friedrichs-Levy stability criterion is satisfied. If the stability criterion is not satisfied,the programmestopsautomatically. Under no circumstancedid the programme stop during the simulationsreported in this study due to failure of the stability criterion. In most of the simulations,the grid sizesusedat the break boundarywere coarser thanthosewhich were usedby Tiley. Theyuseda grid size of 0.0254m for the British Gas testsbut shecould not achievestableresultsfor someof the tests. Tests for which stable resultscould not be achievedarethosefor which the initial pressurewas high and the initial flow ratewaszero. No such problemswere encounteredin this study, even with grid size 253

of up to 0. Im for the British Gastests. The only problem encounteredfor such testswas that the executionspeedof the programmewas very slow and the boundaryconditions, in in is further This the condition explained some cases, were not properly modelled. discussionof the method of characteristics. Tiley (1989) also noticed that the theoretical p-t curves tended to begin their pressure drop too early. This effect was more noticeablethe further the transducerwas from the break. She attributed this phenomenonto the responsetime of pressure transducers,which causesa delay in recording the pressuredrop. In this study, only the transducerswhich are close to the break havebeeninvestigated. The samephenomenon by as was observed They was repeated. However, a comparisonof the pressurewave speedspredictedby the computermodel,with those obtainedfrom the BMI and University of Calgarytestsshow a good agreement. Sincethere were no p-t curvesprovided for the latter two test, it was not possibleto comparethe pattern of the p-t curves: Also since there are no pressurewave curves for the British Gas test which have beenused in this it study, was not possibleto comparethe pressurewave speed. However, the fact that the by favourably from the predicted model compare pressurewaves with experimentalvalues the BMI, Universityof CalgaryandFoothillstests indicate that there is a problem with the British Gas test data. They (1989) did contact the British Gas in order to ascertainthe accuracyof their data. Shereportedthat they could not ascertainthe accuracyof their data, but they indicatedthat they were not entirely satisfiedwith it. Although not conclusively confirmed, they believedthat one of their transducers,PT2, was malfunctioning.

in the final pressurereachedafterthe break. In some Theynoticedinconsistency casesthe pressurewas higher and in somecaseslower. She attributed this condition to inaccuraciesincurred in the calculation of equalizationpressureat the break. No such problems were encountered in this study due to the accurate model developed for calculatingequalizationpressure.It shouldbe noted that the data for the curve PTI in the British gastest BGTI arethoseobtainedfrom predictionsusing the British Gas theoretical model and not from experiments. In simulating the Foothills test, They (1989) useda grid spacingof 0.01m at the break. This spacingis much finer than that which is usedin this study (0.1m to 0.5m, and it would be expectedto produce better results. However, They usedtime stepsvarying from 0.65 to 0.75ms (quoted as s in her thesis). This time step is too big for the value of Ax usedandfails the Courant-Friedrichs-Levystability criterion. The ratio Ax/At usedby 254

Tiley wasbetween13.33and 15.39m/s,which is much lessthan 1000m/swhich is usedin this study. This gross error in the time step is probably the main reasonfor the instability by by decreasing Tiley. The Tiley the time step could that argument problemsencountered havelead to an increasein the accumulativeround-off error in the resultsis not valid. Simulation results produced in this study using the method of characteristics data, the as the p-t curves and the pressurewave experimental compare very well with speeds are modelled correctly. They argued that a possible reason for her model's overestimation of the wave speedswas the second-orderapproximation used near the break. This argument is not correct becausethe second order approximation has an oppositeeffecti.e. producinglower wave speedsthan those producedusing the first order approximation.As wasthe casewith Tiley's simulations,the final pressurescalculatedfor the Foothill tests were slightly higher than the experimentalvalues. The reasonfor this is that the theoreticalmodelsdid not account for the crack propagationalong the length of the pipe. The crack propagation has the effect of moving the point at which the equalizationwould occuralongthe pipe. In this study,this effect was more severewith test NABTF7 thanwith testNABTFI becauseof the long crack propagationin the former test (18.28min NABTF7 comparedwith 4.21m in NABTF 1). The crossingof the p-t curves, which was attributed to temperaturerelated zero drift during testing [Foothills Pipeline (Yukon) Ltd. (1981)] does not exist in prediction obtained using the method of characteristics andthe Warming-Kutter-Lomaxmethods. However, this phenomenonwas observed when the MacCormack method, indicating that the method is not suitable for analysisof transientflow following a break in a high-pressuregas pipelines. The two other setsof data,namelythe SNGSOand the Alberta PetroleumIndustry data, were not usedby Tiley. Both setsinvolved relatively long pipesand the data seem to be consistent and satisfactory. A good agreement was obtained between the experimentaldata and the prediction results,evenwith the big grid size used. The major weaknessof the data is that they do not contain sufficient information about the gas used, some specificationsof the testing systemand accuracyof measurementsrecorded. The intact end pressurewhich was calculatedusing the Alberta PetroleumIndustry data differ from the experimentaland their model prediction data considerably. The second-order methodresultsshowthe biggestdiscrepancy.Also the massflow ratescalculatedusing the secondorder-methodarehigher. The valuesproducedwith the first-order method compare muchbetterwith the experimentaldatathan thoseobtainedusing the second-ordermethod.

255

The valuesmassflow ratesproducedby the first-order method of characteristicsare lower thanthe experimentalvalues. The reasonsfor this discrepancyare errors in calculatingthe° gas density and the big grid spacingused. The results obtainedfrom this simulationgive an indicationthat the first-order method of characteristicsis the most suitablein situations big is where a grid spacing used. A comparisonof the model prediction and experimentaldata for the SNGSO test show a reasonablygood agreement. The pressurestarts to drop at the sametime (Figs. 6.23 and6.24), which indicate that the wave speedspredicted by the computer model are correct. However,the magnitudesof the pressuredrop and flow velocity predictedby the computer model are less than experimental values. This discrepancy is caused by interpolationerror andthe big grid spacingused. Data which was obtainedusing a smaller grid spacing(Ax=1m) and the second-ordermethod of characteristicscould not be used for validationbecausethe interpolation error and round-off error of the computeraffected the resultsvery much. This effect is explainedin detail in the grid size discussion.

THE QUANT SOFTWARE PROPERTIES OF FLUIDS

FOR

THERMODYNAMIC

AND

TRANSPORT

The fluid propertydata,which are usedby the programmes,are calculatedusing the QUANT Software. Heattransferis calculatedusingthe recovery factor and adiabaticwall temperaturemethodand therefore the valuesof Pr, thermal conductivity of the fluid and CARare not required. However, the three valueslisted aboveare usedin other variations of the computermodel. Three main limitations of the QUANT software are unavailability of Pr, kf and t at somelow temperatures,unavailability of all output valuesat much lower temperaturesandtemperaturesbelow 200K and lack of all outputs at high pressures.The composition of the gas mixture used in this study is given in Table 6.2. The highest pressurefor which output is available,but only for part of the temperaturerangerequired, is around7MPa. This limitationhasbeenovercomeby using the fluid propertieswhich are availableat the closeststateof the fluid andadditionaldata for the gas mixture which was provided by the suppliersof the QUANT software on request [Silberring (1995)]. The additional data was calculated using a later version of QUANT which is still being developed. The later version of QUANT is expectedto be able to produce output at higher pressures,similar to those encounteredin high-pressure natural gas pipelines.

256

It is observedthat the temperatureof the gas near the break boundary drops to values below the 200K minimum limit of the QUANT. This situation occurs when the initial pressureof thegasin the pipelineis sufficientlyhigh and the initial temperatureof the gasis relativelylow. This meansthat in order for QUANT software to be able to cover all linebreakeventsin high-pressurenatural gas fully, its rangeof temperaturecoveragehas to be extendedto temperatureslower than the presentminimum limit of 200K. Silberring (1993-95) and Flatt (1993-96) arguedthat the temperaturescould never drop below the 200K limit. The routine FLDPROPV calculates the fluid properties from the data file FLUID. DAT, using an interpolation procedure. When used with the transient analysis programmesproblemswere encountereddue to the limited computer memory available. An alternative procedurein which the closest properties are used, produced satisfactory results. However the alternative procedure requires that a close spacing of the input parameterbe used. A spacingof 0.1MPa and 5K was usedfor pressureand temperature respectively.A comparisonof executiontimes when constantfluid propertieswere used and when the proceduredescribed abovewas used, revealedthat the biggest proportion is the time of execution spentin calculatingthe fluid propertiesfrom the data file generated by QUANT. As explainedin Section4.5, the procedurewhich is used in this study is the fastestoption available. INITIAL CONDITIONS BEFORE THE BREAK In all the analysesreported in this chapter,calculationof the initial conditions of the gas is greatly simplified by the fact that the initial flow velocity is zero. Any of the different steady state analysismodelspresentedin Section 4.2.2 would have produced the same results. However, the non-adiabaticnon-isothermalsteadystate analysismodel is used. Also since the flow velocity is zero, transient analysisbefore the break would producethe sameresultsasthoseproducedby steadystate analysis. An application of the transientanalysisbefore the break are presentedin section7.4. CONDITIONS AT THE BREAK BOUNDARY The methodof characteristics is usedfor solutionat the boundary nodes,in all the transient analysis programmes. In the caseofthe break boundary the initial break conditions are calculatedasdescribedin Section4.2.3. The initial conditions calculated,enabledresults 257

t=0 be Using time the this at pressure procedure, attained. stable which are numerically drops to its equalisationvalue which is calculated as describedin Section 4.2.3. The density Also falling the but keeps first of is time drop passes. as on temperature small at break the The decompression remains drops during at sound speed of the process. the gas break drop drops. the The drops but at pressure the sameat t=Os, as the temperature boundary, below the equalisationpressureis natural. There is no time dependentcurve broken The drop the by (1989) They end. to the at pressure model such as the one used time dependentcurve may lead to in incorrect results becausethe exact rate of pressure drop is not known. When using the MacCormacksecond-ordermethod, solution at the node next to by been has in This boundary the controlled velocity. the node produced an overshoot limiting the magnitudeof flow velocity from exceedingits correspondingspeedof sound.

NUMERICAL METHODS OF SOLUTION (i) The Method of Characteristics Tiley (1989)usedthe second-ordermethodbecausethe first-order method did not meetthe required accuracy and stability criteria.

Even with the second-order method, she

For instability certain grid and accuracyof results. encounteredproblemsof numerical became initial the solution unstableat random points along the conditions, points and by be She totally that the using an problem could alleviated recommended pipeline. bishnoi by Picard The and sameproblemswere encountered alternativenumericalmethod. (1989). Flatt (1989) encounteredsimilar problemswhich he called singularity. Problems suchasthoseencounteredby Flatt (1986) and Tiley (1989), do not exist with this model. A convergence tolerance of ± 5% and ±1% is used for the first and second-order calculations, respectively. The number of iterations required is very small. In most calculations,oneor two iterations for eachorder are sufficient. The fluid propertiesused for calculationof the coefficientsof the threesimultaneous equations,in the first-order step, are those at position M in Fig. 4.4. The use of the fluid properties at position M rather thanthoseat positionsQ, R andS was madein order to reducethe size of the programme. However,a significanterror could be introduceddue to this simplification, especiallyif the fluid propertiesvary considerablybetweentwo grid points, such as is the casein ruptured high-pressurenatural gas pipelines. Sincein this casethe first-order method is usedas a 258

first-order introduced by for the the the second-ordercalculation, error rough estimate first-order by is In the the the that second-ordermethod. case only calculation corrected fluid is for then the the propertiesat respective complete solution, values of method used positionsQ, R and S should be used in the first iteration, and their averageswith valuesat iterations. For iteration be in for the P the the should used subsequent previous position fluid the properties used to calculate the coefficients of the second-order method, Q, between R and the of position newly established simultaneousequationsare averages S and those previously calculatedat position P. It was observedin this study that in situationswere there are sharpchangesin the fluid properties, such as during the first few At's in the region around the break, the is fails first-order For the numerically. such cases, method used second-order step throughout the calculation. When the method of characteristicswas used to model the flow reversal in the sectionof the pipelinedownstreamof the break, it producedresultswhich tend to lean on the valuesat the intact end of the pipeline section. This directional bias results in a very in drop the broken section of the pipeline,including the broken boundary. slow pressure The problem of directional bias does not exist either in the programme based on the MacCormackmethod, nor the one basedon the method of characteristics for the section break. investigation into the problem did not reveal any An the the of pipelineupstreamof error in the calculation procedureor computer coding. In fact the valuescalculated at variousstagesarecomparablein magnitudewith thosecalculatedby the programmefor the pipeline section upstream of the break and the characteristic curves are positioned correctly. The problemof directionalbias with the method of characteristicswas observed only whenthe second-orderapproximationwas used. Both the upstreamand downstream modelsproduce comparableresults,when the first-order method is used. It was observedthat the second-ordermethod of characteristicsis more accurate than the first-order method, especiallyif the simulationinvolves the areain the vicinity of the breakandvery small grid size. The second-ordermethod producesoutputs which are beyondthe rangeof maximumand minimum output around the leadingedgeof the of the decompressionwave. This causessevereproblemswith the numericalalgorithm. The situation is controlled by using the first-order approximation. The computer programme includesa provision to check if the condition is encounteredand calls the programmefor the first-orderapproximationto correct the situation. It was also observedthat the use of 259

a smaller grid size helps to minimize the problem to some extent. Computation speed of the first-order method is over twice as fast as that of the second-order method. For some boundary conditions and initial conditions of the gas, the second-order method becomes slower.

For example, when modelling the British Gas test BGT3 (initial pressure and

temperature are 12.4MPa and 273K respectively and Ax=O.OIm), the second-order method failed first-order Also to that the the ten times method method. second-order slower was handle the choking boundary condition well. The first-order method was used in this case and produced satisfactory results. When a larger grid size (Ax=O. lm) was used, the problem with the second-order method was minimal and computational speed was the same as with the other Brirish Gas test data. The problem of modelling the choking condition was caused by the failure of the numerical procedure to model properly the temperature drop and therefore increase in density, which follows the initial decompression. The situation was corrected by tuning in the temperature drop and also increasing the convergence interval of the second-order method.

The major differencebetweenthe first- and second-ordermethodsis that the pressuredrops predictedby the second-ordermethodare smallerthan those predictedby the first-order method. Also the speedof propagation of the pressurewaves is slightly fasterwith the first-ordermethod. However, both resultscomparewell with experimental data. The main reasonfor the discrepancyin the first- and second-ordermethodsis the different convergencecriteria use (1% and 5% for the first- and second-ordermethods respectively).The two valueswere selected,firstly becausethey are thought to correspond with relative accuraciesof the two methods,and secondlybecausethe first-order method is alsousedasa roughapproximationof the second-ordermethod. Ideally, a more accurate criterion should be used if the first-order method is used on its own to calculatethe final 'solution. However, it is thought that a value which is smallerthan 5% would not too fine in relative to the overall accuracy of the first-order approximation. Unfortunately, the pressurewave propagationspeedsof the two methodscould not be comparedat further positions from the broken boundary and for a longer run time. A comparison of the predictedmassflow ratesand pressureat the intact end, using the first- and second-order methodsandthe AlbertaPetroleumIndustry data, showsa better agreementwith the firstorder method than with the second-ordermethod. This result was achievedeven though a finer grid spacing was used with the secondorder method. An explanation for the discrepancyis given under the generaldiscussion. 260

(ii) The MacCormack Second-order Method

Beaucheminand Marche (1992) concludedthat the MacCormack method is superior to the method of characteristicswhen Cn differs appreciablyfrom unity. When C. is much smaller than unity, the MacCormackmethod producesresultswith a precisionthat could numberof computationnodes,when the method of the not be attainedwith anyreasonable 1 2 is They that the the and which are alternatives concluded use of also characteristic used. described in Section 3.31, in successionon time steps could introduce significant oscillationsin the solution, especiallywhere the basicequationsare poorly approximated. Directional bias could be avoided by using exclusively one of the two calculation alternatives.

The directional bias is important only when working with two space

dimensions, in which case it

both that the methods average of was recommended

(alternative 3) be used. Beaucheminand Marche (1992) claimedthat doing so did seem to smearthe shock slightly and the computationtime was doubled. In this study the computer programmeis written in a way that any of the three is be is despite it fact This the that recommendedin Section3.3.1 alternativescould used. that alternativeI be used exclusively. Resultsfrom the three alternativesare presentedin Fig. 6.12,6.13 and6.14,for British GasBGT2. Due to limited computer memory,only the first-order calculation of the method of characteristicscould be used at the boundary points, when using alternative 3. The MacCormack method is extremely simple to programme, comparedwith the method of characteristics. The execution speedof the MacCormackmethodis fasterthanthat of the method of characteristics. But in this case, fluid from to the the time the QUANT software constitutes calculate properties used where the biggestproportionof the CPU time, the two methodshaveexecution speedswhich do not differ much. The sameargumentappliesfor the differencewhich is to be expected betweenalternative3 and the other two alternativesof the MacCormack method. It is statedin Section3.2.2.2that in the presenceof shocks,explicit finite-difference methodsof higherthanfirst order produceconsiderableovershoot and oscillatory systems. Resultsobtainedfrom the MacCormackmethod are oscillatory especiallynear the broken end. The oscillations are more severefor smaller L/D values. Also an overshoot was observed in the flow velocity, at the node next to the boundary node at the break. The overshoot was controlled by limiting the magnitudeof the flow velocity to that of the correspondingspeedof sound. Theseproblemswere not encounteredwith the method of

261

characteristics. With the MacCormack method the problem of directional bias which was encounteredwith the method of characteristics, in the section of the pipeline down stream the break, does not exist.

The threealternativesof using the MacCormack method were comparedusing the British Gas test BGT2 data. Results are presented in Figs. 6.12,6.13 and 6.14. Alternatives I and 3 producedsimilar results. Alternative 2 producedthe worst results, with much bigger oscillations and pressuresfalling fastest. The computation speedof alternative 1 is higher than that of alternative3. (iii) The Warming -Kutler - Lomax Method

The computer programme for the Warming-Kutler-Lomax method was written and but it could not run fully becauseof the limited computer memory. compiledsuccessfully, Even in the main frame computer, the memory of 20MB which was allocated is not enough. In order to simplifythe programme,the method of characteristicswas used both for the boundary node and the node next to the boundary. Bhallamudi and Chaudhry (1990) recommendedthat the MacCormack second-ordermethod be usedfor solution at the node next to the boundary node. Even this did not make it possible to run the programme on a pc. However, by using constant values for the thermodynamicand transport propertiesof the fluid, it was possibleto run the programmefor the third-order Warming-Kutler-Lomax methodon a PENTIUM P75 pc.

The methodfailedto producegood resultsin the regionof sharpchangesin fluid in Results properties. produced this region included negativevelocities (while all other velocitieswere positive)andpressureswhich are higher than the initial pressurebefore the break. The first-order method of characteristicswas used to smoothenthe sharpedgeof the decompression wave,andalsothe flow velocity andpressurewere limited to abovezero andwithin the initial pressurebeforethe break. The British Gas test BGT3, whose run time is 35ms,wasusedto test the Warming-Kutter-Lomaxmethod. The first-order method of characteristics was used during the first lOms, but this was not enoughto smoothenthe sharp edge of the pressureprofile so that the Warming-Kutler-Lomax method could be applied. Finallythe initial 20msof the run time werefound to be sufficient, for "smoothing" by the methodof characteristicsand the Warming-Kutler-Lomax method was used for the remaining 15ms. Under this condition, good and stable results were obtained. No oscillations or overshoot, such as those observedwith the MacCormack method were 262

Warming-Kutler-Lomax Results the method are compared with with produced present. MacCormack by the method the of characteristics, those produced second-order method 6.19. in 6.16 Figs. data. to The results are presented and experimental

(iv) Comparison of the Different Numerical Methods The transient analysismodelsbasedon the method of characteristics,the MacCormack based Warming-Kutler-Lomax on accuracy,stability compared the are method methodand MacCormack including for The the main reason of results and computationaleconomy. in two the Warming-Kutler-Lomax to methods this whether study confirm was methods and high-pressure in break following flow for gas the transient a modelling are suitable in during literature this study and summarized reviewwhich was conducted pipelines.The Chapter 3, indicated that these methods could be more suitable for high-pressuregas linebreak applicationsthan the method of characteristics. Tiley (1989) used the secondby be better but that using could obtained results order methodof characteristics, concluded They [Thorley In and publication solution. a previous an alternativenumericalmethod of (1987)], the MacCormack method was recommendedas the most suitable for linebreak problems. Based on the comparisonmadein this study, it is concludedthat the method of linebreak investigate, for best is the three the applications. methods of all characteristics The criteriausedin comparingthe threemodelis accuracyand stability of results,computer in CPU The Warming-Kutler-Lomax time third-order method, requirements. memoryand combinationwith the methodof characteristics,produced resultswhich are closeto those produced by the second-order method of characteristics. The computer memory be found big. latter MacCormack The to too the method was methodwas requirementof flow following linebreak in high-pressure for transient natural gas a unsuitable modelling in in low It the p-t results resulted oscillating pressure region, which pipelines. produced in low It the the pressure other. predicts wave speeds reasonably well curvescrossingeach it in the high pressureregion. The magnitudeof equalization region,but it underestimates pressure is slightly higher than that calculatedwith the method of characteristics. The computationspeedof the MacCormackmethodis oneanda half time faster than that of the first the or secondalternativesare used. second-order method of characteristics,when When the third alternativeis used,the computationspeedof the MacCormack method is the sameas that of the secondorder-methodof characteristics 263

GRID SIZE

The factor of two which is usedbetweenone grid spacingand the next, and the secondA Taylor's interpolation theorem, results. the produce good using order polynomial in long flow the problemwas encounteredwith the variablegrid method,when modelling in time The next the whose points node results some computer error of round-off pipelines. level for calculationhasnot beenreached,being includedin the calculations. The error occursbecauseof the big differencebetweenthe grid spacingat the break boundaryandnodeswhich are nearerto the intact end. Even when a value of Ax=1m is intact is long 512m for (Ax broken 11.8km the end), the at = end, a pipelinewhich usedat the error still affects the programme. Although this error does not seemto affect the flow it has the velocity on a significant negative effect pressure predictions much, boundaries internal In somecasesnegativevelocities are obtainedat some predictions. betweendifferent grid sizes,while the rest of the velocities are positive. The net effect is flow An the the of magnitudes of exampleof a casewhere this velocity. an underestimation involving in is Fig 6.23. For shorter the simulations error affected predictionresults shown pipelines,this error doesnot affect the programmeand the variablegrid model produces good results. A bigger grid spacing, of Ax=10m at the break, was used for the SNGSO data. Simulation results are presented in Figs. 6.23 to 6.25. The computer round-off error did not cause serious problems in the solution.

The first-order method of characteristics was

intact biggest (at The the grid size end) was ix=320m. used.

Better results than those

With break, Ax=1m, the the obtained. method and at were second-order produced with the first-order method and the bigger grid spacing, the equalization pressure reached the ambient value within a few tens of milliseconds, while with the second-order method and the smaller grid spacing, ambient pressure was not reached during the 25s run time. This illustrated the significance of the grid size on the predicted decompression behaviour. This is one of the parameters which were investigated in this study. It was observed that for tests in which the positions of the transducers were very close to the broken end, such as the British Gas tests, a finer grid size was required. For tests in which the transducers are positioned further away from the broken end, such as the Foothills tests a coarser grid size produced sufficiently accurate results. A grid size of 0.1 m at the broken end was adequate for the Foothills tests, where most transducers were positioned tens of metres away from the break. For the British Gas tests, where the transducers were positioned within two metres away from the broken end, a grid spacing of 0.1m was not adequate. As seen in

264

Fig. 6.21, the wave speeds produced by the grid spacings of 0.01 and 0.1 m are almost the falling for 6.19 6.17 to Figs. However, the that the time taken to start show pressure same. is longer when a smaller grid spacing is used. With test where the results are required broken 10m, from break, Ox the end, to the at thousands of metres away values of up Industry Petroleum for Alberta This the the case even was produced satisfactory results. from is Im the flow away tests, where the mass rate was calculated at a position which

break. A simplerule of thumb is recommendedfor selectionof grid sizes. If the result is length from break is few the the of the and/or metresaway requiredat a position which a be between 0.1 long, 0.01 is Ax less 100m should than and a value of ranging pipesection break from is is the If tens the away of metres result required at a position which used. l. Om between long, 0.1 is hundreds Ax and the metres a ranging of value of pipe and/or further 100m if is longer be For that the sections and/or required pipe result should used. break, be from Ax 10m to the of of up values could used. away EXPERIMENTAL DATA All the experimentaldata usedto validate the computer model were availablein the form of printedgraphs.The datawas convertedinto a numericalform in order to enableplotting data In together the the the order same graph. results with experimental on predicted of to simplifythe conversionof the dataand also to minimize errors, the graphswere scanned and the necessarycoordinates required to reproduce the graphs were read using the Autocadgraphicssoftware. The accuracyof the dataobtainedis just as good as that of the original graphs. PRESSURE Both the method of characteristics and the MacCormack method predict the pressure following a break reasonablywell. With the MacCormack method predictions, the pressure flattened faster than with the method of characteristics. This is because of the higher speeds of sound which are obtained from the MacCormack method predictions. Unlike with the MacCormack method, the method of characteristics results are not oscillatory and a consistent pattern is maintained, whereby the pressure decreasesas the broken boundary is approached and for any particular position in the pipeline it decreaseswith time.

The pressure plateau which was demonstrated by Jones and Gough (1981), signifying the two-phase region, was not observed. The reason is that the QUANT 265

for for liquid does the transition region data the and phase not produce output software from liquid to gas phases. TEMPERATURE

Temperaturevariationsin a gas undergoingdecompressionfrom a high-pressurepipeline, is one of the parameterswhich has not been studied much. Experimental data for for is temperaturevariations very scarceand nonewere obtained validation of the computer model. It is obviousthat the temperaturedrops after the pipeline breaks. There is a debate temperature the magnitude of maximum the relative regarding exists which over uncertainty dropsresultingfrom high-pressurepipeline breaks. Richardson(1993-96) arguesthat the temperaturedrop could be as high as 50K at the break. In this study temperaturesfell to below the 200K minimumlimit of output for the QUANT software (drops of up to 100K) in the region near the break, when the method of characteristicswas used. With the MacCormack method, the temperature drops were much lower. Flatt (1993-96) and Silberring(1993-95)arguethat the temperaturedrop could not be that big and should not fall below the minimum limit of the QUANT software. The temperaturedrops predicted using the method of characteristicsare much higher than those predicted using the MacCormack method.

In some cases the

temperaturesfell below 200K when using the method of characteristics,but stayedwell initial MacCormack For 200K the the conditions of the method. same using above when in Ud Also drop higher the temperature some smaller pipes of were used. when was gas, fell broken boundary Foothills NABTF7, for temperature the test the the at cases, example below 200K just when the break condition was introduced i.e. t= Os. It was not possible to know exactly what the temperaturewas, in caseswhere it fell below 200K. For such for 200K temperatureand the valuesat the 200K temperature the used cases, valueof was fluid from for QUANT An the the properties required error other software. used were would be introducedin the calculationespeciallyin the heat transfer. However, the error introducedis smallbecausethe pressurein such situationsis low and therefore variations in temperature would cause a small variations in the thermodynamic and transport propertiesof the gas.

266

FLOW VELOCITY The flow velocity at the broken end at time t= Osafter the break is equal to the speedof falls before break. flow The break falls the the the sound sound velocity at speedof as Jul flow is (M, 1) the to the when choked = according choking condition = a. When the pressureat the brokenendfalls to ambientpressure,the choking condition no longer exists in flow be lower flow The the than the the sound. speedof and magnitudeof velocity will this low pressureregime is not investigatedbecauseof the long time required for the condition to be attainedand the limited time which is availableto test the model. The flow velocity predicted using the MacCormack method is higher than that flow Also former the the the characteristics. predictedusing methodof whenusing method, fasterandspreadsout fasterthan with the latter method. As a result, the velocity increases velocity gradient is much less with the MacCormack method than with the method of characteristics. FLUID DENSITY After the break, the density of the fluid falls in a similar manneras the pressurebut by a lesserproportion. With the MacCormackmethoddensity predictionsare lower than those predictedusing the method of characteristics. PRESSURE WAVES PROPAGATION SPEED The valueof the speedof soundin the fluid (a), at the break, remainsthe sameat t= Osas it was before the break (Refer to Fig. 4.3). Due to the temperaturedrop which continue to take placethereafter,the densityof the fluid at the break increasesslightly. The increase in densityresultsin a drop on the speedof sound,accordingto equation(A-11), and hence a drop in the flow velocity accordingto the chocking condition Jul= a. Sincethe drop in pressureis higher in proportion than the drop in density, the speedof sound drops as the pressuredrops. This patternis transmittedtowards the intact end of the pipe as time goes and is lessrapid away from the break and as time passes. The predicted output of "a" differs significantly from the pressurewave speeds presentedin experimentaldata. However,the predictedoutput valuesof "a" comparewell with thosepresentedin other datasuchasthe additional dataprovided by Silberring (1995) and the data presentedby Straty (1974) and Tsumuraand Straty (1977). The reasonfor

267

the discrepancybetweenexperimentalwave speedsand predictedvaluesof "a" is that the valuesof the speedof soundmeasuredin the experimentsis not surprising. The value given in experimentaldata is the speedat which a pressurewave propagatesin the fluid i.e. it is from it for by travel to the time takes one a particular pressure wave calculated measuring known position to another known position in the pipe. The values of "a", which are front. leading This (A-11), to the the pressurewave refer edgeof calculatedusingequation from be faster the p-t curves. The pressurewave speed than value calculated a mean will is influencedby both the valuesof a and u as explainedby the theory of characteristicsin Section4.3.1. The value of the pressurewave propagationspeedcalculatedfrom the p-t curvesresulting from the computer simulation were comparedwith experimentaldata for the FoothillsNABTF 1 andNABTF7 tests,the BMI testsandthe University of Calgary test. The results are presented in Figs. 6.5,6.7,6.20

and 6.21.

Only the method of

characteristics andthe MacCormackmethodswere used. It could be seenfrom the graphs that apartfrom the Universityof Calgarytest,the methodof characteristicsproducesresults which are very close to the experimentaldata. The wave speedscalculated using the second-ordermethod of characteristicsare slightly lower than those calculatedusing the first-order method. The MacCormack method underestimatesthe wave speedsat higher pressuresand producesvalueswhich are slightly higher than experimentalvalue at lower pressures. Data which was obtainedusing the second-ordermethod of characteristicsfor the Foothills NABTF7 test could not be used for validation becausea wrong value was used for the initial density. The wave speeds obtained were very much lower than the experimentalvalues. A finer grid spacingAx=0.01mwasusedfor the University of Calgary test,in an attemptto obtaina betteragreementwith experimentaldata. The result obtained is almostthe sameaswith the grid spacingof Ax--0.1m. It has beenconfirmed that further reduction of the grid spacingswould not produce more accurateresults. Also sincethe computer model has produced good results for all the other tests, it is believedthat the discrepancybetweenthe two sets of data in the University of Calgary test is causedby experimentalerror. COMPUTING RESOURCES The biggestconstraint,asfar asthe computermemoryrequirementand computation speed are concerned,was introducedwhen non-constantfluid property data were used. Initially, 268

the programmeswere run with constantfluid property data, which were obtainedfrom the QUANT software. Even though the basic equations are used almost without any simplifications, the computation was very fast. The introduction of varying fluid data greatly reducedthe computation speed. In situations where the initial pressurebefore the break is low and temperature increase fluid data is the the of constant property greatly use will output not required, computation speedwhile maintaininga reasonablygood accuracyof results. When the it in this programmeswere run manner, was possibleto use a 486 pc with a RAM memory of 8MB. When varying fluid property data was used the 486 pc could not cope and a PENTIUM P75 was sufficient to run all the programmes based on the method of characteristicand the MacCormack method. The programmesbasedon the third-order Warming-Kutler-Lomax method was successfullycompiled but could not run on the PENTIUM pc becauseof insufficient memory. The problemof insufficient memory could be solved by either using a bigger pc or a mainframecomputer. The UNIX basedC++ compilers run the programmeswhich are usingpc basedBorlandC++ (in this caseVersion 2.1) compilerswith very few alterations, namely adding an "include file" stdlib.h and changing the commandsin the "system()" statementsto correspondwith those usedby the UNIX operating system. A user memory of 20MB on the mainframecomputer could only perform as well asthe 486 pc. Evenby usingthe "temporary" directory in the mainframecomputer, which shouldhavemorecapacityavailable,did not improve the situation. The remainingoptions are therefore to use either a more powerful pc than PENTIUM P75 or the mainframe computerwith morememoryallocationthan the present20MB. Computation speedof the pc could be increasedby installinga chip similar to the Microway-i860 which is being used at presenton a 286 pc, at the Imperial College [Richardson(1993-96)]. DISCUSSION OF ERRORS Three main categoriesof error in validating computer modelsfor linebreakanalysis,with experimentaldata were discussedby Tiley (1989). The three categoriesare

(i) Calibration,measurement data. andrecordingof experimental (ii) Assumptionsand simplificationsmadein the basictheoretical equations. (iii) Errors inherentin the numericalmodelling procedure.

269

The error in experimentaldata also includesthe error in converting graphicaldata into numerical data. Regarding the error in preforming the experiments and recording experimental data, the best one could do is obtain an estimate of the error so that it could be accounted for in the simulation No error estimate was provided with the experimental data. Even the rupture time was not provided and in some casesnot all parameters required by the computer model e.g. gas composition, pipe materia etc were provided. Through informed Gas, (1989) British Tiley that the pressure was private communication with measuring system as a whole, in the British Gas tests was believed to be within 5%. The data in which is provided in graphical form into numerical error converting experimental data has been reduced significantly by employing the scanning techniques and Autocad. With this procedure, the accuracy of the results is as good as those of the computer is better than using a manual technique such as the one which much software used, which was used by Tiley. Tiley (1989) quoted an accuracy of±0.01

and 5m/s for pressure ratio

and wave speed respectively.

The error dueto assumptionsandsimplificationsin the theoreticalmodel hasgreatly been minimized in this model, compared with Tiley's. This model still contains some simplificationssuchas one-dimensionalflow, single-phaseflow, non-elasticpipe, no fluid structureinteractionand neglectingminor losses. However, the error introduced by these assumptions should be minimal since in this study, only straight horizontal pipes with have beenused. The error due to inaccuratefluid properties constant cross-sectionarea hasalsobeenminimizedby using the QUANT software for thermodynamicand transport properties of fluids. There still remains some possible error in calculating the fluid propertiesdueto the limitationsof the QUANT software. Thesehavebeendiscussedunder the discussionof the QUANT software. The error is significant in simulationswhere the initial gaspressureexceeds5 to 6MPa andtemperaturesafter the break fall below the 200K minimum limit of QUANT. The errors in estimatingthe friction factor and heat transfer havebeenminimizedby usingflow dependentvalues,which are specificfor eachcalculation step. They assumedthat the friction factor and Stanton number(used to calculateheat transfer) were constantalong the length of the pipe. The errors inherentin the numericalmodelling procedureinclude smearing,when a fixed grid is used; round-off error incaseof iterative methods;and computer round-off error. Tiley (1989)arguedthat the smearingerror is minimizedby using a smallergrid size, while the round-off error is increasedby reducing the grid size. There seemsto be a contradictionin this argumentand the whole idea of using a fine grid spacing. The roundoff error of iterativemethodsdependson a predeterminedaccuracycriteria (as long as the 270

solution converges) rather than the grid size. In the model developed in this study, an accuracy of 5% and 1% was specified for the first- and second-order methods of in is been interpolation has It this that study not necessarily observed error characteristics., reduced by using higher-order approximations. In actual fact, approximations of higher order than one produce wrong results for rapidly varying flows. first-order approximation

hasprovedto bethemostpopularin modellingtransientflow followinglinebreakin highbecause Tiley's model contained an additional of the method error pressure gaspipelines. length final the to the the required result at positions used obtain along of the pipe. The methodusedby Tiley wasto approximatethe result to the valuesat the nearestgrid point. Such approximation could introduce a very big error in the solution, especially if the positionconcernedis close to the broken boundary and in the early time regimewhere the fluid propertiesvary rapidly. The approximation could also introduce significant error at positionsfar from the breakbecauseof the biggergrid size. In the model developedduring this study,linearinterpolation is used to calculatethe final result at the required positions. It is not possibleto establishthe magnitudeof the accuracyof the computer model developedin this study,with certainty,becauseof the poor quality of the experimentaldata which has been used for validation. However, in most casesthe predicted results are in good agreementwith the experimentaldata used.

271

CHAPTER

7

CASE STUDY: THE SONGO SONGO-DAR ES SALAAM NATURAL GAS PIPELINE 7.1 A BRIEF OVERVIEW OF TANZANIA'S ENERGY SECTOR The United Republic of Tanzaniais located on the EasternCoast of Africa, just south of the Equator. The country hasa total area of 945000 squarekilometres, 6% of which is half land is by by Nearly forest reserves,which explainsthe the of covered covered water. great overdependenceof Tanzania on woodfuel as a source of energy. The present populationof Tanzaniais about29 million, with a population growth of 2.55% per annum. Dar es Salaamis the capital of Tanzania and its population is about 5% of the total population. The economy of Tanzaniais predominantly agricultural. Other economic include industries, and processing activities manufacturing mining and tourism. Most of the industriesare concentratedin urban centresand in particular Dar es Salaam. Tanzania'sindigenousenergyresourcesare large and diverse,although they have not yet beenexploited to an appreciableextent. At present,the major energysourcesare hydroelectricity imported petroleum, and indigenous coal. woodfuel,

Tanzania's

hydroelectricpower potentialis in excessof 4700MW of installedcapacity,but only about 10%of it hasactually beendeveloped. Coal reservesare estimatedat about 1900 million tonnes, of which 304 million tonnes are consideredproven. Exploration activities for petroleumbasedfuels havetaken placein Tanzaniafor severaldecades,with the resulting discovery of natural gas. No oil discoverieshave been madeyet and exploration is still underway, with prospectsof findingmorenaturalgas and possiblyoil. The proven natural gasresourceis 20.7 million tonnesoil equivalent(toe), but it is estimatedto increaseto 47 million toe. The total energyconsumptionin Tanzaniain 1990wasjust over 17 million tonnes of oil equivalent (toe) which came from woodfuel (88.5%), biomass (6%), electricity (0.6%), imported crude oil and refined petroleum products (5%) and coal (0.1%). More than 70% of Tanzania'selectricity is generatedfrom hydropower. In the past six years, Tanzania'selectricity demandhas beengrowing at an averageof 12% per annum. The presentelectricitydemandis estimatedat 400MW, but the demandis suppresseddue to the shortageof generationcapacity. The demandis estimatedto double by the turn of this 272

h

P.

ý

ýD" CM vt-n

S.O..... e

J

Map of Africa Showing the Location of Tanzania

Fig. 7.1

I=

ýur. oy

i KENYA

r_ ty

., r

ý'ýf""ý

ZAINF

ý-l `ý\

ZA

,'M

61A

""

M OWA

raw

Fig. 7.2

AJI BIQUE

......

ý,

Map of TanzaniaIndicatingthe Location of the SongoSongo GasDevelopmentProject 273

investment development (1985-2010) five A twenty and programme was set. year century. The programme is aimed at reducing dependenceon hydropower and cutting the oil import bill that consumesabout 60% of the country's foreign currency earnings each year. During

the early years of the 1990's,Tanzaniasuffered serious power problems,as a result of hydropower in levels drought the plants. of reservoirs which reduced water prolonged Power rationinghadto be imposedand as a result industrial production fell drastically and lifestyles had to adjust to the situation. their to change people The Tanzaniangovernmentis implementingan EmergencyPower Project (EPP). A 40MW dieselpoweredplant hasbeencommissionedin Dar es Salaam. Electricity from the plant, which could also use natural gas as a fuel, is linked to the national power grid thereby boosting its output. On the longer term, there are other projects including the development hydropower Kihansi, 200MW plant of the coal resource . at constructionof a is focus Gas Songo Natural Development Project, Songo the the of this study. which and The project is describedin detail in Section 7.2.

7.2 DESCRIPTION OF THE SONGO SONGO GAS DEVELOPMENT PROJECT 7.2.1 BACKGROUND INFORMATION Naturalgaswas discoveredat SongoSongo,which is a tiny island off the Tanzaniancoast, in the Indian oceansouth of Dar es Salaam,in 1974. Apart from the Songo Songo gas including discovery Mnazi The Bay. there at was are other smaller reserves, one reserve, followed by proving between 1976 and 1983 and subsequentpre-feasibilityand feasibility studies. Severaloptions for developingthe resourcewere considered,including export of the gasto earnforeignexchange,productionof ammoniaand urea fertilizer and usesin the domesticmarkets. The latter option includespower generation;pipeline to local markets; and fertilizer, methanol and compressednatural gas facilities. It was consideredthat, althoughthe gas pool was large, it was not large enoughto warrant the developmentof a major compressed natural gas export facility.

The remaining options are being

implemented and in addition, there are plans to export some of the gas to neighbouring Kenya. The SongoSongoGas DevelopmentProject featuresdevelopmentof five existing gas wells, of which three are offshore and two are onshore;tying thesewells to a central 274

island; (25km Songo in Songo marine and a pipeline constructing gas processingplant in building 60MW Dar Salaam; to the power plant 207km onshore)to transport a gas es Dar esSalaamto generateelectricity,andsupplyinggasto the 40MW dual fuel power plant describedin Section 7.1, which is alreadyoperational in Dar es Salaamand also to other Bay for Mnazi implemented is being the but different, reserve. Another, project users. Water jointly by Ministry implemented is being the Songo Songo of The project Energy and Minerals, which plays a coordinating role for all stake holders in the project; Ocelot Tanzania Inc.; and TransCanada Pipelines Limited.

The Tanzania Petroleum

Development Corporation (TPDC)is a partner in the project and the Tanzania Electric Supply Company (TANESCO) will be the major consumer of the gas. The project is in in be 1997. to service anticipated

7.2.2 GAS PROPERTIES Natural gas is found naturally in rock reservoirs below the ground. In its pure state (methane)it is a non-toxic,colourlessand odourlessgas. Methaneis lighter than air. The Songo Songo gas is almost pure natural gas, containing no sulphur compoundsand only is hydrocarbons. Only heavy minimalprocessing required to make the gas smallamountsof

by for the endusers. transport application and subsequent suitable pipeline The Songo Songogas field underliesSongo Songo island. The proximity of the field Fig. 7.3 to the islandis very convenientin that the shorebasedfacilities can be located been have in Three drilled island the wells the onshoreor very shallow water. on and wells drilledonshoreand a further six havebeendrilled in the shallow waters around the island. The ninewells arenamedasSS1 to SS9. Five of the wells have proven reservesin various four have disappointing field Of three the the wells, provided remaining structure. partsof resultsin satellitestructuresand one did not reachthe main reservoir objective as a result of blowout and fire. The main areaof the field hasbeenappraisedby the five wells which have proven the gas volumesin the main block i. e. wells SS3, SS4, SS5, SS7 and SS9. The reservoir fluid propertieswere determinedin a study by Scientific Software-Intercomp(1990), for simulationpurposes. The data includesthe relevantphysicalpropertiesof gas and water. A total of2l analysesof samplesobtainedfrom wells SS5, SS7 and SS9 were conducted. The sampleswere very consistentin composition, exhibiting a high methanecontent of 275

ýDýA

E9 SAIAAM

b4u1-^pl 0

10

20

30

a-W

40

60 K

(k+1 Ip Wmiwon

MAFJ& ISLAND

hot

e"s -

Nye """ Scm. n9s s«pos«wo

@(Hnio ýý i

c

ati. w. o.o

es Salaam Gas Pipeline Location

Fig. 7.3

Kilwa Kivinje-Dar

Fig. 7.4

Pipeline Route from Songo Songo Island to Dar es Salaam 276

Map

fraction low heavy hydrocarbons. 97% a of condensate and and about

The average

is in Table 7.1. the shown gas compositions of The initial reservoir pressureand temperature were specified in order to process the

fluid propertiesdata. Given the averagegas composition data in Table 7.1 and reservoir pressureof 2755psi (19MPa) and temperatureof 203°F (368K), the gas propertieswere determined. In calculatingthe gas properties Scientific Software-Intercomp(1990) used their own equation of state. The gas properties which are relevant to this study are presentedin Table7.2 andcomparedwith valuescalculatedusing the QUANT software for thermodynamicand transport properties,which is usedin this study.

7.2.3 PRELIMINARY

PIPELINE DESIGN

The pipelinesystemincludesthe transmissionpipelinefrom Songo Songo to Dar es Salaam and the distribution network in Dar es Salaam. However, it should be noted that the pipelinedesigndescribedin this sectionis basedon the preliminary designdone by Hardy BBT Ltd. (1989). The gasdistribution systemconsistsof a pipeline network which would deliver gas from the city gate station to each individual consumer. The design is in accordancewith CanadianStandardCSA Z184.M Code. Three demandscenariosnamely low, medium and high were considered. The three scenariosare basedon a 20 year life time of the project and projections for power generationby TANESCO. The 20 year life time is based on the design of a fertilizer plant at Kilwa Masoko (refer to Fig. 7.3). Accordingto the design,the proven gas reserveof 20.53bcm(726bcf) was to be drawn at a rate of 0.56bcm/year(19.75bcf/year), by the fertilizer plant, so that its allocation of 11.18bcm(395bco is consumedin 20 years. The operation of the gas gathering facilities, dehydrationplant and marine pipeline to Kilwa Kivinje was to be under one management to supplygasto the fertilizerplant andthe Dar esSalaampipeline. Therefore, sincethe Dar es Salaam pipeline was to be dependenton the fertilizer plant operating, all design and economicstudieswerebasedon a 20 yearsof operation only. However, in the likely event that the proven reserveswere to be enlargedas production proceeded,the pipeline was sizedbasedon growth through the year 2016. The designbasedconditionswere a pressureof 101kPa(14.65psig), a temperature of 15°C anda specificgravity of 0.60. The maximumallowable operating pressureswere set at 7.93MPa for the low demandscenarioand 9.65MPa for both the mediumand high demandscenarios.The maximumallowablegas temperatureis 49°C, while the maximum designgas temperaturewas set at 38°C. The designgas volume took into consideration 277

the averageandthe peakflow requirementsto the year 2016. The projected gas demands for 2011 in 2016 Table 7.3. the the three years and are presented under scenarios GAS

SS9

SS5

SS7

MEAN

ISM NO.

S VAR

2 02

0.860 0.47

0.68 0.350

0.600 0.290

0.713 0.370

0.000 0.000

0.000 0.000

H4

96.820

97.193

97.440

97.151

0.000

0.000

1.050 0.320 0.070 0.090 0.030 0.030 0.040

1.100 0.300 0.070 0.089 0.025 0.026 0.030

0.940 0.310 0.073 0.088 0.028 0.030 0.025

1.030 0.310 0.070 0.089 0.027 0.029 0.032

0.000 0.000 1.000 2.000 1.000 2.000 1.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000

7H16

0.150

0.075

0.100

0.109

1.000

0.000

8H18 91-120

0.060 0.010

0.044 0.018

0.053 0.023

0.052 0.017

1.000 1.000

0.000 0.000

100.000

100.000

100.000

100.000

2H6 3H8 IO -C4H -C4HIO -C51112 -C5H12 6H14

OTAL

Sourceof GasCompositionData: Scientific Software-Intercomp[1990] Table 7.1 Average [Mol. %] Compositionof Songo Songo Gas Pipe sizesand flow capacity of the pipelinewere determinedusing the Panhandle formula, which considers gas pressure, temperature and velocity. Steel, plastic and aluminium were all consideredas possible materials for construction of the pipeline. However,it was laterdecidedthat only a steelpipe would be able to carry the volumesand pressuresrequired for the transmissionline. The gas velocities in the pipeline, for each scenario and pipe size, were estimatedas shown in Table 7.4. The pipe sizes allow fluctuationswhich mayoccurbetweenthe peak day demandand the averageflow, and the expected velocities are well below the target velocity for the pipeline design. In this particular case flow velocities of 457m/min and 610m/min were given as target and maximumvelocities respectively. Initial Pressure [psia] Reservoir Temperature [°F] [K] Molecular weight Specific Gravity Density (Reservoir) [Ib/ft3] [Kg/m3] Density (Surface) [Ib/ft3] [Kg/m3] Viscosity (reservoir) [cp] -1

2755 203 16.746 0.579 7.1007 0.0441 0.018 4

Source:Scientific Software-Intercomp(1990) Table 7.2 SongosongoGasProperties 278

Pipe sizes were chosen for the three scenarios,basedon the initial pressureof 4.80MPa availableat Kilwa Kivinje and a minimum terminal pressure of 2.07MPa (300psig) demand high both in Salaam. For Dar the scenarios, and medium es at the city gate station Kivinje. be Kilwa at added compression must

The selected pipe sizes, compressor and

in 7.5. Table summarized resulting capacityare DEMAND LOW MEDIUM HIGH

YEAR201 0.56 1.21 2-31-

2016 0.68 1.33 *243

YEAR 2011 0.7 2.37 3-47

YEAR 2016 0.85 2.49 5-59

Source:Hardy BBT Ltd (1989) Table 7.3 Projectedgas demandfor three scenariosfor year 2011 and 2016 Predictedpressureprofiles for the line for all the three scenariosare presentedin Figs. 7.5,7.6 and 7.7. For each scenario there are three curves representing maximum

design flow discharge (300psig), 2.07MPa flow pressureof consistentwith a sustainable flow for daily demands the three and rate of to scenarios rate similar maximumaverage 70% of the abovedesignflow rate. The distribution systemis comparativelysmall so that be in it demand, local it to consideredas a major can not variations while serves smoothout in is by in the The this storage case provided on-line major storage componentof storage. transmissionpipeline. SCENARIO LOW MEDIUM IHIGH

VELOCITY [mix/min] 198 235 256.

NOMINAL PIPF SIZF [fi/min] [M4 650 254 770 254 840,305 .

[ial

10 10 12.

Source:Hardy BBT Ltd (1989) Table 7.4: EstimatedGas Velocities and Pipe Sizes Threeblock valvelocationswere proposedin the preliminarypipelinedesign. These would permit eachsection of the transmissionline betweenthe valves to be isolatedand blown down during repair, maintenanceor in an emergencysituation. Operation of the block valves would be in the manual mode. Automatic linebreak control was not consideredin the preliminarydesign. It was suggestedthat additional block valves could be added to the systemduring the final stage of the project dependingon operational

considerations. The three block valve locations are presentedin Fig. 7.7. The 279

line. full-bore in facilitate be The block to the pigging of valves order will sectionalizing blowdown valveswill be sizedto depletethe line section in lessthan 4 hours. The block demand blowdown for time the three scenarios are of each valve parameters and it However, in 7.6. Table was not specifiedto which of the pipeline sections summarized in Fig. 7.7 or demandscenariosthe blowdown times in Table 7.6 refer. DemandScenario Low Medium -High

PressureCompressorSize CapacityFlow Pipe size Actual sizexWall thickness IMC M/da M MP mm 0.78 273.1x4.78 1.83 2.16 9.65 273.1x5.20 3-39 285 9-65 323-9x6-35

Source: Hardy BBT Ltd (1989) Table 7.5 Pipe Parameters The original designof the pipeline, as presentedby Hardy BBT Ltd. (1989), was for the pipeline to passthrough Kilwa Kivinje and also to be tied to the project for the fertilizer plantat Kilwa Masoko. Recentinformation,including a brochure publishedby the Ministry of Water, Energy and Minerals and the two Canadiancompaniesinvolved in the implementation of the project, indicate that a new pipeline route is being followed. The new route is shown in Fig. 7.4. According to the new route, the pipeline is from Songo Songodirectlyto Dar esSalaam,throughSomangaFungaand not Kilwa Kivinje. The new

pipelineroutedoesnot seemto be tied to anyotherprojectapartfrom the powerstation is it design details Although for the new in Salaam. Dar that expected and other users es route havebeenfinalised,it was not possibleto get them for this study. However, sincethe in from length differ does the one the proposedold route, and assuming pipeline much not

in Dar esSalaamremainthe same,the samedatawhichwas thatthegasdemand scenarios proposedfor the old pipeline designis used. DEMAND SCENARIO Block Valve Size Blowdown Stack LOW MEDIUM IHIGH

254 254 305

102 102 152

Source:Hardy BBT Ltd (198 Table 7.6 Block Valves

280

Blowdown Time 3.5 4.0 13

7.3

COMPUTER SIMULATION OF A LINEBREAK SONGO SONGO DAR ES SALAAM GAS PIPELINE

IN THE

In this section,the information and data presentedin the precedingsectionsof this report and the computer model developedin this study are used to simulate a linebreak in the SongoSongoDar esSalaamnaturalgaspipeline. It is known that a more recent study was conducted by Hardy BBT Limited in 1994, to assessthe environmental impact of the pipeline. However, the study report could not be madeavailablefor this study. Also as pointed out in the previous section, no further details on the new pipeline route could be obtained. Due to theseconstraints,it was decidedthat the data given in the preliminary pipelinedesignby Hardy BBT Ltd. (1989), which are presentedin Section 7.2.3, be used for this simulationexercise. The main reasonsfor this analysisare firstly, to simulate the flow following a linebreakin the Songo SongoDar es Salaampipeline and secondlyto test if the computer programme is capableof handling such a long pipeline as the pipe sectionsbetweenthe block valvelocationshownin Fig. 7.7 (73.2km). Pipelinesare equippedwith block valves different locations along their lengths. The valves enableparts of the pipeline to be at isolatedfor maintenanceand in the event of a leak, blowdown or accidentalrupture. The valvelocationsserveas boundarypoints when performing the computer analysis. Even in situations where such valvesare not provided and the pipeline is considerablylong, it is convenient to imposeboundary conditions in such a way that the length of the pipeline betweenthe boundary points could be handledwith the computer programme. Two main parametersare of interest,namelythe time it takes for a pressurewave to be transmitted from the break point to the nearestblock valve location and the time takenfor the sectionof the pipelinebetweenthe two valveson eachside of the break to be emptied. The former parameter determines how quickly the valve operation would be initiated. The limited time available for this study does not permit simulation of the flow until all the gas has escapedfrom the pipeline, because of the very long CPU time required to carry out the simulation. The simulation is performed only for the early time regime. Referring to Fig. 7.7, it is assumed that a rupture occurs at Marendego, just next to the block valve, on the side of Kilwa Kivinje. The interest is to determine how long it takes for the initial pressure wave to reach Kilwa Kivinje, which is 68.0km away, and initiate the operation of a block valve situated at the exit of the compressor station. The high demand scenario and maximum sustainable flow (the bottom curve in Fig. 7.8) are used. 281

, V-

1.6 1.5 1.4 1.3 1.2 1.1

IV^

1

W

0.9

u ,y y

0.8

ýo a -Uýq

0.7 0.6

0.5 0.4 0.0 0.2 0.1

0 00

40

0

flow Rntc fý1MCFD

120

distance

0

+

10.13

200

160

from source km 25.90

0

210

27.99

Source:Hardy BBT Ltd. (1989) Line PressureProfile for Low Scenario

Fig. 7.5

1.8

1.4 1.3 1.2. 1.1

b^ äG v b Ho ILI

1 0.9 0.8 0.7 0.8 0.5 0.4 0.3 0.2 0.1 0 0

00

40

Flow Rate NINICFD

O

35.00

distance +

120 from source 50.00

160

200

km 0

04.73

Source:Hardy BBT Ltd. (1989) Fig. 7.6 Line PressureProfile for Medium Scenario 282

240

The averagecompositionof the SongoSongo gas, which is presentedin Table 7.1, is usedasinput datato the QUANT software. A data file to be usedby the CFD model is generatedfrom the QUANT software. The initial pressureand massflow rate of the gas at Kilwa Kivinje are as presentedin Fig. 7.8 i.e. 141Opsi(9.65MPa) and 100.79MMCFD (25kg/s). The maximumoperating temperatureof the gas is 311K. The rest of the gas properties at the initial conditions at Kilwa Kivinje are calculated using the QUANT software. A minimumgrid spacingof Ax= 10mand At= 0.01s at the broken end and the variablegrid modelareused. This rathercoarsegrid spacingis chosenin order to minimise the CPU time required for the simulation and also to minimisethe effect of the computer round off error, which is explainedin Section 6.3, when applying the variablegrid model on long pipelines. Resultsfrom the computer simulation are presentedin Figs. 7.9 and 7.10.

Mar4--

Ki1wa 1C1v111, 0 krra

A

681cm

radego

Buingo

Mbezi 199.8km

141.2km

Dar

ea

Salasm

Fig. 7.7 Block Valve Location

The Hardy BBT prediction data which is presentedin Figs. 7.8 are used for comparisonwith the data calculatedusing the non-isothermalnon-adiabaticsteady state analysismodelpresentedin Section4.2.2 and also the transientanalysis,before the break, using the method of characteristics. The prediction resultsobtainedare included in Figs. 7.8. All the rupturetestswhosedataareused,in Chapter6, to validate the linebreakmodel were carriedout usingpipelinesin which the gaswas initially stationary. As a result of this situation, transientanalysisbefore the break was not necessary.In the caseof the Songo Songo-Dar es Salaampipeline, the initial flow velocity before the break is not zero. Transientanalysisbeforethe break was performedin order to establishthe initial unsteady flow condition. Resultsare presentedin Fig. 7.8. In both the steadystate and transient analysesbefore the rupture, the variable grid spacing programmesare used. The heat transfer model which was used is that for buried pipeline.

283

O

a O

p

ü c

C p

c

ýQ CL

E 0 :.) O a

tu c d

Co Z 8 I-

ö

8

0 a

ci

N

NI 93 c3 o o

oo°

a U

98

gý8

_

88

8

TTýýý I2

O Of

Co

1-

(0

U)

lt

V)

[edw] amssaJd

284

fV

. -0

LL

10--

7-CL 6

-100.79MMCFD SSA --111--100.79MMCFD TA MOC-1 (dX=10m) BREAK MOC-1 (dX=10m) -ýk-100.79MMCFD

2

4--

a 3 2 1 -I

0 30

20

10

0

40

60

50

Distance tkm]

Fig. 7.9 Line Pressure Profile for Simulated Break in the Proposed Songo Songo-Dar es Salaam Pipeline

2500--

2000-

1500.. ö yN

1000

0

10

20

30

40

50

60

Time (s]

Fig. 7.10 Mass Flow Rate for Simulated Break in the Proposed Songo Songo Dar es Salaam Natural Gas Pipeline

285

70

7.4 DISCUSSION OF SIMULATION

RESULTS

Simulationresultsfor the pressureprofile usingthe steadystateanalysis(SSA) and transient analysis(TA) modelsbefore the break are presentedin Fig. 7.8, together with the original data provided by Hardy BBT Ltd. (1989). The pressureprofiles for steadystate analysis and transient analysisbefore the break are also presentedin Fig. 7.9, together with the break. The massflow rate variation is presentedin Fig.7.10. 60s the after pressureprofile In calculatingthe pressureprofile using the steadystate analysismodel problems were encounteredwhen the pipe diameter of 0.305m, which was, specified was used. Usingthis pipediameterandgasflow rate of 25kg/s, which were given, resultedin a much bigger pressuredrop in the pipeline. The pressurepredicted at Marendego (68km from Kilwa Kivinje) was around 4MPa, which is less than a half of the pressurepresentedby Hardy BBT Ltd. Transientanalysisusing the second-ordermethod of characteristicswas performed, in order to confirm whether the problem was with the steady state analysis model or the design data provided by Hardy BBT Ltd. When a boundary condition correspondingwith the pressurepresentedin the Hardy BBT curves(8.2MPa) was imposed negativevelocitieswereobtainedat the low pressureend. The negativevelocities indicate that the pressureimposed at the boundary was higher than the actual pressurefor the parametersgiven. The massflow ratewas reducedby a half in order to seeif the pressure drop would decreaseto the required value, but the pressuredrop remainedhigher. When the pipe diameterwas increasedto 0.5m, the samepressureprofile as that presentedby Hardy BBT wasobtained. Transientanalysiswas then performed and producedthe same pressureprofile asthe steadystateanalysis. However, after the transientanalysis,the flow velocity was slightly higher than that obtainedusing the steadystate analysismodel, but within the limit specifiedby Hardy BBT. This result indicatesthat the pipe diameter of 0.305m which was specifiedis too small. The first-order method of characteristicsfailed to produce good results in the transient analysisbefore the break. Consequently,the second-ordermethodwas used. There was no output from the QUANT software for Pr, p and kr for all the range of pressureandtemperatureused. Of the threemissingproperties,only p is required by the computer programme. A constantvalue was used throughout the range of pressureand temperatureencountered. In the fluid data usedfor validation of the model in Chapter 6, a maximumvariationin p was 36.5%within the rangeof p and T used. The constantvalue was chosen such that it is in the middle of the range, therefore reducing the error in estimatingp to a maximumof 18.25%. The error introduced is not significant.

286

The first-order methodof characteristicswas used in order to minimize on CPU time. The pressure profile after the break, which is shown in Fig. 7.9, indicates that the leading edge of the pressurewave had not reached the compressor station at Kilwa Kivinje 60s after the break. It had reached a position which is 43.52km from Kilwa Kivinje (24.48km from the break). Due to the limited time available for this study, no further simulations could be performed. The mass flow rate variation is presented in Fig. 7.10. The initial mass flow rate at time tos after the break is less than the peak flow rate. The increase in mass flow rate is caused by increase in the gas density, which occurs after the break because of drop in temperature. During the period of 60s which the flow was simulated, most of the rapid variations in gas outflow took place. The mass flow rate would continue to decay exponentially. The flow rate of around 300kg/s, which was reached after the 60s is comparable in magnitude with the initial flow rate in Alberta Petroleum Industry test APIT3. Assuming that the time taken for the pipe to be emptied is directly proportional to the length of the pipeline and that the flow rate would continue to fall at the same rate as in test APIT3, the time required for the section of the pipe to be emptied is 2.6hours. However, sincethe flow rate of 300kg/s, in this case is reached when the decay rate is much smaller than in the initial flow after the break in test APIT3, a blowdown time which is longer than 2.6hours is expected. Therefore the 3.3hours blowdown time presented by Hardy BBT Ltd is probably correct.

287

CHAPTER 8

RECOMMENDATIONS

FOR FURTHER WORK

In the previous study by Tiley (1989), further work was recommended in four areas namely investigation of the wave speed error, further testing of the model, improvement of the different Through further the the approach model. refinement of stability of solution and by investigation Tiley four has in been the this recommended of areas model, all which used (1989), have been successfully covered. The method of characteristics model which has been developed in this study is sufficiently accurate and reliable for simulating linebreak is bias directional in full-scale The the the model only weakness of situations pipelines. downstream in flow the the the of pipe section of reversal which occurs when modelling the break. This is the only area recommended for further investigation in the model based in for flow In the the the the analysing meantime, programme on method characteristics. the section of the pipeline upstream of the rupture, could be used for the downstream section after the necessary adjustment

of the input data. Alternatively, either the

programmesbasedon the MacCormack or Warming-Kutler-Lomax methods could be used.

Two numericalmodelsin additionto the methodof characteristics, which is the only one covered by Tiley(1989), havebeenused in this study. Both the MacCormack and Warming-Kutler-Lomaxmodels,havesuccessfullybeendeveloped. What remainsto be doneis to test further the models. Specifically,the problemsof overshootand oscillating results, which are associatedwith the MacCormack numerical method, need to be investigated. Another numericalproblemwhich needsto be solved,is that of round-off error, when modelling flow in long pipelinesusing the variablegrid spacingmethodand with very small Ax at the broken end. The QUANT software, in its present version and with all its limitations, provides a sufficiently accurate and optimum means of incorporating real gas properties into the model. However, better results are expected if the current limitation of the software is overcome. The limitations are lack of output for some values of input parameters and for some gases and lack of output at high pressure and temperatures lower than 200K. The former limitation of lack of output at high pressures has already been overcome [Silberring (1993-95)]. But the version of the updated software is not yet available commercially. It is important to follow up on further updates of the software, which would improve the performance of the model.

288

A great deal of time was devoted to exploring, acquiring and developing an optimum method of incorporating the data produced by the software into the numerical programmes.

The method adopted

in this study [refer to Section 4.5] provides a

compromise as far as accuracy and economy in computational resources are concerned. Two areasare recommended for further work in the application of the QUANT software. The first area is in relation to the routine DFGEN for generating the fluid property data file FLUID. DAT. A programme for automatic execution of the QDB programme of the QUANT

software is required, in order to make the process simpler and less time

consuming. A detailed explanation of this procedure is given in Section 6.3. The second area, which is also described in section 6.3, concerns the retrieval of data from the file FLUID. DAT. An interpolation procedure has proved to be impossible due to the limited computer memory, when running the transient analysis programmes. The alternative method i. e. generating fluid property data in which the input parameters (p & T) are so close such that the values which are closest to the parameters for which fluid data is required are used. The latter method is lesspractical with the presently available computing resources. In this study a spacing of 0. IMPa in pressureand 5K in temperature was used. A finer spacing could be used if the domain for the particular analysis is smaller. Further investigation of the possibility of using interpolation procedures is required especially if bigger computing resources are available. The biggest challenge in developing and also applying the computer programmes in this study, has been that of limited computing resources available. It was decided from the beginning(refer to Section 4.6) that the computer codes developed in this study should be capable of being run on personal computers. The best computer which could be used was the best that was availableat the time. But even this requirement could not be met due to limited financial resources which were available for acquiring such computers and also the high speed at which computer technology advances. At first it was thought that a 486DX personal computer would be sufficient. This was true even for the transient analysis including

the variable grid

method before varying fluid property data was

introduced. A PENTIUM P75 is sufficient to run the complete programmes for transient analysis at a speed which is much faster than the one which could be achieved using the 486DX personalcomputer for the simpler version for the transient analysis programme i. e. with constantfluid property data. Better performance is expected if a more powerful pc such as a PENTIUM P166, is used. It is also thought that such a powerful pc's could

1 289

handle a larger number of grid points on one section of pipeline being modelled, than the presentmaximum number of grid points of around 120. This could eliminate the problem of round-ofd' error experienced when modelling long pipelines with the variable grid programmes.It is therefore recommendedthat a more powerful pc than the PENTIUM P75 be used to test the programmes. Also the possibility of incorporating a chip similar to the Microway-i860

[refer to Section 4.6], into the present series of personal computers in

order to increase their computational speed is worth exploring.

A detailed presentationof experimentaldata which exist for validating linebreak modelsis given in ChapterS. Sufficient experimentaldata,which covers both short and long pipelinessectionswereacquiredandused in this study. However, it was not possible to get datawhicharerecentandconsequentlyrather old data were used. The data which areusedin this studyweregathered using testswhich were conducted during the past 16 to 20 years. It is expectedthat in such a long time, there would be a considerable improvementin instrumentationand experimentaltechniqueswhich could result in more accurate experimentaldata. Someof the tests reported in Chapter5, were carried out within the pasttwo years.Data from suchtestsarethereforeexpectedto be more accurate than the one usedin this study and should if possiblebe used. Unfortunately,it was not possibleto haveaccessto this data becauseof commercialsensitivity. The model developedin this study assumesthat the fluid is a homogeneousgas mixture and that the flow is one-dimensional. These assumptions are adequate for practical linebreak modelling. Better results could be obtained if the simplifications of homogenous and one-dimensionalflow are not used. In terms of priority, investigation into these aspects should come as the last step. In order of preference multiphase effects including frictional , forces are probably the most significant, and should be investigated first. The gamdeleps method [Flatt and Trichet (1995)] which is briefly describedin section 2.7, is recommended.

290

CHAPTER 9

CONCLUSIONS After a critical review of the previousstudy by Tiley (1989), a theoreticalmodel hasbeen developedfor analysisof the transientflow following linebreakin high-pressurenaturalgas in has been improvement An made manyaspectsof the model which was used pipelines. by Tiley (1989). Theseincludederivation of the basicequationsof flow and subsequent force, frictional thermodynamicand transport propertiesof simplifications;calculationof the fluid, heattransferto the fluid and numericalmethodof solution. The basicequations developed by flow delta Flatt (1989). The based the methodwhich was of are on gamma threepartialdifferentialequationsof flow are derived for unsteadyquasi-one-dimensional flow of a realgasthrougha non-rigid non-constantcross-sectionalareapipe. Two further in before basic flow them the the simplifications are made on equations of applying developmentof the computermodel for analysisof the flow following a linebreak. The QUANT software for thermodynamicand transport propertiesof the fluids is used. The flow dependentexplicitequationwhichwasdevelopedby Chen(1979) is used to calculate frictional force. The heattransferthrough the pipe is calculatedusing a formula which is basedon the adiabaticwall temperatureand recoveryfactor. The heattransferis also flow dependent;and the calculationprocedureincludesboth pipesexposedto the atmosphere andburiedpipes. An improvedandmoreaccurateequationfor calculatingthe initial break conditions at the break is used. This improved equation avoids the problem of underestimatingthe flow velocity at the outlet end and overestimatingthe temperature drop at the break. A non-uniform grid spacingwhich is very similar to that usedby Tiley (1989) is used, in order to be able to handle long pipelines, to produce stable results and to also adequately model the physical behaviour of the gas, following a rupture. This method involves the use of a fine grid spacing in the vicinity of the break and a coarser spacing as one moves away from the break. A possibility of modelling the flow reversal in the section of the pipeline downstream of the break, similar to the one developed by Tiley (1989) is provided. Four different models for calculating the initial steady-state conditions, before the break have beendeveloped. The four models, all of which are viscous flow models are based on the assumptions of incompressible flow and compressible isothermal, adiabatic or non-isothermal non-adiabaticflow. Three different numerical methods for solution of 291

the theoretical models have beendeveloped. The three numerical methods are a first-order backward method, a first-order forward method and a second-ordermethod. The numerical and theoretical models have beencomparedwith predictions made for the proposed Songo Songo Dar es Salaam Natural Gas Pipeline [Hardy BBT Ltd. (1989)].

The theoreticaltransientanalysismodelis developedinto a pc basedcomputercode usingthe C programminglanguageandthree different numericalmethodsof solution. The three numerical methodsare the method of characteristics,which was used by Tiley (1989); the two-step second-orderexplicit finite-differencemethodof MacCormack;and the third-orderWarming-Kutler-Lomax explicitfinite-differencemethod. Predictionsusing the three numericalmethodsof solution havebeencomparedwith data obtainedthrough both full-scalepipelineand laboratoryexperiments.The test sectionsvary from 6.1m to 11.8kmin lengthand 0.1 to 1.4min diameter. The transient analysismodels,basedon the method of characteristics, produce results which are in agreementwith experimentaldata. Also the steady-stateanalysis modelspredictionare in good agreementwith the predictionsmadefor the SongoSongo Dar es Salaampipeline. A PENTIUM P75 is just adequateto run the transientanalysis programmeswhich are basedon the method of characteristicsand the second-order MacCormackmethod.The programmes basedon the third-order Warming-Kutler-Lomax methodrequirea computerwith a greatercapacity. They could be run on the PENTIUM P75pc available,when constantfluid propertieswere used. The Warming-Kutler-Lomax method producedgood results, when it was used in cooperationwith the method of characteristicsin order to smoothenthe rapid variationsthe flow. The predictionswhich are madeusingthe MacCormackmethodare not satisfactory(refer to section6.3). The QUANT software for thermodynamic and transport properties of fluids containsenough substancesto cover all natural gas mixtures at pressure of up to around 6MPa and temperatures of 200K and above. In most high-pressure natural gas pipelines the pressuresare much higher than 6MPa and during a linebreak the temperature of the gas near the break falls below the 200K lower limit of QUANT.

Additional data were

provided at high pressures,for the particular natural gas mixture used for validation. The present version of QUANT together with the additional data have successfully been incorporatedwith the numericalprogrammes. An updated version of the QUANT which covers the high-pressure zone was not available in time for this study. The optimum way of incorporating the QUANT software into the numerical programmes is by producing 292

data files for each particular gas mixture and range of parameter and using an interpolation procedure to obtain the fluid properties at each required state of the gas. Due to limitations in computer memory the interpolation procedure could not be used and instead the nearest data are used.

Sufficient and suitableexperimentaldata were securedand used to validate the model. Both full-scalepipelineexperimentaldata and data obtainedfrom experimentson shortpipelinesectionsof up to 243mhavebeenused. All the data usedwere obtainedfrom experimentscarried out within the past 16 to 26 years. There is plenty of more recent experimentaldata but its availability for application in this study has proved to be too difficult. It was stated in Section 1.3 that one of the objectives of this study was to investigatefurther the linebreakphenomenon,especiallyat the brokenboundary in order to understand it better. The other objective was to develop a computer model for analysisof linebreakin high-pressurenatural gas pipeline,basedon the previousmodel which was developedby Tiley (1989). Both these objectiveshave been implemented satisfactorily. The theoreticaldecompressiona-p curve Fig 4.3 at the brokenboundary, hasbeenconfirmedby computerthe predictions. The computer model, which is based on the numerical method of characteristics, is very stablenumerically,unlike the Tiley (1989) and Picard and Bishnoi (1989) model and also does not suffer singularity problems such as those experienced by Flatt (1986). Predictions from this model compare better with experimental results than the Tiley's model and the wave speedsare correctly estimated. The model developed in this study contains the additional feature of being able to model heat transfer for caseswhere the pipeline is buried under water, ground or any other medium whose thermal conductivity, and also the depth of the pipe in the medium are known. In addition to the method of characteristicsmodel, computer models based on two other numerical methods have been developed. These two methods were not used by Tiley (1989). The model based on the MacCormack method does not seemto be suitablefor high-pressure linebreak applications. The Warming-Kutter-Lomax method produces results which are comparable with the method of characteristicsprediction. However, the former methods dependstoo much on the latter method, in order to be able to produces those results. Further testing of the models is necessary before making a definite judgement on the suitability and merits of each of these models. 293

It has been confirmed in this study that for typical high-pressure gas pipeline ruptures, the lowest gas temperature at the break could sometimes drop below the 200K minimum limit of the QUANT software. This is contrary to the argument by Silberring (1993-95) and Flatt (1993-96) that the temperature can never drop to such low values. The model based on the MacCormack method predicts lower temperature drops than the method of characteristicsmodel. The former model also predicts higher flow velocities than the latter.

The comparisonwhichwasmadebetweenthe first-order methods and the secondorder methodfor solutionof the steady-stateflow equationsrevealsthat the second-order methodproducesresultswhich aresignificantlydifferent and more accuratethan the firstorder methods. In this study the steady-state analysisresults are usedjust as initial estimates,and are later improved by a transient analysis. Therefore the first-order backwardmethodis sufficient.In caseswhere more accuratesteady-stateanalysisresults arerequired,the second-order methodshouldbe used. Comparisonbetweenthe first-order and the second-ordermethod of characteristicsresults indicate that the second-order methodis significantlymoreaccuratethanthe first-ordermethod,only when the flow in the vicinity of thebreakis considered.However, the second-order methodfails to handlethe rapidvariationsof fluid propertiesat the initial At's nearthe break. In suchsituationsthe first-order methodis used. A major weakness,which is yet to be explained, was discovered when the method of characteristics was used to model the flow reversal on the downstream section of the broken pipeline (negativeflow velocity and x increasing away from the break). The model calculates accurately the fluid properties and the positions on the t-x plane for the intersections of the characteristics curves with the time level t line (Refer to Fig. 4.4). However it producesresults which always tend to be biased to the values at the intact end. This directional bias has the overall effect of predicting a decompression process which is extremelyslow. No such problems were experienced with the other numerical methods. The computer models for transient analysis have been validated with experimental data involving linebreak of pipelines which are as long as 11.8km. It has also been used to simulate the flow following a rupture in a pipeline which is 68km long. Due to the big differencebetweenthe grid spacingsat the broken and the intact ends, the round off error of the computer resultsin a failure of the variable grid method to function properly. This error does not result in failure of the programme to run. But the model produces results which are wrong at some grid points. The net effect is to underestimate the flow velocity. For shorter pipe sections, such a problem does not exist. 294

Sincethe methodof characteristicsis iterative while the explicit finite-difference method are not, it is expectedthat the computation time of the methodof characteristics would be much longer than those of the explicit finite-differencemethods. In this study the methodof characteristic and explicit finite-difference have computationspeedswhich do not differ much.The reasonsfor this situation are the extremelyquick convergenceof the method of characteristicsmodel and comparatively long time taken to retrieve the variable thermodynamicand transport fluid data from the data file produced using the QUANT software. The main sourcesof error in the model prediction and validation results include round off error of the computer and its related errors especiallythose causedon the variablegrid method, whichwere explainedearlier. More errors result from interpolation and also conversion of experimentaldata from graphicalinto numericalform. The latter error hasbeen minimisedby usingthe scanningtechniqueand graphicsAutocad software. The other errors are causedby the use of estimatedfluid properties,and some other parameterssuchasgas composition,initial temperatureetc. Estimatedvalueswere used becauseof lack of output from the QUANT software and also becausesomeparameters necessaryfor the computermodel are not provided with experimentaldata. The previous model by They (1989) had two major weaknesses,namely have the wavespeedsandinstabilityif the solution. Both theseweaknesses overestimating beeneliminatedin this study. The methodof characteristics,which was also usedby Tiley has proved to be the most suitablenumericalof solution for analysisof transient flow following linebreakin high-pressuregas pipelines. The first- and second-ordermethods produce results that are very close together, with the second-ordermethod producing slightlylowerwavespeedsthanthe first-ordermethod. The first-order methodis over two times fasterthan the second-ordermethod in computation speedand in somecasesthe formermethodhandlestheboundaryconditionsbetterthanthe latter method. At positions which are further away from the break, the first order method seemsto producebetter resultsthanthe second-ordermethod. The first order methodis thereforepreferredif the resultsare requiredfor positionswhich are further away from the break, and the analysis involveslongpipelines.The first-orderapproximationhasprovedto a very important tool in modellingof transientflow following linebreakin high-pressuregaspipelines. It is the only numericalmethod,out of thosestudied,which can handleall the initial gas conditions, flow variation andgrid sizesbest. It producessufficientlyaccurateresults,with the least 295

computational resources. It therefore provides a practical method for simulating the flow in long pipelines. Except for the problem of directional bias in the second-order method of characteristics when modelling the flow in the section of the pipeline downstream of the break, this study has produced an accurate computer programme, based on the method of characteristics, for analysis of the flow following a break in high-pressure gas pipelines.

The MacCormack second-ordermethod was previously thought of as being potentiallybetterfor linebreakproblemsthanthe methodof characteristics.In contrast,this studyhasconfirmedthat the MacCormackmethodis completelyunsuitablefor modelling of the flow following linebreakin high-pressuregaspipelines. The third-order WarmingKutler-Lomaxmethodproducesresultswhich are comparablewith those producedby the secondordermethodof characteristics.However, it dependstoo much on the methodof characteristics,to smoothenthe rapid flow following a break. This makesone wonder whetherthe resultsareproducedby the Warming-Kutler-Lomaxmethodor the methodof characteristics.Due to this over-dependence on the methodof characteristics,it was not possibleto determinethe execution speedof the Warming-Kutler-Lomax method. It requiresa lot more computermemorythan the second-ordermethodof characteristics.

296

REFERENCES AND BIBLIOGRAPHY Adkins, C. J. (1983), "Equilibrium Thermodynamics", 3`d Edition, Cambridge University Press.

Alberta PetroleumIndustry GovernmentEnvironmentalCommittee (1979), "Hydrogen SulphideIsopleth Prediction-PhaseII: Pipe Burst Study". Alberta PetroleumIndustry GovernmentEnvironmentalCommittee (1978), "Hydrogen SulphideIsopleth Prediction-PhaseI: Model SensitivityStudy". Allaire,P. E. (1985),"Basicsof the FiniteElementMethod: SolidMechanics,Heat Transfer and Fluid Mechanics",Wm. C. Brown Publishers,Dubuque,Iowa. Ames, W. F. (1977), "Numerical Methods for Partial Differential Equations", Second Edition, Academic Press, New York. Ansari, J. S. (1972), "Influence of Including Radial Flow on the Solution of Unsteady Pipe Flow Equations", Transactions of the ASME, Paper No. 72-FE-28. Arrison, N. L., Hancox, W. T., Sulatisky, M. T. (1977), "Blowdown of a Recirculating Loop With Heat Addition", IMECHE, C202/77, pp. 77-82.

Bannister,F. K. and Mucklow, G. F. (1948), "Wave Action Following SuddenReleaseof CompressedGasfro a Cylinder", Proc. IMECHE, 159,pp. 269. Batchelor, G. K. (1992), "An Introduction to Fluid Dynamics", Cambridge University Press. Beauchemin,P. and Marche, C. (1992), "Transients in Complex Closed-Conduit Systems: Second-Order Explicit Finite Difference Schemes", in Bettess, R. and Watts, J. (Ed. ), "Unsteady Flow and Fluid Transients", A. A. Balkema, Roterdam, pp. 31-39.

Bhallamudi,M. S. and Chaudhry,M. H. (1990), "Analysisof Transientsin Homogeneous Gas-LiquidMixtures", in Thorley, A. R. D. (Ed.), "PressureSurges: Proceedingsof the 6"' InternationalConference",BHRA, pp. 343-348. Bisgaard, C. Sorensen, H. H. and Spangenberg, S. (1987), "A Finite Element Method for Transient Compressible Flow in Pipelines", Int. J. Num. Meth. Fluids, 7, pp. 291-303. Bornstein, R. and Roth, W. A. (1920), "Landolt-Bornstein Physikalisch-Chemische Tabellen", Verlag von Julius Springer, Berlin. Botros, K. K., Jungowski, W. M. and Weiss, M. H. (1989), "Models and Methods of Simulating Gas Pipeline Blowdown", Can. J. Chem. Eng., 67, pp. 529-539. Boulos, P. F., Wood, D. J. and Funk, J. E. (1990), "A Comparison Numerical of and Exact Solution of Pressure Surge Analysis", in Thorley, A. R. D., (Ed. ), "Pressure SurgesProceedings of the 6th International Conference", BHRA, pp. 149-159.

297

Bratland, 0. and Jensen,P. (1990), "Simulation Models of Pressure-Steps in Very Long Hoses or Pipes With Visco-Elastic Walls", in Thorley, A. R. D. (Ed. ), "Pressure Surges: Proceedings of the 6" International Conference, BHRA, pp. 59-68.

Brown, F. T. (1969)"A QuasiMethod of CharacteristicsWith Application to Fluid Lines With Frequency-Dependent Wall ShearandHeatTransfer",Journalof Basic of Engineering Transactionsof the ASME SeriesD, 91(2), pp. 217-227. Budny,D. D., Wiggert,D. C. andHatfield,F. J. (1990),"Energy Dissipationin the AxiallyCoupled Model for Transient Flow", in Thorley, A. R. D. (Ed.), "PressureSurgesProceedingsof the 6' InternationalConference",BHRA, pp. 15-26. Chen, N. H. (1979), "An Explicit Equation for Friction Factor in Pipe", Ind. Eng. Chem. Fund. 18(3), pp. 296-297. Chen, J. R., Richardson,S. M. and Saville, G. (1992), "Numerical Simulation of Full-Bore Ruptures of Pipelines Containing Perfect Gases", Trans. IChemE, 70(B), pp. 59-69. Chen, J. R., Richardson,S. M. and Saville, G. (1993), "A Simplified Numerical Method for Two-Phase Pipe Flow", Trans. IChemE, 71(A), pp. 304-306.

Cheng, L. C. and Bowyer, J. M. (1978), "Tube-TransientCompressibleFlow Code", in Papadakis, C. and Scarton, H. (Eds.), "Fluid Transientsand Acoustics in the Power Industry", Proceedingsof the Winter Annual Meeting of the American Society of MechanicalEngineers,pp. 31-35. Courant, R. Isaacson,E. and Rees, M. (1952), "Communs. Pure Appl. Maths., 5, pp. 243. Courant, R. and Friedrichs, K. 0. (1948), "Supersonic Flow and Shock-Waves", Interscience Publishers. Cronje, J. S., Bishnoi, P. R. and Svrcek, W. K. (1980), "The Application of the CharacteristicMethod to Shock Tube Data that Simulate a Gas Pipeline Rupture", Can. J. Chem. Eng., 58, pp. 289-294. De Bakker, A. G. (1993), GASUNIE, The Netherlands, "Private Communication". Dzung, L. S. (1944), "Beitrage zur Thermodynamik Realer Gase", Schweiz, Archiv Bd. 10, No. 11, pp. 305-313 (German). Earnshaw,S. (1860), "On the MathematicalTheory of Sound", Phil. Trans., Royal Society, 150, pp. 133. Eichinger, P. and Lein, G. (1992), "The Influence of Friction on Unsteady Pipe Flow", in Bettess,J. and Watts, J., "Unsteady Flow and Fluid Transients", Balkema, Roterdam, pp. 41-50.

Fan,D. andTijsseling,A. (1990),"Fluid StructureInteractionWith Cavitationin Transient PipeFlow", ASME J. Fl. Eng., 114(2),pp. 268-274. 298

Fannelop,T. K. and Ryhming, I. L. (1982), "Massive Releaseof Gas From Long Pipelines", Journal of Energy, 6(2), pp. 132-140. Fashbaugh,R. H. and Widawsky, A. (1972), "On the Difference Methods for Shock Wave Propagation in Variable Ducts Including Wall Friction", Naval Civil Engineering Laboratory, Port Heneme, California/ASME.

Feuer,A. (1982),'The PuzzleBook: Puzzlesfor the C ProgrammingLanguage",PrenticeHall. Flatt,R (1985a),"A Singly-IterativeSecondOrderMethod of Characteristicsfor Unsteady Compressible One-DimensionalFlows", Communicationsin Applied NumericalMethods, 1, pp. 269-274. Flatt, R. (1985b), "On the Application of Numerical Methods of Fluid Mechanics to the Dynamics of Real Gases", Forschung im Ingenieurwesen, 51(2), pp. 41-52, in German (English translation: British Gas Translation No. T07823/BG/LRS/LRST817/86). Flatt, R. (1986), "Unsteady Compressible Flow in Long Pipelines Following a Rupture", Int. J. Num. Meth. Fluids, 6, pp. 83-100. Flatt, R. (1989), "Unified Approach of Calculation - of Monophase and Biphase/Homogeneous Flows of Condensable Pure Fluids", Entropie, '149, pp. 48-55 (French).

Flats, R. (1993a), "Derivation of the ConservationEquationsfor UnsteadyQuasi-OneDimensionalFlow of Real GasesThrough Non-Rigid Non-Constant-AreaDucts", Swiss FederalInstitute of Technology,Lausanne(unpublished). Flatt, R. (I993b), "The Basic Relation of the yS Method", SwissFederalInstitute of Technology,Lausanne(unpublished). Flatt, R. (1993-1996), Communication".

Swiss Federal Institute of Technology, Lausanne, "Private

Flatt, R. and Trichet, J. C. (1995), "CFD Oriented Thermodynamics for Single and 2Phase/EquilibriumFlows of Real Gases",EUROMECH Colloquium 331: Flow With Phase Transition, Gottingen. Foothills Pipe Lines (Yukon)Ltd (1981), "Final Report on the Test Program at the Northern Alberta Burst Test Facility", Report No. 18031. Gathmann,IL J., Hebeker, F. K. and Schöffel, S. (1991), "On the Numerical Simulation of Shock Waves in an Annular Crevice and its Implementation on the IBM ES/3090 with Vector Facility", in Brebbia, C. A., Howard, D. and Peters, A. (Eds.), "Application of Supercomputers in Engineering II", Computational Mechanics Publications and Elsevier Applied Science, London, pp. 319-330.

299

Glaister, P. (1994), "An Efficient Shock-CapturingAlgorithm for CompressibleFlows in a Duct of Variable Cross-Section",Int. J. Num. Meth. Fluids, Vol. 18, pp. 107-122. Glass,I. I. (1958),"ShockTubes-PartI: Theory andPerformanceof Simple Shock Tubes". UTIA Review, 2(I). Goldwater M. H. and Fincham, A. E (1981), "Modelling of Gas Supply Systems" in Nicholson,H. (Ed.), "Modelling of Dynamical Systems,Vol. 2", IEE-Control Engineering Series13. Govier, G. W. and Aziz, K. (1972), "Flow of Complex Mixtures in Pipes", van Nostrand Reinhold. Gradle, R. J. (1984) "Design of GasPipelineBlowdowns", Energy Processing,Jan-Feb. Groves, T. K., Bishnoi, P. R. and Wallbridge, J. E. (1978), "Decompression Wave Velocities in Natural Gasesin Pipelines",Can. J. Chem. Eng., 56, pp. 664-668. Groves, T. K. (1976), "MeasuredVelocity of DecompressionWaves in Natural Gas in Pipelines",Foothills Pipelines(Yukon) Ltd. Guy, J. J. (1967), "Computation of Unsteady Gas Flow in Pipe Networks", IChemE SymposiumSeries,23, pp. 139-145. Hall, J. G. (1958),"ShockTubes-PartII: Production of Strong Shock Waves; Shock Tube Applications,Design and Instrumentation",UTIA Review, 2(11). Haque, M. A., Richardson,S. M., Saville, G., Chamberlain,G. and Shirvill, L. (1992), "Slowdown of PressureVessels:H. ExperimentalValidationof Computer Model and Case Studies",Trans. IChemE, 70(B), pp. 10-17. Haque, M. A., Richardson, S. M. and Saville, G. (1992), "Blowdown of Pressure Vessels: 1. Computer Model", Trans. IChemE, Vol. 70(B), pp. 3-9.

Haque, M. A., Richardson, S. M., Saville, G. and Chamberlain, G. (1990), "Rapid Depressurizationof PressureVessels",J. Loss Prev. ProcessInd., Vol. 3, pp. 4-7. Hardy BBT Ltd. (1989), "Kilwa Kivinje-Dar es SalaamGas Pipeline Feasibility Study", TanzaniaPetroleumDevelopmentCorporation, Project No. CG10305 (R). Hartree, D. R. (1964), "Numerical Analysis", SecondEdition, Oxford University Press. Hirose,M. (1971),"Frequency-Dependent Wall Shearin Transient Fluid Flow: Simulation of UnsteadyTurbulent Flow", MSc. Dissertation,MassachusettsInstitute of Technology. Hoffmann, K. A. (1989), "Computational Fluid Dynamics for Engineers", Engineering Education Systems(EES), Texas. Holman, J. P. (1980), "Thermodynamic",Third Edition, McGraw-Hill, New York. 300

Holt, M. (1984), "Numerical Methods in Fluid Dynamics", 2ndEdition, Springer-Verlag, Berlin. Huang, W. and Sloan, D. M. (1993), "A New PseudospectralMethod With Upwind Features",IMA Journal of Numerical Analysis, 13, pp. 413-430. Ilic, V. (1987), "Two-Phase Blowdown Through a Short Tube", in Dodge, F. T. and Moody, F. J., "Fluid Transientsin Fluid-StructureInteraction-1987", FED, 56, ASME, pp. 87-92. IPS (1995.a), Dar es Salaam,20'hJuly. IPS (1995.b), Dar es Salaam,27thSept. UnsteadyCompressibleFlow With Friction and Heat Issa,R. I. (1970), "One-Dimensional Transfer", MSc Dissertation,Imperial College of Scienceand Technology, University of London. Jones,D. G. andGough,D. W. (1981),"Rich GasDecompressionBehaviour in Pipelines", AGA-EPRG Linepipe ResearchSeminarIV, Duisburg, British Gas, Report No. ERS E 293. Kagawa, T., et al. (1983) " High Speedand Accurate Computing Method of FrequencyDependentFriction in LaminarPipeFlow for CharacteristicsMethod", JSME, 49(447), pp. 26- 38 (Japanese). Kawabe, R. (1982), "An Analytical Method for Thermal-Hydraulic Transientsin Piping Networks", Nucl. Eng. and Design, 73, pp. 441-446. Kay, W. B. (1936), "Densityof HydrocarbonGasesand Vapours", Ind. Eng. Chem.,28(9), pp. 1014-1019. Kernighan,B. W. and Ritche, D. (1978), "The C ProgrammingLanguage",Prentice-Hall. Kernighan,B. W. and Pike, R. (1984), "The Unix ProgrammingEnvironment", PrenticeHall. Kimambo,C. Z. M. (1987),"Energyin Tanzania",PaperPresentedat a SeminarOrganised by DevelopmentEducation Centre,Energy for All Project, Cumbria, Feb. 1987. Kimambo,C. Z. M. andThorley,A. (1995), "Computer Modelling of TransientsRuptured High PressureNatural GasPipelines:A Review of Experimentaland Numerical Studies", C502/003, IMECHE. Kiuchi, T. (1993), "An Implicit Method for Transient Flow in Pipe Networks", Paper Submittedfor Publication.

301

Knox, H. W., Atwell, L. D., Angle, R., Willoughby, R. and Dielwart, J. (1980), "Field Verification of Assumptions for Modelling Sour Gas Pipelines and Well Blowouts". PetroleumSociety of CIM, PaperNo. 80-31-46. Kunsch, J. P., Sjoen, K. and Fannelop,T. K. (1991), "Loss Control by Subsea Barrier Valve for Offshore Gas Pipeline", 23'dAnnual Offshore Technology Conference, Houston, Texas, OTC 6616.

Kunsch, J. P., Sjoen, K. and Fannelop, T. K. (1995), "Simulation of Unsteady Flow, Massive Releasesand Blowdown of Long, High-PressurePipelines", 1995 International Gas ResearchConference,Cannes,France. Kuo, S. S. (1972), "Computer Applications of Numerical Methods", Addison-Wesley, Reading. Lang, E. and Fannelop, T. K. (1987), "Efficient Computation of the Pipeline Break Problem", in Dodge, F. T. and Moody, F. J. (Eds.), "Fluid Transient in Fluid-Structure Interaction", FED 56, ASME, pp. 115-123. Lang, E. (1991), "Gas Flow in Pipelines Following a Rupture Computed by a Spectral Method", J. Appl. Maths. and Physics (ZAMP), 42, pp. 183-197.

Lapidus,L. andPinder,G. F. (1982), "Numerical Solution of Partial Differential Equations in Scienceand Engineering",John Wiley and Sons,New York. Lavooij, C. S. W and Tijsseling, A. S. (1990), "Fluid Structure Interaction in Compliant Piping Systems", in Thorley, A. R. D. (Ed), "Pressure Surges-Proceedingsof the 6" International Conference,BHRA, pp.85-100. Lax, P. and Wendroff, B. (1960), "Systemsof ConservationLaws", Comm. Pure Appl. Math., 13, pp. 217-237. Lister,M. (1960),"The NumericalSolutionof Hyperbolic Partial Differential Equationsby the Method of Characteristics",in Ralston, A. and Wilf, H. S. (Eds.), "Mathematical Methods for Digital Computers",John Wiley & Sons,pp. 165-179. Lyczkowski,R W., Grimesey,R A. andSolbrig,C. W. (1978), "Pipe Blowdown Analyses Using Explicit Numerical Schemes",AICHE SymposiumSeries,74 (14), pp. 129-140. MacCormack,R. (1971), "Numerical Solution of the Interaction of a Shock Wave with a LaminarBoundary Layer", in Holt, M. (Ed.), Lecture Notes in Physics;Vol. 8, SpringerVerlag, New York, pp. 151-163. MacCormick,S. F. [Ed.] (1987), "Multigrid Methods: Frontiers in Applied Mathematics", Society for Industrial and Applied Mathematics,Philadelphia,Pennsylvania. Mallinson, J. (1996), British Gas, Gas Research Centre, Loughborough, "Private Communication".

302

Martin, D. J. (1967), "Equations of State", Ind. Eng. Chem., 59(12), pp. 34-52.

Mathers, W. G., Zuzak, W. W., McDonald, B. H. and Hancox, W. T. (1976), "On Finite Difference Solutions to the TransientFlow-Boiling Equations", Committee on Safety of NuclearInstallationsSpecialistsMeetingon Transient Two-PhaseFlow, Toronto, Ontario, PergamonPress,pp. 278-315. Meissner, L. P. and Organick, E. I. (1980), "FORTRAN 77: Featuring Structured Programming",Addison-WesleyPublishingCompany,Reading. Ministry of Water Energy and Minerals, Ocelot TanzaniaInc. and TransCanadaPipelines Ltd. (1994), "Songo Songo Gas DevelopmentProject", Brochure. Ministry of WaterEnergy and Minerals (1990), "Energy Balance:Tanzania, 1990" Dar es Salaam. Moe, R and Bendiksen, K. H. (1993), "Transient Simulation of 2D and 3D Stratified and Intermittent Two-Phase Flows. PART I: Theory", Int. J. Num. Meth. Fluids, 16, pp. 461487.

Monro, D. M. (1985), "FORTRAN 77", Edward Arnold Ltd, London. Morrow, T. (1996), Southwest Research Institute, San Antonio, Texas, "Private Communication". Münster, A. (1969), "Statistical Thermodynamics: Volume I", Springer-Verlag,Berlin. Niessner, H. (1980), "Comparisonof Different Numerical Methods for Calculating One Dimensional Unsteady Flows", Lecture Series 1980-1, Lecture No. 16, von Karman Institute for Fluid Dynamics,Sweden. Ohmi, M., Kyomen, S., and Usui, T. (1985), "Numerical Analysisof TransientTurbulent Flow in a Liquid Line", Bulletin of JSME, Vol. 28, No. 239, pp. 799-806. Oliemans, R. V. A. (1987), "Modelling of Two-Phase Flow of Gas and Condensatein Horizontal and Inclined Pipes", Paper for Presentationat the 6d ASME Annual Pipeline EngineeringSymposium,Dallas. Olorunmaiye,J. A. andImide,N. E. (1993), "Computationof Natural Gas Pipeline Rupture ProblemsUsing the Method of Characteristics",J. HazardousMaterials, 34, pp. 81-98. Osiadacz,A. J. andYedroudj,M. (1989), "A Comparisonof a Finite ElementMethod and a Finite Difference Method for Transient Simulation of a Gas Pipeline", Appl. Math. Modelling, 13, pp 79-85. Osiadacz,A. J (1987), "Simulation of Analysisof GasNetworks", Spon, London. Otwell, R. S., Wiggert,D. C. andHatfield, F. J. (1985), "The Effect of Elbow Restraint on PressureTransients",Journal of Fluids Engineering,Vol. 107, No. 3, ASME. 303

Padmanabhan,M., Ames, W. F. and Martin, C. S. (1978), "Numerical Analysis of Pressure Transients in Bubbly Two-Phase Mixtures by Explicit-Implicit Methods", J. Eng. Maths., 12(1), pp. 83-93.

Picard, D. J. and Bishnoi, P. R. (1988), "The Importance of Real-Fluid Behaviour and Nonisentripic Effects in Modelling DecompressionCharacteristicsof Pipeline Fluids for Application in Ductile FracturePropagationAnalysis", Can. J. Chem. Eng., 66, pp. 3-12. Picard,D. J. andBishnoi,P. R (1987);"Calculationof the ThermodynamicSound Velocity in Two-PhaseMulticomponent Fluids", Int. J. Multiphase Flow, 13(3), pp. 295-308. Picard, D. J. and Bishnoi, P. R. (1989), "The Importance of Real-Fluid Behaviour in PredictingReleaseRatesResultingFrom High-PressureSour-GasPipelineRuptures", Can. J. ChemEng., 67, pp. 3-9. Press,W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery,B. P. (1992), "Numerical Recipesin C: The Art of Scientific Computing", CambridgeUniversity Press,2ndEdition. Rachford,H. H. Jr. and Dupont, T. (1974), "A Fast Highly Accurate Meansof Modelling TransientFlow in GasPipelineSystemsby Variational Methods", Soc. Pet. Eng. J., 14, pp. 165-178. Rachford,H. H. Jr. (1972),"TransientPipelineCalculationsImprove Design", Oil and Gas J., 70(44), pp. 54-57. Rachid, F. B. F. and Stuckenbruck, S. (1990), "Transients in Liquid and Structure in Viscoelastic Pipes", in Thorley, A. R. D. (Ed.), "PressureSurges:Proceedingsof the 61 InternationalConference,BHRA, pp. 69-84. Ralston,A. and Wilf, H. S. [Ed.] (1964), "MathematicalMethods for Digital computers", John Wiley & Sons,Inc., New York. Richardson,S. M. andSaville,G. (1991), "Blowdown of Pipelines",Society of Petroleum Engineers,SPE23070,pp. 369-377. Richardson, S. M. (1993-1996), Department of Chemical Engineering and Chemical Technology,Imperial College,London, "Private Communication". Richtmyer, R. D. (1957), "Difference Methods for Initial-Value Problems", Interscience Tracts in Pure and Applied Mathematics,Number 4, IntersciencePublishers,Inc., New York, ChapterX. Rowlinson,J. S. andWatson, I. D. (1969), "The Prediction of ThermodynamicProperties of FluidsandFluid mixtures-I: The Principle of CorrespondingStatesand its Extensions", Chem.Eng. Sc., 24, pp. 1565-1574. Ryhming I. L. (1987), "On the Expansion Wave Problem in a Long Pipe with Wall Friction", ZAMP, 38.

304

Saad, M. A., (1985) "Compressible Fluid Flow", Prentice-Hall, Inc., Englewood Cliffs, New Jersey.

Saville, G. and Szczepanski,R. (1982), "Methane-Based Equations of State for a CorrespondingStatesReferenceSubstance",Chem. Eng. Sc., 37(5), pp. 719-725. Schildt, H. (1989), "ANSI C Made Easy", McGraw-Hill. Schmidt, F. W., Henderson, R. E. and Wolgemuth, C. H. (1984), "Introduction to Thermal Sciences: Thermodynamics, Fluid Dynamics, Heat Transfer", John Wiley & Sons. Scientific Software-Intercomp (1990), "Songo Songo Reservoir Study", Tanzania Petroleum Development Corporation.

Sens, M., Jouve, P. and Pelletier, R. (1970), "Detection of a Break in a Gas Line", 110' International Gas Conference,Moscow, Paper No. IGU/C37-70, in French, (English Translation:British Gas, TranslationNo. BG. TRANS 4910). Shin, Y. W. (1978), "Two-Dimensional Fluid Transient Analysis by the Method of in PapadakisC. and ScartonH. (Ed.), "Fluid Transientsand Acoustics in Characteristics", the Power Industry", ASME, New York, pp. 179-185. Silberring, L. (1990), "Thermodynamicand Transport Propertiesof Gasesfor Computer Aided Engineering", Paper Presentedat the 21" International CODATA Conference, Columbus,Ohio. Silberring,L. (1991), "A TemperatureFunction for the Virial Coefficients", Entropie, No. 164/165,pp. 139-140. Silberring, L. (1991), "A Computer Program Systemfor Thermodynamicand Transport Propertiesof Real Gases",SilberringEngineeringLtd., Zurich. Silberring,L. (1993-1996),SilberringEngineeringLtd., Zurich, "Private Communication". Silberring,L. (1995), "Thermodynamicand Transport Property Data for a Typical Natural GasMixture at High Pressures",Suppliedin a Diskette on Request(Unpublished). Silvester,P. P. (1984), "The Unix SystemGuidebook: An Introductory Guide for Serious Users", Springer-Verlag. Stittgen,M. andZielke,W. (1990), "Fluid Structure Interaction in Flexible Curved Pipes", in Thorley, A. K D. (Ed.), "Pressure Surges-Proceedingsof the 6' International Conference,BHRA, pp. 101-120. Straty,G. C. (1974), "Velocity of Soundin DenseFluid Methane", Cryogenics,July 1974. Suwan, K. and Anderson, A. (1992), "Method of Lines Applied to Hyperbolic Fluid TransientEquations", Int. J. Num. Meth. Eng., 33, pp. 1501-1511.

305

Swaffield, J. A. (1967), "A Study of the Influence of Bends on Fluid Transients in Pipelines",MPhil Thesis,City University, London. Swaffield, J. A. (1968-69), "The Influence of Bends on Fluid Transients Propagated in IncompressiblePipe Flow", Proc. Inst. Mech. Engrs., Vol. 183, Part 1, No. 29, pp. 603-614. Swamee,P. K. and Jain, A. K. (1976), "Explicit Equations for Pipe-Flow Problems", ASCE J. Hydr. Div., pp. 657-664.

Taylor, D. L. (1992), "Computer-Aided Design", Addinson-WesleyPublishingCompany, Reading. The Guardian(1995), TanzanianNewspaper, 13thOct Thorley, A. R. D. (1991), "Fluid Transientsin PipelineSystems",D. & L. GeorgeLtd. Thorley,A. R. D. andTiley, C. H. (1987), "Unsteadyand TransientFlow of Compressible Fluids in Pipelines:A Review of Theoretical and Some ExperimentalStudies", Int. J. of Heat and Fluid Flow, 8(1), pp. 3-15. Tijsseling, A. S. and Lavooij, C. S. W (1990), "Fluid Structure Interaction and Column Separation in a Straight Elastic Pipe", in Thorley, A. R. D. (Ed.), "PressureSurges: Proceedingsof the 6" InternationalConference",BHRA, pp. 27-42. Tijsseling, A. S. (1993), "Fluid-Structure Interaction in the Caseof WaterhammerWith Cavitation", Communicationson Hydraulic and GeotechnicalEngineeringof the Faculty of Civil Engineering,University of Dundee,Report No. 93-6, ISSN. 0169-6548. Tiley, C. H. (1989), "PressureTransientsin a Ruptured Gas Pipeline With Friction and Thermal Effects Included", PhD Thesis,City University, London. Trikha, A. K. (1975), "An EfficcientMethod for SimulatingFrequencyDependentFriction in TransientLiquid Flow", Tranc. ASME, J. Fluid Eng., 977, pp. 97-105. Tsumura,R andStraty,G. C. (1977),"Speedof Soundin Saturatedand CompressedFluid Ethane", Cryogenics,April 1977. Van Reet, J. D. and Skogman, K. D. (1987), "The Effect of Measurement Uncertainty on Real Time Pipeline Modelling Applications", 1987 ASME Pipeline Engineering Symposium-ETCE, Dallas, Tx. Van Deen, J. K. and Reintsema, S. R. (1983), "Modelling of High-Pressure Gas Transmission Lines", Appl. Math. Modelling, 7, pp. 268-273.

Vardy, A. E.(1993), Faculty of Civil Engineering, University of Dundee, "Private Communication". Vardy, A. E., Hwang, K. L. and Brown, J. M. B. (1993) "AWeighting Function Model of TransientTurbulent Pipe Friction", J. Hyd. Res., 31(4), pp. 533-548. 306

Approximating Unsteady Friction at High Reynolds Numbers", in Vardy, A (1992) ," Bettess, R. and Watts, J. (Eds. ), "Unsteady Flow and Fluid Transients", Balkema, Roterdam, pp. 21-29. Vardy, A. E. and Fan, D. (1990), "Flexural Waves in a Closed Tube", in Thorley, A. R. D. (Ed. ), "PressureSurges: Proceedings of the 6" International Conference", BHRA, pp. 4358.

Vetterling, W. T., Teukosky, S. A., Press,W. H. and Flannery,B. P. (1992), "Numerical Recipes:ExampleBook (C)", 2"aEdition CambridgeUniversity Press. Warming,R F., Kutler, P. and Lomax, H. (1973), "Second-andThird-Order Noncentered DifferenceSchemesfor NonlinearHyperbolicEquations", AIAA Journal, 11, pp. 189-196. Weber SystemsInc. Staff (1984), "C Users' Handbook", Addison-Wesley,Reading. Weimann,A. (1978),"GasDistribution Network Dynamic Modelling and Simulationwith Respectto Network Control and Monitoring", Thesis,Technical University of Munich, in German(English Translation:British Gas Translation No. LRS T 435). Welty, J. R (1978), "EngineeringHeat Transfer", John Wiley & Sons. White F. M. (1988), "Fluid Mechanics",SecondEdition, McGraw-Hill, New York. Wilson,D. J. (1981), "Expansionand Plume Rise of Gas Jets from High PressurePipeline Ruptures", Alberta Environment, Pollution Control Division. Wood, D. J., Dorsch, R. G. and Lightner, C. (1966), "Wave-Plan Analysis of Unsteady Flow in ClosedConduits", ASCE-Journalof Hydraulics Division, 92(HY2), pp. 83-92. Wood, D. J., Funk, J. E. and Boulos, P. F. (1990), "Pipe Network Transients-Distributed and Lumped Parameter Modelling", in Thorley, A. R. D. (Ed.), "Pressure Conference", Surges-Proceedings 6th International BHRA, pp. 131-142. the of Wortman,L. A. and Sidebottom,T. 0. (1984), "The C ProgrammingTutor", Prentice-Hall. Wu, H. L., Pots, B. F. M., Hollenberg, J. F. and Meerhoff, R. (1987), "Two-Phase, Gas/CondensateFlow at High Pressure in an 8-Inch Horizontal Pipe", Paper for Presentationat the 3rdInternationalConferenceon Multiphase Flow, The Hague. Wylie, E. B. and Streeter,V. L. (1978), "Fluid Transients",McGraw-Hill. Wylie, E. B. (1982), "Fluid Transient in Two Spatial Dimensions", Paper for Presentation at the 1982 Pressure Vessel and Piping Conference, ASME, Orlando, Florida.

Yigang,C. andTing-Chao,S. (1990), "An Efficient Approximate Expressionfor Transient Flow of High Viscous Fluid in Hydraulic Pipelines",in Thorley, A. R. D. (Ed.) "Pressure Surges-Proceedings of the 6`hInternationalConference",BHRA, pp. 15-26.

307

Zemansky, M. W. and Dittman, R. H. (1983), "Heat and Thermodynamics",G"'Edition, McGraw-Hill, New York, Chapter 13. Zha, G. C. andBilgen,E. (1993), "NumericalSolution of Euler Equationsby Using a New Flux Vector Splitting Scheme",Int. J. Num. Meth. Fluids, 17, pp. 115-144. Zielke, W. (1968), "Frequency-DependentFriction in Pipe Flow", Journal of Basic Engineering,March, pp. 109-115. Zigrang, D. J. and Sylvester, N. D. (1985), "A Review of Explicit Friction Factor Equations",J. Energy ResourcesTechnology, Trans. ASME, 107(2), pp. 280-283. Zucrow, M. J. and Hoffman, J. D. (1977), "Gas Dynamics:Multidimensional Flow, Vol. 2 ", John Wiley & Sons,Chapter 19, pp. 295-399.

308

APPENDICES EQUATIONS

A GENERAL

Cross-sectionareaof the pipe: (A-I)

Ad

Flow velocity of fluid: A

(A-2)

pA

Friction factor calculatedusing the Chen's explicit equation: 1" fD

c-5.0452 lo g -2.0 3.7065 d

5.8506

IE1.1098 lo g(2.82577 )D

Re

+

Re 0.8911

Frictional force per unit length of the pipe: AdP

u2 fo

(A-4)

Reynold's number: Pud Re µ

(A-5)

Heat transfer acrossthe pipe calculatedby the Stantonnumbermethod: 0=

(A-6)

nd pu CpSt(T,v- T)

Stantonnumber: St .

Nu

h

Pr Re

puCp

St - f(Pr

,

(A-7)

Re)

(A-8)

Staticpressure:

(A-9)

p= Zips Speedof sound/wavespeedfor a perfect gas:

(A-10)

apKKZRT P

andfor a realgas: a"y,

P Y,

ZR

(A-11)

T

P I

309

(A -3 )

Prandtl number Pr

(A-12)

-

k

Nussetl number: Nit

"hd k

Grasshoffnumber: Gr

(A-13)

P2 d3 &T

.ßg

(A-14)

µ2

Mach number: Ma

(A-15)

-u

a

B EXPRESSIONS FOR EQUATION OF STATE The equationof state for a perfect gas is (B-1)

pRT

p"

The general equation of state for a real-gas is P-

(B-2)

ZpRT

Van Der Waal equation of state (1873) p

NRT

NSA

V-NB

V2

(B-3)

Where N Is the number of molts

Dietrici equationof state (1899)

(-. RT A CXP V-B RTV

jr

(B-4)

Berthelotequationof state (1903) P

RT V-B

-A

(B-5)

yy ý

where A and B in equation(B-3) to (B-5) are constants. Redlich-Kwong equation of state (1949) TA PR 0

ZC(VR-0.08664 /Z, )

0.42748 _

$Z. 2 TR . VR(VR. 0.08664 /Zý2

310

(B-6)

Soave-Redlich-Kwong(SRK) equation of state (1972) p.

RT

Jr

ä-b

v"(v".b )

(B-7)

where

RT b"0.08664 PC ä(T) " a(T,,)ä(Trw) R 2Tc 2 J(T) - 0.4247747 PC s)]2 ä(T,. w) " [1"rn(1-T;

1.574W-0.176W2

m-0.480

Peng-Robinson(PR) equationof state (1976) f,

ä(T)

RT -

v-,

_

(B-g)

b). b(v'-b) v"(v". where RT

b"0.077780

c PC

J(7) - ä(Tc)&(Tr.w) R2T2 r ä(TC) " 0.45724 PC ä(T,, w) " [1. n(1-Ti's)]

ni " 0.377464 . 1.54226w-0.26992 w2

Martin equationof state(1949) A2'B2T. C2eý

RT p.

v-b

...

(v-b)2

As'b5T.Cse-xr (v-b)s

A3'B3T'C3e

A4"B4T"C4e

(v-b)4

(v-b)3

A6.B6T.C6e-kr A7.B7T.C7e.kr ov e

(B-9)

e -2ov

where a, b, K, A, to A,, B, to B7 and Cl to C7are constants. The two parameter thermal equation of state of Martin correspondingstatesis Pf "

T,

9(4-T)

Z, ' v, 0.085

64(Z, v, 0.04)2

311

based on the principle of

where P, .LT, pc

"T

VP "

T`

PQ

Ze " 0.28366

p

Note Z. is the pseudo-criticalreal gas factor The Martin equation solved explicitly with respect to T, is ] (Z, * v, - 0.085)[36#64(Z." v, 0.04)2P, Tr -

(B-11)

9(Z, * v, 0.085). 64(ZcO yr 0.04)2

Van Reet-Skogmanequation of state (1987) 1-

1r-

z

P

(B-12)

rY where y Is a constant .

Clausiusequationof state (1930) p(v-b)

(B-13)

RT -

Benedict-Webb-Rubin equation of state CO RTB ö Aö RD T2 P-

RT

0

C

to )e (1" ;2 ; 3'. 2

as

-a

;6

;3

b, Co, Aa, B0, c. a, ca are empirical a, where

(B-14)

constants

Redlich-Kwong equationof state RT

P

(B-15)

-a v-b v(v. b)TO*5 where

(64 RTC)

(48R Tc512 b"0.086

a"0.427

PC

PC

Beattie-Bridgemanequation of state (1928) pwhere

A"

RT

(-)

Ao(1-a), vv

(B-16)

(v.

) -v

B"

Bo(1-b),

c"C vT3

Thevaluesof the five constantsA a, Bo,b, arein literaturefor somegases. Kamerlingh-Onnesequationof state (1902), in the virial form: pv -NRT

(1 "B.

C........

v

)

(B-17)

v2

in termsof powersof the pressuremaybeusedasshownbelow; An expression Bp Gip2 NRT pv " " . ".........)

312

(B..18)

Saville-Szczepanskimethane-basedequationsof state pRT . p24 1T . A2T'12 . A3 . AIT

p-

. p3(f6T.

A7T"n.

As:

ASIT2) .

AylT2)

All A121T) p4(A10T . . " ' pSA13 " p6 ,1T.

A15/T2).

' p'cif , /T

. AI=/T2) (A2jT2

" A211T3)

p442/T2.

AdT3)

" P9A19IT2 . P3 e'YP2[ A2IT4). P2(A2JT2 " . " P642/T2

p7A1JT

4) A2/T /T2 Pe42, . "

" A2dT3)

" P1003dT2 . A21/T3 . A3JT4)1

p-

pRT.

p2(61T. B2Tin.

B3.

(A2.19)

B41T. BSIT2

3) 3) B7IT BVT B101T " p3(B6T . . p4(8g . " ps(11T2

. B12/T2

` p6B131T2

` p8B141T ' p9B15/T ' p11B1a/T2 p3C. YP2E(17/T2

'p

4B20/T2.

4) B11/T ' . p2B19/T4

6B21/T2"p oý22/T2. B, /T 4)J p

313

(A2.20)

C EXPRESSIONS FOR FRICTION FACTOR C-1 LAMINAR FLOW:

Re < 2100

Hagen-Poiseuilleequation 64

(C-1)

ID - Re

C-2 TRANSITION ZONE C-2.1 Implicit Equations Colebrooke equation(1938) d Valid up

1"2.0

E"0.01 (Re f)

E+2.5226 log 3.7065 d Re ID

fD

(C-2)

Oliemansexpressionfor two-phaseflow (1976) 1" FD

18.7 2c log * -2.0 deff Reff p

(C-3)

. 1.74

Re is a two-phase Reynolds number and doffis the effective diameter for the two, - phase mixture.

(1938-39)or Prandtl-Colebrooke Colebrooke-White equation 1FD

1"2 fD

log -2.0

E. log d

2.53 3.7d

1.14 -2

(C-4)

Re lv

(I . 9.35d) log eRe fD

C-2.2 Explicit Equations

Chen equation(1979 : All valuesof Re and e/d 1e_5.0452

fD

log - -2.0 3.7065d

Re

log

1.1098 5.8506 E

2.82577 D

Re0.2911

Churchill equation (1977) : For all Re and e/d /1 AD-8 1Rel

12

11

12 3

(A. B) 2

314

(C-6)

(C-5)

where A-

16

1

2.4547 In

) 0.9

C Re

B-

r 0.277 Id

r 37530 1 16 J Re

Wood equation (1966): Re > 10,000and 105< e/d < 0.04 bRe'

AD -a.

where

(C-7)

0.225 E+0.53 dddd

a-0.094

£

(.

0.44 Ec-1.62 88.0

b-

1)0.134

Swamee-Jain equation (1976): 5000 s Re s 10' and 10-6s e/d s 10'2 5.74 1=2 e log (C-8) . 0.9 3.7d Re fD 1.325

AD . In

e

5.74

37d

ReO9

2

Generalised Haaland equation (1983): n=3 recommended for gas transmission lines 9RfE1.11N 1.1.8 C. 9) lo -6 d 3.7d Re ng fD

Shachamequation(1980) 1" ID

14.5

e log -2.0 3.7d

(C-10)

Ra

Simple Zigrang-Sylverster equatior (1982) 1" ID

13 Re

e log -2.0 3.7d

(C-11)

Moody equation(1976) Rc

() fD " 0.0055 1.20.000

6

3

.

(C-12)

Jain equation (1947) 1-1.14 Fd D

E, log -2

11,25

(C-13 )

o. Re 9

Simpleequation of Chen(1979) IE1.11

ID

log " -2.0

7.15

2.549d

Re

315

0.9

(C-14 )

Simple equationof Serghides(1984) t-

12

c. log -2.0 3.7d

fD

Re

(C-15)

I

IntermediateZigrang-Sylversterequation(1982) 1-2.0 jD

13

c. log 3.7d

e-5.02 log 3.7d Re

(C-16)

Re

IntermediateZigrang-Sylversterequations(1985) 1 fD 1.

E-4.518 lo g -2.0 3.7d Re

FD

e 1.14 - 21og d

e-2.51 log -2.0 3.7d Re

lo8

21.25 o"9 Re

1) (C-17)

11p 1_5__r.

6_9

iC"18)

3.7d

Re

Zigrang-Sylversterequationof highestprecision(1982) 1 -"-

ID

C 5.02 e_ e log

Olog -2.

3.7d Re

5.02

13

loge

(C-19)

3.7d Re

3.7d Re

Serghidesequationsof highestprecision(1984) 1 XD

(3.7d Olog _L_. -2.

5.02 E log log __E_1.02 log 3.7d Re 3.7d Re Re

1.02

13 E 3.7d Re

Zigrang-Sylverster equation of highest precision (1985) fD "A-

(CA)2

(C-21)

C-2B. A

and (A - 4.781)2 4.781 B- 2A . 4.781

AD "

.z (C-22)

where A"

B"

C"

12)

log .R -2 3.7d e IE*2.51 log -2 ` 3.7d log -2

E 3.7d

(C-22a) A1 J Re

) 2.51 Bl Re

316

(C-22 b)

(C-22c)

(C-20)

C-2.3 PARTIALLY

DEVELOPED TURBULENCE

Chaudhryequation (1979) : IR,> 2000 (C-23)

JD - 0.046 Re -0.2

Panhandle'A' equation 1-3.39

E Re0.0773

(C-24)

ID for is is E the the the allows an adjustable parameter which systemand efficiencyof where in losses the and variation pipe roughness. effects of minor

Uhl equation(1965) 1" fD

log F 4 -2

fD (C-25) Re

Where F is the drag factor to account for the effect of bendsand fittings. Smith equation(1956) 1"2

log Re

2D

lD

0.3 "

(C-26)

Blasius equation(1911) : Re> 10'

fD - 0.3160 Re ,0'25

(C-27)

C-2.4 FULLY DEVELOPED TURBULENCE Smith equation(1956) 1"2

log

FD

3.7d E

2.2773 .

(C-28)

Nikuradse equation(1932) : Smooth pipe and 3000 < Re < 3.4x 106 1-2.0 tog(RelD . 0.8 (C-29) AD Von Karman equation (1932): Rough pipe and (d/e)/(Ref1) > 0.01 1-2 E. logt 1.74 (C-30)

lo

`

317

D LIST OF PROGRAMMES AND SUB-ROUTINES BREAK

driver routine for transientanalysisafter the break

BRINC

driver routine for initiating the break and calculatinginitial fluid properties at the break

BRINCU

both in initial break for upstream and conditions calculating routine downstreamsectionswith uniform grid size.

BRINCV

in both initial break for upstream and conditions calculating routine

downstream sectionswith variablegrid size. BRMCCAUD transientanalysisdownstreamafter the breakwith uniform grid size and pipe exposedto the atmosphere BRMCCAUU transientanalysisupstreamafter the break with uniform grid size and pipe exposedto the atmosphere BRMCCAVD transientanalysisdownstreamafter the breakwith variablegrid size and pipe exposedto the atmosphere BRMCCAVU transientanalysisupstreamafter the break with variablegrid size and pipe

exposedto the atmosphere break downstream BRMOCAUDtransientanalysis the with uniformgrid sizeandpipe after exposedto the atmosphere BRMOCAUU transientanalysisupstreamafter the break with uniform grid size and pipe

exposedto the atmosphere BRMOCAVD transientanalysisdownstreamafter the breakwith variablegrid size and pipe exposedto the atmosphere

BRMOCAVUtransientanalysisupstreamafterthe breakwith variablegrid sizeandpipe exposedto the atmosphere BRWKLAUD transientanalysisdownstreamafter the breakwith uniform grid size and pipe exposedto the atmosphere BRWKLAUU transientanalysisupstreamafter the break with uniform grid size and pipe exposedto the atmosphere BRWKLAVD transientanalysisdownstreamafter the breakwith variablegrid size and pipe exposedto the atmosphere BRWKLAVU transientanalysisupstreamafter the break with variablegrid size and pipe

exposedto the atmosphere

318

CALCQDB

QDB MS-DOS for of of execution execution automatic commands in programme the QUANT software

CLEARSCR

clearing the screen

CTWORK

promptingthe userto indicatewhetheror not he wishesto continue running the programme

DFGEN

Creating ASCII data books from the QUANT programmefor the domain fluid parameter of

FLDPARAS

for in file input binary paras. usr running randomaccess parameters writing the QUANT programme

FLDPROPB Readingoutput datafrom binary random accessfile props.dat producedby the QUANT programme FLDPROPS Readingoutput data from ASCII file producedby the FLDPROPB FLDPROPV Reading output data from ASCII data book produced by the subroutine interpolating for the required input parameters and GRIDGEN

sub-programmefor grid generation

HEATATM

Calculation of heat transferthrough a pipe exposedto the atmosphere

HEATGR

Calculationof heat transferthrough a buried pipe

INIQDB

execution of MS-DOS commandsto call the programme QDB in the QUANT software, for manual inputting of the input parametersand executionof calculation

INITIAL

initiating the executionof the main programme

MCC

sub-programmefor second-orderanalysisusing the MacCormack method

MOCD

sub-programme for second-order method of characteristics analysis downstreamthe break

MOCU

sub-programme for second-order method of characteristics analysis upstreamthe break

MTRANS

main/driver routine

SBMOCAU1 first-order calculationof new dependentparameters p, u and p using the method of characteristics,with uniform grid size and pipe exposedto the atmosphere SBMOCAU2 second-ordercalculationof new dependentparametersp, u and p using the method of characteristicsand with uniform grid size

319

SBMOCAV 1 first-order calculation of new dependentparameters p, u and p using the to the exposed pipe size and grid variable method of characteristics, with atmosphere SBMOCAV2 second-ordercalculationof new dependentparametersp, u and p using the to the exposed pipe and size grid variable with characteristics, method of atmosphere

SBMOCGUI first-ordercalculationof newdependentparametersp, u and p usingthe buried pipe and size method of characteristics,with uniform grid SBMOCGU2 second-ordercalculationof new dependentparametersp, u and p using the buried pipe and size grid method of characteristics,with uniform SBMOCGV 1 first-order calculation of new dependentparameters p, u and p using the buried pipe and method of characteristics,with variablegrid size SBMOCGV2 second-ordercalculationof new dependentparametersp, u and p using the

buried pipe size and with variablegrid methodof characteristics, SBMOCVUI first-order calculation of new dependentparameters p, u and p using the method of characteristicsand with variablegrid size and without interpolation SDMOCAV2 second-ordercalculation of new dependentparametersp, u and p for the

downstream with variablegrid section,usingthe methodof characteristics, sizeandpipeexposedto the atmosphere for dependent the SDMOCGV2 second-order u and p, p parameters calculationof new downstreamsection,usingthe methodof characteristics,with variable grid size and pipe exposedto the atmosphere

downstream STADIBUD adiabatic compressible with backwarddifrerencinganduniform grid size STADIBUU

dif%rencing backward and uniform adiabatic compressibleupstreamwith grid size

STADIBVD

downstreamwith backwarddifferencingand variable adiabaticcompressible grid size

STADIBVU

adiabaticcompressibleupstreamwith backward differencingand variable grid size

STADIFUD adiabatic downstream with forwarddifFerencing compressible anduniform grid size 320

STADIFUU

dilTerencing forward and uniform grid adiabatic compressibleupstream with size

STAD 1FVD

differencing forward downstream and variable with compressible adiabatic grid size

STAD 1FVU

difTerencing forward and variable grid upstream with adiabaticcompressible size

STAD2UD

difi'erencing downstream and with second-order adiabatic compressible uniform grid size

STAD2UU

difTerencing and uniform second-order upstream with adiabaticcompressible grid size

STAD2VD

ditTerencing downstream and compressible with second-order adiabatic variablegrid size

STAD2VU

dif erencingand variable compressible upstream second-order adiabatic with grid size

STEAD

driver routine for steadystateanalysisroutines

STINADUD

break incompressible downstream the of with uniform grid size adiabatic

STINADUU

break incompressible the with uniform grid size upstreamof adiabatic

STINADVD

break incompressible downstream the of with variablegrid size adiabatic

STINADVU

break incompressible the upstream with variable grid size adiabatic

STINISUD

isothermalincompressibledownstreamof the break with uniform grid size

STINISUU

isothermalincompressibleupstreamof the break with uniform grid size

STINISVD

isothermalincompressibledownstreamof the break with variablegrid size

STINISVU

isothermalincompressibleupstreamthe break with variablegrid size

STISOUD

isothermalcompressibledownstreamwith uniform grid size

STISOUU

isothermalcompressibleupstreamwith uniform grid size

STISOVD

isothermalcompressibledownstreamwith variable grid size

STISOVU

isothermalcompressibleupstreamwith variablegrid size

STNONAUD non-adiabaticnon-isothermalcompressibledownstreamwith uniform grid size and pipe exposedto the atmosphere STNONAUU non-adiabaticnon-isothermalcompressibleupstreamwith uniform grid size and pipe exposedto the atmosphere STNONAVD non-adiabaticnon-isothermalcompressibledownstreamwith variable grid size and pipe exposedto the atmosphere 321

STNONAVU

non-adiabatic non-isothermal compressible upstream with variable grid size and pipe exposed to the atmosphere

SUBHTA

Calculation of frictional force and heat transferthrough a pipe exposedto the atmosphere

SUBHTG

Calculation of frictional force and heat transfer through a buried pipe

SYSDATA

inputting generalsystemand gas data

TRANS

driver routine for transientanalysisbefore the break

WKL

Warming-Kutter-Lomax for the third-order analysis using sub-programme method

XSVBKSB

driver routine for solution of linear simultaneousequations

XZRHQR4

driver routine for solution of fourth-order polynomial.

YESNO

prompting the user to answerYES (Y) or NO (N)

NOTE: (i) The namesof the subroutinesusedfor first- andsecond-orderapproximationsin the break by for before the transient analysis are obtained method of characteristics, break, i. letters letters SB SN. first those the two the after of used e. with replacing (ii) The namesof routinesusedfor transientanalysisare obtainedby replacingthe first two letters of those usedafter the break, i.e. BR with letters TR. (iii) The namesof sub-routinesand routines used in the case of buried pipeline are obtainedby replacingthe letter A in the namesof thoseused for pipeline which are by letter G. to the the atmosphere, exposed

322

Smile Life

When life gives you a hundred reasons to cry, show life that you have a thousand reasons to smile

Get in touch

© Copyright 2015 - 2024 PDFFOX.COM - All rights reserved.