The deflagration pressure analysis code (DPAC) has been upgraded for use in modeling hydrogen deflagration transients. The upgraded code is benchmarked using data from vented hydrogen deflagration tests conducted at the HYDRO-SC Test Facility at the University of Pisa. DPAC originally was written to calculate peak pressures for deflagrations in radioactive waste storage tanks and process facilities at the Savannah River Site. Upgrades include the addition of a laminar flame speed correlation for hydrogen deflagrations and a mechanistic model for turbulent flame propagation, incorporation of inertial effects during venting, and inclusion of the effect of water vapor condensation on vessel walls. In addition, DPAC has been coupled with chemical equilibrium with applications (CEA), a NASA combustion chemistry code. The deflagration tests are modeled as end-to-end deflagrations. The improved DPAC code successfully predicts both the peak pressures during the deflagration tests and the times at which the pressure peaks.

## Introduction

The deflagration pressure analysis code (DPAC) code originally was written to evaluate peak pressures for deflagrations in radioactive waste storage and process facilities at the Savannah River Site [1]. The primary purpose for the code was to provide a safety analysis method to address risks from hypothetical deflagrations. The original version of the code coupled a thermodynamic model of adiabatic isochoric complete combustion (AICC) for hydrogen/air and/or benzene/air mixtures with a simple flame propagation calculation. An option was added for flame propagation in one, two, or three dimensions. In lieu of a detailed thermal hydraulic model, the code simply divided the deflagration vessel into two volumes on either side of the flame front, one for unburned gas and the other for burned gas. The code included calculations for venting, radiative heat losses, and structural deformation due to pressurization from the deflagration. DPAC was written in standard Fortran and can be compiled to run on a personal computer.

Recently, there has been renewed interest in modeling hydrogen deflagrations. This interest, combined with the availability of additional data from vented hydrogen deflagration tests, provided motivation to modify DPAC.

The modified version of DPAC is limited to linear, end-to-end propagation of hydrogen/air deflagrations. This combination covers the majority of hypothetical deflagrations encountered in safety analyses. Upgrades to DPAC include the addition of a concentration-dependent laminar flame speed correlation, a model for turbulent enhancement of the flame speed, and the incorporation of inertial effects during the flame propagation. In the modified code, the thermodynamic calculations are performed by the NASA combustion code chemical equilibrium with applications (CEA) [2,3], which is invoked within the Fortran executable. Structural analyses have been deleted, and a model for water vapor condensation on the vessel walls has been added.

Previous work on DPAC and the relationship of the code to other deflagration modeling work has been reviewed previously [1], so the following discussion focuses on the benchmark tests, the changes to DPAC, and a comparison of DPAC calculations with the benchmark test results.

## Description of Benchmark Tests

A series of tests conducted using the HYDRO-SC Test Facility at the University of Pisa was selected for benchmarking of the deflagration calculations [4,5]. Carcassi and Fineschi conducted a series of unvented and vented tests in this facility, using 14% hydrogen in air. The HYDRO-SC test vessel was an upright cylindrical steel tank with a diameter of 0.650 m and a total height of 1.628 m. The test vessel had top and bottom end caps, each with a radius of 0.520 m. Tests were conducted with an unvented tank and with 0.030-m, 0.05-m, 0.07-m, and 0.10-m diameter vents. The benchmarking calculations do not include any losses for these vents other than a 0.5 velocity head loss for a sudden flow entry. During the tests, the gas mixture was ignited at the bottom of the tank, and the vent was located at the top. Consequently, these tests can be modeled using a model based on end-to-end flame propagation.

The Carcassi and Fineschi tests were chosen for benchmarking, because they: (1) used hydrogen gas in air, (2) were conducted using a tank that, although somewhat smaller than typical processing vessels or storage tanks, was representative in size, and (3) produced data for a range of vent sizes. In addition, an examination of the data indicates that the pressure transducers that were used had a sufficiently rapid response time to generate an accurate pressure profile. Figure 1 depicts the pressure transients for the Carcassi and Fineschi tests.

## Overview of Changes to Deflagration Pressure Analysis Code

