The Blast Os The Explosion

Broad Lane, Sheffield S3 7HQ Telephone: 0114 289 2000 Facsimile: 0114 289 2500

Explosion Hazard Assessment: A Study of the Feasibility and Benefits of Extending Current HSE Methodology to take Account of Blast Sheltering HSL/2001/04 Project Leader: Dr M Ivings Authors: Dr C Catlin, Dr M Ivings, Mr S Myatt Dr D Ingram, Prof D Causon, Dr L Qian Fire and Explosion Group

© Crown copyright 2001

Summary

Objectives / Background This study concerns the methodology used by MSDU, HSE in assessing Land Use Planning cases (LUP) near Hazardous Installations storing LPG and presenting a Vapour Cloud Explosion (VCE) hazard. The work is also relevant to all risk assessment and consequence models that are used to assess the hazard posed by blast. The methodology currently used by MSDU assumes that the blast propagates without interaction with any buildings or terrain features that may lie between the explosion source and the proposed site. Thus the estimate of the incident peak positive overpressure is likely to be conservative. The probability of fatality is then inferred from the peak positive overpressure using war-time data. This allows the vulnerability of the public to be assessed from which consultation distance and hazard zone boundaries are determined. The leading shock wave is the highest frequency component of the blast wave and will, therefore, be the most affected by interaction with obstacles. A revised methodology which takes account of obstacle interaction will, therefore, exhibit sheltering effects. Thus there are situations where buildings that lie between the explosion source and the proposed site will provide substantial sheltering in the region of the proposed construction. In such situations it is conceivable that the existing methodology places the hazard zone boundaries at substantially larger distances than is required to maintain an acceptable probability of fatality. A revised methodology, therefore, which takes sheltering into account could potentially free-up substantial areas of land around existing and future hazardous installations.

Main Findings Predictions of Computational Fluid Dynamic (CFD) simulations have been compared against a limited number of field-scale experiments which investigate the effects on the blast overpressure of the intervening obstacles. This has established that sufficient accuracy can be obtained using CFD methods without a prohibitive requirement of staff or computing resources. It has been shown that for simple linear rows of square cross-section buildings that three-dimensional CFD simulations provide a very accurate predictive method but require substantial computing and staff-time resources. However the application of a two-dimensional axi-symmetric approximation can also be used in the computations to yield adequate accuracy using contemporary workstations and a simulation turnover of less than a day. This has enabled the two-dimensional CFD to be used to simulate a range of cases for different building layouts and sizes. The peak positive overpressures determined from these simulations have provided sufficient data with which to propose a preliminary methodology for incorporating the effects of sheltering into existing explosion hazard assessment models. The methodology amounts simply to determining the TNT mass-scaled height of the buildings, and their intervening distances, and referring to look up tables to determine the relevant sheltering factors. These can then be used

to determine the inward shifts of the hazard zone boundaries. The preliminary methodology has been applied to two representative scenarios, the larger being on the scale of the Flixborough disaster, and has shown that the outer and inner boundaries of the Outer Zone can be moved in by as much as 80 m and 50 m respectively. These distances represent substantial free-up of land and therefore suggest that a revised methodology could significantly increase the land available for development around a major hazard site and/or lead to less restraints on the site owner.

Main Recommendations A preliminary methodology for incorporating the effects of sheltering into existing explosion hazard assessment models has been proposed. MSDU can now decide on the appropriateness of including the sheltering effects of buildings in their assessment of LUP cases. This will crucially depend on the fact that intervening buildings are not necessarily permanent and their removal will change the consultation distance. The further development of the methodology described in this report will lead to a useful tool that could be used in the preparation or assessment of COMAH safety reports where blast is a significant hazard. The development of the methodology should be based on a programme of two and three dimensional CFD simulations validated against further blast experiments.

Contents

1. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2. BLAST PREDICTION METHODS . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1. TNT Equivalence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2. TNO Multi-Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3. Consequence models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4. LPG-RISKAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3 3 3 5 6

3. METHODOLOGY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.1. Nominal experimental programme . . . . . . . . . . . . . . . . . . . . 7 3.2. Experimental method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.3. Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.4. Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.5. Numerical method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4. VALIDATION OF THE METHODOLOGY . . . . . . . . . . . . . . . . . . . 4.1. The appropriateness of axi-symmetry . . . . . . . . . . . . . . . . 4.2. Mesh refinement study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3. 2D Axisymmetric CFD calculations . . . . . . . . . . . . . . . . . . . 4.4. Parametric study CFD calculations . . . . . . . . . . . . . . . . . . .

24 24 27 29 40

5. SHELTERING MODEL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1. Extended TNT Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2. Two-dimensional CFD parameter variation study . . . . . . 5.3. Applicability of revised methodology . . . . . . . . . . . . . . . . .

41 41 42 53

6. CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 6.1. Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 7. REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

1.

INTRODUCTION

This study concerns the methodology used by MSDU, HSE in assessing Land Use Planning cases (LUP) near Hazardous Installations storing LPG and presenting a Vapour Cloud Explosion (VCE) hazard. This work is also relevant to all risk and consequence models that are used to assess the hazard posed by blast. Necessarily the methodology must be quick to use yet conservative. Current methods assume unobstructed terrain and yet the fatality statistics are largely based on war time data which would have involved the additional effects of blast structure interaction, structural damage and projectiles. The current methodology, therefore, can be considered as conservative. A revised methodology, which results in a substantial reduction in the consultation distances, would have large financial implications if blast waves were the only or dominant hazard at that distance. For example the sheltering effects of buildings sited on the hazardous installation could be taken into account and free up land space for development. For the reasons outlined above the methodology must be simple but versatile enough to apply to the many potential scenarios. The current methodology achieves this by assuming a flat unobstructed terrain between the blast source and the site. Thus the complexities of the effects of the intervening buildings are overlooked for the sake of versatility but with the assurance that the blast prediction is conservative. This is because for those instances when the blast front comprises a shock wave its interaction with buildings before arrival at the site will tend to remove energy from the blast wave and lessen the peak shock overpressure. The exceptional case might for example occur when the initial blast wave is preferentially aligned with the buildings so as to amplify the incident shock overpressure relative to the standard methodology. Such exceptions are presently addressed by basing the fatality data upon worst case data. The time is also opportune for the development of more sophisticated blast prediction methodologies since they at present lag behind those used for gas dispersion calculations. In many instances gas dispersion calculations are performed prior to a blast calculation in order to estimate the explosion inventory in the event of a flammable vapour release. Such dispersion methods cater for the improved mixing caused when a vapour disperses over built up terrain. These methods have been implemented into PC software in a form that they can be easily and quickly used by MSDU. This study necessarily considers the extension of the current TNT and TNO blast prediction methods. Both of these methods cater for ideal blast waves for which sheltering or amplification are most likely to occur. The shorter frequency component of the blast wave, namely the leading shock wave and immediate later times, is the most likely to be affected through interaction with buildings. The lower frequency components will be least affected. This study will be primarily concerned with peak incident overpressure in view of it being the principal measure of blast severity in the RISKAT methodology. One of the key issues concerns the representation of the complex variety of buildings and layout patterns. In the case of gas dispersion these can be adequately represented by roughness parameters that describe mean separation and building heights. In the case of explosions the scope for a simple approach is less obvious since we are studying the interaction of the higher frequency component of the blast wave which, therefore, is sensitive to the details of the geometry. Thus

1

the development of a comprehensive methodology necessarily requires the investigation of a large number of obstacle layouts and geometries and also a range of incident blast wave lengths. The effects of blast-structure interaction will be most pronounced when the blast wavelengths are similar in size to that of the building whereas the effects will be negligible for much larger wavelengths. In general there are many three dimensional layouts to be studied and in each case there can be a large surface area of terrain to search over to find the maximum incident overpressures. Thus an experimental programme alone has limitations in view of the need for the large number of blast measurements for each geometry studied as well as the large number of geometries. Computational Fluid Dynamic methods have an excellent track history for simulating ideal blast wave propagation and obstacle interaction but their accuracy can be questioned when applied for the first time to complex geometries. This question principally concerns the ability to capture the important physical features without the need for prohibitive computer and staff time. While in principle the methods can predict the near exact behaviour, in practice there are substantial constraints placed by the size and speed of the computer. A limited experimental programme was envisaged, therefore, to establish the sufficiency of the combined numerical and computational resources. The primary objective was to assess whether there was scope for developing a simple methodology. It was decided to investigate where significant effects of sheltering might be observed. Linear rows of buildings without gaps were therefore studied in the absence of obstructions to the down wind side: Depending upon their proximity to the hazardous installation, downwind structures could reflect and amplify the shock pressure incident upon the site. If, therefore, a valuable methodology could be distilled from these cases then this would provide a sound financial incentive for future more exhaustive investigation. This approach has the additional benefits that the obstacle layouts used in the experiments were simple and needed few changes. It also allowed the numerical representation to be axi-symmetric thus making it feasible to service the larger part of the study using standard workstations except for some isolated three-dimensional runs which were used to validate the axi-symmetric approximation.

2

2.

BLAST PREDICTION METHODS

2.1.

TNT Equivalence

The longest established methods for estimating incident blast overpressure are based upon high explosive (HE) blast. This approach is directly compatible with the estimates of damage which, having incorporated war-time data, are generally characteristic of high explosives. Because the source of energy is concentrated into a small volume, ideally in a mathematical sense to an infinitesimal point, the corresponding blast waves are considered ideal. The lack of geometrical scales results in the blast wave being characterised by only one effective length scale, namely the cube root of the charge mass. Thus the blast waves from different masses of TNT are all similar if the distance and time scales are normalized by the cube root of mass. This is the basis of Hopkinson scaling [1]. In practice the blast waves from a wide range of energetic explosions can be adequately characterised by that of the detonation of a suitable mass of TNT provided that the corresponding conversion efficiency of chemical energy into the blast wave energy is known. It is noted that many of the blast parameters, e.g. peak negative overpressure, positive and negative overpressure phase impulses and durations can be normalised using Hopkinson scaling and represented by one family of graphs [1]. It is important to note that in this study the TNT blast wave parameters correspond to ground reflected overpressure, namely that achieved when the given mass of TNT is placed on the ground and the overpressures measured at ground level. In many instances overpressures are quoted relative to free-air bursts, namely those not influenced by ground reflection, for which twice the mass of TNT is required to achieve the ground reflected overpressure. Clearly ground reflected overpressures are the more relevant to this application. In general the blast waves generated by explosions are not the same as those from HE. This is particularly true of Vapour Cloud Explosions (VCEs) in which the flame speeds are not sufficiently high to ensure a shock wave at the blast front. This limitation is addressed in the TNO Multi-Energy method described below. However for many vapour cloud explosions with high flame speeds [2] the blast from a VCE can be adequately determined by assuming a suitable conversion efficiency between the heat of combustion of the gas mixture and an equivalent mass of TNT [2]. 2.2.

TNO Multi-Energy

