Ind. Eng. Chem. Res. 2001, 40, 5847-5857
5847
Emissions from Fixed-Roof Storage Tanks: Modeling and Experiments Renato Rota,*,† Simone Frattini,† Sabrina Astori,‡ and Renato Paludetto§ Dipartimento di Chimica Fisica Applicata, Politecnico di Milano, Via Mancinelli 7, 20131 Milano, Italy, Centro Ricerche Enichem, Via della Chimica, 30175 Marghera, Italy, and Divisione Poliuretani e Cloro Enichem, Piazza Boldrini 1, 20097 S. Donato Milanese, Italy
The two significant types of emissions from fixed-roof storage tanks are “storage” and “working” losses. While the former arise from changes in environmental parameters, the latter result from liquid level changes in the tank. Optimizing the storage policy requires the prediction of the influence of the main operating parameters on the overall emissions. The aim of this work was to develop a mathematical model able to represent the main phenomena involved in such emissions and to validate it by comparison with both laboratory-scale and real-size experimental data. Moreover, because semiempirical relations are available to estimate monthly or annual emissions, a comparison with these relations has been carried out both to further validate the model and to discuss the range of reliability of such relations. Introduction Emissions from fixed-roof storage tanks are responsible not only for a depletion of the product supply but also for contributions to atmospheric pollution. Moreover, like any other atmospheric emission, they can also be regulated by state or regional agencies. For several decades, semiempirical relations have been proposed to estimate the overall amount of emissions from fixedroof storage tanks1-3 in terms of monthly or annual amounts. However, while the overall emissions from an industrial site can give a broad idea of the extension of the problem, predicting the influence of the main operating parameters on the emission rate allows for optimization of the storage policy, thus reducing the emissions. Previous contributions concerning the mathematical simulation of tank emissions have been presented in the literature.4-8 However, they usually consider separately the simulation of the so-called “breathing” (or “standing”) losses from that of the “working” ones. While the former are related to the expansion of the vapor phase due to changes in environmental parameters (e.g., solar radiation, barometric pressure, or day/night temperature variations) and occur without any liquid level change in the tank, the latter arise during the filling phase when the liquid level rises, the pressure inside the tank exceeds the relief pressure, and vapors exit the tank through the top vent. Separating the simulation of these phenomena is an approximation because they obviously proceed simultaneously. Moreover, the simultaneous simulation of both phenomena is important because the initial conditions inside the tank vapor phase (that is, the temperature and concentration profiles) play a role for estimating the emissions in the following working cycle (that is, filling, emptying, and standing). In other words, the prediction of the overall emissions in a given period, say 1 year, requires a series * To whom correspondence should be addressed. Tel: +39 0223993154. Fax: +39 0223993180. E-mail: renato.rota@ polimi.it. † Politecnico di Milano. ‡ Centro Ricerche Enichem. Fax: +39 0412912530. E-mail:
[email protected]. § Divisione Poliuretani e Cloro Enichem. Fax: +39 0252032884. E-mail:
[email protected].
of computations where the final conditions of a working cycle become the initial conditions of the following one. Moreover, the available mathematical models are usually validated by comparison with the overall emission data from fixed-roof storage tanks. However, validating complex models by comparison with integral data (such as monthly or annual emissions) can be misleading because of the cancellation of errors. As a consequence, the main aim of this work has been to develop a mathematical model representing the main phenomena involved in emissions from fixed-roof storage tanks (that is, those referring to the breathing and working losses). This allows one to discuss the influence of the main operating parameters on the emissions and to optimize the storage policy. In particular, because the aim of this work has been to develop a daily working tool for evaluating the storage policy with respect to environmental problems, mathematical models developed in the frame of computational fluid dynamics have not been considered because they still require dedicated resources and cannot be considered daily working tools to be used by plant people. The model predictions have been validated by comparison not only with the overall emission data but also with concentration values measured in both laboratoryscale and real-size fixed-roof storage tanks. Then, model results have also been compared with the predictions of well-established semiempirical relations, thus allowing one to discuss the range of reliability of such relations. Finally, it should be mentioned that the developed model can also be useful for setting up the design basis of the tank emission treatment unit before discharge to the atmosphere. Mathematical Model The vapor inside a fixed-roof tank is stratified, thus leading to a clear vertical gradient of both composition and temperature.5 Such a gradient depends on both operating and environmental parameters and strongly influences the emissions through the upper layer near the tank top. Consequently, it is not possible to utilize a simple lumped parameter model to predict correctly the influence of various operating parameters. However, considering the tank as a vertical cylinder, it is reasonable to assume no radial (that is, horizontal) gradients,
10.1021/ie010111m CCC: $20.00 © 2001 American Chemical Society Published on Web 11/02/2001
5848
Ind. Eng. Chem. Res., Vol. 40, No. 24, 2001
with the boundary condition ν|z)νjt ) ν*. Here ν* accounts for both the vapor-liquid interface movement and the mass transport between liquid and vapor phases due to evaporation or condensation: ν* ) νj - (D/F)(∂C/ ∂z)|z)νjt. However, when νj * 0, it usually prevails and νj = ν*, while the transport term can play a role during the breathing phase, particularly for highly volatile compounds. The mass transport equation is
∂C ∂2C ∂(Cν) -E 2 + )0 ∂t ∂z ∂z
Figure 1. System coordinate and integration domain (dashed region) as a function of time.
Also in this case, E represents the axial dispersion coefficient, which equals the diffusion coefficient in motionless conditions, and it has been computed through standard correlations for empty tubes9 during filling and emptying. The boundary condition at the liquid surface is
at least as far as convective motions inside the tank are negligible. The validity of this assumption will be verified by comparing model predictions with experimental data, as discussed in the next section. Accordingly, the main phenomena involved in both breathing and working losses can be described by the continuity equation and the mass and energy balance equations together with suitable boundary conditions to represent the tank volume filled by vapors. Note that the liquid phase has not been included in the integration domain, as shown in Figure 1, because its thermal inertia is much higher than that of the vapor phase. Moreover, the liquid temperature can also be controlled by heat exchangers when compounds requiring cooling or heating to maintain the desired characteristics are stored. Using these approximations, the energy balance equation is then represented by (at constant pressure)
∂(FcPT) ∂2T ∂(FνcPT) - ET 2 + )0 ∂t ∂z ∂z
C|z)νjt ) Csat
E
ET is the thermal dispersion coefficient, which equals the thermal conductivity when the vapor is at rest, while during vapor motion it can be computed using standard correlations for empty tubes.9 The boundary conditions for this equation are
(2)
{
|
z)H
∂C ∂z
Here, TILG is the vapor-liquid interface temperature and Tt is the tank internal roof temperature. Note that because during filling or emptying the liquid level changes, the lower boundary condition is imposed on a boundary moving at a velocity equal to νj, the vaporliquid interface velocity, which is an arbitrary function of time. In particular, it can be positive (filling phase), negative (emptying phase), or even zero (breathing phase). For the sake of simplicity, in the following, only step functions will be considered, as shown in Figure 2, which is reasonable when pumping is carried out at a constant rate. The continuity equation is given by
) Cν|z)H
inhale (ν|z)H < 0)
|
)0
z)H
exhale (ν|z)H > 0)
ν ) ν* +
( |
(
| )
∂T R ∂T TILG ∂z z ∂z
(7)
(8)
z)νjt
where R ) ET/FILGcP. When this relation is introduced in both the continuity and mass-transfer equations, the system reduces to the two following equations:
| )
∂T ∂T 2 ∂T ∂T ∂T R ∂2T + ν* T 2+ )0 ∂t ∂z TILG ∂z ∂z ∂z ∂z z)νjt ∂C ∂C ∂2C R ∂2T R ∂T ∂T - E 2 + ν* + +C )0 ∂t TILG ∂z ∂z z)νjt ∂z TILG ∂z2 ∂z
(
(6)
Using both the ideal-gas and the constant-pressure approximations, the accumulation term in the energy transport equation disappears because ∂(PMPcPT/RT)/ ∂t ) 0. Consequently, the same equation can be integrated to deduce an explicit relation for the vapor velocity:
(3)
(
∂C ∂z
stating that diffusive flux equals the convective one. During exhale conditions (ν|z)H > 0, typically during tank filling or heating due to environmental conditions), the mass flux of pollutant is assumed to be due solely to the convective flow, that is
T|z)νjt ) TILG
∂F ∂(Fν) + )0 ∂t ∂z
(5)
where Csat is the saturation concentration at the interface conditions. At the top dome there are two different boundary conditions according to whether the tank is in inhale or exhale conditions, that is, air enters the tank or vapors exit, respectively. During inhale conditions (that is, ν|z)H < 0, which are typical, for instance, of a tank emptying or cooling due to environmental temperature decreases), at the top dome the mass flux of pollutant is zero and there is no emission from the tank. In this case, the boundary condition has been set as
(1)
T|z)H ) Tt
(4)
( )
| ))
(9)
Ind. Eng. Chem. Res., Vol. 40, No. 24, 2001 5849
This reduced system enhances the robustness of the numerical method with respect to the original one. A peculiar characteristic of this problem is that it involves a moving boundary. In other words, the mesh point in contact with the liquid surface is moving with a velocity equal to νj, while the position of the mesh point at the tank top is fixed. This means that a classic integration technique, such as the Lagrangian or the Eulerian one, cannot be used. An arbitrary Lagrangian Eulerian (ALE) technique11-13 has been used to handle this problem. This method involves a finite difference mesh whose vertexes can move with the same velocity as the fluid (in this case the Lagrangian approach is recovered), can be held fixed (such as in the Eulerian method), or can move in any other prescribed way. In other words, the ALE approach consists of recasting the equations from the physical domain to a new one. This is done by moving any cell vertex, with a different velocity allowing one to keep the boundaries fixed in the new domain. In this case, the expression used for the mesh vertexes velocity is
z-H w ) νj νjt - H
(10)
so that the velocity of the cell in contact with the liquid surface (that is, at z ) νjt) is w ) νj, while that of the cell close to the tank top (where z ) H) is zero. Because all of the cells are moving at a velocity equal to w, a new dimensionless coordinate within the new domain can be defined as follows:
y)
νjt - z z - wt ) H νjt - H
(11)
The peculiar feature of this dimensionless coordinate is that at any time the cell at the liquid surface (characterized by a value of z ) νjt) has a new coordinate y ) 0, while the new coordinate of the cell at the tank top (characterized by a value of z ) H) is y ) 1. In other words, using the ALE technique allows one to remove the movement of the boundary because in the new domain the boundary conditions are enforced at y ) 0 (the liquid surface) and y ) 1 (the tank top). As a consequence, it is possible to solve the problem in the new domain with a classic Eulerian technique. Equation 9 can be recast in the new domain by the following coordinate transformation:14
νj(1 - y) ∂ ∂ ∂y ∂ ∂ ∂ ) + ) + ∂t ∂y ∂t ∂t νjt - H ∂y ∂t ∂ ∂ ∂ ∂y 1 ) )∂z ∂y ∂z νjt - H ∂y ∂2 1 ∂2 ) ∂z2 (νjt - H)2 ∂y2
(12)
Constant boundary velocity (νj) has been assumed to deduce the last equalities. This means that these equations are correct inside each region of the step function representing νj (see Figure 2). Introducing eq 12 into eq 9 leads to the following new system of partial differential equations:
{
(
Figure 2. Vapor-liquid interface velocity as a function of time.
with the boundary conditions (in the new domain)
T|y)1 ) Tt T|y)0 ) TILG
{
C|y)0 ) Csat
|
∂C )0 exhale ∂y y)1;ν(y)1)>0 ∂C E ) (Cν)|y)1(H - νjt) inhale ∂y y)1;ν(y)1)