To obtain a better model to benchmark the pressure transients from the Carcassi and Fineschi tests, the DPAC source code was extensively revised. The most basic change was to the structure of the transient calculations within DPAC. In the preceding version, an energy balance was applied to simultaneously calculate the burning velocity and the rates of pressure and gas temperature changes due to burning, venting, and radiation heat transfer. In the modified version, these changes are calculated separately and consecutively at each calculation time step. This approach is valid if the time step is sufficiently small.

The flame speed is estimated by calculating the extent to which a burning flame compresses the unburned gas in the vessel. It is assumed that AICC takes place at the burning surface, and that the freshly burned gas then expands and simultaneously compresses both the unburned and previously burned gas in the vessel. Compression of the unburned and previously burned gases is assumed to occur adiabatically and isentropically, until the pressure of the freshly burned gas equals the vessel pressure.

A summary of other code modifications follows.

## Changes to Combustion Calculations

The combustion calculation in DPAC is replaced by a call to an executable version of the NASA combustion code CEA [2,3]. CEA is used to compute gas properties both before and after combustion. The combustion is modeled as a constant volume combustion, so CEA gives the burned gas pressure and temperature prior to expansion into the volume occupied by the unburned gas. This expansion provides an acceleration factor to convert the turbulent burning velocity into an actual flame speed.

## Changes to Laminar Burning Velocity

The previous version of DPAC assigned the laminar burning velocity and the turbulent velocity multiplier as input items. In the modified version, an empirical laminar burning velocity correlation is added so that the laminar burning velocity can be calculated as a function of the temperature, the pressure, and the unburned gas composition. The correlation is fit to data compiled by Dahoe at varying hydrogen/air ratios for the unburned gas [6], as shown by Fig. 2.

The empirical curve fit constants $c1$, $c2$, $c3$, $c4$, and $c5$ have the values 4.216679, −0.2023, 2.217257, 1.425145, and −4.25637, respectively.

## Changes to Turbulent Velocity Multiplier

The modified version of DPAC also includes a turbulence model to estimate the ratio of the turbulent burning velocity to the laminar burning velocity. The burning velocity prior to the expansion of the burned gas is estimated by combining an analysis conducted by Chen et al. [8], which gives the turbulent burning velocity as a function of the turbulence strength for the flame, with an estimate of the turbulence strength at the flame front. The Chen et al. correlation is an asymptotic combination of a turbulent diffusion model, originally proposed by Damköhler [9], and a limiting value for the turbulent flame speed obtained from experiments. The Damköhler model states that, for low level of turbulence, the flame front is thin. For thin flames, the turbulent flame speed is modeled as a function of the turbulent diffusivity, which in turn is proportional to the turbulence strength. Experiments with highly turbulent flames indicate that for thick flames, the increase in the turbulent burning velocity is due to the increase in the flame surface area caused by “wrinkling” of the flame front. This increase is limited to a value of twice the turbulence strength.

## Changes to Flame Speed Calculation

The acceleration of the flame by expansion of the combustion gases is modeled by calculating the volume of unburned gas displaced as the burning gases expand from their AICC pressure to the pressure inside the vessel. The degree of acceleration depends on whether the flame behaves as a “confined” or a freely expanding flame. The rate of pressure increase for the Carcassi and Fineschi benchmark tests initially is high, indicating that the flame is “confined” and later drops to a level consistent with that of a freely expanding flame. Separate flame acceleration factors are derived for “confined” and free flames.

*D*, this equation is

where *P*_{tr} is the pressure at the transition from a confined to an expanding flame. A comparison of Eqs. (8) and (20) shows that the rate of pressure increase for an expanding flame is half that for a confined flame.

The transition from a confined flame to an expanding flame is applied when the resistance to expansion into the volume of unburned gas in the vessel becomes equal to the resistance to expansion into the burned gas. This is assumed to occur when the numbers of moles of unburned and burned gas remaining in the vessel are equal.

