Fast Numerical Simulation for Full Bore Rupture

of Pressurized Pipelines
Haroun Mahgerefteh, Pratik Saha, and Ioannis G. Economou
Dept. of Chemical Engineering, University College London, London WC1E 7JE, United Kingdom

An efficient numerical simulation (CNGS-MOC), based on the method of character-

istics for simulating full bore rupture of long pipelines containing two-phase hydrocar-
bons, was de®eloped. The use of cur®ed characteristics, in conjunction with a com-
pound nested grid system, as well as a fast mathematical algorithm, lead to a significant
reduction of CPU time, while impro®ing accuracy. The model is ®alidated extensi®ely
against field data including those obtained during the Piper Alpha tragedy, as well as the
Isle of Grain depressurization tests. Its predictions are compared with those based on
other mathematical models including PLAC, META-HEM, MSM-CS, as well as
BLOWDOWN. Both CNGS-MOC and META-HEM produce reasonably accurate pre-
dictions with the remaining models assessed performing relati®ely poorly.

Introduction
The most common way of exporting oil and gas from off- temperature drops in the fluid and the pipeline are of pri-
shore installations is through long pipelines. The combina- mary importance in order to avoid pipeline brittle fracture or
tion of the high pressure needed to induce flow, as well as the formation of solid hydrates, which can in turn lead to
the massive amount of inventory present in the pipeline, pose blockage.
a risk of full bore rupture ŽFBR., which, in the event of oc- The depressurization process is not amenable to simple
curring, can lead to one of the most catastrophic types of analysis due to its highly transient unsteady-state nature. De-
accidents possible in the offshore industry. During the Piper spite this, in recent years, several workers have proposed em-
Alpha tragedy ŽCullen, 1990., for example, the rupture of the pirical correlations based on experimental data ŽWeiss et al.,
three gas export pipelines, although not the primary cause of 1988.. These are, however, very limited in the range of appli-
the accident, resulted in the major escalation of the incident cability and particularly inappropriate for pipelines contain-
and the eventual collapse of the platform. There have also ing two-phase or condensable mixtures. Also, apart from pro-
been several other documented cases relating to FBR, which viding an estimate of the release rate with time, no other
have resulted in numerous loss of life and significant damage information on the variation of the fluid dynamics within the
to the environment ŽBond, 1996.. pipeline during depressurization is provided. The critical role
Consequently, an important part of the safety assessment of this type of data is already demonstrated in governing the
case ŽHMSO, 1989. for all offshore platforms utilizing subsea dynamic response, as well as the integrity of emergency shut-
pipelines greater than 40 cm in diameter is the prediction of down valves following FBR of gas pipelines ŽMahgerefteh et
the amount of inventory released and its variation with time al., 1997, 1998..
in the event of FBR. Such information has a direct impact on The partial differential equations pertaining to conserva-
almost every safety aspect of the platform including the sur- tion of mass, momentum, and energy for an element of the
vival time of the temporary safe refuge. fluid medium within the pipeline constitute a system of es-
Another area of concern is during controlled release or sentially Euler equations. These incorporate a stiff source
blowdown. In such cases, prior estimations of the resulting term due to the friction in the momentum equation. Heat-
transfer effects, on the other hand, are taken into account in
the energy equation. These non-negligible nonisentropic ef-
Correspondence concerning this article should be addressed to H. Mahgerefteh. fects together with the choked flow boundary condition ren-
Current address of I. G. Economou: Molecular Modelling of Materials Lab., Insti- der modeling of depressurization as a result of FBR different
tute of Physical Chemistry, National Research Centre for Physical Sciences ‘‘De-
mokritos,’’ GR-15310 Aghia Paraskevi Attikis, Greece. from flow in a shock tube ŽLyczkowski et al., 1978 a,b..

AIChE Journal June 1999 Vol. 45, No. 6 1191

 The inviscid part of the conservation equations can be given by ŽChen et al., 1993.
shown to be hyperbolic ŽZucrow and Hoffmann, 1976. with
three real eigenvalues. The resulting Euler equations can be rt q r u x q u r x s 0
solved using a variety of numerical techniques with varying r u t q r uu x q Px y b s 0
degrees of success. Previous attempts include a finite differ-
ence method ŽFDM. ŽChen et al., 1993, 1995a,b; Bendiksen Pt q uPx y a2 Ž r t q u r x . y c s 0 Ž1.
et al., 1991., a finite-element method ŽFEM. ŽLang, 1991;
Bisgaard et al., 1987. and the method of characteristics where r , u, P, b , and a are the density, velocity, pressure,
ŽMOC. ŽChen et al., 1992; Flatt, 1986; Olorunmaiye and friction force term, and acoustic velocity of the fluid, respec-
Imide, 1993.. As all of these techniques in essence involve tively. Subscripts t and x denote property derivative with re-
the numerical discretization of the pipeline into a large num- spect to time and space Žin the x-direction., respectively. Fur-
ber of elements, a solution invariably requires very long CPU thermore, c is the nonisentropic term incorporating heat
time Ž40 days on a 300 MHz Pentium II processor can be transfer and frictional effects given by
typical.. This problem is exceptionally acute in the case of qh y u b
MOC particularly when simulating FBR of long pipelines Ž ) c sw Ž2.
100 km., typical of North Sea operations. rT
In this study, the MOC is adopted to simulate FBR or
blowdown of long pipelines containing condensable or two- where qh is the heat transfer from the pipe wall to the fluid.
phase hydrocarbon mixtures. This technique is employed as The thermodynamic function w is given by
opposed to FDM and FEM, as both have difficulty in han- ­P
dling the choking condition at the ruptured end. The MOC ws ž / Ž3.
handles choked flow intrinsically via the Mach line character- ­s r

istics. Moreover, MOC is considered to be more accurate than

FDM, as it is based on the characteristics of wave propaga- The MOC involves the replacement of the three conservation
tion. Hence, numerical diffusion associated with a finite dif- equations by the following characteristic and compatibility
ference approximation of partial derivatives is avoided ŽChen equations
et al., 1992.. The effects of valve closure during emergency dt 1
shutdown on the fluid dynamics within the pipeline are sim-
ply and conveniently handled by specifying the appropriate
ž /
dx p
s
u
Ž path line characteristic . Ž4.
boundary conditions across the valve at selected time and dt 1
space grids. The number of iterations involved are dramati-
cally reduced in MOC by setting up and solving a system of
ž /
dx "
s
u" a
Ž Mach line characteristics . Ž5.