The Multi-Energy methodology [3] was developed in recognition that many explosions, particularly VCE’s where flame acceleration isinvolved, can vary substantially in the efficiency with which the chemical energy released by the combustion is converted into mechanical energy of the blast wave. Detonation, whether of gas or HE is the most efficient conversion process and the blast waves from detonations can therefore be accurately related to those of a suitable equivalent mass of TNT. In fact the Multi-Energy method categorises the explosion strength into ten levels, the peak incident overpressure and positive overpressure phase duration being given in Figures 2.1 and 2.2. Normalization of distance and time scales is here performed in terms of the combustion energy (E/P0)1/3 to allow one set of graphs to be used for different inventories of vapour cloud. HE blast corresponds to explosion strength 10. In fact it can be seen that in the range 70 mbar < Ps < 140 mbar there is effectively no difference in the peak incident

3

overpressure from explosion strengths 6 through 10. In the range 70 mbar < Ps < 140 mbar there is at most 20% underestimate in the positive phase duration. Explosion strengths 1 through 3 on the other hand are not sufficient to generate the lower peak overpressure threshold (70 mbar) of interest in this study. This leaves explosion strengths 4 and 5 neither of which is severe enough to produce shocks at the high level of relevance to this study, namely 140 mbar. Explosion strength 5 only acquires the ideal blast wave-form at peak positive overpressures substantially below 70 mbar. Outside the ideal blast wave region the Multi-Energy method assumes that the blast front has a finite rise time with total positive phase duration between 2 to 3 times longer that the equivalent HE blast wave. In view of the absence of the shock and the longer duration such blast waves will be far less susceptible to sheltering or amplification effects and have, therefore, not been considered in this study.

Dimensionless maximum side-on overpressure (Ps )

100.000

1 2 3 4 5 6 7 8 9 10 600 mbar 140 mbar 70 mbar

10.000

1.000

0.100

0.010

Ideal blast wave region

0.001 0.1

1

10

100

Combustion energy-scaled distance (R) Figure 2.1 TNO Multi-Energy methodology : Decay with distance of peak positive overpressure for ten explosion strengths

4

Dimensionless positive phase duration (Ts )

10

1 2 3 4 5 6

1

Ideal blast wave region

7 8 9 10 600 mbar 140 mbar 70 mbar

0 0.1

1 10 Combustion energy-scaled distance (R)

100

Figure 2.2 TNO Multi-Energy methodology : Variation with distance of positive overpressure phase duration for ten explosion strengths 2.3.

Consequence models

Consequence models are used by MSDU to assess Land Use Planning cases where sites are storing flammable substances and the dominant hazard is that of blast. These models may also be used by MSDU when assessing the predictive elements of COMAH safety cases where a ‘snapshot’ of the blast hazard is required. In Land Use Planning cases decisions are made based upon where the proposed site lies relative to the source of the explosion. The development is categorised according to the vulnerability of those persons that will be at the development. The advice is determined by a matrix of overpressure versus distance (Table 2.1). The matrix is divided into three zones (Inner, Middle and Outer) according to the peak incident shock overpressure (Ps). Beyond the radius of the Outer Zone, the so called Consultation Distance (CD) where Ps < 70 mbar, fatalities are unlikely. The Outer Zone corresponds to the overpressure range 70 mbar < Ps < 140 mbar where a 1% chance of fatality amongst the general public is expected. The Middle Zone is defined by the region 140 mbar < Ps < 600 mbar . It is assumed that there is a near certain chance of a fatality in the Inner zone, where Ps > 600 mbar .

5

Table 2.1 Decision function used in the consequence model for assessing Land Use Planning cases Inner Zone Middle Zone Outer Zone Ps > 600 140 < Ps < 600 70 < Ps < 140 Advise Against Category D Consult Consult Advise Against Don’t Advise Against Category A Consult Don’t Advise Against Category C Consult Consult Category B Don’t Advise Against Don’t Advise Against Don’t Advise Against The vulnerability axis is split into four categories: A - residential such as houses and hotels, B - industrial namely factories and offices, C - leisure/retail which includes shops, restaurants, sports centres and D - institutional namely hospitals and schools. A revised methodology which accounts for sheltering, for example, might propose a shift of the zone boundaries toward the source thus freeing up land for development. 2.4.

LPG-RISKAT

LPG-RISKAT [4] is used by MSDU when assessing applications for the development of new major hazard sites. If a risk based decision matrix is used then the matrix given in Table 2.2 is adopted. It is similar to that for the consequence model above and uses the same categories for vulnerability. Table 2.2 Decision matrix used in LPG-RISKAT (where DD is dangerous dose or worse and cpm is chances per million) Inner Zone Middle Zone Outer Zone risk > 10 cpm DD 1< risk <10 cpm DD risk < 1/3 cpm DD Category D Advise Against Consult Consult Category A Advise Against Don’t Advise Against Consult Category C Don’t Advise Against Consult Consult Category B Don’t Advise Against Don’t Advise Against Don’t Advise Against In the present case, where we are interested in blast effects, a dangerous dose is an overpressure of 140 mbar. In the case of fireball, the dangerous dose is 1000 tdu.

6

3.

METHODOLOGY

3.1.

Nominal experimental programme

As explained previously the terrain has been assumed flat in view of the relatively small distances (within a 500 m radius) between the centre of the blast source and the site in question. For example, if the TNT equivalent to the Flixborough disaster of 16 Te is assumed, the distances to the 600 mbar, 140 mbar and 70 mbar are given respectively by 105 m, 254 m and 413 m. The 600 mbar contour, therefore, is very close to the explosion source and the benefits of sheltering are far less in terms of the distances by which the zone boundaries could be moved closer. Thus the greatest significance of a revised methodology would be on the inner and outer boundaries of the Outer Zone since the distances involved are larger and, therefore, the area of land affected larger. As explained previously the generic scenario to be studied involves the placement, between the blast source and the site of interest, of unbroken rows of square cross-section obstacles. Obstacles on the downstream side of the site have been omitted since they could potentially counter any beneficial effect of sheltering. The generic geometries chosen were based upon terraced housing, either separated by back-toback gardens, or else an access street, and modern industrial layouts also separated by access routes. A survey of several industrial sites was conducted: Older industrial sites were commonly housed in square cross-section brick buildings typically cramped together. The layout was seldom alike the tidy linear rows in the tests but the gaps between buildings were on average one building height. The newer sites appeared to be constructed from generic building pods whose base width to height ratio can vary between 2 to 4. These were typically constructed on brown-field sites where the uncluttered start had enabled a regular building pattern, alike that of the tests. The layout on a site depends upon the floor space required by the businesses housed in each building and also the type of vehicle access required. Many small businesses in consecutive buildings all requiring articulated lorry access would result in the widest gaps between buildings (typically equal to the base width of the generic building pod). When a single business is housed and the vehicle access is from the narrower end of the building then the land space between buildings is less wide than the articulated lorry access described above. The limiting cases correspond to a large business requiring larger than normal floor space. In this case there is no gap and the building pods are attached to one another to form a building with a many to one base width to height ratio. The obstacles layout used in the study are shown in Figure 3.1. The proposed configurations, therefore, approximate the older industrial sites and also help to investigate the question of how a relatively narrow gap between rows affects the energy reflection. This configuration then approximates the newer industrial layouts of pods with the lowest base width to height ratio.

7

Figure 3.1 Layouts of obstacle rows used in field-scale experiments The generic obstacle layout is shown in Figure 3.2 and is characterised by the dimensions H, W, D, Xs and G.

W

G D

H

Xs

Figure 3.2 Generic geometry and definition of parameters It was important to specify the test conditions so as to identify a significant sheltering effect in the pressure range of relevance to the consequence model, i.e. 70 mbar to 600 mbar, and which could also be effectively accommodated on the HSL test site. Too large a charge could result in ground shock with disturbance to the neighbouring public. Too small a charge could result in

8

asymmetry of the blast wave and significant dissipation of the blast front energy through its interaction with the ground surface. Thus simple criteria were applied to size the charge weights and obstacles. Namely, that the wave length of the positive overpressure phase of the free air blast wave should be of a similar size to the height of the obstacles. Thus a range of HE charge weights were estimated, based upon the available obstacles (namely 0.6 m × 0.6 m square crosssection concrete blocks) that resulted in the appropriate blast wavelength at the 140 mbar and 70 mbar distances. The smallest physical height possible therefore was 0.6 m and multiples of 0.6 m achieved by placing the blocks next to and above one another. As explained previously this linear arrangement also had the advantage of simplifying the grid geometry used in the numerical calculations. The parameters used in the test programme are given in Table 3.1.

Test no. 1 2 3 4 5 6 3.2.

Table 3.1 Configuration parameters used in experimental programme Ps* H W/H Xs/M1/3 M Xs H/M1/3 (m) (m/kg1/3) (mbar) (kg) (m) (m/kg1/3) 0.46 7.7 9.97 140 0.6 1 0.78 0.31 11 16.29 70 0.6 1 0.89 1.09 10.3 10 140 1.2 1 1.17 0.73 14.6 16.22 70 1.2 1 1.33 0.73 14.6 16.22 70 0.6 2 0.67 1.09 10.3 10 140 0.6 2 0.58

Experimental method

The experiments were performed on the Blast Range facility of the Health and Safety Laboratory (HSL) at Buxton, Derbyshire by Explosives Section staff. The Section has great experience in measuring overpressures from a wide variety of explosive events[5,6]. The range, 11 m wide and 112 m long, had a rolled hardcore surface which provides a flat area on which the concrete pendine block obstructions could be built, and which reduced the likelihood of spurious shockwave reflections from reflecting surfaces. PE4, a plastic explosive with a TNT equivalence of approximately 1.3 [7], was used for all the tests. The explosive was formed into a hemisphere, placed on a mild steel plate (500x500x28 mm), and detonated remotely using an L2A1 detonator inserted into the top of the explosive (Figure 3.3).

9

Figure 3.3 Arrangement of explosive system The pendine blocks (600x600x1800 mm) were arranged into obstructions as shown in Figure 3.1. Each obstruction extending over the middle 9 m of the range width for all tests. Blast waves were measured using up to ten Meclec FQ11C piezo-electric gauges (resonant frequency 80 kHz) mounted in B12 baffles (Figure 3.4), positioned as shown in Figures 3.1 and 3.5. The gauges were mounted on steel supports with bases that were resting on 30 mm thick expanded foam sheets in order to reduce the effects of ground shockwaves which could interfere with shockwaves propagated through the air. The explosives were positioned centrally with the obstructions on one side. This allowed comparison of the interaction of the blast wave with the obstructions and with those on the ‘clear field’ side while using the same explosive source. A test using only a detonator was also performed to determine its contribution to the shockwave; this was found to be small.

10

Figure 3.4: Photograph of Meclec FQ11C transducer in B12 baffle

Figure 3.5 File allocations for gauges used during tests Generally, the gauges were positioned along the central axis of the range at the half height of the obstructions. However, a gauge at each end of the measured range was also positioned at 1.5 times the height of the obstructions along the central axis, and two other gauges, one at each end, were positioned in an offset position 2.125 m from the central axis. In some instances 11 transducers were required when only 10 were available. In such cases the offset transducer at the furthest point from the explosion on the clear field side was removed to a more useful position. Since complex blast wave interactions were not expected on the clear field side of the explosion, and 2 transducers would remain at that distance (on the central axis), it was

11

considered that this would not compromise the results. A photograph of a typical experimental set-up (Test 1) is shown in Figure 3.6.

