Calculating The Release and Dispersion of Gaseous, Liquid and Supercritical CO
Calculating The Release and Dispersion of Gaseous, Liquid and Supercritical CO
Calculating The Release and Dispersion of Gaseous, Liquid and Supercritical CO
Calculating the Release and Dispersion of Gaseous, Liquid and Supercritical CO2
CHRIS DIXON
AND
MICHAEL HASSON
MMI Engineering Ltd., The Brew House, Wilderspool Park, Greenalls Avenue, Warrington, WA4 6HL
1. Introduction
The importance of reducing emissions of carbon dioxide is now widely accepted, with a number of industrial scale power generation projects now underway featuring carbon capture and storage (CCS). A number of the proposed plants utilize pre-combustion carbon capture, in which the fuel stream is rst decarbonised, with the subsequent carbon dioxide stream being transmitted elsewhere for sequestration in underground structures, such as depleting oil and gas reservoirs or saline aquifers. Inevitably, the large scale precombustion decarbonisation of fuel streams results in a corresponding requirement to process and transport safely bulk quantities of carbon dioxide. Within the UK health and safety legislative framework, quantitative risk assessment (QRA) is commonly used to demonstrate that the risks to plant personnel and wider communities are acceptable, or as low as reasonably practicable (ALARP). Central to the QRA process is the integration of the consequences of a range of postulated accident scenarios and corresponding event frequencies, from which a picture of the risk posed by a facility emerges. It can be seen that a reliable understanding of risk is heavily dependent on the quality of the assessment of accident consequences. As a family of technologies supporting low carbon power generation and carbon capture emerges, it is especially important that the risks posed by such new facilities are robustly understood and managed. This need is further reinforced by the fact that many major CO2 sources, which may be candidates for carbon capture, are located in areas of higher population density. This paper describes some of the challenges associated with predicting the consequences of accidental releases of carbon dioxide, under a range of thermodynamic conditions representative of those encountered in facilities designed to extract, process and transport large volumes of this potentially toxic, asphyxiating material. A number of techniques to address these challenges are presented, along with eorts made to validate the proposed methods against experimental data. The release scenarios addressed by these techniques include releases from gaseous, sub-cooled liquid, saturated liquid and supercritical vessels and the blow-down of pipelines tranporting liquid and supercritical CO2 . Methods to calculate the release rates from these scenarios will be presented, together with methods for calculating the subsequent gas dispersion using both integral methods for sub-sea and atmospheric dispersion and CFD modelling for the detailed investigation of dispersion around process plant.
2. Overview of Modelling
The principal concerns of this paper are how to model the release and subsequent dispersion of CO2 as an aid to understanding the hazards which may be presented by operational, emergency and accidental releases. Carbon dioxide is well known to be an
asphixiant, but is also a toxic substance with the toxic eect depending on the concentration. The concentrations can extend from levels which are widely accepted as posing an immediate risk of death, through to a level which is published as acceptable for exposure throughout a working day. Table 1 shows the broad levels of concentration with corresponding hazard levels. These hazards are generally presented in the form of contours or isosurfaces of CO2 concentration at values which relate to diering levels of hazard, such as long and short term exposure limits. In general it is ideal for models to be conservative, but not overly so, since highly conservative models can severely limit design constraints without need. Modelling for consequence analysis typically takes place in two stages. In the rst stage the release rate from a given inventory or scenario is calculated. In the second, the dispersion of the released CO2 is calculated. The level of complexity or diculty associated with each of these tasks depends on the particular scenario. For release rate calculation, the thermodynamic state of the inventory is a key parameter in determining this. This is also true for dispersion, where additional parameters which need to be considered in choosing a methodology include whether or not surrounding buildings need to accounted for and the length scales of interest. For example, if the requirement is to calculate low concentration level contours from a warm, vertical, low speed gaseous release in to a moderate wind with no inuence from surrounding buildings, then a simple Gaussian plume model may suce. On the other hand, for the release of a very cold gas and its interaction with the modules of an o-shore platform in a low wind, such a model would be entirely inappropriate and a computational uid dynamics (CFD) model would be needed. Cases intermediate between these two extremes should be treated individually. In some cases, such as the release of CO2 from a liquid inventory a high proportion of solids can be produced. In these cases, o the shelf codes may not be able to handle the physics so that either special procedures must be carried out (examples are given later) or a particular code version obtained.
allowing background plant operation to continue, whilst the trip or unavailability is addressed. During such operational or planned venting, the CO2 will normally be released in the gaseous phase. Usually, the vent will be located at a high point such as a stack at an on-shore site, or potentially on the are tower of an oshore CO2 injection, hydrocarbon production installation. Planned releases can occur at a variety of temperatures from close to the sublimation temperature up to over 100 C. Oshore, where the CO2 may be used for enhanced oil recovery (EOR) the gas may not be pure CO2 , but may contain hydrocarbon components. Venting can be designed to give either a sub-sonic release or a release which is sonic at the exit. In some cases the venting can occur over a long period of time in which case a steadystate calculation can be carried out. In other cases, though the ow rate may fall over the period during which the venting takes place, the rate at which it falls may be suciently small to allow a pseudo-transient calculation to be used. By this it is meant that a number of separate steady-state calculations are carried out, each corresponding to ow rates at xed points during the transient. For some calculation methods, such as CFD, this technique can substantially reduce the required computing time and should be used when practical. 3.2. Accidental Releases Accidental releases can occur from a huge variety of process conditions, including both gaseous and liquid inventories together with the supercritical uid state. In line with comparable oil and gas or process industry facilities, these releases will generally be above water either on-shore, or on the topsides of oshore platforms. However, sub-sea releases from pipelines or platform risers are examples of other major accident scenarios which would be considered. Gaseous releases from a hole in a vessel or broken pipe-work are similar to planned releases, though there are notable dierences: In general they will be sonic releases from pressurised vessels, plant or pipework. They will generally be relatively close to ground level (in the on-shore case), or associated with one of the process modules (in o-shore cases), rather than from a deliberately elevated vent stack. They may be of short duration and inherently transient if the vessel or process inventory is small, with factors such as emergency isolation and blowdown further complicating the character of the transient. Releases from vessels with liquid and certain supercritical inventories are signicantly more complex. In these cases, as the CO2 enters the atmosphere it makes a transition from a liquid to a two-phase gas/solid mixture where the solid fraction depends on the upstream conditions. During this transition the uid expands in a characteristic tulip shape. The solid particles which are formed then sublime back to a gas. A second, distinct, class of accidental release is that from a pipeline rupture. Normally, for the purpose of identifying bounding scenarios in a pipeline case for safely, a full-bore rupture is considered as this is the worst event in terms of peak ow rate. Again, while the case of a gaseous pipeline is relatively straightforward with analytical expressions in common use, the case of a pipeline transporting CO2 in supercritical (e.g. immediately downstream of compressors) or liquid (e.g. along the sea bed) states is signicantly more complex to calculate.
Release and Dispersion of CO2 4.1. Release from Vessels with Gaseous Inventory
This method works equally well for any upstream state and can therefore be used for the gaseous releases as well, and in that case is equivalent to the usual method of setting the Mach number at the orice to unity. For a vessel containing either a sub-cooled liquid or a supercritical uid with a thermodynamic state to the liquid side of the saturation dome, a typical thermodynamic path from the vessel conditions to the exit plane is shown by case A in gure 7. In this case the corresponding curve of ow rate against exit pressure has a sharp maximum in the curve which corresponds to the point where the isentrop crosses the saturation curve. The reason for the maximum being in this position is that a further drop in pressure along the isentrop results in a liquid/vapour two-phase state at the exit plane where the density rapidly falls due to the vapour portion. In some cases, a gaseous vessel condition may be close enough to the saturation line that
the isentropic path leads into the two phase region, Such a thermodynamic path shown by case B in gure 7. In this case, the reason for the path terminating on the saturation dome is more subtle, being due to a discontinuity in the gradient of the isentrops across the saturation line on the ph plane. It can be noted that this discontinuity is dicult to see in typical ph diagrams due to the logarithmic scale used for pressure. It can be noted that for a saturated liquid, where the thermodynamic state lies on the saturation dome (for example the bottom point of line A in gure 7), the modied Bernoulli equation is not valid, and indeed as this point is approached it must begin to break down. The normal form of the Bernoulli equation can still be used, but as shown in the example below it is very pessimistic. There is some uncertainty to the discharge coecient for the liquid case. For liquid discharge, a value of 0.62 is commonly used (Lees 1996). Since there is additional physics taking place here with ashing at the exit, it is not clear that this is the best value. Cumber (2002) suggests a value of 0.8 for both gaseous and two phase releases. For the releases considered in this section, the ow will generally be a saturated two phase mixture at, or close to, the outlet. Hence, if the thermodynamic state of the uid within the vessel is a sub-cooled liquid, a two-phase liquid/vapour mixture or a supercritical uid or a gas, it is suggested that the discharge coecient is chosen to be 0.8. Intermediate between the simple approach of using the Bernoulli equation with the exit pressure equal to the saturation pressure and the approach of maximising the ow rate with respect to the exit pressure is the -method (Fauske & Epstein 1988). In this method a correlating factor () for compressibility is used with an assumed equation of state. The advantages of the -method are that only upstream conditions are required and the process of varying the exit pressure required in the full method described above is not needed. In addition, because the thermodynamic data required by the method is less than in the full method, tabulated data can more readily be used. Furthermore, the method has several extensions, for example to allow it to deal with a release which takes place through a length of pipework. The main disadvantage of the -method is that it is only valid for pressures below the critical pressure. In order to compare the various methods described for sub-cooled and supercritical liquids, their predicitons are presented in gure 3 over a range of pressures. The graph shows mass ow rate per unit area as a function of pressure, without any discharge coecient applied. It can be seen that the straightforward application of the Bernoulli equation gives a signicantly higher ow rate than the other methods. On the other hand, all three of the other methods are in good agreement. It can be noted that these three methods are in good agreement for a reason, this being that they are doing essentially the same thing in the particular region of thermodynamic states under consideration well above the saturation dome. Figure 3 also shows a similar graph for the case of a saturated liquid inventory. In this case it is not possible to use a modied form of the Bernoulli equation as this would give a ow rate of zero. The remaining three methods are compared, and again the -method and the general method are found to be in good agreement, but the standard Bernoulli equation is higher by up to a factor of four. Hence it is clear that the Bernoulli equation should not be used in this case. 4.3. Transient Releases In order to calculate the transient owrate from a vessel, in addition to an outow model, a thermodynamics model and a vessel model are required. The overall method for the calculation is based on that described by Speranza & Terenzi (2005) though they consider only releases which are gaseous at the exit. In this work, the thermodynamics of the gas
has been represented using the Peng-Robinson equation of state (Poling et al. 2001). The eectiveness of this equation of state for carbon dioxide can be assessed from gure 4 which compares the output from this equation to data published by a major industrial gases supplier. The agreement is reasonable, and is certainly considered adequate for the purpose of major accident consequence modelling. In the vessel model, although time has no thermodynamic signicance and pressure stepping is often used (Mahgerefteh & Wong 1999), a time-stepping method is used here for convenience and though a small time-step is required run times are still short. Heat transfer is accounted for once the uid reaches a saturation condition where it is assumed that it occurs by boiling heat transfer. As a test of the method, calculated results are compared with the experiments of Gebbeken & Eggers (1996) in which pure carbon dioxide at initially supercritical conditions was released from a top vented lab-scale vessel. In these experiments, the test vessel was cylindrical with a volume of 0.05m3 , an internal diameter of 0.242m and a height-to-diameter ratio of 4.5. A venting pipe was connected to the top of the vessel with a cross-sectional area of 50mm2 and an orice was installed for varying the minimum cross-sectional area of the release. Releases were undertaken from vessel pressures of 150, 200, 250 and 300bar with initial uid temperatures of 298, 313 and 323K. The minimum release area was set to either 17mm2 or 50mm2 (i.e. the pipe area). Measurements are presented for vessel pressure, as well as vessel temperature and void fraction at various vertical heights within the vessel. Results are not presented for all combinations of initial pressure, temperature and orice size with pressure transients presented for a range of combinations, but temperature at only one. Since both pressure and temperature data are available, this is used as one validation test for the code. The comparison between the calculation and experiment for the vessel pressure and temperature for this case, which has initial conditions of 150bar and 313K and a 17mm2 orice are shown in gure 5. For this case the agreement is seen to be extremely good. A pressure trace is also provided by (Gebbeken & Eggers 1996) for a case with a 17mm2 area starting at 200bar and 313K and the same level of agreement is observed in that case (not shown here). The results for the cases with a 17mm2 orice are very good, but the results are not so good for the other cases presented by (Gebbeken & Eggers 1996), all of which have a minimum area of 50mm2 . This is the same area as the pipework, the length of which is not given, so that there is no orice as such in these cases. Hence, the discharge coecient is not neccessarily the same in these cases as in those with a 17mm2 orice. With this in mind, the pressure curves for the case already presented, together with two additional cases which have a 50mm2 mimimum area are shown in gure 6. The discharge coecient in these two cases is set to 0.6 and the results are reasonable. Without knowing the details of the pipework, it is dicult to know if this adjustment is justied. For calculations aimed at supporting cases for safety, it would typically be assumed that there was no connecting pipework as such information is often unavialble at the stage in a project when much consequence modelling needs to be undertaken. In that case the default discharge coecient of 0.8 would be used and the release ow would be pessimistic at least in terms of rate, if not duration.
5. Dispersion Calculation
In this section, examples are provided showing dispersion results for gaseous and liquid releases. These dispersion cases could be calculated in various ways, but in the cases shown here a computational uid dynamics (CFD) code has been used. There are advan-
tages and disadvantages to using a CFD code for this purpose. Advantages include the ability to deal with arbitrary terrain and buildings, and an ability to handle cases where the wind and the release direction are not aligned. The main drawback is that CFD is more numerically demanding than integral methods and therefore, in general, it takes longer to produce results. 5.1. Dispersion From Gaseous Releases For a rather hot gaseous release, such as may occur downstream of a rst compression stage, the CO2 can be a gas within the plant pressure boundary and remain a gas throughout the subsequent expansion and dispersion. In this case there is an underexpanded jet with a shock structure within it. This shock structure is at a very small scale compared to that of the overall dispersion, so that it cannot be resolved within a dispersion calculation which may require a domain of tens or hundreds of metres. The way to overcome this is to use a pseudo-source, some way downstream of the actual release point. One way do this is to carry out a separate CFD calculation to attempt to resolve the local structures within the under-expanded region, thus supplying source terms for a subsequent dispersion analysis. Alternatively, one can use a simpler method published by the HSE based on overall conservation laws (Ivings 2004). This latter approach constitutes an essentially standard approach, but also provides a prototype for the more complex case of a liquid release discussed in the next section. Figure 9 shows some typical results for a hypothetical site with a sloping terrain and a number of buildings and plant items on the site. The two rectangular blocks in the middle foreground represent items of plant which are relatively uncongested, which is to say that they consist of, for example, steel structure and pipework but are largely open. Hence the jet passes through this region relatively undisturbed while the other buildings are solid and are seen to strongly interact with the jet. In this case the wind direction is from the bottom right of the gure towards the top left and the plume is carried o in this direction. The results shown here are essentially illustrative, but the actual purpose of this particular calculation was to ivestigate the addition of a small quantity of easily detectable gas to the CO2 stream to study whether this could be used for detection beyond the region where CO2 concentration was dangerously high. The isosurface of larger volume in the gure shows the detectable limit of this gas while the smaller one shows a CO2 isosurface. 5.2. Dispersion From Liquid Releases In this section, an example is provided showing dispersion results for a liquid release. The dispersion could be calculated in various ways, but in this case a computational uid dynamics (CFD) code has been used. For a liquid release, the expected path on the thermodynamic pressure-enthalpy diagram from stagnation conditions to the exit plane is shown schematically by path A in gure 7. At the exit plane a saturated liquid is expected. Hence, beyond the exit the liquid will ash to a vapour whilst simultaneously cooling such that particles of solid carbon dioxide are formed. These particles will then sublime external to the vessel. This process is expected to be rather complex with three phases present, gas/liquid, liquid/solid and solid/gas phase changes occurring and very large gradients. In addition, the length and time scales are separated from those of the far-eld dispersion, so that including the detail of the region immediately adjacent to the point of release in a CFD model for the far eld dispersion is not considered viable. Hence, a simplied method for dealing with the region immediately adjacent to the point of release has been employed,
based on the work of Fauske & Epstein (1988). The rst stage is to nd the ow rate and the exit plane parameters as described in section 4.2. The depressurisation zone beyond the nozzle (an idealised point of release) is shown schematically in gure 7, where the jet is divided into two regions: a depressurisation zone and a two-phase entrainment zone. The inlet to the CFD model is taken to be at the end of the depressurisation zone where the pressure in the jet has reduced to atmospheric pressure. The following assumptions are made: It is assumed that there is no entrainment into the jet in the depressurisation zone, due to the pressure being above ambient in this region. Friction and heat transfer are neglected in the depressurisation zone. Conservation principles are then used to nd the jet properties at the end of the depressurisation zone, specically conservation of momentum ux, conservation of energy and conservation of mass. The following points should be noted: The sharp distinction between the depressurisation and entrainment zones in gure 7 is an idealisation, but in the absence of detailed information, this simplication is made. The method does not attempt to calculate the distance between the exit plane of the nozzle and the inlet plane to the CFD model: it is assumed that this distance is small compared to the scales over which dispersion is occurring. In order to dene a geometry within the CFD domain, it is assumed that the depressurisation tulip takes the shape of part of a sphere. In liquid release scenarios, a substantial fraction of the ejected CO2 can be predicted to be in the form of solids at the end of the depressurisation zone, and this should be accounted for by some means within the analysis. Several options can be considered to do this; the method chosen here was to track a scalar equation for the solids concentration. For example, if the total mass ow rate of gaseous CO2 into the domain amounted to, say, 70% of the mass ow rate the solids concentration will account for the rest, and this allows the scalar representing solids concentration to be set at the inlet. The equation for solids concentration is coupled to both the energy equation and the CO2 mass fraction equation via source terms, since as the solids concentration reduces due to sublimation, CO2 gas is produced and energy is absorbed. It can be noted that other approaches are possible. To demonstrate the level to which the modeling approach has been successful for a liquid release, some comparisons with large scale experiments, recently conducted by BP, are now presented. For reasons of commercial condentiality it is not possible to allow access to the experimental data at this time, so scales have been removed from the results. Although this may reduce the usefulness of the results, it nonetheless allows an assessment of the agreement between experiment and calculation to be made. The results presented are for the release of sub-cooled liquid CO2 at high pressure. In order to perform a comparison between the measured and calculated jet/plume shapes, a still photograph taken from a steady-state test showing the edge of the visible plume in side elevation for a particular test was digitized. It is not clear a priori which, if any, concentration level the visible plume corresponds to. In fact it may be a better match to some particular temperature level rather than a concentration level. It is not attempted here to make such a detailed comparison, but rather a pair of CO2 concentration levels are plotted, superimposed on the game graph as the visible plume and an overall comparison made.
It can be noted that the method employed to produce an inlet diameter for the effective release does not calculate the distance from the nozzle exit plane to this point. Hence, in order to allow the comparison of jet/plume shape to be made, the CFD and measured jet shapes are aligned in the horizontal direction such that their diameters match at the start of the predicted jet. The CFD calculation can then be considered to be accurate if the plume shape downstream of this point is well represented. The results of the comparison between CFD and measurement are shown in gure 8 which showsthe jet shape in side-elevation. The contour levels which have been extracted from the CFD results are at 0.1% and 1%, thus covering an order of magnitude in concentration level. Close to the release these levels are in very good agreement, while as one moves further downstream they part as would be expected. In the authors opinion, at the top of the plume the comparison is extremely favourable, whilst at the bottom, the comparison is not as favourable, but is still reasonable. It can be noted that the plume edge is less distinct at the bottom; i.e. if an attempt was made to place error-bars on the measured plume position they would be larger at the bottom. In addition, the distinct increase in spreading rate at the bottom of the jet which seems to occur downstream of the inlet may be caused by a measurement stand at this position interacting with the bottom of the plume. In addition to the measured shape of the plume close to the release point, concentration measurements were made along the centreline of the plume to a considerable downstream distance. A comparison between the measured and calculated concentration is shown in 10 and the agreement is found to be good being within the error bounds, apart from very close to the release. Finally a comparison can be made between the measured and calculated temperature values close to the release point which is done in gure 11. Here again the agreement is found to be reasonable.
7. Acknowledgments
The authors would like to thank BP for being allowed to take part in their recent CO2 release and dispersion validation exercise, thus gaining access to the experimental results presented here. The authors would also like to thank their collegues Dr Oliver Heynes and Dr Chris Robinson for their assistance.
10
Cumber, P. 2002 Modelling top venting vessels undergoing level swell. J. Hazardous Materials 89, 109125. Fauske, H.K. & Epstein, M. 1988 Source term considerations in connection with chemical accidents and vapour cloud modelling. J. Loss Prev. Process Ind. 1, 7583. Gebbeken, B. & Eggers, R. 1996 Blowdown of carbon dioxide from initially supercritical conditions. J. Loss Prev. Process Ind. 9 (4), 285293. Ivings, M. 2004 Outstanding safety questions concerning the use of gas turbines for power generation: Final report on the CFD modelling programme of work, Report CM/03/08 . HSE. Lees, F.P. 1996 Loss Prevention in the Process Industries, 2nd edn. Butterworth-Heinemann. Mahgerefteh, H. & Wong, S.M.A. 1999 A numerical blowdown simulation incorporating cubic equations of state. Computers and Chemical Engineering 23 (1999), 13091317. Poling, B.E., Prausnitz, J.M. & OConnell, J.P. 2001 The Properties of Gases and Liquids, 5th edn. McGraw-Hill. Speranza, A. & Terenzi, A. 2005 Blowdown of hydrocarbons pressure vessels with partial phase separation. Series on Advances in Mathematics for Applied Sciences 69 (2005), 508 521.
11
Figure 1. Typical thermodynamic paths from vessel (stagnation) to exit plane conditions for three vessel thermodynamic states. [htbp] le=,scale=0.8 Figure 2. default Figure 3. Mass ow rate per unit area as a function of pressure for sub-cooled liquid and supercritical inventories (left) and saturated liquid inventories (right). Flow rates per unit area are calculated using the Bernoulli equation, the Bernoulli equation assuming the saturated vapour pressure at the exit, the general method, and the omega method.
Figure 4. Applicability of the Peng-Robinson equation of state to CO2 by comparison with manufacturers data; lines are from the equation of state while points are manufacturers data. The left hand side shows isotherms on a pressure-enthalpy diagram while the left hand side shows saturated liquid and vapour densities.
Figure 5. Comparison between calculation and experiment for blowdown of a lab-scale vessel containing CO2 initially at 150bar and 313K with a 17mm2 orice. The left shows the pressure transient while the right shows the corresponding temperature; experiment shown by points, calculation by lines.
Figure 6. Comparison between calculation and experiment for blowdown of a lab-scale vessel containing CO2 initially at 150bar. Here only the pressure is shown for three cases. It should be noted that for the 7.98mm diameter (50mm2 area) cases the discharge coecient has been tuned (to 0.6 in both cases) as the pipe length is not known. The 4.65mm case has an orice and the discharge coecient is set to 0.8. Figure 7. Schematic diagram showing the expansion process during a ashing release. Figure 8. Comparison between measured and calculated jet/plume shape close to the release point for a liquid release. the edge of the visible plume is shown by the solid line. The calculated concentration contours at 0.1 and 1% are showns by the dashed lines (the 1% level forms the narrower band.) Figure 9. CO2 isosurfaces for a horizontal, sonic gaseous release on a hypothetical site. Figure 10. Comparison between calculated and measured CO2 concentrations downstream of a horizontal high pressure liquid release. Experimental data shown by points; calculated using CFD shown by line. Figure 11. Comparison between calculated and measured temperature downstream of a horizontal high pressure liquid release. Experimental data shown by points; calculated using CFD shown by line.