simultaneous equations. The above together with the use of The characteristic Eqs. 4 and 5 stipulate the way in which
second-order or curved characteristics, and a nested grid sys- information is propagated through a flow field. In regions
tem, result in a significant reduction in CPU time. where fluid properties change dramatically, in a nonlinear
This procedure allowed for the first time extensive evalua- manner, the characteristics are curved ŽFlatt, 1985.. This is
tion of the performance of MOC in simulating FBR through- usually the case for two-phase mixtures or condensable gases.
out the discharge process by comparison with field data, as It is assumed that straight line characteristics can lead to sig-
well as other solution techniques including META-HEM nificant global errors. Figure 1 shows the characteristics
ŽChen et al., 1995a,b; Chen, 1993., MSM-CS ŽChen et al., curves.
1995a,b; Chen 1993., BLOWDOWN ŽRichardson and Saville
1991. and PLAC ŽHall et al., 1993.. The field data were from
pipeline depressurization tests carried out in the Isle of Grain
ŽRichardson and Saville, 1996., as well as those recorded dur-
ing the night of the Piper Alpha tragedy ŽChen, 1993.. The
simulations were performed on the basis of an homogeneous
equilibrium model ŽHEM. in which all phases are assumed to
be at thermal and phase equilibrium. Comparison with exper-
imental data ŽPicard and Bishnoi, 1989; Chen, 1993; Chen et
al., 1995. has shown that this assumption is applicable to FBR
of long Ž )100 m. pipelines.

Model Development
Characteristic and compatibility equations
For one-dimensional Ž1-D. multiphase homogeneous flow
where thermodynamic and phase equilibrium exist simultane-
ously and that fluidrstructure interactions ŽLawooij and Tijs-
selling, 1990; Stittgen and Zielke, 1990. are negligible, the
continuity, momentum, and energy conservation equations for
an element of the fluid within the pipeline are respectively Figure 1. Curved (second-order) characteristics.

1192 June 1999 Vol. 45, No. 6 AIChE Journal

 The pathline compatibility C o is where the subscripts g and l denote the gas and liquid phase,
respectively, and x refers to the mass of vapor per unit mass
c of bulk fluid. The speed of sound a and w , for real multicom-
dP y a2 d r s dx s c dt Ž6.
u ponent multiphase fluids, are written as ŽGroves et al., 1978;
Picard and Bishnoi, 1988. follows. For the case of single phase
The positive Cq and negative Cy compatibility equations can mixtures, it is
be written as
g
d " P " r ad " us w c " a b x dt Ž7. a2 s Ž 10.
kr
The values of P, r , u, and a as a function of time and rj Ta2
distance along the pipeline Žsay, points, p, o and n, in Figure ws Ž 11.
1. may be obtained by the inverse marching method of char- cp
acteristics. This involves dividing the pipeline into a large
number of distance Ž D x . and time elements Ž D t ., expressing where g is the ratio of specific heats, and k and j are the
the compatibility equations in finite difference form, and isothermal and isobaric coefficients of volumetric expansion,
solving them at the intersection of the characteristics curves respectively, and c p is the specific heat capacity at constant
with the spatial axis Žpoints p, o, n, and j.. The loci of points, pressure. All these parameters are easily obtained from the
p, o, n, and j are found by accounting for the curvature of the Peng-Robinson EoS. For two-phase flows, the analytical de-
characteristics by considering them as arcs of parabolas and termination of g and c p becomes complex. Hence, the pa-
making use of the geometrical fact that the tangents at point rameters a and w are evaluated numerically at a given tem-
j and p meet halfway along the time axis at point q. The perature and pressure. The sound velocity at a certain tem-
previous time line properties are calculated directly from so- perature T and pressure P can be expressed as
lution of quadratic interpolation formulas for spatial dis-
cretization. The effect of accounting for curved characteris- DP
a2 s Ž 12.
tics as opposed to the extensively used linear approximation r ŽT , P . s y r ŽT U , P y D P . s
ŽChen et al., 1992; Flatt, 1986; Picard and Bishnoi, 1989. on
the CPU time required and on the accuracy is quantified later. where the subscript s denotes a constant entropy condition.
Numerical stability requires the fulfillment of the Courant- To solve for T U , a Newton-Raphson iteration is performed
Friedrich-Lewy Criterion ŽCourant et al., 1926; Zucrow and where the objective function is written as
Hoffman, 1976. given by

Dx v Ž n. s s Ž T , P . y s Ž T U Ž n. , P y D P . Ž 13.
Dts Ž8.
Ž < uq a< max . t s 0 where superscript Ž n. denotes the iteration level.
The coupling of the flow equations using a compound The numerical expression for the parameter w can be ob-
nested grid system ŽCNGS, see later. together with the ther- tained from the expression
modynamic and phase equilibrium equations based on HEM
form the basis of the CNGS employed in this study. Before DP DT
ws sr2
the calculation procedure for this model is presented, the im-
portant thermophysical and hydrodynamic constitutive rela-
ž /
Ds r ž /
Dr s
Ž 14.

