Fuites dans recipients et conductions: emission en phase liquide
Leakages in vessels and pipes: emission in liquid phase
Redactor:
Emilio Turmo Sierra
Ingeniero Industrial
CENTRO NACIONAL DE CONDICIONES DE TRABAJO
Tres peligros causantes de accidentes mayores en actividades industriales-incendio, explosión y fuga tóxica tienen normalmente su origen en la emisión de una sustancia peligrosa de un recipiente, equipo o conducción, la cual si es en fase gaseosa o en forma de vapor, va directamente a dispersarse en la atmósfera y si es en fase líquida y no es posible su recogida en un lugar seguro, estará sometida a la consiguiente vaporización y dispersión en el aire atmosférico.
Para los cálculos del caudal de emisión, se debe hacer la distinción entre gases, líquidos y productos que son gaseosos a presión atmosférica y temperatura ambiente, pero que se encuentran licuados a presión. En esta Nota Técnica, se considera el fenómeno físico de la emisión en fase líquida de un producto inflamable o tóxico por un orificio de un determinado tamaño, existente u ocasionado de forma accidental.
La existencia de recipientes de almacenamiento y de procesos presenta situaciones de riesgo por escape de una sustancia inflamable o tóxica. En el caso de una sustancia inflamable, al mezclarse la fase gaseosa con el aire dará lugar a la formación de una nube de vapor inflamable que al desplazarse puede encontrar fácilmente una fuente de ignición, ocasionando un incendio u explosión de una nube de vapor no confinada, lo que afectará a la zona del siniestro y posiblemente a zonas habitadas en un entorno próximo. En el caso de fuga de una sustancia tóxica, se formará una nube de vapor tóxico con una determinada concentración en función de la distancia a la fuente de emisión, que también afectará a la planta de proceso o almacenamiento y a su entorno, pudiendo generar efectos nocivos a una distancia considerable del punto de emisión.
La previsión de estas situaciones obliga a realizar estimaciones de las posibles emisiones que pueden tener lugar en diferentes instalaciones, de forma que pueda conocerse de antemano la magnitud aproximada del caudal de emisión por un supuesto orificio.
Entre las muchas circunstancias que pueden ser origen de emisiones peligrosas, aparece frecuentemente el fallo del propio equipo contenedor de la sustancia.
También es importante considerar otras situaciones de escapes por válvulas que se quedan abiertas o por venteos forzados en emergencias. Un ejemplo de esta última situación seria el fallo de la refrigeración en un recipiente de almacenamiento de un gas licuado inflamable refrigerado a baja presión, que daría como resultado un venteo forzado con una gran liberación de vapor.
Las situaciones que dan origen a la emisión de contaminantes se pueden clasificar de la forma siguiente:
Según el fluido:
gas/vapor.
líquido.
mezcla de vapor y líquido.
Según el equipo afectado:
recipientes.
otros equipos.
conducciones de tuberías.
Según la abertura:
rotura completa.
abertura limitada (válvula de alivio, disco de rotura, orificio, grieta, conexión, purga, toma de muestras, cierres de bombas, bridas, extremos o rotura de tuberías, etc.).
Según el recinto:
dentro de un edificio.
al aire libre.
Según la altura de la emisión:
a nivel inferior al suelo.
a nivel del suelo.
a nivel superior al suelo.
Según el impulso del fluido:
bajo impulso.
gran impulso.
Para estudiar la emisión debe conocerse la fase en que sale del recipiente. Como norma general puede adoptarse la siguiente: si el escape procede de un recipiente que contiene líquido a presión, normalmente saldrá líquido si la abertura está por debajo del nivel de líquido y vapor o mezcla de vapor y líquido si está por encima del nivel del líquido. Además debe tenerse en cuenta que para una diferencia de presión dada, el caudal másico de emisión es normalmente mucho mayor para un escape en fase líquida o mezcla de vapor y líquido que para gas o vapor. En la figura 1 se muestran los tres tipos de fugas concurrentes en un almacenamiento de gas licuado.
Fig. 1: Tres tipos de fugas en almacenamiento de gases licuados
Al estudiar una emisión deben considerarse diversos aspectos, entre los que cabe destacar:
La distancia entre válvulas de cierre y el tiempo de respuesta para determinar la duración de la fuga. Las Normas de Seguridad Federales para Gas Natural Licuado del Departamento de Transporte de EE.UU. especifican una duración de fuga de 10 minutos. Otras instituciones como la holandesa Rijnmond Public Authority, utilizan un tiempo de 3 minutos, si se dispone de un sistema de detección de fugas combinado con válvulas de aislamiento accionadas a distancia.
La fiabilidad de los sistemas técnicos de protección, tales como válvulas de retención, válvulas de exceso de flujo, etc.
La dependencia del tiempo en las fugas transitorias: los caudales de fuga son decrecientes conforme disminuye la presión aguas arriba del punto de emisión.
La reducción del flujo: las válvulas, bombas u otras restricciones en las tuberías que pudieran reducir el caudal por debajo del estimado a partir de la caída de presión y el área de la sección de descarga.
Un dato importante de partida para los cálculos es el tamaño del orificio. Si es una emisión prevista técnicamente se emplea la dimensión real de la válvula o tubería correspondiente. Cuando la emisión es accidental se debe hacer una estimación, que puede estar orientada mediante identificación de riesgos, considerando los posibles escenarios de incidentes.
Cabe destacar que los puntos más débiles de las instalaciones expuestos a fugas y que generan frecuentes incidentes son las bombas de impulsión de fluidos por fallos de estanqueidad en el sellado de las mismas, así como las bridas de unión entre las propias bombas y las tuberías. Adicionalmente y de mayor gravedad es la rotura accidental de una tubería por el impacto de un vehículo u otro elemento mecánico.
Para el estudio del movimiento de fluidos se deben recordar algunos conceptos fundamentales.
El caudal volumétrico Qv es el producto de la velocidad v del fluido por el área A de la sección transversal del flujo.
El caudal másico Qm es el producto de la velocidad v del fluido por el área A de la sección transversal y por la densidad r del fluido.
Recibe el nombre de velocidad másica G el producto de la velocidad v del fluido por la densidad r del mismo.
Para el estudio del escape de un líquido por un orificio en un recipiente, se parte de la ecuación de Bernouilli:
p + rgh + (rv2/2) = Cte.
p = presión absoluta (Pa)
r = densidad del líquido (kg/m3)
h = altura respecto a un nivel de referencia (m)
g = aceleración de la gravedad (m/s2)
v = velocidad del líquido (m/s)
Esta ecuación es válida para un flujo estacionario, no viscoso e incompresible y se cumple en dos puntos cualesquiera de una línea de corriente. Si consideramos como línea de corriente la vena del líquido que pasa por el orificio del recipiente y aplicamos la ecuación anterior a un punto 1 situado en la superficie de nivel del líquido y a otro punto 2 situado exactamente en el centro de la sección del orificio de salida, se tendrá:
p1 + rgh1 + (rv12/2) = p2 + rgh2 + (rv22/2)
Si el diámetro del recipiente es grande respecto al diámetro del orificio, la velocidad v1 de descenso de la superficie del líquido se puede considerar nula.
p1 = presión absoluta en el recipiente a nivel del líquido. Si el recipiente está abierto o en comunicación con la atmósfera, esa presión p1 será la atmosférica pa y si está cerrado y sometido a otra presión, p1 será la presión de saturación pv o cualquier presión genérica Ph.
P2 = pa (presión atmosférica)
La expresión anterior se podrá escribir de la siguiente forma:
Ph - Pa + rg (h1 - h2) = rv22/2
h1 - h2 = h = diferencia de altura entre el nivel del líquido y el orificio de salida.
y la velocidad teórica de salida del líquido será:
v2 = [2 (ph - pa) r -1 + 2 gh]1/2
El caudal másico de emisión Qm se define como la masa de fluido que sale por unidad de tiempo y su expresión es:
Qm = área de sección efectiva · velocidad · densidad
El área de la sección efectiva es el producto del área A de la sección transversal del orificio por un coeficiente de contracción Cc. La velocidad real es menor que la teórica calculada debido al rozamiento y a la viscosidad, por lo que la velocidad teórica biene afectada por un coeficiente de velocidad Cv. El coeficiente de descarga Cd se define como la relación entre el caudal real que sale por un orificio y el teórico. Puede expresarse también en función de Cv y de Cc, es decir Cd = Cv x Cc con lo cual se tiene en cuenta la influencia del rozamiento y la contracción de la vena de corriente por el orificio de salida. El coeficiente de descarga Cd es aproximadamente 0,62 para un orificio de borde afilado y la unidad para una boquilla de borde redondeado. Para un tramo muy corto de tubería o una boquilla externa conectadas a un recipiente, este coeficiente es aproximadamente 0,81.
El caudal másico Qm será:
Qm = Cd Ar [2 (ph - pa) r -1 + 2gh]1/2
Si la presión de almacenamiento es la atmosférica, se tendrá ph = pa y la expresión anterior se transforma en:
Qm = Cd Ar (2gh)1/2
En caso de escape de un gas licuado, el caudal de emisión disminuye rápidamente debido a la caída de presión y temperatura conforme va escapando (proceso adiabático) para lo cual se necesita su estudio particular, aunque en primera aproximación y si se considera una emisión de corta duración en un recipiente a presión, pueden aplicarse las fórmulas presentadas.
Durante la emisión la altura h varía con el tiempo y depende de la forma y posición del recipiente. Para un recipiente cilíndrico vertical con una sección transversal Ag (área proyectada sobre el suelo) y con una altura inicial h0 de líquido, se demuestra matemáticamente que la masa m escapada en el tiempo t es:
y el caudal másico Q en función del tiempo t es:
Cuando el escape se produce en un orificio de una tubería, se puede considerar como un escape de un recipiente al cual esté conectada. Si la tubería es corta y no existe una bomba de trasiego, el caudal másico se puede calcular con las fórmulas anteriores.
En el caso de conducciones largas, hay que tener en cuenta la fricción del líquido en las paredes de la tubería. Los codos, reducciones, válvulas y otros accesorios, aumentan las resistencias al paso del fluido y se contabilizan en forma de longitud equivalente de tubería o según la pérdida de carga que ocasionan. La pérdida de carga se da en alturas de columna del líquido (m) y la pérdida equivalente de presión en unidades de presión (Pa).
La pérdida de carga en un flujo estacionario a través de una tubería de sección circular se calcula mediante la expresión de Darcy - Weisbach.
hf = pérdida de carga (m)
f = factor de fricción de Moody
l = longitud de tubería (m)
d = diámetro de tubería (m)
v = velocidad del líquido (m/s)
g = aceleración de la gravedad (9,81 m/s2)
El factor de fricción de Moody f depende de la rugosidad relativa de la pared interna de la tubería (razón entre la rugosidad absoluta e y el diámetro d de la tubería) y del número de Reynolds Re.
Re = número de Reynolds (adimensional)
p = densidad del líquido (kg/m3)
µ = viscosidad dinámica o absoluta (N s/m2)
n = viscosidad cinemática (m2/s)
El número de Reynolds Re es un parámetro adimensional de gran importancia en el estudio del flujo de fluidos y con la experimentación se ha encontrado que para valores de Re menores que aproximadamente 2000 el flujo es generalmente laminar y para valores superiores a 4000 el flujo es generalmente turbulento. No existe un cambio brusco de flujo laminar a turbulento, sino que la transición se produce gradualmente a medida que aumenta la velocidad. El valor concreto de la velocidad dependerá de la geometría y de la rugosidad de la tubería. Cuando el número de Reynolds Re está comprendido entre 2000 y 4000 el flujo está cambiando rápidamente de laminar a turbulento y la frecuencia con que se presenta el flujo turbulento aumenta a medida que el Re se acerca a 4000.
El factor o coeficiente de fricción de Moody puede deducirse matemáticamente y su valor es f = 64/Re en el caso de régimen laminar. Para régimen turbulento no se dispone de relaciones matemáticas sencillas para obtener el valor de f en función del número de Reynolds y además los investigadores de este campo han encontrado que sobre el valor f también influye la rugosidad relativa de la tubería. La fórmula más aceptada para todas las tuberías es la ecuación de Colebrook:
En esta ecuación la incógnita f está en forma implícita y su resolución debe hacerse por iteración partiendo de un valor aproximado. Otra forma de resolverla es mediante un diagrama que da las relaciones existentes entre el coeficiente de fricción f, el número de Reynolds Re y la rugosidad relativa e/d.
Puede observarse que para tuberías lisas el valor e/d es muy pequeño y puede despreciarse el primer término entre paréntesis de la ecuación anterior quedando en la forma:
que se transforma en:
En la ecuación de Colebrook, para números de Reynolds Re muy grandes, el segundo término del paréntesis es despreciable y en tales casos la viscosidad no influye prácticamente y el coeficiente de fricción f tan solo depende de la rugosidad relativa e/d de la tubería. Este hecho se pone de manifiesto en el diagrama ya que las curvas se vuelven horizontales para números de Reynolds elevados.
La rugosidad de la superficie de una tubería tiene una importancia considerable para valores elevados del número de Reynolds como puede observarse en el diagrama. En la tabla 1 se dan valores de la rugosidad absoluta para diversos materiales. Cuando las tuberías están corroídas, el valor de la rugosidad aumenta y puede ser hasta diez veces mayor.
Tabla 1: Valores de la rugosidad absoluta e |
Para el cálculo de la pérdida de carga es necesario partir de una velocidad v del líquido, que en el caso de una tubería será inferior al calculado en el escape directo de un recipiente unido a ella. En el caso corriente de la existencia de una bomba impulsora se conocería el caudal en función dela altura manométrica, con lo cual se podría calcular la velocidad del líquido para el diámetro existente de la tubería. Conociendo la velocidad, se podrá calcular el número de Reynolds y mediante el diagrama de la fig. 2 se obtendrá el factor de fricción de Moody y en consecuencia se podrá calcular la pérdida de carga hf del tramo de tubería hasta el punto en que se considere el posible escape.
Fig. 2
Para el flujo de un líquido a través de una tubería, las pérdidas de carga por fricción están dadas por
Shf = hft + hfa + hfc
Shf = pérdida por fricción total
hft = pérdida de carga por fricción debida a la tubería
hfa = pérdida de carga por fricción debida a accesorios
hfc = pérdida de carga por fricción debida a aumento o contracción brusco de la vena líquida
Para la pérdida de carga por fricción en la tubería ya se ha indicado el proceso a seguir.
La pérdida de carga por fricción debida a cada accesorio instalado y a contracciones bruscas de la conducción se expresa mediante la fórmula:
hf = kv2/2g
Una lista de valores del coeficiente k se expone en las tablas 2, 3 y 4.
Tabla 2: Valores del coeficiente k para el cálculo de la pérdida de carga |
Tabla 3: Coeficiente k para contracción brusca del diámetro de tubería. (Ref. 7) |
Tabla 4: Valores del coeficiente k para tuberías con contracciones y ensanchamientos. (Según King, Handbook of Hydraulics. McGraw-Hill Book Co.) |
A veces interesa determinar la distancia que alcanza el escape de un chorro de líquido desde una abertura horizontal de un depósito para determinar si puede sobrepasar un cubeto de retención y la distancia hasta incidir sobre el suelo. Para una velocidad de descarga v (previamente calculada con lo indicado en el apartado 3), la distancia horizontal x recorrida en el tiempo t es:
x = vt
y la distancia vertical y, diferencia de cotas entre orificio y suelo es:
y = gt2/2
Eliminando t entre estas dos ecuaciones, la distancia x alcanzada es:
x = v (2y/g)1/2
Calcular el caudal de emisión en una esfera de almacenamiento de propano licuado a través de un orificio de 10 mm de diámetro. La presión de almacenamiento es la presión de vapor del propano a la temperatura de 25 ºC y el punto de descarga está situado a 2 m por debajo del nivel del gas licuado.
Datos
Densidad del propano licuado = 490 kg/m3
Presión de vapor del propano a 25 ºC = 9,3 bar abs.
Coeficiente de descarga Cd = 0,61
Solución: La presión a emplear en la expresión del caudal másico es la presión absoluta interna de la fase gaseosa, que será la presión de vapor saturado del propano en equilibrio con la fase líquida (Ph = 9,3 bar = 9,3.105 Pa).
Qm = Cd Ar {[2 (ph - pa)/r] + 2 gh}1/2 = 1,37 kg/s
El cálculo realizado es válido como valor inicial de la emisión, ya que el caudal de descarga variará con el tiempo, debido a la disminución de la presión absoluta interna p y la altura h de columna de liquido, conforme vaya progresando el derrame hasta su control. Para una mayor aproximación a una emisión real también se debería considerar que el proceso se acerca más a un proceso adiabático que isotérmico.
Un tanque cilíndrico vertical de almacenamiento de acrilonitrilo a presión atmosférica, de 5 m de diámetro y 4 m de altura y con nivel de líquido 3 m descarga a través de una tubería horizontal de acero estirado de 50 mm de diámetro, conectada a la base del tanque según el esquema de la figura 3. Calcular la cantidad de acrilonitrilo derramada en un accidente en que resulta seccionada esa tubería a 100 m de su conexión al tanque y suponiendo que se tarda 3 minutos en accionar la válvula de cierre más próxima. Los accesorios instalados entre el tanque y el punto de derrame son: 3 válvulas de compuerta, 1 válvula de retención, 1 bomba centrífuga, una válvula de control de asiento y 4 codos de 90º. Entre el punto de derrame y el tanque receptor hay 1 válvula de retención y 2 codos de 90º.
Fig. 3: Croquis de la conducción de distribución de acrilonitrilo en la que se produce una rotura
Datos
Viscosidad a 25 ºC = 0,34 cP (1 cP = 1 mPa.s)
Densidad = 0,8 kg/dm3
Caudal másico trasegado por la bomba = 0,750 kg/s
Rugosidad absoluta de la tubería = 0,00024 cm
Coeficiente de válvula de compuerta: k = 0,25
Coeficiente de codo 90º: k = 0,75
Coeficiente de válvula de retención: k = 2
Coeficiente de válvula de control abierta: k = 3
La curva característica de la bomba se da en la figura 4.
Solución:
Caudal másico Qm = 0,750 kg/s
Caudal volumétrico Qv = Qm/r = 9,375.10-4 m3/s
Aplicando el teorema de Bernouilli, queda la expresión en la que la altura manométrica HD que disponemos en el tanque de almacenamiento más la altura proporcionada por la bomba HB debe equilibrar la altura hDF a que se debe elevar el líquido más la pérdida de carga total shf.
HD + HB = HDF + Shf HD = 3 m HDF = 5 m
La pérdida de carga total se calcula a partir del tramo de tubería más los accesorios existentes entre los tramos considerados. Siguiendo el croquis e incluyendo el factor k correspondiente a cada accesorio, obtenido de la tabla 2 tenemos:
Shf = pérdidas de carga de [conexión del tanque de almacenamiento a tubería a ras de la pared + 3 válvulas de compuerta + 2 válvulas de retención + 1 válvula de control + 6 codos 90º + conexión de tubería a tanque receptor (tubería entrante ) + 100 m de tubería de d = 50 mm]. El tramo de tubería entre el seccionamiento y el tanque receptor se ha despreciado en comparación con los accesorios instalados, debido a su corta longitud.
Shf = (0,50 + 3 · 0,25 + 2 · 2 + 1 · 3 + 6 · 0,75 + 1) v2/ 2g + (f · l / d) · (v2/ 2g)
sabiendo que
v = Qv / A = Qv/ pd2 / 4; v2 = 16Qv2 / p2d4
y sustituyendo
Shf = [0,50 + 3 · 0,25 + 2 · 2 + 1 · 3 + 6 · 0,75 + 1 + f 100/ 0,05] 8 Qv2/ (p2 d4g)
El valor del coeficiente de fricción de Moody depende del número de Reynolds Re y de la rugosidad relativa e/d. A su vez Re depende de la velocidad de circulación del fluido, por lo que partiremos de la velocidad estacionaria antes de ocurrir el accidente.
Área de la sección de tubería
A = p d2 /4 = p 0,052 /4 = 1,963 · 10-3 m2
Velocidad
v = Qv / A = 0,477m/s
Re= vdr/µ = 5,61 · 104
e/d = 4,8 - 10-5
Mediante estos dos valores en el diagrama de la figura 2 se obtiene un coeficiente de fricción f = 0,021 con lo que
Shf = 0,648 m
Al ocurrir el accidente con seccionamiento de la conducción cambian las condiciones de trabajo de la bomba centrífuga, para lo cual se debe establecer la nueva ecuación de equilibrio:
HD + HB = HDF + shf HD= 3m HDF = 0 m
Shf = (0,50 + 3 · 0,25 + 1 · 2 + 1 · 3 + 4 · 0,75) v2 / 2g + (f l v2 ) / (d2g)
La pérdida de carga total se puede expresar en función del caudal volumétrico Qv sabiendo que v = Qv / A.
Shf = (0,50 + 3 · 0,25 + 1 · 2 + 1 · 3 + 4 · 0,75 + f · 100 / 0,05) (Qv / A)2 / 2g=
= (9,25 + f · 2 000) 8 Qv2 / (p2 d4g)
Aunque el factor de fricción de Moody también depende del caudal (y por consiguiente de la velocidad ) y del número de Reynolds, su variación es muy inferior a la variación del caudal al cuadrado, por lo que la pérdida de carga se puede expresar aproximadamente en función de Qv2 manteniendo el mismo factor de fricción de Moody f obtenido anteriormente o uno ligeramente inferior ( p. e. 0,016 ). Utilizando el valor f = 0,021 se tiene:
Shf = (9,25 + 0,021 · 2 000 ) 8Qv2 / (p2 · 0,054 · 9,8) = 678 231,6Qv2
La ecuación de equilibrio pasaría a ser:
3 + HB = 0 + 678 231,6Qv2 HB = - 3 + 678 231,6Qv2
Supongamos un Qv = Qv1 = 1,5 I/s = 1,5 · 10-3 m3 / s
La curva característica de la bomba centrífuga nos da para ese caudal una altura manométrica de 6 m.
Sustituyendo el valor Qv1 = 1,5 · 10-3 m3 /s en la ecuación tenemos
HB = - 3 + 678 231,6 (1,5 · 10-3 )2 = - 1,47 distinto de 6, por lo que no se da por bueno el caudal supuesto.
Elegimos otro valor Qv = Qv2 = 2 l / s = 2 · 10-3 m3/s
La curva característica nos da una altura manométrica de 5 m.
HB = - 3 + 678 231,6 (2. 10-3)2 = - 0,28 distinto de 5. Tampoco es válido el valor supuesto.
Elegimos otro valor Qv = Qv3 = 3 I/s = 3 · 10-3 m3/s
La curva característica nos da una altura manométrica de 3 m.
HB = - 3 + 678 231,6 (3 · 10-3)2 = 3,10 semejante a 3
con lo que se puede dar por finalizada la iteración, ya que con un gráfico es difícil una mayor aproximación.
Debe observarse que al ocurrir el seccionamiento de la conducción, durante unos segundos se tendrá un transitorio mientras el caudal volumétrico que proporciona la bomba va pasando del caudal inicial estacionario de 0, 9375 I/s al estabilizado a unos 3 I/s, hasta su control. En este documento no se ha considerado el estudio de estos transitorios.
Se han supuesto tres minutos hasta el corte de la emisión, por lo que aproximadamente se derramarán
3 I/s · 3 min · 60 s/min = 540 litros de acrilonitrilo.
En el ejercicio no se ha planteado la consiguiente evaporación del acrilonitrilo derramado, producto muy tóxico y que requeriría severas medidas de intervención.
Fig. 4: Curva característica de la bomba centrífuga instalada en la red de distribución del acrilonitrilo
(1) AMERICAN INSTITUTE OF CHEMICAL ENGINEERS
Guidelines for Chemical Process Quantitative Risk Analysis
Nueva York, American Institute of Chemical Engineers, 1989, 585 págs.
(2) CÁTEDRA DE MECÁNICA DE FLUIDOS. UNIVERSIDAD POLITÉCNICA DE
VALENCIA
Curso de Ingeniería Hidráulica
Madrid, Instituto de Estudios de Administración Local, 1987, 664 págs.
(3) COULSON, J.M. y RICHARDSON, J.F.
Chemical Engineering (Vol. 1)
Oxford, Pergamon Press, 1990, 6 vols.
(4) GILES, R.V.
Mecánica de los Fluidos e Hidráulica
Nueva York, Mc Graw-Hill, 1969, 273 págs.
(5) KING, H.W.
Manual de hidráulica
México D.F., UTEHA, 1992
(6) LEES, F.P.
Loss Prevention in the Process Industries
Londres, Butterwoth & Co. (Publishers) Ltd., 1980, 2 vols., 1316 págs.
(7) PERRY, R.H. y GREEN, D.
Perry's Chemical Engineers' Hanclibook. 6ª edición
Nueva York, Mc Graw-Hill, 1984
(8) TNO.
Methods for the Calculation of Physical Effects (Yellow Book). 2ª edición
Voorburg (Holanda), TNO, 1992