Figure 3.6 Photograph of typical experimental arrangement (Test 1) Because the number of gauges varied depending on the test being performed, the allocated gauge numbers were not the same in all cases. Figure 3.5 shows the relative positions of the gauges and their allocated numbers (their file identification codes are also provided). The gauge closest to the explosion (Gauge 5) was positioned perpendicular to the central axis and was used to trigger the data logging system. The distance of the gauge from the centre of the explosive charge varied between tests as shown in Table 1. Table 3.2 Distance of trigger gauge from centre of explosion Test No Distance (m) 1&2 + Detonator test 1.25 3&4 1.00 5&6 2.00 Electrical signals from the gauges were amplified using Kistler 5011 charge amplifiers (having a 200kHz frequency range), and recorded on a Nicolet 500 series datalogger. 12 bit samples were taken at a sampling rate of 500kHz. The mean sensitivity of each transducer was programmed into its associated charge amplifier so that the voltage readings obtained from all the transducers were directly comparable. All of the equipment was calibrated and quoted results are traceable to national standards where appropriate.

12

Each test was performed in triplicate (except the detonator-only test) and the data was processed using FAMOS software, a package designed to manipulate waveforms to enable quantitative data to be obtained. Where possible the following parameters were obtained for each trace: é Shock wave rise time é Peak pressure é Impulse é Time of peak pressure é Positive phase duration é Negative phase duration é Peak negative pressure The traces from the obstructed side of the explosion did not mirror the clear field traces due to complex wave interactions. This made it more difficult to obtain accurate measurements of the parameters listed above. Generally, where broad pressure pulses were observed with sharper pulses superimposed on the broader peak, the peak pressure and its time were assigned to the largest sharp peak. 3.3.

Experimental results

Examples of shock wave traces from the clear field and obstructed sides of the explosion are shown in Figures 3.7 and 3.8, respectively (1 mV = 0.33 mbar). The repeatability of the three firings for each of the tests was found to be high. Mean values averaged over the three firings for each of the parameters measured in the tests are summarised in Tables 3.3 to 3.9. Aaa005a Aaa001a

Aaa002a

Aaa003a

Aaa004a

V 7.0

TEST 1. Firing 1. (17/3/00) Clear field side of explosion

Gauge 5 1.25m from charge (30cm high) Rise Time: 67microseconds Max pressure: 2278mbar Time of max pressure: 0.060ms

6.5 6.0

320.7g PE4 initiated by L2A1 detonator Triggered by arrival of blast wave at first transducer (1.25m from explosion)

5.5 5.0 4.5 Gauge 3 4.7m from charge (30cm high) Rise Time: 67microseconds Max pressure: 395.0mbar Time of max pressure: 7.99ms

4.0 3.5 Gauge 4 2.9m from charge (30cm high) Rise Time: 60microseconds Max pressure: 672.2mbar Time of max pressure: 3.26ms

3.0 2.5

Gauge 1 7.7m from charge (90cm high) Rise Time: 70microseconds Max pressure: 164.7mbar Time of max pressure: 16.1ms

2.0 Gauge 2 7.7m from charge (30cm high) Rise Time: 63microseconds Max pressure: 157.1mbar Time of max pressure: 16.0ms

1.5 1.0 0.5 0.0 -0.5 -1.0 -5

0

5

10

15

20

25 ms

Figure 3.7 Typical traces from clear field side of explosion (Test 1)

13

Aaa005b Aaa001b

Aaa002b

Aaa003b

Aaa004b

mV 700

Gauge 6 2.9m from charge (30cm high) Max pressure: 218.5mbar Time of max pressure: 3.96ms

650 600

TEST 1. Firing 1. (17/3/00) Obstructed side of explosion 320.7g PE4 initiated by L2A1 detonator

550 500

Triggered by arrival of blast wave at first transducer (1.25m from explosion)

Gauge 7 4.7m from charge (30cm high) Max pressure: 100.7mbar Time of max pressure: 12.2ms

450 400

Gauge 8 7.7m from charge (90cm high) Max pressure: 78.21mbar Time of max pressure: 17.4ms

350 300 250 200 150