tions are to be addressed. All hydrodynamic constitutive The latter equation is solved numerically by performing an
relations pertaining to momentum exchange between two isentropic flash calculation using Eq. 12.
phases are nonexistent as we assume HEM.
Frictional force effects
Thermophysical properties
At the fluid-wall interface, the presence of the viscous drag
Phase equilibrium calculations are performed with the force can be modeled as a linear combination of unsteady
widely used Peng-Robinson equation-of-state EoS ŽPeng and and steady wall drag or friction ŽChen, 1993.. From an earlier
Robinson, 1976., which is particularly useful for handling review ŽSaha, 1997., it seems that there is still some uncer-
multicomponent hydrocarbon mixtures. The number and the tainty in the accurate prediction of an unsteady friction fac-
appropriate fluid phaseŽs. present are obtained using the sta- tor in rough pipes where highly turbulent flows prevail.
bility test based on the Gibbs tangent plane criterion devel- Moreover, from the few studies performed to date ŽZucrow
oped by Michelsen Ž1982a,b, 1987.. For unstable systems, the and Hoffman, 1976; Wood and Funk, 1970., the effect of the
technique also provides the composition of a new phase which unsteady component in wall shear stress calculations seems
can be split off to decrease the Gibbs energy of the mixture. to diminish in magnitude as the turbulence of flow increases.
Pseudo-fluid properties for mixtures are calculated on the In the absence of any theoretically and experimentally justi-
basis of pure liquid and gas properties obtained from the EoS. fied transient turbulent friction calculation for rough pipes, it
For example, the mixture density r is given by was decided to ignore this component in the model proposed
rg rl in the present study.
rs Ž9. The steady flow friction factor was estimated using the
r g Ž 1q x . q r l x
Moody ŽMassey, 1983. approximation to Colebrooke’s equa-

AIChE Journal June 1999 Vol. 45, No. 6 1193

tion. It is the most accurate expression available, represent- The change in velocity as a result of this pressure drop is
ing results within 5% of experimental data, and is given by
DP
e 10 6
1r3
ž
u i s u iy1 1y
Ž r ZRT . iy1 / Ž 20.
f w s 0.001375 1q 20,000
ž D
q
Re / Ž 15.
where the subscript i denotes the grid points at which the
fluid conditions are being calculated, and iy1 refers to the
which is valid for ReG 2,000. fluid conditions at the previous grid point.
For ReF 2,000, the following well established laminar flow Rupture Plane Calculation at t s 0. For a change in pres-
correlation is employed sure from initial pressure Pin to a given pressure P following
FBR, the fluid discharge velocity at the rupture plane is given
16 by
fw s Ž 16.
Re
P dP
usy HP Ž 21.
Two-phase mixtures are simply handled by replacing sin- in
ra
gle-phase properties by two-phase mixture properties. There-
fore, for a two-phase mixture, For a condensable gas mixture, the thermodynamic relation-
ship between r , a, and P becomes highly nonlinear and Eq.
f wm 21 can only be solved numerically. Expressing this equation
b sy2 rm u < u < Ž 17. in finite difference form results in
D

where f wm is the two-phase mixture friction factor and can Prq1 y Pr

u rq1 s u r y Ž 22.
be determined either from Eq. 15 or 16. Ž r a. r
The mixture viscosity is given by
where r shows the current level of decompression and r q1
1 x Ž 1y x . denotes the new level resulting from an incremental decom-
s q Ž 18.
mm mg ml pression step. In this study the new pressure is taken as

Prq1 s 0.95Pr Ž 23.

The gas and liquid viscosities are calculated according to the
Ely and Hanley scheme for nonpolar gaseous mixtures, and
the Dymond and Assael scheme for liquid mixtures ŽMassey, The calculation procedure for the boundary condition at the
1983.. rupture plane at t s 0 assuming isentropic expansion is as
follows:
v Start the decompression process at the initial pressure

Numerical discretization of the pipeline Pin and temperature Tin and calculate the bulk density and
The grid structure for numerical discretization of the speed of sound by performing an isothermal flash Žthe en-
pipeline length is based on a Compound Nested Grid System tropy is kept constant throughout the decompression process..
ŽCNGS. described in earlier publications ŽMahgerefteh et al., v Calculate the new fluid pressure and velocity using Eqs.

1997, 1998.. Briefly, it involves the use of increasingly finer 23 and 22, respectively. Perform an isentropic pressure-
grids near the rupture plane where the transients are most entropy flash based on the initial entropy, and the new pres-
rapid. This procedure results in a significant reduction of CPU sure to obtain a new density and speed of sound.
v Compare the new speed of sound a
time as coarser grids can in turn be used at locations away rq1 and velocity u rq1 .
from the rupture plane. Since the smaller grids are geometri- If velocity is less than speed of sound, then repeat calculation
cally similar, and contained within the large normal mesh, a until u rq1 s a rq1 , when the fluid conditions at the rupture
consistent Courant number ŽCourant et al., 1926. is main- plane are known based on an isentropic decompression.
tained throughout the discharge process and numerical insta- Interior Point Calculation at t ) 0. The methodology for the
bility is avoided. calculation of the conditions at a solution point Žsuch as point
j in Figure 1. was described earlier. Calculation of the corre-
Initial condition sponding temperature involves the solution of the following
equation
Steady-State Flow Prior to Rupture. It is assumed that the
flow in the pipeline prior to rupture can be taken to be r j y r Ž Tjr , Pj . s 0 Ž 24.
isothermal steady state. It can be shown that the correspond-
ing pressure drop across an element D x is given by
based on an iterative numerical scheme. The subscript j de-
notes conditions at the solution point Žsee Figure 1., and the
biy1
Pi s Piy1 q Dx Ž 19. superscript r is related to the unknown temperature. The so-
u2 lution of Eq. 24 becomes a root finding problem where a
1y temperature is sought to match the density obtained from the
ZRT iy1 compatibility equations to that calculated from an isothermal

1194 June 1999 Vol. 45, No. 6 AIChE Journal

 and 7, respectively. The calculations of bulk density and speed
of sound are obtained from Eqs. 9 and 10 or 12, respectively.
Calculation at the Boundary between Grids of Different Size.
The use of a nested grid system creates boundaries between
different size grids. To deal with these regions, a direct solu-
tion technique has been developed which is particularly suited
to highly transient two-phase flows. In Figure 3, this tech-
nique is illustrated for the calculation at solution point jq4.
The path line compatibility is given by

do t 1 4D t 1
lo s s s ´ x o s x i y4D t 1 u o Ž 27.
do x uo xi y xo

The positive and negative characteristics are respectively

given by

dq t 1 4D t 1
s s ´ x p s x i y4D t 1 Ž u p q a p . Ž 28 .
Figure 2. Boundary cell layout. dq x u p q ap xi y x p