The flame speed calculation sets the flame speed (with respect to a fixed location) equal to the product of the burning velocity (the velocity at which the unburned gas flows into the flame) and a volume ratio for expansion of the gas that has burned during the current time step. It is assumed that the gas expands until its pressure equilibrates with the pressure in the vessel. The starting pressure prior to expansion is set equal to the AICC pressure. To ensure that the final vessel pressure for complete combustion will reach the AICC pressure (not accounting for heat losses), the volumetric expansion calculations are based on simple, constant temperature volumetric displacement rather than isentropic expansion and compression. It is recognized that this simplistic approach may lead to errors in temperature calculations, but the approach is deemed adequate for calculating pressure transients.

## Changes to Venting Analysis

An exact, iterative calculation of the venting flow rate replaces the semi-empirical estimate of venting flow rates for nonchoking flow. As in the previous version of DPAC, the flow out the vents is assumed to be at steady-state with respect to the relatively slower combustion transient. The vent flow is modeled using equations for isentropic flow [11], which are solved iteratively. The relative venting rates of unburned and burned gases are proportioned by the relative volumes of unburned and burned gases remaining in the vessel; this apportioning was also used in the previous versions of DPAC.

where *v*_{out} is the velocity out the vent, at the inlet temperature and pressure.

It may be noted that the mass balance is expressed in terms of the substantial derivative. This is done to include inertial nonequilibrium terms, so that the analysis is applied to a shock wave that propagates back into the vessel from the vent.

where $((1/P)(dP/dt))out$ represents the normalized rate of pressure change due to flow out the vent.

For vent areas less than the area calculated by this expression, a transient shock will not form at the vent, and venting will start immediately after ignition. For larger areas, a shock will develop, and venting will be delayed until the pressure increases to the value specified by Eq. (21).

## Changes to Heat Transfer Analysis

The emissivity for the vessel surface is assumed to be that for an oxidized steel surface at ambient temperature, which is approximately 0.8.

where $Peff$ is the effective pressure representing water vapor radiation absorption band broadening by air, and $L\epsilon $ is the characteristic length for thermal radiation in vessel. The characteristic radiation length is assumed to be equal to 3.5 times the gas space volume, divided by the interior surface area [14].

## Comparison of Deflagration Pressure Analysis Code Calculations With Benchmark Deflagration Transients

Figures 3–5 compare deflagration transients calculated by the modified version of DPAC with the measurements of Carcassi and Fineschi. Transients are compared for the no venting case (Fig. 3), the smallest vent (Fig. 4), and the largest vent (Fig. 5). The predicted deflagration transients shown in these figures fairly closely match the measured transients, with respect to both the peak pressure and the rate of pressure change during the deflagration and subsequent venting. Although DPAC provides an accurate estimate of the overall rate for the deflagration transient, the initial rates of pressure increase calculated by DPAC exceed the measured rates. This probably is due to the fact that the DPAC model assumes that the flame front starts as a planar front that covers the entire bottom surface of the vessel, whereas an actual flame probably ignites over a small surface area. This implies that the initial rate of pressure should be lower than calculated, until the flame spreads to cover the entire cross section of the vessel.

The comparison for the unvented test provides a check of the accuracy of the thermal radiation and condensation heat transfer models within DPAC. As Fig. 3 shows, without any heat losses the pressure in the vessel increases until it reaches the AICC pressure of 5.448 bar. With heat losses included, the predicted peak pressure is 4.702 bar, compared to a measured peak pressure of 4.507 bar. This comparison indicates that the DPAC predictions are slightly conservative in that they underestimate the amount of heat loss.

For the vented tests, the comparisons also indicate how well DPAC models the effect of venting. DPAC accurately predicts the pressure transients for the two smaller vents, as illustrated by Fig. 4. The predicted pressure transients for the two larger vents did not agree with the measured transients until inertial effects are included, as Fig. 5 demonstrates. As stated previously, one inertial effect is the development of a pressure shock that propagates backward from the vent into the vessel for sufficiently large vents, and a second effect is the acceleration of the flame front by venting, again for large vents. The DPAC code models the first effect by suppressing venting until the vessel pressure reaches the minimum value for sonic flow through the vents. DPAC accounts for the second effect by allowing complete venting of the unburned gases when the difference between the vent flow rates, scaled to the flame front area, first exceeds the flame speed. With these changes, DPAC generates a better estimate of the peak pressure and a better estimate of the rate of venting pressure losses for the 100 mm vent.