Gauge 9 7.7m from charge (30cm high) Max pressure: 50.82mbar Time of max pressure: 18.2ms Gauge 10 7.7m from charge (30cm high, Offset Max pressure: 50.82mbar Time of max pressure: 18.7ms

100 50 0 -50 -100 -150 -200 -250 -300 -350 -400 0

5

10

15

20

25

30 ms

Figure 3.8 Typical traces from obstructed side of explosion (Test 1) The blast traces for some of the gauges had perturbations in the baseline prior to the arrival of the air blast. These perturbations were not always observed on the clear field side and are often more prominent in data from the obstructed side of the explosion, Figure 3.9 shows a typical example. These perturbations are likely to be due to ground shockwaves passing through the ground at a faster speed than the air blast.

14

Aal003a

Aal004b

mV 280

mV 280

260

260

240

240

Gauge 3 Clear field trace

220

220

200

200

180

180

160

160

Gauge 9 Obstructed field trace

140

140

120

120

100

100

80

80

60

60

40

40

20

20

0

0

-20

-20

-40

-40

-60

-60

-80

-80

-100

-100

-120

-120 0

10

20

30

40

50

60 ms

Figure 3.9 Comparison of blast traces from the clear and obstructed side of Test 4. (Gauges 14.6 m from the explosion) 3.4.

Discussion

The data obtained indicates that complex interactions occur between the air shockwave and the obstructions. Gauges positioned at the half height of the pendine block obstructions tended to record much lower pressure values than gauges which were positioned above them. The signals from the obstructed side were also much smaller than those recorded on the clear field side of the tests. These results indicate that the obstacles create a significant sheltering effect. The presence of signals before the arrival of the air blast indicates that shockwaves are being transmitted through the ground. This occurs at a faster rate than shockwave propagation in the air. The ground shock signals are small in comparison to the air shockwaves on the clear field side of the tests, particularly close to the explosion, but can become significant when the air blast is small i.e. when ‘sheltering’ of the gauge has occurred. Where comparable ground shock data are available on the clear and obstructed side of a test, it appears that the baseline perturbations due to ground shock are decaying by the time that the air blast arrives at the gauge (Figure 3.10). The clear field trace is typical of results likely to be produced by an air blast i.e. the shock front gives the characteristic step rise in pressure followed by a sharply decaying back slope. These observations indicate that data on the interaction of the air blast with the obstructions is unlikely to have been compromised by the ground shock effect acting directly on the gauges.

15

Aal004a mV 340 320 300 280 260 240 220 200 180 160 140 120 100 80 60 40 20 0 -20 -40 -60 -80 -100 -120 -140 -160

Aal001b mV 340 320 300 280 260 240 220 200 180 160 140 120 100 80 60 40 20 0 -20 -40 -60 -80 -100 -120 -140 -160

Gauge 4 Clear field trace Gauge 6 Obstructed field trace

0

10

20

30

40

50

60

70

80 ms

Figure 3.10 Comparison of blast traces from clear and obstructed sides of Test 4. (Gauges 11.0 m from explosion point) Since direct ground shock effects do not appear to affect the data significantly in the region of interest (ie. after the arrival of the air blast), the more complex traces from the obstructed side of the tests are likely to be caused by the interaction of the air blast with the obstructions ie. reflections between blocks and diffraction from the block edges. These interactions are likely to cause some retardation of the blast propagation which is seen as a later arrival time of the air blast on the obstructed side compared to that on the clear field side. There is a possibility that ground shock could transmit into the pendine blocks which could propagate the shock wave back into the air at their surfaces. The waves generated in this manner are likely to occur in the area of interest but the magnitude of their contribution is not known. The large difference between the density of air and the concrete pendine blocks suggests that energy transfer between the two would be inefficient and would result in small shock waves. However, the air waves generated by this process will be produced over a large flat surface which may have the effect of amplifying the wave amplitude over the distances involved in the tests. No experimental work has been done to quantify these effects.

16

Gauge

10 1 8 2 9 3 7 4 6 5

Gauge

10 1 8 2 9 3 7 4 6 5

Table 3.3 Results from Detonator test (n/m = not measurable) Positive phase parameters Negative phase parameters Rise time Peak Time to Impulse Duration Peak Duration (ms) pressure peak (kPa.s) (ms) pressure (ms) (kPa) pressure (kPa) (ms) n/m n/m n/m n/m n/m n/m n/m 48 0.59 102.70 n/m n/m n/m n/m n/m 0.26 103.50 n/m n/m n/m n/m 45 0.63 102.70 n/m n/m n/m n/m n/m n/m n/m n/m n/m n/m n/m 63 1.35 90.10 n/m n/m n/m n/m n/m n/m 94.80 n/m n/m n/m n/m 49 1.91 84.80 n/m 0.67 0.86 1.90 213 0.43 85.90 n/m n/m n/m n/m 72 8.45 0.07 n/m 0.51 3.68 3.78

Table 3.4 Mean results from Test 1; Mean mass of PE4 = 320.2g Positive phase parameters Negative phase parameters Rise time Peak Time to Impulse Duration Peak Duration (ms) pressure peak (kPa.s) (ms) pressure (ms) (kPa) pressure (kPa) (ms) 816 5.15 18.70 0.02 6.56 3.41 53.56 63 15.31 15.98 0.02 3.07 4.51 13.57 69 7.85 17.38 0.01 6.38 2.42 29.04 69 15.71 16.12 0.02 3.16 4.40 13.00 752 5.08 17.97 0.02 6.60 2.53 40.95 67 38.18 8.00 0.03 2.44 9.01 14.42 3452 9.97 12.24 0.03 5.93 6.65 41.44 59 72.80 3.28 0.05 2.00 12.32 12.93 142 21.22 4.01 0.05 5.17 12.65 39.39 57 296.80 0.05 0.07 0.96 23.10 12.15

17

Gauge

10 1 8 2 9 3 7 4 6 5

Gauge

1 7 2 8 3 9 4 6 5

Table 3.5 Mean results from Test 2 ; Mean mass of PE4 =182.5g Positive phase parameters Negative phase parameters Rise time Peak Time to Impulse Duration Peak Duration (ms) pressure peak (kPa.s) (ms) pressure (ms) (kPa) pressure (kPa) (ms) 599 2.41 28.25 0.01 5.66 1.54 32.70 61 7.06 25.60 0.01 3.06 2.20 11.61 65 3.04 27.39 0.01 5.71 1.32 40.67 65 7.33 25.55 0.01 2.60 2.86 10.08 575 2.67 28.05 0.01 5.74 1.32 25.47 357 15.44 13.74 0.02 2.51 4.51 10.24 3566 4.85 18.76 0.01 6.32 3.45 58.57 62 22.47 8.53 0.02 2.20 4.84 9.43 3467 8.51 13.53 0.02 5.15 6.27 54.33 36 236.48 0.03 0.06 0.97 17.93 11.37

Table 3.6 Mean results from Test 3 ; Mean mass of PE4 = 757.0g Positive phase parameters Negative phase parameters Rise time Peak Time to Impulse Duration Peak Duration (ms) pressure peak (kPa.s) (ms) pressure (ms) (kPa) pressure (kPa) (ms) 63 15.53 22.49 0.03 4.35 6.50 12.36 3 3.78 28.98 0.02 9.16 2.97 14.84 67 15.14 21.88 0.03 4.23 6.70 12.00 3 3.69 26.35 0.02 10.85 3.08 12.91 17 17.85 21.93 0.03 4.37 6.75 11.62 65 10.66 23.54 0.02 10.25 3.98 14.09 60 29.40 12.09 0.04 3.50 8.36 9.41 21 9.79 16.96 0.04 8.75 11.69 13.69 5 913.33 0.03 0.16 1.27 28.38 18.63

18

Gauge

1 7 2 8 3 9 4 6 5

Gauge

1 7 2 8 3 9 4 6 5

Table 3.7 Mean results from Test 4 ; Mean mass of PE4 = 431.0g Positive phase parameters Negative phase parameters Rise time Peak Time to Impulse Duration Peak Duration (ms) pressure peak (kPa.s) (ms) pressure (ms) (kPa) pressure (kPa) (ms) 71 7.47 35.24 0.01 3.85 1.71 8.42 3 2.43 39.83 0.01 8.24 1.73 15.32 69 7.61 35.12 0.01 3.27 2.81 8.20 3 2.23 39.16 0.01 11.25 1.12 10.28 14 8.98 35.20 0.01 4.03 3.91 9.27 104 5.16 36.24 0.01 6.01 1.49 12.87 61 11.11 25.78 0.02 4.27 3.18 8.56 3 4.24 29.70 0.02 8.77 4.52 11.77 7 713.90 0.00 0.12 1.10 29.26 10.87

Table 3.8 Mean results from Test 5 ; Mean mass of PE4 = 431.0g Positive phase parameters Negative phase parameters Rise time Peak Time to Impulse Duration Peak Duration (ms) pressure peak (kPa.s) (ms) pressure (ms) (kPa) pressure (kPa) (ms) 66 6.79 34.32 0.01 3.49 1.97 12.29 3 3.38 35.92 0.01 6.28 2.27 12.25 68 7.65 34.31 0.01 3.34 2.60 8.38 3 0.40 10.52 <0.01 4.38 1.39 13.56 35 8.01 34.29 0.01 4.03 3.85 9.82 59 5.03 34.46 0.01 5.89 2.35 12.09 72 10.10 24.92 0.02 4.18 3.63 9.13 34 6.59 27.56 0.02 6.17 4.85 9.55 21 212.73 0.05 0.07 1.59 23.03 13.05

19

Table 3.9 Mean results from Test 6 ; Mean mass of PE4 = 757.0g (n/m = not measurable) Positive phase parameters Negative phase parameters Gauge Rise time Peak Time to Impulse Duration Peak Duration (ms) pressure peak (kPa.s) (ms) pressure (ms) (kPa) pressure (kPa) (ms) 1 74 18.21 21.11 0.03 4.31 4.76 14.69 7 3 4.74 24.03 0.02 7.19 3.91 13.77 2 69 16.34 20.58 0.03 4.30 7.01 11.68 8 3 1.08 50.74 n/m n/m 1.39 n/m 3 47 16.85 20.56 0.03 3.96 6.01 12.68 9 65 9.82 22.16 0.02 6.42 5.52 35.38 4 71 29.72 11.01 0.04 3.43 8.15 13.88 6 60 13.57 15.53 0.04 5.85 11.12 11.48 5 43 250.31 0.06 0.10 1.48 24.62 12.77

3.5.

Numerical method

The hydrocode in this study has been under continuous development at Manchester Metropolitan University (MMU) since the early 1990’s and is based on modern high resolution, finite volume, time marching methods. The hydrocode was originally developed using total variation diminishing (TVD) schemes developed for aeronautical applications [8]. The code has been updated to utilise solution methods based on MUSCL reconstruction [9] and approximate Riemann solvers [10,11] which more accurately represent the physics of wave propagation. The fully three dimensional/axi-symmetric hydrocode uses a structured body fitted grid topology, with a sub grid scale modelling capability, models chemistry using reduced reaction kinetics mechanisms [12,13] and solves either the inviscid Euler equations or the compressible Favre averaged form of the Navier-Stokes equations [14,13]. The code has been used to study a number of industrial blast wave and explosion hazard problems including the impact of the detonation of an explosive charge on a French blast shelter [15,8], the effect of blast wave impact on a Tank farm [16] and blast diffraction in a twin side-by-side jet engine configuration [17]. In addition the code has been used for contract and consultancy work for the Defence Evaluation and Research Agency (DERA) and British Aerospace, plc and has attracted funding from the Research Councils [18]. The code currently reflects the state of the art in blast and shock wave capturing methods and has been proven to be accurate and robust. Because of the simplicity of the geometry in the present study a rectangular, uniformly spaced, Cartesian grid has been used throughout and in most cases axi-symmetry has been assumed to minimise the computational costs. Furthermore, the inviscid, non-reactive, Euler equations have been solved since viscous effects are expected to be negligible in the region of interest and since high explosive charges are being used, gaseous detonation models are not required. The initial conditions for the blast wave have been obtained using the balloon analogue [19,20]. This approach has been chosen since a complete blast calculation requires a reactive-flow computation for the detonation process within the solid charge, an approach which is both computationally expensive and only warranted where detailed data is required in the near field region. In the balloon analogue, initial flow conditions for a statically pressurised gas are specified within a

20

balloon of arbitrary shape. The blast propagation is initiated by the rupture of the balloon, causing the pressurised gas to expand outwards forming a shock wave while a rarefaction travels inwards giving a blast wave like decay in the pressure behind the initial shock, see Figure 3.11. Control of the shape of the balloon and the internal conditions allows adjustment of the initial conditions to match a wide range of actual blasts waves in the mid to far-field. Figure 3.12 shows the typical blast wave profile generated by the balloon analogue. The blast wave and its decay are generated by the shock and rarefaction waves which form when the balloon ruptures, while the subsequent reshock is generated by the reflection of the rarefaction waves from the centre of the balloon. The magnitude of the blast wave, the positive and negative phase durations, and the magnitude of the reshock are all determined by appropriate choices of parameters. The energy potential of a balloon filled with ideal gas is given by [20]

E=

(P−P 0 )V c−1

(3.1)

where E is the total blast energy of the charge (in J), P is the initial static pressure within the balloon (in Pa), Po is the ambient pressure (in Pa), V is the volume of the balloon (in m3), and γ is the ratio of specific heats for the gas inside the balloon. The values of the free parameters P, V and γ, at t=0, together with the density of the gas inside the balloon, ρ (in kgm−3) are selected to give the best representation of the real blast wave. It is important to note that although the shape of the balloon should be similar to the shape of the explosive charge, there is no requirement for the volume of the balloon, V, to be equal to the charge volume. To simplify parameter selection in the present study, the balloon pressure (P) is set to the detonation pressure of TNT (21.0 GPa) and the ratio of specific heats γ is set to that for air, (1.4). The volume of the charge can then be calculated using equation (3.1), leaving either the density (ρ) or temperature (T) inside the balloon as the sole free parameter. In order to select an appropriate density for the gas inside the balloon, a series of numerical simulations starting with varying initial density has been compared with the overpressure measured on the unobstructed side of the experiments. Figure 3.13 shows the pressure traces obtained at Gauge 4 (2.9 m from the charge) with ρ = 8, 11 and 16 kg m−3 respectively for a charge weight of 0.461 Kg (TNT equivalent). All three values of density give good agreement with both the peak overpressure and the decay rate in the post shock expansion. The main effect of varying density is to alter the magnitude and arrival time of the reshock. Higher air density in the balloon results in lower pressures in the negative phase and a stronger, more delayed, reshock. As a result of these comparisons, an initial balloon ρ = 11 kg m−3 has been adopted as the optimised value for the following calculations. It is important to note that unless the model includes a detailed gas dynamic representation of the flow conditions and chemical composition of detonation products within the balloon, accurate prediction of the arrival time and magnitude of the reshock is impossible. The reshock is a small feature in comparison to the main blast wave and will have a small effect on the overpressure and impulse loading at the points of interest.

21

Initial Balloon Shock Wave

Rarefaction

Figure 3.11 Balloon analogue: Waves generated after the initial rupturing of the balloon

Initial Shock P

Reshock

Positive Phase Duration

Negative Phase Duration

t

Figure 3.12 Balloon analogue: Typical blast wave profile generated by the balloon rupture

22

Gauge 4, DX=DY=0.01 2.5 Experiment 8.0 11.0 16.0

Overpressure(V)

2

1.5

1

0.5

0

-0.5 0

0.005

0.01

0.015 Time(s)

0.02

0.025

0.03

Figure 3.13 Balloon analogue: Effect of the initial density of the balloon gas, ρ=8, 11 and 16 kg m−3

23

4.

VALIDATION OF THE METHODOLOGY

4.1.

The appropriateness of axi-symmetry

In order to be able to simulate the large number of cases required by the parametric study it is necessary for the computation to be as efficient as possible. This efficiency would be easily obtainable if 2D axi-symmetric calculations can be shown to give sufficient agreement for the centreline pressure measurements. In order to ascertain whether the axi-symmetric assumption is appropriate a three dimensional calculation has been performed and comparisons made both with experimental measurements and a two dimensional axi-symmetric simulation of the same case. The configuration used is that from Test 1 (see Figure 4.1). In the 3D simulation, the symmetry along the centreline of the experiment has been exploited to reduce the required mesh size. The charge is placed on the symmetry (Z=0) plane. In the simulation each of the Pendine blocks is 4.25 m long and 0.6 m wide with a square cross-section. The inter-block spacing is 1.2 m. Pressure monitoring points are placed 0.3 m above the floor at the mid-points of the cavities (0.6 m from each wall) on the symmetry plane (Z=0.0 m) and off the axis (Z=2.125 m). Comparisons between the centre line and off axis predictions and the experimental measurements have been made. The computational mesh (8 m long, 4 m high and 6 m wide) contained 3 million grid cells with equal mesh spacing, ∆X=∆Y=∆Z=0.04 m. Because of the computational cost associated with the 3D calculations, it has been necessary to obtain temporary access to a high performance α-powered workstation. The simulation was run on a 600 MHz DEC-α RISC processor with 512 Mbytes of memory and took 23 hours 30 minutes to compute. In addition to gaining access to the α workstation additional staff effort was required to port and optimise the code for the α processor. Pressure histories at six measuring points as shown in Figure 4.1 are given in Figure 4.2. These Figures show a great deal of similarity between the centreline (Z=0.0 m) and off axis (Z=2.125 m) histories. The only difference is a slight time lag (around 1.2 ms), which results from the differences in distance from the centre of the charge and the monitoring point. The numerical predictions for peak overpressure and positive phase impulse agree to approximately 2%. Gauges 8 and 9, in the experiment, are at the same locations as the last two monitoring points and a similar time lag can be observed. The agreement between the experimental measurements and simulated pressure histories is encouraging. It is worth noting that the oscillations which occur prior to the arrival of the incident shock in Gauges 8 and 9 are due to the detection of the ground wave, which is not present in the simulations. A further comparison between the 3D results, already discussed, and a 2D axi-symmetric calculation (on an equivalent mesh) has been performed in order to establish equivalence between the centre line and axi-symmetric monitoring points. Table 4.1 shows the peak overpressure and impulse from both the 2D and 3D calculations together with percentage errors. This shows that the errors are acceptable. Figure 4.3 compares the pressure histories at Gauges 6, 7 and 9 from both the 2D and 3D simulations. This Figure shows good agreement between the 2D and 3D simulations during the positive phase of the blast wave, while in the negative phase some discrepancies start to become apparent. Since the main focus of this study is the effects of the initial blast wave, these discrepancies are acceptable

24

Table 4.1: Comparison between 2-D and 3-D calculated results 2-D 3-D Gauge No. Pmax(Err%) TPmax(Err%) Imp.(Err%) Pmax TPmax Imp. mbar 6 7 9

ms

mbar s

mbar ms mbar s

146.2(1.7%) 5.41(0.0%) 0.42(2.4%) 148.8 5.41 0.41 73.9(3.4%) 12.76(0.7%) 0.26(3.7%) 71.3 12.66 0.27 42.0(8.4%) 18.10(0.2%) 0.13(7.1%) 45.9 18.06 0.14

y

z

Pendine blocks

4.25m charge 2.125m 0.6m x 0.6m

monitoring points

Figure 4.1 3D Test: Geometrical Layout

25

3-D, test 1, gauge 6, 200x100x150

3-D, test 1, gauge 7, 200x100x150

0.8

0.4 Z=0.000m Z=2.125m Experiment Overpressure(V)

Overpressure(V)

0.6

Z=0.000m Z=2.125m Experiment

0.3

0.4

0.2

0

-0.2

0.2 0.1 0 -0.1 -0.2

-0.4

-0.3 0

0.005

0.01

0.015 0.02 Time(s)

0.025

0.03

0.035

0

0.005

0.01

0.015 Time(s)

0.02

0.025

0.03

3-D, test 1, Gauge 9 & 10, 200x100x150 0.2 Z=0.000m Z=2.125m Gauge 9 Gauge 10

Overpressure(V)

0.15 0.1 0.05 0 -0.05 -0.1 -0.15 0

0.005

0.01

0.015 0.02 Time(s)

0.025

0.03

0.035

Figure 4.2 3D Simulation: Comparison between centre line and off axis pressure histories

26

Test 1, Gauge 6, DX=DY=0.04

Test 1, Gauge 7, DX=DY=0.04

0.5

0.25 2-D 3-D

0.4

0.15 Overpressure(V)

0.3 Overpressure(V)

2-D 3-D

0.2

0.2 0.1 0 -0.1

0.1 0.05 0 -0.05 -0.1

-0.2

-0.15

-0.3

-0.2

-0.4

-0.25 0

0.005

0.01

0.015 Time(s)

0.02

0.025

0.03

0

0.005

0.01

0.015 Time(s)

0.02

0.025

0.03

Test 1, Gauge 9, DX=DY=0.04 0.15 2-D 3-D

Overpressure(V)

0.1

0.05

0

-0.05

-0.1

-0.15 0.005

0.01

0.015 0.02 Time(s)

0.025

0.03

Figure 4.3 Comparison between 2-D and 3-D pressure histories 4.2.

Mesh refinement study

In order to examine the grid convergence of the numerical scheme a systematic mesh refinement study has been conducted using four grids, with mesh spacings of ∆X=∆Y=0.04 m, ∆X=∆Y=0.02 m, ∆X=∆Y=0.01 m and ∆X=∆Y=0.005 m. Each of the grids was used to simulate the unobstructed side of test 1. Figure 4.4 shows the pressure history obtained on each grid at the monitoring point furthest from the point of detonation. These results show the expected steepening and gain in amplitude of the incident blast wave with reducing mesh spacing. This convergence can be better demonstrated by using the grid convergence index (GCI). If we choose the maximum positive overpressure Ps as the required indicator of the grid convergence, then given the Ps values at the same gauge location on a fine grid (PS1) and a coarse grid (PS2), the GCI can be obtained from, e

GCI 21 = 3 r p21−1

(4.1)

27

where p is the formal order of accuracy of the numerical method (2 in the present case), r is the mesh refinement factor

r=

DX 2 DX1

(4.2)

and ε21 is a measure of the percentage relative error between the two solutions at the given location, i.e.

e 21 = 100

Ps 1 −Ps 2 Ps 1

(4.3)

Table 4.2 shows the GCI for the three gauge positions used in the Test 1. The GCI estimates an error of about 5% in the peak overpressure between the finest two grids, thus indicating that acceptable accuracy can be obtained on mesh ∆X=∆Y=0.01 m.

Table 4.2: Maximum overpressure values (in mbar) at three gauge points on four grids with the Grid Convergence Index (GCI) at each location Gauge No. Ps4 Ps3 Ps2 Ps1 GCI43 GCI32 GCI21 2 121.17 144.82 158.82 167.67 16.3 8.82 5.27 3 199.74 239.00 261.09 276.73 16.42 8.46 5.65 4 461.31 556.80 619.44 656.70 17.15 10.112 5.67 ∆X, ∆Y (m)

0.04

0.02

0.01

0.005

In order to generalise these findings it is useful to non-dimensionalise the mesh spacing with respect to the nominal blast wave length, λ. This length scale is related to the TNT equivalent charge weight, W through a 1/3rd power scaling, i.e. λ = W1/3 Figure 4.5, shows the peak overpressure, Ps, at the three gauge locations plotted against the nondimensionalised grid spacing λ/∆x. This graph indicates that acceptable ( ≈ 5%) accuracy will be obtained for 100 < λ/∆x < 150. Finally, it is important to note that increasing the spatial accuracy of the grid will also increase the temporal accuracy of the integration scheme. This is a direct consequence of the CourantFriedrichs-Lewy (CFL) stability condition which states that  ∆X ∆Y ∆t < min ,  max ij uij + cij max ij vij + cij 

   

(4.4)

where ∆t is the time step, ∆X and ∆Y are the grid density, uij and vij are the velocity components in the x and y coordinate directions and cij is the local sound speed.

28

Grid Convergence, Test 1, Gauge 3 300 0.04 0.02 0.01 .005

Overpressure(mbar)

250 200 150 100 50 0 -50 -100 -4

-2

0

2 4 Time(ms)

6

8

10

Figure 4.4 Mesh refinement: pressure history at Gauge 3 (Test 1) on meshes with spacing ∆X=∆Y=0.005, ∆X=∆Y=0.01, ∆X=∆Y=0.02 and ∆X=∆Y=0.04

Grid Convergence Study 700 Gauge 2 Gauge 3 Gauge 4

Overpressure(mbar)

600

500

400

300

200

100 0

50

100 150 lambda/DX

200

250

Figure 4.5 Mesh refinement: Peak overpressure vs. nondimensionalised mesh spacing 4.3.

2D Axisymmetric CFD calculations

All the two dimensional calculations have been carried out on a SGI workstation and the calculations require about 0.11 ms per time step per control volume. Thus, for a typical case using 800×400 grid points and 2000 time steps, the computational time would be around 22 hours. During the calculation, the time history of the overpressure at selected gauge points are recorded

29

at each time step. The flow variables such as pressure P, density ρ, temperature T, and total energy E etc. for the full domain are also recorded at a regular interval and post-processing software developed at MMU is used to draw the contour plots. This provides a clear image of the position and strength of the shock waves and a useful mechanism to analyse the underlying physical processes of the blast-obstacle interactions. In all the two dimensional calculations the flow is assumed to be axi-symmetric about the vertical axis and comparisons are made with experimental measurements from the central plane of the experiment (Figure 4.6). y Upper boundary (transmissive)

Symmetry Axis

Downrange boundary (transmissive) Pendine blocks

Charge

Ground Surface

r(x)

Figure 4.6 Set-up of the numerical modelling

Table 4.3: Experimental test configurations Test No. M Xs No. of H W G (kg TNT) (m) obstacles (m) (m) (m) 1 0.461 7.7 3 0.6 0.6 1.2 2 0.308 11.0 3 0.6 0.6 1.2 3 1.094 10.3 2 1.2 1.2 1.2 4 0.729 14.6 2 1.2 1.2 1.2 5 0.729 14.6 2 0.6 1.2 1.2 6 1.094 10.3 2 0.6 1.2 1.2 Detailed comparisons have been made between experiments and simulations for six test cases (Table 4.3 shows the main parameters of each test). In all the experiments pressure histories were recorded at various locations on both the obstructed and unobstructed sides of the charge. Figures 4.7-4.12 show simulated and experimental pressure histories at various gauges. Tables 4.4-4.9 compare the peak pressure and positive phase impulse recorded at the same gauges. These tables also show the time of arrival of the pressure wave and the percentage errors between the simulated and experimental peak pressures and impulses. In examining these results it is useful to note that Gauges 2, 3, 4 and 5 are on the unobstructed side of the charge while

30

Gauges 6, 7, 8 and 9 are on the obstructed side. Difficulties with Gauge 8 during Tests 5 and 6 prevents meaningful comparisons to be made at this station for these cases. One trend that can be seen in all the comparisons is that the errors observed at Gauge 5 tend to be quite high, this is a direct result of the proximity of Gauge 5 to the charge and the consequent strength of the pressure pulse. Two problems arise - firstly, the blast wave at this location has a very small wavelength and is not well resolved by the grid. Increased mesh resolution would improve this measurement but would dramatically increase the computational costs without any significant improvement of the results in the medium to far field. Secondly, the pressure transducers used in the experiment have a finite response time and are unlikely to measure the peak overpressure accurately in this case. Due to the short wavelength of the blast wave, missing the peak by a fraction of a millisecond can result in a dramatic reduction in pressure. In all the Figures the experimental and numerical pressure histories closely agree in terms of underlying trends. This is particularly clear in Tests 1 and 2 (Figures 4.10 and 4.11) where Gauges 6 and 7 show quite complex patterns of multiple spikes during the initial blast wave. These spikes are caused by the blast wave reverberating in the street canyons between the pendine blocks. The reverberation effects can lead (see Gauge 7) to peak pressure being observed after the initial blast wave. The peak pressures in this case are predicted to around 10% by the numerical method (see Table 4.4 and 4.5). In all the experiments, but particularly in the latter cases, the measurements from the pressure transducers have been contaminated by ground shocks. This causes oscillations in the pressure trace prior to the arrival of the blast wave. In these cases comparisons between the CFD predictions (which assumes a perfectly rigid ground) and the experiments must be treated with caution. However, there is a high level of agreement between the CFD and the experiment in many cases on both the clear field and obstructed sides of the explosion, for example see the Gauge 6 results in Figures 4.10 and 4.11. These results give strong support to the earlier argument that the effect of the ground shock has decayed before the arrival of the air blast. Furthermore, there is sufficient experimental data, where ground shocks do not interfere with the experimental data, that are in good agreement with the CFD predictions to give confidence in the numerical simulations. Thus indicating that the 2D axi-symmetric simulations will provide a reasonable basis for the development of a sheltering methodology.

31

Test 1, Gauge 2

Test 1, Gauge 3

0.5

1.2 Simulation Experiment

Simulation Experiment

1 Overpressure(V)

Overpressure(V)

0.4 0.3 0.2 0.1 0 -0.1

0.8 0.6 0.4 0.2 0 -0.2

-0.2

-0.4 0

0.005

0.01

0.015 0.02 Time(s)

0.025

0.03

0

0.005

0.025

0.03

14 Simulation Experiment

Simulation Experiment

12 Overpressure(V)

2 Overpressure(V)

0.015 0.02 Time(s)

Test 1, Gauge 5

Test 1, Gauge 4 2.5

1.5 1 0.5 0

10 8 6 4 2 0 -2 -0.005

-0.5 0

0.005

0.01

0.015 0.02 Time(s)

0.025

0.03

0

Tset 1, Gauge 6

0.005

0.01 0.015 Time(s)

0.02

0.025

Test1, Gauge7

0.8

0.4 Simulation Experiment

Simulation Experiment

0.3 Overpressure(V)

0.6 0.4 0.2 0 -0.2 -0.4

0.2 0.1 0 -0.1 -0.2

-0.6

-0.3 0

0.005

0.01

0.015 0.02 Time(s)

0.025

0.03

0

0.005

0.01

0.015 0.02 Time(s)

0.025

Test 1, Gauge 9 0.2 Simulation Experiment

0.15 Overpressure(V)

Overpressure(V)

0.01

0.1 0.05 0 -0.05 -0.1 0

0.005

0.01

0.015 0.02 Time(s)

0.025

0.03

Figure 4.7 Test 1: Computational and Experimental pressure histories.

32

0.03

Test 2, Gauge2

Test 2, Gauge 3

0.25

0.5 Simulation Experiment

Simulation Experiment

0.4 Overpressure(V)

Overpressure(V)

0.2 0.15 0.1 0.05 0 -0.05

0.3 0.2 0.1 0 -0.1

-0.1 0.01

-0.2 0.015

0.02 0.025 Time(s)

0.03

0.035

0

0.005 0.01 0.015 0.02 0.025 0.03 0.035 Time(s)

Test 2, Gauge 4

Test 2, Gauge 5

0.8

9 Simulation Experiment

0.7

7 Overpressure(V)

0.6 Overpressure(V)

Simulation Experiment

8

0.5 0.4 0.3 0.2 0.1 0

6 5 4 3 2 1

-0.1

0

-0.2 0

-1 -0.01-0.005

0.005 0.01 0.015 0.02 0.025 0.03 0.035 Time(s) Test 2, Gauge 6

0.005 0.01 0.015 0.02 0.025 Time(s) Test 2, Gauge 7

0.3

0.15 Simulation Experiment

0.25

Simulation Experiment

0.1 Overpressure(V)

0.2 0.15 0.1 0.05 0 -0.05 -0.1 -0.15

0.05 0 -0.05 -0.1 -0.15

-0.2 -0.25

-0.2 0

0.005 0.01 0.015 0.02 0.025 0.03 0.035 Time(s)

0

0.005 0.01 0.015 0.02 0.025 0.03 0.035 Time(s)

Test 2, Gauge 9 0.1 Simulation Experiment

0.08 Overpressure(V)

Overpressure(V)

0

0.06 0.04 0.02 0 -0.02 -0.04 0

0.005 0.01 0.015 0.02 0.025 0.03 0.035 Time(s)

Figure 4.8 Test 2: Computational and Experimental pressure histories

33

Test 3, Gauge 2

Test 3, Gauge 3

0.5

0.6 Simulation Experiment

Simulation Experiment

0.5 Overpressure(V)

Overpressure(V)

0.4 0.3 0.2 0.1 0 -0.1

0.4 0.3 0.2 0.1 0 -0.1

-0.2

-0.2 0.0160.018 0.02 0.0220.0240.0260.028 0.03 Time(s)

0.0160.018 0.02 0.0220.0240.0260.028 0.03 Time(s)

Test 3, Gauge 4

Test 3, Gauge 5

1

40 Simulation Experiment

Simulation Experiment

35 Overpressure(V)

Overpressure(V)

0.8 0.6 0.4 0.2 0

30 25 20 15 10 5

-0.2

0

-0.4 0.005

0.01

0.015 0.02 Time(s)

0.025

-5 -0.01

0.03

0

0.01 0.02 Time(s)

Test 3, Gauge 6

0.04

Test 3, Gauge 8

0.4

0.15 Simulation Experiment

0.3

Simulation Experiment

0.1 Overpressure(V)

0.2 0.1 0 -0.1 -0.2

0.05 0 -0.05 -0.1

-0.3 -0.4

-0.15 0

0.01

0.02 0.03 Time(s)

0.04

0.05

0

0.01

0.02 0.03 Time(s)

0.04

Test 3, Gauge 9 0.4 Simulation Experiment

0.35 0.3 Overpressure(V)

Overpressure(V)

0.03

0.25 0.2 0.15 0.1 0.05 0 -0.05 -0.1 -0.15 0

0.01

0.02 0.03 Time(s)

0.04

0.05

Figure 4.9 Test 3: Computational and Experimental pressure histories.

34

0.05

Test 4, Gauge 2

Test 4, Gauge 3

0.25

0.3 Simulation Experiment

Simulation Experiment

0.25 Overpressure(V)

Overpressure(V)

0.2 0.15 0.1 0.05 0 -0.05

0.2 0.15 0.1 0.05 0 -0.05 -0.1

-0.1

-0.15 0

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s)