dy t 1 D t1
s s ´ x n s x i y D t 1 Ž u n y a n . Ž 29 .
pressure-temperature flash. For this purpose, the method of dy x u n y an xi y xn
Brent ŽPress et al., 1994. was employed according to which
the bracketing of the root is performed, and, for our calcula- Using Eqs. 28]29 and linear interpolation, the simultaneous
tions, it produces rapid convergence. Once the temperature equations for the solution at points u p and a p are
is obtained at the solution point, the speed of sound and the
parameter w are found by a flash calculation. The above steps u i y u iy1 u i y u iy1
are repeated until convergence is achieved for the dependent
variables P, r and u.
ž
u p 1q
D x2
4D t 1 q
/ D x2
4D t 1 a p s u i Ž 30.
Rupture Plane Calculation at t ) 0. The methodology used
a i y a iy1 a i y a iy1
for the curved characteristics in the boundary cell involves
the addition of an extra ghost cell adjacent to the boundary
ž
a p 1q
D x2
4D t 1 q
/ D x2
4D t 1 u p s a i Ž 31.
cell as shown in Figure 2. This permits quadratic interpola-
tion for points p and o. The conditions at the node iq1 are Similarly, a 2=2 system of equations can be set up for ob-
the same as at the node i, that is taining u n and a n. These are given by

V iq1 s V i Ž 25. u iq1 y u i u iq1 y u i

D t1 y D t1 an s u i
where
u n 1q
ž D x1 / D x1
Ž 32.

V s P , u, r , T , or a. a iq1 y a i a iq1 y a i
D t1 y D t1 u n s ai
The rest of the calculations are exactly the same as those
a n 1q
ž D x1 / D x1
Ž 33.

for the interior point with the exception that simultaneous

solution of the positive and pathline compatibility equations
is not possible, since no algebraic relationship exists for the
speed of sound of a two-phase mixture. The negative or left
running characteristic at the rupture plane will be vertical
wthe gradient, 1rŽ uya. is infinityx, and, hence, perpendicular
to the x-axis, as shown in Figure 2. Since the negative Mach
line is vertical, the conditions at the previous time level x i
are already known, so that the negative compatibility equa-
tion expressed in a finite difference form becomes

1 1
Pj y Pi y
2
Ž Ž r a. i q Ž r a. j . Ž u j y u i . s 2 Ž Ž c y a b . i
qŽ c y ab . j . Ž t j y ti . Ž 26.

This compatibility equation can be solved together with the

pathline and positive compatibility equations given by Eqs. 6 Figure 3. Boundary between fine and coarse meshes.

AIChE Journal June 1999 Vol. 45, No. 6 1195