Table 1 summarizes the comparison of measured and calculated peak pressures and venting rates for the Carcassi and Fineschi tests. It may be seen that the modified version of DPAC accurately predicts the peak pressures and the deflagration times required to reach the peaks, when inertial effects are included to model the two tests with the largest vents.

Vent diameter (m) | Measured peak press. (0.1 MPa, bar) | Measured time to peak (s) | Calculated peak press.^{a} (0.1 MPa, bar) | Calculated time to peak (s) |
---|---|---|---|---|

No vent | 4.507 | 0.278 | 4.702 | 0.300 |

0.03 | 3.882 | 0.271 | 4.053 | 0.307 |

0.05 | 3.334 | 0.264 | 3.081 | 0.320 |

0.07 | 2.670 | 0.145 | 2.761 | 0.131 |

0.10 | 2.187 | 0.119 | 2.382 | 0.138 |

Vent diameter (m) | Measured peak press. (0.1 MPa, bar) | Measured time to peak (s) | Calculated peak press.^{a} (0.1 MPa, bar) | Calculated time to peak (s) |
---|---|---|---|---|

No vent | 4.507 | 0.278 | 4.702 | 0.300 |

0.03 | 3.882 | 0.271 | 4.053 | 0.307 |

0.05 | 3.334 | 0.264 | 3.081 | 0.320 |

0.07 | 2.670 | 0.145 | 2.761 | 0.131 |

0.10 | 2.187 | 0.119 | 2.382 | 0.138 |

The calculated peak pressures for the 0.070-m and 0.10-m vent tests account for inertial effects associated with propagation of a pressure shock back from the vent during the initial stages of the transient and with subsequent acceleration of the flame front due to venting.

## Summary and Conclusions

The DPAC transient deflagration code has been modified and successfully benchmarked using a series of vented deflagration tests performed at the University of Pisa. The modified version of DPAC serves as a tool to predict maximum pressures for safety analyses of hypothetical hydrogen deflagrations at the Savannah River Site. Modification to DPAC include the addition of a laminar burning velocity correlation and a mechanistic model for the turbulent burning velocity multiplier, improvement of the modeling of radiative heat transfer, incorporation of a model for condensate heat transfer, and the inclusion of inertial effects in the calculation of the flame speed. The benchmarking studies demonstrated that the inertial effects play a significant role in determining the rate of pressure rise and are needed to successfully model well-vented deflagrations.

Several features present in the previous version of DPAC were judged to be of minor use in addressing current safety analysis concerns and therefore were deleted. The modified version of DPAC no longer performs combustion chemistry calculations for gas mixtures other than hydrogen and air and does not include any structural analyses or wall temperature calculations. In addition, the code is restricted to the modeling of one-dimensional, end-to-end flame propagation. Even with these limitations, it is anticipated that DPAC will be able to address most safety analysis needs.

Deflagration pressure analysis code is configured as a Fortran executable on a personal computer platform. A code user is required to use an input text file. Use of the code requires the presence of an executable copy of Version 2 of the NASA combustion code CEA on the same platform or at an addressable location.

## Acknowledgment

This manuscript has been authored by Savannah River Nuclear Solutions, LLC under contract no. DE-AC09-08SR22470 with the U.S. Department of Energy.

## Funding Data

U.S. Department of Energy (Contract No. DE-AC09-08SR22470).

## Nomenclature

- $A$ =
vessel interior surface area, m

^{2}*a*_{4}=turbulent diffusivity coefficient, 0.78

- $AF$ =
flame front area, m

^{2}- $Av$ =
cross-sectional flow area of the vent, m

^{2}- $B$ =
arbitrary proportionality constant

*b*_{1}=limiting ratio of turbulent burning velocity to turbulence intensity, 2.0

*b*_{2}=proportionality constant, 1.78

*b*_{3}=coefficient for the relative increase in the ratio of the turbulent and laminar burning velocities as a function of the square root of the ratio of the turbulent and laminar diffusivities, 1.0