0

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(S)

Test 4, Gauge 4

Test 4, Gauge 5

0.35

25 Simulation Experiment

Simulation Experiment

20

0.25

Overpressure(V)

Overpressure(V)

0.3

0.2 0.15 0.1 0.05 0

15 10 5 0

-0.05 -0.1 0

-5 -0.01

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s)

0

Test 4, Gauge 6

0.03

0.04

Test 4, Gauge 8

0.2

0.08 Simulation Experiment

0.15

Simulation Experiment

0.06 Overpressure(V)

0.1 0.05 0 -0.05 -0.1

0.04 0.02 0 -0.02 -0.04

-0.15 -0.2

-0.06 0

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s)

0

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s)

Test 4, Gauge 9 0.2 Simulation Experiment 0.15 Overpressure(V)

Overpressure(V)

0.01 0.02 Time(s)

0.1

0.05

0

-0.05 0

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s)

Figure 4.10 Test 4: Computational and Experimental pressure histories.

35

Test 5, Gauge 2

Test 5, Gauge 3

0.25

0.25 Simulation Experiment

Simulation Experiment

0.2 Overpressure(V)

Overpressure(V)

0.2 0.15 0.1 0.05 0 -0.05

0.15 0.1 0.05 0 -0.05 -0.1

-0.1