The solution for u o depends on the direction of fluid flow. Table 1. Piper-Alpha to MCP-01 Pipeline and Isle of
For l o ) 0 then Grain Test DataU
Isle of Grain
ui
uo s Ž 34. Piper Alpha Test P40 Test P42
u i y u iy1
ž 1q
D x2
4D t 1
/ Pipeline length, L Žkm.
Inner dia., D Žm.
Initial pres., Pin Žbar.
54
0.419
117
0.1
0.15
21.6
0.1
0.15
11.3
Initial temp., Tin ŽK. 283 290.95 293.15
For l o - 0 Heat-transfer coeff. Uh 5.0 ŽChen 100 ŽChen, 100
ŽWrm2 ? K. et al., 1992. 1993.
ui Pipe roughness factor Žmm. 0.05 0.05 0.05
uo s Ž 35. Ambient temp., T` ŽK. 283 292.25 291.75
u iq1 y u i
D t1
ž 1q
D x1 / U
Piper Alpha inventory Žmolar % .; CH 4 Ž73.6 ., C 2 H 6 Ž13.4 ., C 3 H 8 Ž7.4 .,
i-C 4 H 10 Ž0.4 ., n-C 4 H 10 Ž1.0 ., i-C 5 H 12 Ž0.08 ., n-C 5 H 12 Ž0.07 ., n-C 6 H 14
Ž0.02 ., N 2 Ž4.03 .. Isle of Grain inventory Žmolar % .; C 3 H 8 Ž95 ., n-C 4 H 10
Ž5 ..
The locations of x o , x p , and x n can now be obtained di-
rectly from Eqs. 27]29 by substituting the calculated values
for u p , a p , u n , a n , and u o from the above equations. The
values of P and T at the three initial points are then calcu- Case Study
lated using linear interpolation, and the corresponding den- For obvious reasons, the amount of experimental data
sity r and speed of sound a are obtained by performing an available relating to FBR is very scarce. The limited data re-
isothermal flash. ported are mainly confined to short Žca. 100 m. pipelines of
All the initial point flow variables are now available to small diameter Žca. 0.1 m. containing single component flu-
compute the flow conditions at the solution point j. For the ids. Following an extensive literature search, two sets of field
positive compatibility, it is data were identified, which are suitable for modeling pur-
poses. The first set includes the intact end pressure data
Pj y Pp q Ž r a . p Ž u j y u p . s Ž c q a b . p 4D t 1 s K 1 Ž 36 . recorded during FBR of a long gas line, containing mainly
methane connecting the Piper Alpha and MCP-01 platforms
while for the negative compatibility, it is ŽChen, 1993.. The other set contains field data related to a
series of pipeline depressurization tests carried out by Shell
and British Petroleum on the Isle of Grain ŽRichardson and
Pj y Pn y Ž r a . n Ž u j y u n . s Ž c y a b . n D t 1 s K 2 Ž 37 .
Saville, 1996.. In these tests, extensively instrumented
pipelines containing commercial liquefied petroleum gas
Furthermore, the pathline compatibility for l0 ) 0 is ŽLPG. and incorporating pressure transducers and thermo-
couples were used. Inventory and holdup were also recorded
Pj y Po y Ž a2 . o Ž r j y r o . s co4D t 1 Ž 38. using load-cells and neutron back scattering. Table 1 shows
the prevailing experimental conditions, as well as the fluid
compositions for the above cases.
and the pathline compatibility for l o - 0 is

Pj y Po y Ž a2 . o Ž r j y r o . s co D t 1 Ž 39. Results and Discussion

Optimization of CNGS-MOC; effect of using a second-order
Solving Eqs. 36 and 37 simultaneously for u j one can write or cur©ed characteristic on the simulation accuracy and on
CPU time
K 1 y K 2 q Ž r a . p u p q Ž r a . n u n q Pp y Pn Figure 4 demonstrates the effect of using second-order
uj s Ž 40. characteristics on the CPU time when simulating FBR Žsee
Ž r a. n q Ž r a. p Table 1.. The data show the variation of predicted release
rate with time over the first 450 s following FBR using differ-
The pressure may be calculated either from Eq. 36 or Eq. 37. ent numerical discretization techniques. All simulations were
The density at the solution point can now be obtained from performed on a DEC Alpha server 8400 5r440. Curve A shows
the pathline compatibility. For l o ) 0 the results using a simple grid system ŽSGS. with a uniform
D x of only 10 m. This calculation is used as a reference for
Ž Pj y Po . q a2o r 0 y co4D t1 accuracy as it uses the finest grid system compared to the
rj s Ž 41. other simulations, and, hence, errors introduced as a result
a2o of numerical discretization are minimized. Curve B shows the
results using the second-order Žcurved characteristics . solu-
The above procedures, also called the predictor steps, pro- tion for CNGS with a coarse grid D x of 500 m. Curve C
duce first estimates of the various flow parameters. Better shows the CNGS results for the second-order solution with
estimations for these variables are obtained using an iterative D x of 250 m and, finally, curve D shows the first-order Žlin-
procedure known as the corrector steps ŽZucrow and Hoff- ear approximation of characteristics . CNGS solution with D x
man, 1976.. of 250 m. From these calculations, it is clear that the first-

1196 June 1999 Vol. 45, No. 6 AIChE Journal

 Figure 5. Refraction of characteristic lines at the single-
phaser
rtwo-phase interface.

were to cross such an interface, significant refraction would

occur. Figure 5 shows this type of phenomenon. To avoid
instability, a substantially smaller value of D t compared to
that dictated by the CFL stability criterion is required. Expe-
rience has shown that a typical D t value required in order to
avoid instability is approximately 10% of the maximum value
dictated by the CFL stability criterion.
Interestingly, it is found that even such a significant reduc-
Figure 4. Piper Alpha: MCP line release rate profiles for tion in D t does not result in an appreciable increase in the
different grid settings. CPU time. This is primarily a consequence of the smaller
Curve A: First-order simple grid system ŽSGS ., D x s 10 m, number of iterations required in the corrector step in order
CPU time s 250 h; Curve B: Second-order ŽCNGS-MOC ., to reach a solution.
D x s 500 m, CPU time s 3.75 h; Curve C: Second-order
ŽCNGS-MOC ., D x s 250 m, CPU time s 12.2 h; Curve D:
First-order ŽCNGS-MOC ., D x s 250 m, CPU time s 12.1 h.
Validation of the CNGS-MOC
Figure 6 shows the measured intact end pressurertime his-
order CNGS solution Žcurve D. consistently underestimates tory following the FBR of Piper Alpha to MCP-01 subsea
the release rate, whereas both second-order solutions pro-
duce good agreements. For the second-order solutions Žcurves
B and C., an increase in D x does not lead to an appreciable
loss of accuracy, but the reduction of the CPU time is signifi-
cant Žcf. 12.2 h with 3.7 h.. It is noteworthy that the use of
the second-order characteristics in conjunction with the
CNGS results in approximately 75 fold reduction Žcf. 250 h;
curve A with 3.75 h curve B. in CPU time. Based on the
above, the second-order CNGS with a grid spacing of 500 m
is chosen for the remaining of the validation tests.
Picard and Bishnoi Ž1989. reported diverging oscillations in
pressure and density values with rapid termination of their
FBR simulation when using a first-order method of charac-
teristics. The authors were thus unable to report any data
beyond the first few seconds of the modeling of a hypotheti-
cal break to a sour gas pipeline. The investigations show that
this is mainly a consequence of the following:
Ž1. Use of linear interpolation for calculating bulk mixture
density at intermediate points along the grid system. This
presupposes a linear variation of fluid density with distance
along the pipeline, thus giving rise to escalating errors. These
errors become exacerbated in the case of two-phase mixtures
particularly in the vicinity of the rupture plane. This problem
is avoided by directly obtaining the fluid density from the
EoS via an isothermal flash calculation at the specified pres- Figure 6. Intact end pressure vs. time profiles for Piper
sure. Alpha-MCP pipeline.
Ž2. The abrupt change in the speed of sound on crossing Curve A: Field data; Curve B: CNGS-MOC, CPU times 6
the gasrliquid interface means that if the characteristic lines days; Curve C: CNGS-MOC ideal gas, CPU times 1.5 min.

AIChE Journal June 1999 Vol. 45, No. 6 1197

 possibly due to the reflection of expansion waves in the com-
pressed liquid off the closed end of the pipe. The theoretical
predictions represented by curves C and D are in fair agree-
ment with the experimental data. Nevertheless, since imme-
diate transition to the saturation conditions, is assumed one
is unable to predict the initial undershoot due to the transi-
tion from liquid to two-phase flow. The finite discrepancies
between theory and experiment may be attributed to a num-
ber of factors including the uncertainty associated with the
measurement of pressure, the inaccuracies associated with the
prediction of VLE data, the uncertainty in the selected pipe
roughness, as well as the lack of accurate information on the
fluid composition. The fluid composition chosen in this study
is an approximate one for commercial LPG.
It is noteworthy that the measured data in Figure 8, as well
as those presented in Figures 9]12 exhibit finite oscillations,
although none of the simulations reviewed in this study, in-
cluding ours, exhibit such behavior. It is strongly believed that
this is a consequence of ignoring fluid-structure interaction
ŽLawooij and Tijsselling, 1990; Stittgen and Zielke, 1990. in
the modeling. This is essentially a dynamic phenomenon, with
the interaction being caused by dynamic forces which act
conversely on fluid and pipe.
Figure 7. Condensation volume fraction at closed end The simulations consider only rupture in straight, horizon-
for Piper Alpha-MCP line. tal well-anchored pipelines in which the fluid compressibility
is by far smaller than the pipe wall elasticity. Fluid-structure
interaction can effectively be ignored.
line. Curve A shows measured data, whereas curve B shows Figure 9 shows the measured Žcurves A and B. and pre-
the predictions using the second-order CNGS-MOC, as de- dicted temperature profiles Žcurves C and D. at the open and
scribed above. This is the first time that the entire release closed pipeline ends for test run P40. The predicted tempera-
process has been simulated using a two-phase MOC model.
Theory and experiment are in excellent agreement. Curve C
shows the corresponding data generated using first-order
characteristics in conjunction with the ideal gas assumption,
CNGS-ideal as reported in our previous publication
ŽMahgerefteh et al., 1997.. The discrepancy between ideal gas
prediction and field data Žcurve A. demonstrates the impor-
tance of taking into account real fluid behavior when requir-
ing accurate predictions.
The results for CNGS-MOC Žcurve B. show a slight deflec-
tion at about 1500 s. This corresponds to the onset of con-
densation at the closed end as indicated in Figure 7. These
data show the variation of volume fraction of condensate with
time at the closed end. The maximum liquid volume fraction
is ca. 0.034 corresponding to a liquid mass fraction of about
15%. As the pipeline continues to depressurize, less inven-
tory stays in the liquid phase and at about 10,000 s following
FBR the closed end is exposed to gas only.

FBR experiments using LPG at the Isle of Grain

Figure 8 shows predictions for the open- and closed-end
pressure-time histories for the LPG mixture, as compared to
experimental data. Curves A and B show the measured data,
whereas curves C and D represent the corresponding simu-
lated data using CNGS-MOC. The extremely rapid initial de-
crease in pressure is due to an almost instantaneous change
Figure 8. Pressure-time profiles at closed and open
of phase from compressed liquid propane to a flashing two-
ends for the P40 test (LPG).
phase mixture. The subsequent slight undershoot measured
Curve A: Field data Žopen end .; Curve B: Field data Žclosed
in pressure Žcurves A and B. is probably attributed to end .; Curve C: Open end, CNGS-MOC; Curve D: Closed
nonequilibrium effects such as delayed bubble nucleation, and end, CNGS-MOC.

1198 June 1999 Vol. 45, No. 6 AIChE Journal

Figure 9. Temperature-time profiles at the open and
closed ends for the P40 test (LPG). Figure 11. Pressure-time profiles at the open end for the
Curve A: Field data Žopen end .; Curve B: Field data Žclosed P42 test (LPG).
end .; Curve C: CNGS-MOC Žopen end .; Curve D: CNGS-
MOC Žclosed end .. Curve A: Field data; Curve B: CNGS-MOC; Curve C:
META-HEM; Curve D: MSM-CS; Curve E: Blowdown;
Curve F: Plac.

ture profiles show similar trends with those observed above

for the pressure profiles when compared to experimental data.
rapid increase in the measured temperature at the open end
The open-end predictions are slightly higher, whereas the Žcurve A. can be observed towards the end of the depressur-
closed-end predictions are slightly lower than measured data.
ization process. This rise corresponds to the cessation of
The maximum discrepancy is about 58C. A significant and

Figure 12. Pressure-time profile predictions at the

Figure 10. Total line inventory predictions for the P40 closed end for the P42 test (LPG).
test (LPG). Curve A: Field data; Curve B: CNGS-MOC; Curve C:
Curve A: Field data; Curve B: META-HEM; Curve C: META-HEM; Curve D: MSM-CS; Curve E: BLOW-
CNGS-MOC. DOWN; Curve F: PLAC.

AIChE Journal June 1999 Vol. 45, No. 6 1199

two-phase flow with the subsequent onset of gas-phase flow namic prediction due to the narrow phase envelop involved.
at the rupture plane, and is probably a reflection of different In PLAC, fluid physical properties are calculated at a given
heat-transfer rates from the wall to the fluid. This trend can- pressure and temperature in a cell by interpolation from
not be reproduced with the model as a constant heat-transfer look-up tables generated using a phase equilibrium package
coefficient is employed. The slight leveling of the tempera- before transient flow calculations begin.
ture immediately before its subsequent rise corresponds to Based on the comparison of total line inventory prediction
the moment when the flow at the rupture plane ceases to be for test P42, Chen Ž1993. showed that the initial inventory
choked. from PLAC starts at a much lower value Ž695 kg. than the
Figure 10 shows a comparison for the total line inventory experimental Ž970 kg., and drops to less than 20% within 5
during depressurization for the test run P40. Curve A shows seconds. The considerable under-estimation can be at-
the load cell values, whereas curves B and C show the predic- tributed to two possible reasons. Firstly, too much gas forma-
tions using META-HEM ŽChen et al., 1993, 1995; Chen, 1995. tion is predicted in the pipeline by PLAC resulting in lower
and the CNGS-MOC model, respectively. META-HEM is a line inventory prediction. Hall et al. Ž1993. take the initial
homogeneous equilibrium model that uses a FDM to resolve state of the fluid as 80% liquid which means that in PLAC
the Euler equations. Both models slightly overestimate line transient flow calculations start at a point that is already quite
inventory, but give very similar results. well into two-phase flow as opposed to near the saturation
In Figure 11, experimental data Žcurve A. and theoretical point. In the CNGS-MOC, calculations start at the saturation
predictions are shown for the open-end pressure for the P42 point producing good initial agreement. A second reason
case using CNGS Žcurve B. and other models including could be due to overestimation of the release rate by the ho-
META-HEM Žcurve C., MSM-CS Žcurve D., BLOWDOWN mogeneous frozen critical flow boundary condition.
Žcurve E., and PLAC Žcurve F. for the LPG mixture. In Fig- Figure 12 shows the same data as those presented in Fig-
ure 12, the corresponding predictions at the closed end are ure 11, but for the closed end. CNGS-MOC, META-HEM,
shown. The BLOWDOWN model was developed by Richard- and BLOWDOWN give very similar predictions with MSM
son and Saville Ž1991., the MSM model by Chen Ž1993., and doing less well and PLAC performing very poorly. As for the
PLAC by Hall et al. Ž1993.. The BLOWDOWN simulation open-end pressure predictions, PLAC shows a very sharp drop
Žcurve E. for pipelines is based on quasi-steady-state, equilib- in pressure at ca. 5 s and continues to drop at a much faster
rium, and homogeneous two-phase flow assumptions. The rate than the other models thereafter. Once again, the
MSM-CS model Žcurve D. takes into account the slip velocity CNGS-MOC provides relatively good predictions of the con-
between the two phases, and, thus, incorporates the relevant ditions following FBR.
hydrodynamic constitutive relations that accompany this ef-
fect. As a result, two momentum equations, one for each
phase, are solved instead of just one in META-HEM. The Conclusions
MSM is based on the same solution techniques as META- In this article, the development of an efficient numerical
HEM. PLAC Žcurve F. is an extension of the nuclear reactor simulation ŽCNGS-MOC. is described that is based on the
safety code TRAC in that it incorporates a flash calculation method of characteristics for simulating the variations in fluid
package MULTIFLASH to predict VLE data for hydrocar- dynamics following full bore rupture of long pipelines con-
bon mixtures. In addition, the original critical flow boundary taining high-pressure two-phase hydrocarbon mixtures. The
condition in TRAC is replaced by a homogeneous frozen flow long CPU time has been largely addressed, and this has been
model that allows different temperatures for each of the two synonymous so far with such types of simulations by using
phases, and, thus, permits thermal nonequilibrium since two curved characteristics in conjunction with a compound nested
energy equations are solved for each phase. This creates a grid system, as well as setting up and solving a system of si-
dilemma as to which temperature should be used for the multaneous equations. As opposed to the commonly used lin-
equilibrium phase behavior, an inevitable problem in deter- ear characteristics, curved characteristics are found to afford
mining the phase behavior of multicomponent mixtures. the use of much larger discretization grids, while at the same
From Figure 11 it is apparent that the best predictions are time improve the global accuracy. Also, the number of itera-
obtained from CNGS-MOC and META-HEM. BLOW- tions involved in arriving at a solution in the corrector step is
DOWN substantially underestimates open-end pressure be- significantly reduced.
cause of the quasi-steady flow assumption and MSM-CS gives It is expected that the application of curved characteristics
almost the same predictions as BLOWDOWN. It is surpris- is particularly pertinent in the case of two-phase flows, espe-
ing that the former cannot provide a better prediction of the cially near the rupture plane where the variations in fluid
pressure profile at the open end. This could be due to the properties are highly nonlinear. Instability problems are over-
use of unrealistic hydrodynamic constitutive relations. An- come, which invariably result in premature termination of
other possible reason is the inaccuracy in predicting the wave such types of simulations by obtaining the relevant fluid
propagation velocities. properties at intermediate points along the grids from the
PLAC also gives rise to poor predictions. It predicts a large equation of state directly. This is an improvement over previ-
drop in pressure of at least 2 bara at about 5 s after rupture ous work, where an interpolation was used which relies on a
and continues to drop thereafter. One reason for this dis- linear variation of fluid properties in between grid points.
crepancy might be due to the fact that no actual flash calcu- Extensive comparison with experimental data revealed that
lations are performed to determine fluid properties during the model produces reliable results in terms of predicting
the depressurization. Depressurization profiles for mixtures various parameters including pressure, temperature, and to-
such as LPG are especially sensitive to accurate thermody- tal line inventory variations with time following FBR. It also

1200 June 1999 Vol. 45, No. 6 AIChE Journal

compares very favorably with META-HEM, while, in com- of Two-phase Hydrocarbon Flows in Pipelines,’’ Proc. Euro. Two-
parison, other models such as BLOWDOWN, MSM-CS, and Phase Flow Group Meeting, Hannover, Germany ŽJune 6]10, 1993..
HMSO, ‘‘The Offshore Installations Emergency Pipeline Valve Reg-
PLAC perform surprisingly poorly. The success of the ulations, Statutory Instruments,’’ Health and Safety Executive,
CNGS-MOC in simulating FBR is indicative of the applica- Bootle, U.K. Ž1989..
bility of HEM in such types of processes. Lawooij, C. S. W., and A. S. Tijsselling, ‘‘Waterhammer with Fluid-
In conclusion, this work shows that FBR can be simulated Structure Interaction,’’ App. Scient. Res., 47, 273 Ž1990..
Lang, E., ‘‘Gas Flow in Pipelines Following a Rupture Computed by
effectively using the method of characteristics in conjunction
a Spectral Method,’’ J. Appl. Math. and Physics (ZAMP), 42, 183
with the homogeneous equilibrium model, given an accurate Ž1991..
thermodynamic and phase behavior model. Apart from its ac- Lyczkowski, R. W., D. Gidaspow, C. W. Solbrig, and E. D. Hughes,
curacy, a particularly attractive feature of MOC is the ease ‘‘Characteristics and Stability Analysis of Transient One-dimen-
with which it may be used for simulating the dynamic re- sional Two-Phase Flow Equations and Their Finite Difference Ap-
proximations,’’ Nucl. Sci. Eng., 66, 378 Ž1978a..
sponse of emergency shutdown valves following FBR by sim- Lyczkowski, R. W., R. A. Grimesey, and C. W. Solbrig, ‘‘Pipe Blow-
ple implementation of the appropriate boundary conditions down Analyses Using Explicit Numerical Schemes,’’ AIChE Symp.
at the location of the valves. Ser., 174, 129 Ž1978b..
So far, all of the existing numerical mathematical treat- Mahgerefteh, H., P. Saha, and I. G. Economou, ‘‘A Study of the
Dynamic Response of Emergency Shutdown Valves Following Full
ments for simulating FBR require the use of powerful work-
Bore Rupture of Gas Pipelines,’’ Trans. Inst. Chem. Eng., Part B,
stations. Given the current rapid pace in computer process- 75, 201 Ž1997..
ing speeds, it is believed that the CNGS-MOC described in Mahgerefteh, H., P. Saha, and I. G. Economou, ‘‘Control Valves for
this study offers the best chance for carrying out such types Pipe Rupture,’’ The Chem. Engineer, 666, 26 Ž1998..
of simulations accurately on personal computers in reason- Massey, B. S., Mechanics of Fluids, Van Nostrand Reinhold, Woking-
ham, U.K. Ž1983..
ably practical time scales. Michelsen, M. L., ‘‘The Isothermal Flash Problem: I. Stability,’’ Fluid
Phase Equil., 9, 1 Ž1982a..
Michelsen, M. L., ‘‘The Isothermal Flash Problem: II. Phase-Split
Literature Cited Calculation,’’ Fluid Phase Equil., 9, 21 Ž1982b..
Assael, M. J., J. P. Martin Trusler, and T. Tsolakis, Thermophysical Michelsen, M. L., ‘‘Multi-Phase Isenthalpic and Isentropic Flash Al-
Properties of Fluids, Imperial College Press, London Ž1996.. gorithms,’’ Fluid Phase Equil., 33, 13 Ž1987..
Bendiksen, K. H., D. Malnes, R. Moe, and S. Nuland, ‘‘The Dynamic Olorunmaiye, J. A., and N. E. Imide, ‘‘Computation of Natural Gas
Two-Fluid Model OLGA: Theory and Application,’’ SPE Produc- Pipeline Rupture Problems Using the Method of Characteristics,’’
tion Eng., 6, 171 Ž1991.. J. Hazardous Mat., 34, 81 Ž1993..
Bisgaard, C., H. H. Sorensen, and S. Spangenberg, ‘‘A Finite Ele- Peng, D. Y., and D. B. Robinson, ‘‘A New Two-constant Equation of
ment Method for Transient Compressible Flow in Pipelines,’’ Int. State,’’ Ind. Eng. Chem. Fund., 15, 59 Ž1976..
J. Num. Meth. Fluids, 7, 291 Ž1987.. Picard, D. J., and P. R. Bishnoi, ‘‘The Importance of Real-fluid Be-
Bond, J., ‘‘IChemE Accidents Database,’’ IChemE, Rugby, U.K. haviour in Predicting Release Rates Resulting from High Pressure
Ž1996.. Sour-gas Pipeline Ruptures,’’ Can. J. Chem. Eng., 67, 3 Ž1989..
Chen, C. R., S. M. Richardson, and G. Saville, ‘‘Numerical Simula- Picard, D. J., and P. R. Bishnoi, ‘‘The Importance of Real-fluid Be-
tion of Full-bore Ruptures of Pipelines Containing Perfect Gases,’’ haviour and Nonisentropic Effects in Modelling Decompression
Trans. Inst. Chem. Eng. Part B, 70, 59 Ž1992.. Characteristics of Pipeline Fluids for Application in Ductile Frac-
Chen, J. R., ‘‘Modelling of Transient Flow in Pipeline Blowdown ture Propagation Analysis,’’ Can. J. Chem. Eng., 66, 3 Ž1988..
Problems,’’ PhD Thesis, Univ. of London, Imperial College Ž1993.. Press, H. W., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,
Chen, C. R., S. M. Richardson, and G. Saville, ‘‘A Simplified Numer- Numerical Recipes in Fortran, 2nd ed., Cambridge Univ. Press,
ical Method for Transient Two-phase Flow,’’ Trans. IChem. E., 71A, Cambridge Ž1994..
304 Ž1993.. Richardson, S. M., and G. Saville, ‘‘Blowdown of Pipelines,’’ paper
Chen, C. R., S. M. Richardson, and G. Saville, ‘‘Modelling of Two- SPE 23070, Society of Petroleum Engineers Offshore Europe 91,
phase Blowdown from Pipelines: I. A Hyperbolic Model Based on Aberdeen, U.K., 369 Ž1991..
Vibrational Principles,’’ Chem. Eng. Sci., 50, 695 Ž1995.. Richardson, S. M., and G. Saville, ‘‘Blowdown of LPG Pipelines,’’
Chen, C. R., S. M. Richardson, and G. Saville, ‘‘Modelling of Two- Trans. Inst. Chem. Eng. Part B, 74, 235 Ž1996..
Phase Blowdown from Pipelines: II. A Simplified Numerical Saha, P., ‘‘Modelling the Dynamic Response of Emergency Shut-
Method for Multi-component Mixtures,’’ Chem. Eng. Sci., 50, 2173 down Valves Following Full Bore Rupture of Long Pipelines,’’ PhD
Ž1995.. Thesis, Univ. College London Ž1997..
Courant, R., K. O. Friedrichs, and H. Lewy, ‘‘Uberdie Partiellen Dif- Stittgen, M., and W. Zielke, ‘‘Fluid Structure Interactions in Flexible
ferenzialgleichungen der Mathematisehen, Physik,’’ Mathematische Curved Pipes,’’ Proc. of 6th Int. Conf. on Pressure Surges, A. R. D.
Annalen, 100, 32 Ž1926.. Thorley, ed., BHRA, 101 Ž1990..
Cullen, W. D., ‘‘The Public Inquiry into the Piper Alpha Disaster,’’ Weiss, M. H., K. K. Botros, and W. M. Jungowski, ‘‘Simple Method
Dept. of Energy, HMSO, London Ž1990.. Predicts Gas Blowdown Time,’’ Oil and Gas. J., 12, 55 Ž1988..
Flatt, R., ‘‘A Singly-iterative Second Order Method of Characteris- Wood, D. J., and J. E. Funk, ‘‘A Boundary Layer Theory for Tran-
tics for Unsteady Compressible One Dimensional Flows,’’ Comm. sient Viscous Losses in Turbulent Flow,’’ ASME J. Basic Eng., 92,
in Appl. Num. Meth., 1, 269 Ž1985.. 865 Ž1970..
Flatt, R., ‘‘Unsteady Compressible Flow in Long Pipelines Following Zucrow, M. J., and J. D. Hoffman, Gas Dynamics, Vols. I and II,
a Rupture,’’ Int. J. Num. Meth. Fluids, 6, 83 Ž1986.. Wiley, New York Ž1976..
Groves, T. K., P. R. Bishnoi, and J. M. E. Wallbridge, ‘‘Decompres-
sion Wave Velocities in Natural Gases in Pipelines,’’ Can. J. Chem.
Eng., 56, 664 Ž1978..
Hall, A. R. W., G. R. Butcher, and C. E. The, ‘‘Transient Simulation Manuscript recei®ed Dec. 7, 1998, and re®ision recei®ed Mar. 29, 1999.

AIChE Journal a990507--KB-155 June 1999 Vol. 45, No. 6 1201