- $ci$ =
empirical curve fit parameters for Dahoe laminar flame speed measurements (

*i*= 1–5)- $dt$ =
vessel diameter, m

- Da
=_{T} Damköhler number for turbulent flames

- $FA$ =
fraction of the vessel surface area that receives radiation, equivalent to a viewing factor

- $hg$ =
height of the total gas space in the vessel, m

- $kP$ =
heat capacity ratio for burned gas, J/mol/K

- $kR$ =
heat capacity ratio for unburned gas, J/mol/K

- $L$ =
length scale for turbulent diffusion across the flame, m

- $LF$ =
flame thickness, m

- $L\epsilon $ =
characteristic length for thermal radiation in the vessel, m

- $mP$ =
mass of burned gas in the vessel, kg

- $mR$ =
mass of unburned gas in the vessel, kg

- $Mg$ =
molecular mass of gas, kg/mol

- $Ms$ =
Mach number at stagnation conditions

- $P$ =
pressure in the vessel, Pa

- $Pa$ =
atmospheric pressure, 0.1013 MPa

- $PAICC$ =
pressure for adiabatic, isochoric complete combustion (AICC), Pa

*P*_{eff}=effective pressure representing water vapor radiation absorption band broadening by air, Pa

- $PH2O$ =
water vapor partial pressure, Pa

- $Pi$ =
initial pressure in the vessel, Pa

- $Ps$ =
standard pressure for the laminar burning velocity correlation, 0.1 MPa (1 bar)

*P*_{tr}=pressure at the transition from a confined to an expanding flame, Pa

*P*_{0}=stagnation pressure, Pa

*P*_{2}=exit pressure, Pa

- $qc$ =
condensation heat transfer rate, W/m

^{2}*q*_{rad}=thermal radiation heat transfer rate, W/m

^{2}- $qrad,c$ =
total heat transfer rate including thermal radiation and condensation, W/m

^{2}- $Rg$ =
ideal gas law constant, 8.314 J/mol/K

- $sL$ =
laminar burning velocity in Lagrangian coordinates (following the flame front), m/s

- $sL,s$ =
laminar burning velocity at standard temperature (300 K) and pressure (0.1 MPa or 1.0 bar), m/s

- $sT$ =
turbulent burning velocity in Lagrangian coordinates (following the flame front), m/s

- $t$ =
time, s

- $TP$ =
product gas temperature for constant pressure complete combustion, K

- $TR,i$ =
initial unburned gas temperature, K

- $Ts$ =
standard temperature for the laminar burning velocity correlation, 300 K

*T*_{0}=temperature of unburned gas, stagnation temperature, K

- $V$ =
vessel volume, m

^{3}- $v\u2032$ =
turbulence strength, m/s

- $vF$ =
flame speed, m/s

*v*_{out}=velocity out the vent, at the inlet temperature and pressure, m/s

- $vT$ =
total velocity at the flame front, m/s

- $VF$ =
fraction of the vessel gas space occupied by burned gas

- $VP$ =
burned gas volume, m

^{3}*V*_{P}_{,}_{tr}=burned gas volume at the transition from a confined to an expanding flame, m

^{3}- $xH2O$ =
mole fraction of water vapor in the burned gas

- $\Delta hvap$ =
molar heat of vaporization for water at the surface temperature, J/mol

- $\Delta P$ =
change in pressure without including the effect of condensation, Pa

- $\Delta Pc,tot$ =
total change in pressure for venting of burned gas, including the effect of condensation on the vessel walls, Pa

- $\Delta t$ =
time step, s

- $\epsilon eff$ =
effective emissivity for the burned gas

- $\epsilon gas$ =
emissivity of the burned gas

- $\epsilon surf$ =
emissivity of the vessel surface

- $\rho P$ =
burned gas density, kg/m

^{3}- $\rho R$ =
unburned gas density, kg/m

^{3}- $\sigma $ =
Stefan–Boltzmann constant, 5.671 × 10

^{−8}W/m^{2}/K^{4}- $\phi $ =
hydrogen equivalency ratio, equal to the hydrogen molar concentration divided by twice the oxygen molar concentration