-0.15 0

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s)

0

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s)

Test 5, Gauge 4

Test 5, Gauge 5

0.3

8 Simulation Experiment

0.25

Overpressure(V)

Overpressure(V)

0.2 0.15 0.1 0.05 0 -0.05

6 5 4 3 2

-0.1

1

-0.15

0

-0.2 0

Simulation Experiment

7

-1 -0.005

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s) Test 5, Gauge 6

0.005 0.01 Time(s)

0.015

0.02

Test 5, Gauge 9

0.2

0.15 Simulation Experiment

0.15

Simulation Experiment

0.1

0.1

Overpressure(V)

Overpressure(V)

0

0.05 0 -0.05 -0.1

0.05 0 -0.05 -0.1

-0.15 -0.2

-0.15 0

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s)

0

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s)

Figure 4.11 Test 5: Computational and Experimental pressure histories.

36

Test 6, Gauge 2

Test 6, Gauge 3

0.5

0.5 Simulation Experiment

Simulation Experiment

0.4 Overpressure(V)

Overpressure(V)

0.4 0.3 0.2 0.1 0 -0.1

0.3 0.2 0.1 0 -0.1

-0.2

-0.2 0

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s)

0

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s)

Test 6, Gauge 4

Test 6, Gauge 5

1

8 Simulation Experiment

Simulation Experiment

7 Overpressure(V)

Overpressure(V)

0.8 0.6 0.4 0.2 0 -0.2

6 5 4 3 2 1 0

-0.4 0

-1 -0.01

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s) Test 6, Gauge 6

0.01

0.02 0.03 Time(s)

0.04

0.05

Test 6, Gauge 9

0.5

0.3 Simulation Experiment

0.4

Simulation Experiment

0.25 0.2

0.3

Overpressure(V)

Overpressure(V)

0

0.2 0.1 0 -0.1 -0.2

0.15 0.1 0.05 0 -0.05 -0.1

-0.3

-0.15

-0.4

-0.2 0

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s)

0

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time(s)

Figure 4.12 Test 6: Computational and Experimental pressure histories.

37

Table 4.4: Comparison Between Computational Results and Experiments measurements for Test 1 Simulated Experimental Error Gauge No. Pmax TPmax Imp. Pmax TPmax Imp. Pmax Imp. (mbar) (ms) (Pa s) (mbar) (ms) (Pa s) (%) (%) 2 138.7 15.88 20.34 157.1 16.12 19.04 11.7 6.8 3 289.9 7.73 33.07 381.8 8.00 34.56 24.0 4.3 4 660.9 3.22 50.43 728.0 3.28 50.07 9.2 0.7 5 4214.1 0.06 94.04 2968.0 0.05 69.77 42.8 34.7 6 188.1 4.03 46.9 212.2 4.01 45.9 11.3 2.1 7 93.3 11.95 30.3 99.7 12.24 29.4 6.4 3.0 9 49.7 17.75 16.1 50.8 17.97 14.2 2.2 11.8 Table 4.5: Comparison Between Computational Results and Experiments measurements for Test 2 Simulated Experimental Error Gauge No. Pmax TPmax Imp. Pmax TPmax Imp. Pmax Imp. (mbar) (ms) (Pa s) (mbar) (ms) (Pa s) (%) (%) 2 65.1 25.97 9.91 70.6 25.60 8.97 7.7 10.4 3 125.6 14.10 15.98 154.4 13.74 17.16 18.6 6.8 4 195.9 9.11 21.58 224.7 8.53 21.24 12.8 1.6 5 2739.6 0.08 67.21 2364.8 0.03 61.70 15.8 8.9 6 74.4 10.92 22.3 85.4 13.53 23.0 12.8 3.0 7 48.6 17.91 12.4 48.5 18.76 16.1 0.2 22.9 9 28.8 27.14 8.23 26.7 26.73 7.7 7.8 6.8 Table 4.6: Comparison Between Computational Results and Experiments measurements for Test 3 Simulated Experimental Error Gauge No. Pmax TPmax Imp. Pmax TPmax Imp. Pmax Imp. (mbar) (ms) (Pa s) (mbar) (ms) (Pa s) (%) (%) 2 133.0 21.68 263.80 151.4 21.88 259.16 12.1 1.7 3 132.0 22.06 255.47 178.5 21.93 284.13 26.0 10.0 4 249.2 11.86 400.73 294.0 12.09 400.95 15.0 0.0 5 11410.0 0.032 2090.3 9133.3 0.03 1630.2 24.9 28.2 6 106.3 16.13 483.41 97.9 16.96 440.55 8.6 9.7 8 40.3 25.02 184.07 36.9 26.35 193.05 9.2 4.6 9 75.9 22.44 131.51 106.6 23.54 231.88 28.7 43.2

38

Table 4.7: Comparison Between Computational Results and Experiments measurements for Test 4 Simulated Experimental Error Gauge No. Pmax TPmax Imp. Pmax TPmax Imp. Pmax Imp. (%) (mbar) (ms) (Pa s) (mbar) (ms) (Pa s) (%) 2 67.0 34.72 141.54 76.1 35.12 119.68 12.0 18.3 3 66.0 34.99 140.46 89.8 35.20 134.09 26.5 4.8 4 98.3 24.50 187.57 111.1 25.78 173.80 11.5 7.9 5 6654.8 0.069 1194.7 7139.0 0.00 1213.3 6.8 1.5 6 50.8 28.54 239.05 42.4 29.70 183.37 19.8 30.4 8 22.8 37.75 104.30 22.3 39.16 139.04 2.2 25.0 9 39.6 35.09 77.60 51.6 36.24 74.91 23.3 3.6

Table 4.8: Comparison Between Computational Results and Experiments measurements for Test 5 Simulated Experimental Error Gauge No. Pmax TPmax Imp. Pmax TPmax Imp. Pmax Imp. (mbar) (ms) (Pa s) (mbar) (ms) (Pa s) (%) (%) 2 66.0 33.72 141.64 76.5 34.31 119.02 13.7 19.0 3 66.0 33.77 141.36 80.1 34.29 135.41 17.6 4.4 4 98.7 23.47 187.81 101.0 24.92 179.96 2.3 4.4 5 1844.7 0.056 835.55 2127.3 0.05 713.35 13.3 17.1 6 57.8 26.85 239.29 65.9 27.56 220.55 12.3 8.5 9 36.3 33.86 131.41 50.3 34.46 102.52 27.8 28.2

Table 4.9: Comparison Between Computational Results and Experiments measurements for Test 6 Simulated Experimental Error Gauge No. Pmax TPmax Imp. Pmax TPmax Imp. Pmax Imp. (mbar) (ms) (Pa s) (mbar) (ms) (Pa s) (%) (%) 2 132.0 20.95 265.01 163.4 20.58 262.24 19.2 1.1 3 132.0 21.04 263.95 168.5 20.56 281.82 21.7 6.3 4 249.8 10.10 402.07 297.2 11.01 403.92 15.9 0.5 5 2607.0 0.126 1059.6 2503.1 0.06 957.22 4.2 9.6 6 136.6 15.24 497.70 135.7 15.53 393.80 0.7 26.4 9 72.6 21.31 237.98 98.2 22.16 183.37 26.1 29.8

39

4.4.

Parametric study CFD calculations

In order to provide data for the parametric study a series of CFD calculations were carried out giving rise to peak pressures of 140 mbar and 70 mbar at the point of interest. For both peak pressures the number of obstacles upstream of the point of interest has been varied between none, one and two. In all cases an array of computational pressure transducers has been used to give measurements over a range of normalized distances (-2/M1/3....+2/M1/3) around the point of interest. The array extends vertically from ground level up to 1 or 2 building heights (0 < Y < 2H) depending on the case. For each of these locations the peak positive and negative over pressures, the positive and negative impulses and the positive and negative phase durations have been computed. The resulting parametric data was post processed and presented in the form of a series of spreadsheets.

40

5.

SHELTERING MODEL

5.1.

Extended TNT Methodology

The approach toward developing a revised model must necessarily be compatible with the existing TNT and TNO models. This warrants representing geometrical features in terms of normalized dimensions, such as the building height, the distance from the obstacles to the 140 mbar or 70 mbar point and, for more than one row of obstacles, the gap distance between the obstacles. Thus dimensions have been normalized against the cube root of the ground reflected TNT charge mass. This is merely a convention which could be readily converted to other forms of normalization such as the combustion energy scaled distance used in the TNO Multi-Energy methodology. These normalized distances are interrelated via appropriate constant multiplying factors. It becomes possible, therefore, to determine the inward displacements of the normalized zone boundaries from the sheltering parameter, here defined as the ratio of the sheltered and free air incident peak positive overpressure (Ps/Ps*). The more effective the sheltering, the smaller Ps and the closer the zone boundary can be moved toward the explosion source. The relationships between the shifts in the 140 mbar and 70 mbar boundaries and the sheltering parameter are given approximately in Figure 5.1. These are calculated by assuming that the overpressure at the shifted zone boundary, when multiplied by the sheltering factor, gives the nominal zone boundary peak overpressure (i.e. 140 mbar or 70 mbar). As previously explained the same degree of sheltering gives rise to a larger inwards shift of the 70 mbar boundary compared to that of the 140 mbar boundary.

1.0

Sheltering factor (Ps /Ps *)

0.9 140 mbar 70 mbar

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0

2

4

6

8

10

12

14

16

18

1/3

Normalized distance (Xs /M )

Figure 5.1 Zone boundary shifts corresponding to sheltering factors

41

20

5.2.

Two-dimensional CFD parameter variation study

The scope of the two-dimensional parameter variation study is shown in Table 5.1. Two generic building geometries have been considered, namely one or two building-rows of square crosssection buildings. In the case of two building-rows they have been assumed to be one buildingheight apart. This corresponds to an intermediate level of separation between buildings as found in both housing and industrial estates (see Section 3.1). Two placements for the site have been considered namely distances of two building-heights (D = 2H) and four building heights (D = 4H) downwind of the nearest building row. This was done to establish to what extent sheltering was confined to the immediate down wind vicinity of the buildings. It was also necessary to investigate the effect of the ratio of building-height to the blast wavelength. This was done by varying the TNT mass-normalized building heights between 0.5 and 1.0. In the limited scope offered by this study these values were chosen to correspond to representative building heights (10 m) for inventories similar to the Flixborough disaster. Calculations were done for all eight of these cases for the two boundaries of the Outer Zone, namely incident peak overpressures of 140 mbar and 70 mbar respectively, thus making sixteen simulations in all. Table 5.1 Configuration parameters for the two-dimensional CFD calculations Rows Ps* W/H H/M1/3 G/H D/H 1 2 140 70 1 0.5 1.0 N/A 1 2 4

The accuracy of the numerical calculations were maintained by selecting the same computational cell dimensions, when normalized against the cube root of charge weight. In addition to maintaining accuracy this had the additional benefit of enabling the pressure monitoring positions to be placed at the same normalized distances relative to the nominal zone boundary. Pressure-time variation was, therefore, calculated in all cases on the same array of forty five cells whose centres extended horizontally ±2 m/kg1/3 to either side of the zone boundary and from ground level to two building heights (0 < Y < 2H).

42

Normalized positive peak overpressure (Ps /Ps *)

2.0 Free air 140 mbar

Obstructed 70 mbar

1.5 Height range {0
1.0

0.5

0.0 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

Figure 5.2(a) Maximum positive overpressure (0
43

1/3

Normalized positive phase impulse (Is /[Ps *M ])

0.003

0.002

0.001

Height range {0
Obstructed 70 mbar

0.000 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

1/3

Normalized positive phase duration (Ts /M )

Figure 5.2(b) Maximum positive phase impulse (0
0.008

0.006

0.004 Height range {0
Obstructed 70 mbar

0.000 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

Figure 5.2(c) Maximum positive phase duration (0
44

has been associated with a reduction in the probability of fatality. It is important to note, however, that should positive phase impulse be used instead there would be a much reduced beneficial effect of sheltering.

Normalized positive peak overpressure (Ps /Ps *)

A repeat of the above analysis was also performed for the case of two building-rows (See Figures 5.4(a-c) and Figures 5.5(a-c). Again, the incident peak overpressure is substantially reduced and the magnitude of sheltering found to be greater than for the case of 1 building-row. Corresponding observations were made for positive phase impulse and duration. 2.0 Free air 140 mbar

Obstructed 70 mbar

1.5 Height range {0
1.0

0.5

0.0 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

1/3

Normalized positive phase impulse (Is /[Ps *M ])

Figure 5.3(a) Maximum positive overpressure (0
0.002

0.001

Height range {0
Obstructed 70 mbar

0.000 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

Figure 5.3(b) Maximum positive phase impulse (0
45

1/3

Normalized positive phase duration (Ts /M )

0.010

0.008

0.006

0.004 Height range {0
Obstructed 70 mbar

0.000 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

Normalized positive peak overpressure (P s /Ps *)

Figure 5.3(c) Maximum positive phase duration (0
2.0 Free air 140 mbar

Obstructed 70 mbar

1.5 Height range {0
1.0

0.5

0.0 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

Figure 5.4(a) Maximum positive overpressure (0
46

1/3

Normalized positive phase impulse (Is /[Ps *M ])

0.003

0.002

0.001

Height range {0
Obstructed 70 mbar

0.000 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

1/3

Normalized positive phase duration (Ts /M )

Figure 5.4(b) Maximum positive phase impulse (0
0.010

0.008

0.006

0.004 Height range {0
Obstructed 70 mbar

0.000 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

Figure 5.4(c) Maximum positive phase duration (0
47

Normalized positive peak overpressure (Ps /Ps *)

2.0 Free air 140 mbar

Obstructed 70 mbar

1.5 Height range {0
1.0

0.5

0.0 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

1/3

Normalized positive phase impulse (Is /[Ps *M ])

Figure 5.5(a) Maximum positive overpressure (0
0.003

0.002

0.001

Height range {0
Obstructed 70 mbar

0.000 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

Figure 5.5(b) Maximum positive phase impulse (0
48

1/3

Normalized positive phase duration (Ts /M )

0.010

0.008

0.006

0.004 Height range {0
Obstructed 70 mbar

0.000 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

Figure 5.5(c) Maximum positive phase duration (0
49

Normalized peak positive overpressure (Ps /Ps *)

D=2H ; Max{0
D=2H ; Max{0
1.50 ______ Free air

1.25 1.00 0.75 0.50 0.25 0.00 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

Normalized peak positive overpressure (Ps /Ps *)

Figure 5.6(a) Sheltering factors for one building row with height (H/M1/3=0.5)

D=2H ; Max{0
D=2H ; Max{0
1.50 ______ Free air

1.25 1.00 0.75 0.50 0.25 0.00 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

Figure 5.6(b) Sheltering factors for one building row with height (H/M1/3=1.0)

50

The above observations also follow for two building-rows as shown in Figure 5.7(a-b) except that the sheltering factors are significantly larger. The observations concerning the blast wave energy diffracted above the building are also equally valid. The larger sized building is responsible for diffracting more energy. Thus sheltering factors for the region H < Y < 2H are smaller than for 0 < Y < H. Note also that the sheltering factors in this region are less at the outer zone boundary There are immediate conclusions from these observations. First of all global sheltering factors for the region 0 < Y < H can be assigned to both zone boundaries and that these can be shown to be a simple functions of the normalized building height and the particular building geometry (in this case either 1 building-row or 2 building-rows with a gap of one building height). Such representation is shown in Figure 5.8 highlighting areas of uncertainty remaining from this analysis.

Normalized peak positive overpressure (Ps /Ps *)

Second the benefits of sheltering are far less if the proposed construction is substantially higher than the average building height on its upwind side. This observation stands to sense because the upwind obstacles tend to deflect blast energy upwards and, therefore, toward the upper level of the proposed building. D=2H ; Max{0
D=2H ; Max{0
1.50 ______ Free air

1.25 1.00 0.75 0.50 0.25 0.00 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

Figure 5.7(a) Sheltering factors for two building rows with height (H/M1/3=0.5)

51

Normalized peak positive overpressure (Ps /Ps *)

D=2H ; Max{0
D=2H ; Max{0
1.50 ______ Free air

1.25 1.00 0.75 0.50 0.25 0.00 8

10

12

14

16

18

20

1/3

Normalized distance (Xs /M )

Normalized peak positive overpressure (Ps /Ps *)

Figure 5.7(b) Sheltering factors for two building rows with height (H/M 1/3=1.0)

1 Building row

2 Building rows

1.00

0.75

?

?

0.50

?

0.25 MAX{0
0.5

1

1.5 1/3

Normalized building height (H/M )

Figure 5.8 Sheltering factors as a function of building height showing areas of uncertainty .

52

5.3.

Applicability of revised methodology

In order to provide an assessment of the applicability of the new methodology two scenarios have been considered in which buildings are assumed to have a nominal height of 10 m. Two scenarios have been considered. The larger inventory scenario (A) concerns one on the scale of the Flixborough disaster, namely a TNT equivalent (in ground burst TNT mass equivalent) of 16 Te (16,000 kg). Comparison with a smaller scenario (B) of 2 Te has also been considered. In the following calculation we shall determine the shift in meters of the Outer Zone boundaries for the two cases of one and two building rows. The calculation values are provided in Table 5.2. Table 5.2 Calculated values for scenarios used in methodology assessment A B Scenario 16 2 TNT mass (Te) 1 2 1 2 no. building rows 0.84 0.74 0.75 0.61 Sheltering factor 70 140 70 140 70 140 70 140 Zone boundary (mbar) 49 31 83 48 39 23 62 35 Boundary shift (m) 2 46 194 69 46 17 68 24 Land free-up (1000 m ) 120 First we determine the TNT equivalent mass-normalised height (H/M1/3) of the building row which yields H/M1/3 = 0.4 and 0.8 for scenarios A and B respectively. Referring to Figure 5.8 we can now calculate the corresponding sheltering factors for the two scenarios for the 1 and 2 building-row cases. These are found to range between 0.74-0.84 for scenario A and 0.61-0.75 for scenario B. These four sheltering factors can then used in conjunction with Figure 5.1 to determine the TNT mass normalized shifts for the 70 mbar and 140 mbar boundaries of the Outer Zone. The absolute boundary shift distances are given for the cases of 1 and 2 buildingrows in Table 5.2. In Case A, which has a comparable scale to the Flixborough disaster, the absolute shifts for the Outer Zone boundaries in the case of 2 building-rows are 48 m and 83 m respectively. For the case of 1 building-row these reduce to 31 m and 49 m . The corresponding zone boundary shift for the smaller inventory scenario (Case B) are 35 m and 62 m for the case of 2 building-rows and 23 m and 39 m for the case of one building-row. Clearly the zone shift distances for the larger inventory scenario represent substantial free-up of land. When the whole perimeter of the inner and outer boundaries of the Outer Zone are taken into account these can be expressed as areas lying between the concentric circles given by the existing and revised zone boundaries (see Table 5.2). These calculations of land free-up are approximate conservative values based upon the smallest sheltering values for a site centred within four building heights of the last building-row. Clearly substantial improvement might be expected if there were more rows of buildings between the proposed site and explosion source.

53

6.

CONCLUSIONS

A review of TNT and TNO Multi-Energy explosion hazard methodologies has established that there is an opportunity to reduce the conservatism of current methods used by MSDU to advise planning authorities. This is because the models apply to the ideal blast wave regime in which the leading shock is most affected by interaction with structures. Comparison between CFD calculations and measurements of overpressure in field-scale experiments has established that sufficient accuracy can be established in the prediction of blast wave parameters. The ability to predict the detailed pressure-time histories in three-dimensions has also been established when suitably powerful (workstation) computers can be dedicated to the study. Simplified geometries however allow 2-dimensional computational representations which have the benefit that individual calculations can be performed within a day using present day workstations. Thus, without prohibitive investment of staff-time or computing resources many calculations can be performed to explore the effects of blast interaction with simplified building geometries in the ranges of existing consultation distances. The scope for field-scale experimental studies at HSL has been found limited by the size of the charge weight necessary. Larger charges generate ground shocks which cause substantial low frequency corruption to the far field overpressure signals and can also propagate off-site to cause interference to the general public. The high quality of agreement between experiment and CFD for the uncorrupted measurements however has established that, provided due consideration is taken to the choice of grid cell sizes, then adequate accuracy can be determined from the CFD predictions. The investigation of more general geometries however, will require the application of three-dimensional methods with a corresponding need for suitably high-performance workstations and proportionally more staff time in servicing the calculations. The parameter variation exercises offered by the computational study have enabled a simple methodology to be proposed which takes account of sheltering caused by buildings on the upwind side of the proposed site. This methodology employs TNT mass scaling as is consistent with existing TNT and TNO Multi-Energy methods. By calculating the mass normalized building height for the given scenario a simple reference to graphs can establish the worst case sheltering factors for the building layouts of relevance. In this study the extended methodology is limited to one and two continuous rows of buildings oriented at right angles to the line between the explosion source and the proposed site. Given the sheltering factors, and the building geometry, the TNT inventory masses are then used to determine the shift in the inner (140 mbar) and outer (70 mbar) boundaries of the Outer Hazard Zone. For a representative building height of 10 m two typical scenarios were selected, the largest of which is on a scale with the Flixborough disaster. The sheltering has the effect of moving the inner and outer boundaries of the Outer Zone up to 50 m and 80 m respectively. These distances represent substantial free-up of land. Buildings that provide a sheltering effect may collapse as a result of the blast overpressure. However, the time-scale for the building collapse is large compared to the time-scale of the blast wave and therefore the building collapse will not decrease the sheltering effect of the building. It is important to note that intervening buildings that provide a sheltering effect may be demolished within the lifetime of the proposed new development. This in turn will change the

54

consultation distance and will clearly have an effect on MSDU’s assessment of the LUP case. Therefore MSDU need to make a policy decision on whether the sheltering effect of buildings should impact on their methodology used for assessing LUP cases. However, in COMAH safety reports, where a ‘snapshot’ of the blast hazard is required, such considerations do not need to be taken into account. Therefore these sheltering factors could be used and potentially lead to significant financial savings for major hazard site owners. 6.1.

Recommendations

The new blast sheltering methodology proposed in this report needs further development before it can be used for realistic scenarios. The scope of the CFD simulation study should be widened to cover more building sites and layouts primarily to establish whether the extended findings support the conclusions of this study, namely that a simple revision can be made to existing methodologies to account for sheltering. A further limited field-scale experimental programme should be considered for more general three-dimensional building layouts, provided the difficulties caused by ground-shock are managed to acceptable levels, to establish whether the CFD methods can yield satisfactory agreement with experiment without the need for excessive computational or staff resources. Suggestions for study concern gaps between buildings, different heights of buildings on the upwind side, orientation of the building rows to the direction of the incident shock and buildings on the downwind side of the proposed site. Following the successful conclusion of the above studies an extensive three-dimensional simulation study should be considered to provide the necessary blast data upon which to base a revised methodology.

55

7.

REFERENCES [1].

W Baker, P Cox, P Westine, J Kulesz and R Strehlow Explosions Hazards and Evaluation, volume 5 of Fundamental Studies in Engineering, Elsevier, 1983

[2].

R J Harris, M J Wickens Understanding of Vapour Cloud Explosions – An Experimental Study, Institute of Gas Engineers, 55th Autumn Meeting, Communication 1408, 1989; Advisory Committee on Major Hazards, 2nd report, Health and Safety Executive, 1979

[3].

J B M M Eggen, GAME: development of guidance for the application of the multi-energy method, TNO Prins Maurits Laboratory, 1995

[4].

N W Hurst, R D Fitzpatrick, G Clay Development of RISKAT for LPG Part I. Whole Tank Failures, Calculation of Overpressure and Radiation, IR/L/HA/89/1, HSE

[5].

S A Formby and R K Wharton. ‘Blast characteristics and TNT equivalence values for some commercial explosives detonated at ground level’, J Haz. Mats., Volume 50, p183 (1996).

[6].

S A Formby. ‘Quantification of the air blast characteristics of commercial explosives’, I Chem E. 2nd Eur. Conf. Major Hazards Onshore and Offshore’, UMIST. p147, (1995).

[7].

S A Formby. ‘Assessment of the explosion hazard of low order detonating explosives - final report’. HSL Report IR/EX/96/01

[8].

DM Ingram. Numerical Prediction of Blast Wave Flows around Rigid Structures. PhD thesis, Department of Mathematics and Physics, Manchester Metropolitan University, November 1992.

[9].

B van Leer. On the relation between the upwind-differencing schemes of Godunov, Engquist-Osher and Roe. SIAM Journal on Scientific and Statistical Computing, 5(1):1-20, 1984.

[10].

A Harten, PD Lax, and B van Leer. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Review, 25(1):35-61, 1983.

[11].

EF Toro, M Spruce, and W Speares. Restoration of the contact surface in the hll Riemann solver. Shock Waves, Vol 4, pages 25-34, 1994.

[12].

B Jiang, DM Ingram, DM Causon, and R Saunders. A global simulation method for obtaining reduced reaction mechanisms for use in reactive blast wave flows. Shock Waves, 5(1/2):81-87, 1995.

[13].

DM Ingram, B Jiang, and DM Causon. On the role of turbulence in detonation induced by mach stem reflection. Shock Waves, 8(6):327-336, 1998.

56

[14].

P Batten, DM Ingram, R Saunders, and DM Causon. An implicit viscous solver for the compressible navier-stokes equations. Computers and Fluids, 25(4):421-341, 1996.

[15].

DM Ingram, N.S. Ellis, and DM Causon. Hydrodynamic code calculations of an airburst. In K. Takayama, editor, Shock Waves: Procceedings of the 18th International Symposium on Shock Waves, Sendai, Japan, 1991, volume II, pages 905-910, Berlin, Germany, 1992. Springer-Verlag.

[16].

DM Ingram, P Batten, C Lambert, and DM Causon. Hydrodynamic code calculations of a blast on a tank farm. In R. Brun and L.Z. Dumitrescu, editors, Shock Waves @ Marseille: Proceedings of the 19th International Symposium on Shock Waves, Marseille, France 1993, volume IV, pages 419-424, Berlin, Germany, 1995. Springer-Verlag.

[17].

DM Causon and DM Ingram. Numerical simulation of unsteady flow in a twin side-by-side intake system. Aeronautical Journal, 101(1008):365-370, 1997.

[18].

DM Causon, R Saunders, DM Ingram, and P Batten. Numerical modelling of blast wave in the process industry. Final report - SERC grant GR/H18654, Centre for Mathematical Modelling and Flow Analysis, Department of Mathematics and Physics, Manchester Metropolitan University, Manchester M1 5GD, December 1993.

[19].

W Baker, P Cox, P Westine, J Kulesz, and R Strehlow. Explosions Hazards and Evaluation, volume 5 of Fundamental Studies in Engineering. Elsevier, 1983.

[20].

DV Ritzel and K Matthews. An adjustable explosion-source model for CFD blast calculations. In AFP Houwing and A Paull, editors, Proceedings of the 21st International Symposium on Shock Waves, page Paper 6590. Panther Publishing, Fyshwick Australia, 1998.

57