A quantitative risk-assessment system (QR-AS) evaluating operation safety of Organic Rankine Cycle using flammable mixture working fluid

Accepted Manuscript Title: A Quantitative Risk-Assessment System (QR-AS) Evaluating Operation Safety of Organic Rankine Cycle using Flammable Mixture Working Fluid Authors: Hua Tian, Xueying Wang, Gequn Shu, Mingqiang Wu, Nanhua Yan, Xiaonan Ma PII: DOI: Reference:

S0304-3894(17)30389-8 http://dx.doi.org/doi:10.1016/j.jhazmat.2017.05.039 HAZMAT 18598

To appear in:

Journal of Hazardous Materials

Received date: Revised date: Accepted date:

9-2-2017 17-5-2017 21-5-2017

Please cite this article as: Hua Tian, Xueying Wang, Gequn Shu, Mingqiang Wu, Nanhua Yan, Xiaonan Ma, A Quantitative Risk-Assessment System (QR-AS) Evaluating Operation Safety of Organic Rankine Cycle using Flammable Mixture Working Fluid, Journal of Hazardous Materialshttp://dx.doi.org/10.1016/j.jhazmat.2017.05.039 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A Quantitative Risk-Assessment System (QR-AS) Evaluating Operation Safety of Organic Rankine Cycle using Flammable Mixture Working Fluid

Hua Tian, Xueying Wang, Gequn Shu*, Mingqiang Wu, Nanhua Yan, Xiaonan Ma State Key Laboratory of Engines, Tianjin University, People’s Republic of China *Corresponding author: Gequn Shu, State Key Laboratory of Engines, Tianjin University, No. 92, Weijin Road, Nankai Region, Tianjin 300072,People’s Republic of China Tel: +086 022-27409558 E-mail address: [email protected]

Highlights:     

A comprehensive Quantitative Risk-Assessment System (QR-AS) evaluating operation safety of Organic Rankine Cycle is proposed. The QR-AS covers analysis of concentration distribution, explosion risk assessment and mitigation measures. Explosive zone, damage range of explosion and risk prevention measures can be figured out by QR-AS. Results of QR-AS can provide valuable guidance for ORC application with safety operation. A case of propane-carbon dioxide leaking from ORC is exemplified applying the QR-AS.

Abstract Mixture of hydrocarbon and carbon dioxide shows excellent cycle performance in Organic Rankine Cycle (ORC) used for engine waste heat recovery, but the unavoidable leakage in practical application is a threat for safety due to its flammability. In this work, a quantitative risk assessment system (QR-AS) is established aiming at providing a general method of risk assessment for flammable working fluid leakage. The QR-AS covers three main aspects: analysis of concentration distribution based on CFD simulations, explosive risk assessment based on the TNT equivalent method and risk mitigation based on evaluation results. A typical case of propane/carbon dioxide mixture leaking from ORC is investigated to illustrate the application of QR-AS. According to the assessment results, proper ventilation speed，safe mixture ratio and location of gas-detecting devices have been

proposed to guarantee the security in case of leakage. The results revealed that this presented QR-AS was reliable for the practical application and the evaluation results could provide valuable guidance for the design of mitigation measures to improve the safe performance of ORC system.

Key words: Organic Rankin Cycle, hydrocarbon/carbon dioxide mixture, leakage, quantitative risk assessment

Nomenclature Abbreviation CPR LFL ORC QR-AS RD UFL Symbols a Cμ C1

C2 

D E ET hj Jj

k keff L p P P1 P1c Pa Δp Qf QTNT qv ref R R0.5 Rd, 0.5

Critical Pressure Ratio lower flammability limit Organic Rankine Cycle quantitative risk assessment system relative difference upper flammability limit efficiency factor k-ε model specific constant k-ε model specific constant k-ε model specific constant coefficient total specific energy the total explosion energy (J) enthalpy of species j (J) diffusion flux of species j (kg/m2·s) turbulent kinetic energy (m2/s2) effective thermal conductivity distance to the explosive source (m) pressure (Pa) Probability (%) the pressure at the leak position (Pa) critical pressure (Pa) ambient pressure (Pa) the overpressure of explosion shock wave (kPa) combustion heat of TNT per kilogram (kJ/kg) combustion heat of releasing gas per kilogram (kJ/kg) leaking rate (m3/s) the value in reference gas constant radius of death zone (m) serious injury radius (m)

Rd, 0.01 u ui ui 'u j '

Wf WTNT Sh t Y Yi Z Greek letters ρ γ ε μ σk σε φ

slight injury radius (m) release rate (m/s) velocity along the xi direction (m/s) “ Reynolds stresses” the total mass of leaked gas (kg) the equivalent mass of TNT source term time (s) probit variable mass fraction of specie i the explosion characteristic length density (kg/m3) adiabatic index turbulent dissipation rate (m2/s3) viscosity of the medium (kg/ms) Prandtl number for k Prandtl numbers for  flow coefficient

Introduction To recover the waste heat of internal combustion engines, Organic Rankine Cycle (ORC) is a promising solution with outstanding system performance. Many studies have illustrated that hydrocarbons are better choice as working fluids of ORC for their good matching with high temperature engine exhaust [1, 2]. However, such compounds are usually flammable, which limits their practical application [3]. Zeotropic mixtures based on hydrocarbons and carbon dioxide may be good alternatives to solve this issue owing to better environmental characteristics and excellent thermo-physical properties [4, 5]. However, the leakage of working fluid from ORC system is unavoidable in practical application, which may cause explosion once the leaked gas reaches to a critical concentration. Thus, in view of security, quantitative risk assessment of such mixtures leakage must be paid full attention. Nowadays, researches on flammable mixtures leakage are mainly focus on the composition change and the effect on cycle performance. Kim and Didion [6] simulated the isothermal and adiabatic leak process of zeotropic mixtures to study the mass fraction changes in refigenrating system. A. Johansson [7] proposed a method to estimate the circulated composition of mixture in heat pump system under different leakage scenarios by measuring temperature and pressure. Ren et al. [8] predicted the influence of leakage on the cycle performance of auto-cascade refrigeration system. Regarding to the safety of hydrocarbon/carbon dioxide mixture, many theoretical models have been developed to predict their flammability limits [9-11]. These studies are focused on mixture safety without combination of practical application and still

have limitation for the lack of risk assessment on ORC safe performance. Recently, risk assessment of hazardous gas leakage has been extensively applied in other fields. Jose Luis et al. [12] developed a methodology for qualitative risk assessment in industrial plants, based on the combination of Hazard and Operability analysis with fuzzy numbers. Fu et.al [13] proposed a framework for assessing the leakage risk of LNG-fueled vessels by using event tree analysis and CFD simulations. Han et.al [14] developed an integrated risk analysis method for natural gas pipeline network. Elena Stefana et. al [15] proposed a risk assessment model for a new Dual Fuel (LNG-Diesel) system of heavy-duty trucks, which combined Reliability Block Diagram, bow-tie analysis, Fault Tree Analysis and risk matrix. Zhang et al. [16] performed numerical simulations to study explosion hazards of liquefied petroleum gas explosion in vented enclosure with obstacles. L. Ehrhardt et al. [17] evaluated the accuracy of two blast threat assessment tools on the basis of experimental data and proposed a new fit model for quick decision aid applications in the battlefield. However, a limitation can be observed regarding to these methods for not comprehensively analyzing the dispersion characteristic of hazardous gas. As a matter of fact, adequate dispersion model is a fundamental part in risk assessment, which enables an accurate estimation of concentration and positive impact on both human and financial resources. Recently, numerical simulation based on computational fluid dynamics (CFD) has been proved as a powerful tool to investigate gas leakage and dispersion [18-20]. Ryuichi [21] simulated the gas spread in residential space for traceable hazard analysis by using an open-source software OpenFOAM. M. Siddiqui [22] proposed a CFD based model for indoor risk assessment considering accidental release of chlorine in an industrial indoor environment. He et al [23] investigated the leakage mechnism in different leaking conditions for long-distance pipelines through numerical simulations and experiments. Zhang et al. [24] simulated the three dimensional dispersion and conversion behaviors of silicon terachloride release in polysilicon industry. For the leakage of flammable gas, flammability limits are the thresholds of explosion. Thus, analysis of temporal and spatial evolution of the concentration field combined with flammability limits must be a critical part of risk assessment. With the analysis above, there is a lack of risk assessment for ORC related to flammable mixture working fluid leakage, and more comprehensive analysis on the dispersion characteristic of hazardous gas is needed concerning the existing risk assessment models. In response to the current gaps, this paper aims at establishing a comprehensive quantitative risk assessment system (QR-AS) evaluating operation safety of Organic Rankine Cycle using flammable mixture working fluid. In QR-AS, numerical simulations of mixture leakage is applied to establish fundamental data for the following risk assessment based on the TNT equivalent method. In addition, specific risk mitigation measures including the design of ventilation speed, optimization of safe mixture ratio and location of gas-detect device are determined on the basis of further simulations.

Description of QR-AS The QR-AS proposed in this paper is comprised by three parts, as the analysis of concentration distribution, explosive risk assessment and risk mitigation measures. They are illustrated in the schematic diagram in Fig.1 and are explained in details as follows. Part 1 of the QR-AS focuses on analysis of the concentration distribution of hazardous gas. For risk assessment of flammable gas leakage, the evaluation of the concentration field under different scenarios is significant for investigating the diffusion patterns of leakage gas. Besides, it is a critical index to determine the explosive risk related to flammability limits. Part 1 is the fundamental part for risk evaluation of leakage which is essential for further investigation. The detailed modeling is presented in Section 2.1. Part 2 of the QR-AS focuses on the evaluation of explosion consequences based on the TNT equivalent method. It is generally agreed that the major consequence of explosion are connected with blast effects, which is usually of the prime interest for determination of the hazard degree [25]. Through the evaluation of concentration in Part 1, the damage degree of explosion can’t be figure out. Thus, Part 2 aims at dividing the explosive zone where explosion will occur in the presence of ignition source. After that, the impact of explosion overpressure on buildings and individuals is evaluated. The evaluation results in Part 2 will provide dependable data to risk analysis. The detailed modeling is presented in Section 2.2. Part 3 of the QR-AS investigates the effectiveness of specific risk mitigation measures through further simulations. It will provide reasonable scientific theory to prevent gas explosion which is of great significance to improve the security of the ORC system and reduce potential hazard caused by leakage. The detailed modeling is presented in Section 2.3. The evaluation procedure is basically shown in the flow diagram of Fig.2. This procedure clearly reflects the methodology of the QR-AS. In addition, the principle is generally applicable to other ORC with any kind of flammable working fluids by adjusting specific parameters correspondingly using QR-AS. 2.1 QR-AS Part1: Analysis of concentration distribution 2.1.1 Release rate calculation The leakage and diffusion of working medium can be simplified as a jet sprayed into a large room, then mixed with the ambient air. The calculation process of release rate includes three steps as described below. Step 1 determines the operating parameters (temperature and pressure) of the leak source by thermodynamic calculations or experimental measurements. Step 2 calculates the Critical Pressure Ratio (CPR) with equation (1). Since the hole release doesn’t have an expansion section, the maximum speed can only reach the sonic speed which is depended on the CPR [26]. 

P  2   1 CPR  a    P1c    1 

(1)

where Pa is the ambient pressure, P1c is the critical pressure at the leakage point, γ is the adiabatic index of leaking gas. Step3 determines the release rate by following equations [27, 28]: When

Pa Pa  , the flow is subsonic and the release rate can be calculated by equation P1 P1c

(2),  1    Pa    2  u  RT 1      1 1   P1    

When

(2)

Pa Pa  , the flow is sonic and the release rate can be calculated by the follow P1 P1c

equation:  1      P 2 a  u  RT 1       1 1   P1 c   

(3)

where R is the gas constant, T1 is the temperature of leakage gas, P1 is the pressure at the leak position,  is the flow coefficient, 0.97~0.98 [29]. 2.1.2 Numerical simulation of working medium leakage and dispersion In this work, numerical simulation tool Fluent 16.0 is used to study the special and temporal evolution of concentration field. The separated solver and the implicit scheme model are chosen. Simulations are performed by SIMPLE algorithm. Convective terms of the transport equation are discretized using second-order upwind scheme. Step 1 builds up physical model and grid model, since they are foundations of subsequent simulations. Step 2 selects the conservation equations of mass, energy, momentum and species transport. Continuity equation:     ui   0 t xi

(4)

where  is the density, t is the time, and ui is the velocity component along xi direction. Energy equation:

     E     v   E  p      keff T   h j J j   eff  v   Sh t j  

(5)

where E is the total specific energy, keff is the effective thermal conductivity, h j is

the enthalpy of species j, J j is the diffusion flux of species j, and S h is the source term. Momentum equation:

   p   u  ui    uiu j       i   ui 'u j '  t xi xi x j  x j 

(6)

where  is the viscosity of the medium，and ui 'u j ' is the “ Reynolds stresses”. The standard k-ε turbulence model is the most commonly used model for engineering applications, and has been widely validated in dispersion simulations [30-33]. Compared with other models, such as large eddy simulation or algebraic stress models, it has a lower computational cell requirement and faster computational run times [34,35]. Hence，the standard k-ε turbulence model is chosen in this work.

    ku j   k   t x j x j





 t    k 

    u j      t x j x j





 u i u j  k     t    x j   x j xi

 t     

 u i     x j

      u i u j    C1  t    x k  x    j   j xi

(7)  u i  2  C    2 k  x j 

(8) The turbulent viscosity is linked to the turbulent kinetic energy and dissipation via the relation,

t   C

k2



(9)

The model constants are given by

C1 =1.44, C2 =1.92, C = 0.09,  k =1.0 and   =1.3 To determine the concentration of the leakage gas, species transport equation need to be solved by treating the leakage gas as a non-reacting species.

  Yi      uYi    J i t

(10)

Step 3 inputs boundary conditions and initial conditions. Boundary conditions include the leaked gas inlet, air inlet, outlet and wall. The initial thermodynamic state should be determined in accordance to the actual ambient temperature and pressure. Prior to the simulation of the working medium leakage, the flow field in the laboratory due to the ventilation flows should be calculated and used as initial state for the leak simulations. 2.1.3 Model validation A validation of the models has been performed first by recurring the existing data reported in the literature [36] before putting them into application. The RD (relative difference) is determined as follows:

RD 

sim  ref ref

100%

(11)

where sim symbol means the simulated value, ref symbol means the value in reference. In ref [36], simulations of the continuous release and dispersion process of carbon dioxide in a ventilated room are carried out by FLUENT. The solution of mathematical model is based on the equations of mass, momentum, energy and species transport conservation as well as the standard k-ε turbulence model. In addition, an experiment about carbon dioxide leaking into the room from a high pressure cylinder has been carried out. The validation of the simulation method is made by comparison of simulation results with experimental data, which proves the correctness of their results. To study the leakage scenario in ref [36], the same model and boundary conditions are used. The relative difference between our simulation results with reference [36] is shown in Fig.3. The maximum relative deviations are 6% at 50s and 5.8% at 100s, while the average relative differences are 2.1% and 3.3% respectively. The differences are very small and the relative deviations are in an acceptable range. So the models are considered to be correct to predict heavy gas leakage and are adopted to carry out the subsequent simulations. 2.2 QR-AS Part2: Explosive risk assessment 2.2.1 Acquire the flammability limits of leaked gas The flammability limits include lower flammability limit (LFL) and upper flammability limit (UFL). It is an essential feature of flammability characteristics which define the range of concentrations for flame propagation to occur. The flammability limits can be obtained from related literature [11]. 2.2.2 Divide explosive zone The explosive zone is defined as the area where the propane concentration is between LFL and UFL. Based on the flammability limits and predicted concentration, the explosive zone can be figured out by post-processing software CFD-POST 16.0. 2.2.3 Explosive risk assessment based on TNT equivalent method TNT equivalent weight method is a widely used model for simulating the explosion overpressure which can obtain the satisfactory results [37]. Overpressure is the pressure caused by a shock wave over normal atmospheric pressure. It is formed by compression of air in front of a blast wave which heats and accelerates air molecules. In order to predict the overpressure at specified distance, the energy available in explosive gas is converted into an effective explosive weight of TNT based on following equations:

WTNT 

aW f Q f QTNT

W f   qvt

(12) (13)

where WTNT is the equivalent mass of TNT, a is an efficiency factor, 0.04 [29]. Wf is

the total mass of explosive gas. The overpressure of explosion shock wave can be calculated by the Henrych empirical equation [29]: 1.379 0.543 0.035 0.0006 (0.05  Z  0.3) p   2  3  (14) Z Z Z Z4

p  p 

0.6076 0.032 0.209  2  3 Z Z Z

( 0 . 3Z 

0.0649 0.3973 0.3226   Z Z2 Z3

( 1 Z  1 0 )

Z

3

L WTNT

1)

(15) (16) (17)

where p is the overpressure, L is the distance to explosive source , Z is the explosion characteristic length in propagation. For the purpose of predicting the influence range of explosion, vicinity of explosion source is divided into four regions: death zone, serious injury zone, slight injury zone and safety zone. The radius of death zone is recorded as R0.5, which can be determined by Van den Berg and Lannoy equation as follows [38]: W  R0.5  13.6  TNT   1000 

0.37

(18)

The scope of serious injury zone is R0.5
p  D  Rd ,0.5  0   ET 

(20)

ET  WTNT  QTNT

(21)

in which ET is the total explosion energy, D is a constant. Rd, 0.01 is the external radius of slight injury zone. With a corresponding shock wave overpressure 17kpa, it can be determined by equations (18)-(20). The scope of safety zone is R> Rd, 0.01. People in this zone will be free of harm, and the death probability is nearly zero. In addition, “Probit analysis” is used for assessing the damage probability of given injuries, since it has been widely applied in many hazard assessments [40,41]. Probit equations allow the correlation of overpressure to percentage of people affected by a certain damage. The probit is a connection between the percentage of people submitted to a certain effect with a given intensity (ΔP). The connection between

probit Y and probability P is given by following equation[42]

P

1 2



Y 5



 u 2  exp  du  2 

(22)

The probability range (0-1) is generally replaced in probit work by a percentage (range 0-100). The relationship between percentages and probit has been usually taken from Table3. The Probit equation for eardrum rupture is: Y  15.6  1.93ln P (23) Similarly, Probit equation for death from lung hemorrahage is Y  77.1  6.91ln P (24) 2.3. QR-AS Part3: Risk mitigation According to characteristic of concentration distribution and risk assessment results, relative mitigation measures should be proposed. The comprehensive protection measures should contain the following aspects: location of gas-detect alarm devices, design of ventilation system based on the effect of different wind speed, optimization of the mixture ratio based on the relative simulations results.

3. A case of QR-AS application 3.1. The physical model and assumption Fig.4 shows the ORC system for engine waste heat recovery in the authors’ laboratory. The physical model (length 4.5m, width 3.4m, height 2.5m) with a simplified ORC system is shown in Fig.5 while the grids model is shown in Fig.6. As working fluid, propane/carbon dioxide (0.7/0.3) leaks from a hole with a diameter of 10mm located at the outlet of the pump. Leakage is most likely to occur at the outlet of pump where the pressure is highest. Besides, the pump is always located near the ground where the dense gas is easier to accumulate to higher concentration. According to the thermal dynamic calculation by MATLAB, the operating parameter is 307K, 5MPa at the outlet of the pump. The release rate calculated by equation (3) is 245m/s. Simple algorithm is adopted in the transient simulations with the time step 0.01s. As it is shown in Fig.6, the hybrid mesh in computational domain is comprised by unstructured tetrahedron grid and hexahedron grid. At the initial time, propane/carbon dioxide (0.7/0.3) is specified for the gas inlet while 100% air is specified for the air inlet. Detailed boundary conditions are listed in Table 1. The following assumptions are made: (1)The released mixture turns into the gas phase completely at the leakage location. (2)The change of the pressure of leakage gas is ignored and there is no chemical reaction. (3) Wind speed and direction do not change with time, location and height. 3.2. Results for Part 1: analysis of concentration distribution The concentration distribution of propane is a significant concern for the reason that if the concentration exceed LFL, explosion may happen in the presence of ignition. In this paper, concentration is expressed in terms of mole fraction. The plane where the leakage hole exists (y=1.23m) is chosen to study the distribution characteristic at different time. As Fig.7 shows, propane is sprayed upwards from the leakage hole at a

high initial velocity and reaches to the roof quickly, and then spreads around the wall. When it touches the side walls, it will flow in the opposite direction. Such gas will blend with the wind and form two main recirculation regions: one on the upper left corner of the room while the other on the right corner near the ground. Fig.8 shows that the concentration on the windward side of the evaporator is lower than the leeward. It indicates that the presence of evaporator prevent the wind flow to the other side to some extent. When the wind comes across obstacles, it may change direction and lead to a lower concentration on the windward side of the evaporator. Besides, taking the recirculation regions into account, the concentration on the leeward side becomes even higher. To account for the gravitational influence on the leaked gas, the concentration distribution at different heights is shown in Fig.9. As we can see, the propane concentration increases with the height reduction. Since the leaked gas is denser than air, this property determines that the leaking gas will accumulate near the ground. Another reason for the distinct effect of gravity is that the initial velocity decreases in the diffusion process, and the kinetic energy is not as high as that near the leakage source. After some distance from the leaking source, the concentration decreases dramatically because the obstacle prevents the gas flow to the other side to some extent. In Fig.10, propane concentration increases suddenly at the distance about x=1.73m, which is in corresponding to the area of gas leakage. As the jet from the leakage has an important effect on the dispersion, the leaked gas diffuses vertically upwards due to the high initial momentum. Because of the impact of wind, the concentration accumulates along the direction of wind. 3.3. Results for Part 2: Explosive risk assessment 3.3.1. Explosive zone division For propane/carbon dioxide (0.7/0.3), the LFL and UFL are 0.0204 and 0.08699 respectively [11]. When the wind speed is 0.5m/s, the scope of explosive zone at different time is shown in Fig.11. It illustrates that the explosive zone first appears above the leaking source at 10s and 50s. Along with the accumulated concentration, the explosive area expands significantly at 100s. The phenomenon illustrates that it is necessary to install alarm device to ensure timely security measures. 3.3.2. Explosive damage degree When conducting explosion risk assessment, two hypothetical explosion scenarios are discussed. Scenario A: The leakage amount is limited, only the leakage gas explodes. Scenario B: Once the explosion occurs, total propane in the system explodes due to the destruction of system. In Scenario B, the worst case with the maximum filling amount of propane/carbon dioxide is investigated to evaluate explosion damage. Based our experimental data, the maximum filling amount is 50kg. Fig.12 shows the radius of dangerous zones at different time in Scenario A. It can be seen that the external radius of death zone is 1.03m at 100s which represents the distance from the explosion site. The range of serious injury zone is 1.03m < R< 3.78m, while that of slight injury zone is 3.78m < R< 6.79m at 100s. The external radius of slight injury zone also represents the range of safety zone. When R > 6.79m,

people in this zone will be free from harm. Fig.13 shows the attenuation of peak overpressure with the increment of distance in Scenario A. When the distance from explosion site is relatively small (L<0.5m), the peak overpressure attenuate very quickly. After a certain distance, the attenuations of peak overpressures become slow. It can be seen in Fig 13 that at the distance of 2m from explosion site, the peak overpressure is 40kPa at 10s. That means heavy damage on reinforced concrete buildings based on the criteria in Table 2 According to the effects of blast overpressure in Table 2, the consequence caused by explosion at different distances can be approximately assessed in this way. In Scenario B, Probit analysis is used for estimating the damage probability of given injuries. With the overpressure estimated previously by equation (14) - (17), the Probit units for eardrum rupture and lung hemorrhage can be calculated respectively. The individual risk descending with the increase of distance away from explosion center can be seen in Fig.14 (a). It is indicated that the individual risk is almost unchangeable from the center to a certain distance 5m, because of the high fatality probability near the leakage source. When the distance is larger than 7m, the probability of lung hemorrhage is zero. In the eardrum rupture case, the probability of eardrum rupture is 100% when the distance is smaller than 6m. For the injury of eardrum rupture and lung hemorrhage, the threshold 22m represents a safe distance for individuals to be free from harm. These threshold values may be used as the significant guidance values for the design of risk prevention system. Fig.14(b) shows the attenuation of peak overpressure with the increment of distance in Scenario B. It is noteworthy that the tendency in Fig.14 (b) is similar to Fig13 while the value of overpressure is much larger. The most serious damage to buildings at specific distance can be figured out according to Table2. Results for Part3：Mitigation measures Location of the gas detectors In view of security，gas-detecting alarm device should be installed at proper location which it is helpful for relevant staff adopt emergency measures in time to mitigate consequences of release. Planes y=2.05m and y=2.55m with different obstacles existing are chosen to investigate the concentration distribution. In Fig.15, there are two regions with higher concentraion, one at the left corner near the roof while the other at the right corner near the floor. Due to the high initial release speed, the leaking gas reaches to the roof quickly and changes direction which lead to the formation of recirculation region at the left corner near the roof. Along with the decrease of the initial velocity in the diffusion process, the gravity and wind direction have an obvious effect on the dispersion tendency. The concentration accumulates near the floor due to the force of gravity wind. Taking previous analysis in section 3.2 into consideration, the higher concentration all occurs at the downwind zone in different cases. According to the relative requirements of “petrochemical flammable gases and toxic gas detection and alarm design specification” (GB50493-2009), gas-detecting alarm device is preferably installed at the downwind position under the leaking hole. 3.4.2. Ventilation method

From Fig.16, we can see that the scope of explosive zone changes dramatically under different wind speed, the explosive zone only appears on the top of the leaking source when the wind speed reaches to 1m/s. Because with the increment of wind speed, the transportation of gas cloud and dilution effect strengthen, so gas rapidly disperses correspondingly. The analysis of effect of wind speed will provide guidance to the design of ventilation system. Since the ventilation conditions play a significant part in the diffusion characters of the leaking gas, vent or air blower should be set. If the wind speed is slower than 1m/s, forced ventilation measures should be taken. 3.4.3. Mixture ratio Regarding to different ratio of propane, the explosive zone is shown in Fig.17. Along with the increment of propane ratio, the scope of explosive zone expended. Even though the operating conditions of ORC system is different under different mixture ratio, which may lead to different release rate. Analysis of the effect of mixture ratio on the explosive zone is helpful to provide guidance for the optimization of the safe ratio used in ORC system. According to Fig.17, propane ratio in the mixture can’t exceed 0.4 in view of explosive risk.

4. Conclusions The paper proposed a comprehensive risk assessment system (QR-AS) to evaluate the operation safety of Organic Rankine Cycle using flammable mixture working fluid, which covers concentration distribution analysis, explosive risk assessment and risk mitigation measures. A typical case of propane/carbon dioxide mixture leakage is investigated to illustrate the applications of the QR-AS. Concentration distribution analysis provides fundamental data for explosive risk assessment and specific mitigation measures. Based on the assessment results, risk mitigation measures have been proposed, including the safe ventilation speed, safe mixture ratio and location of gas-detecting device. The evaluation results can provide valuable guidance for ORC application with safe operation. Although the study has focused on a mixture of hydrocarbon and carbon dioxide, the QR-AS model is generic enough to deal with cases of any kind of flammable gas. The model in this paper can also be widely applied to refrigeration and pump systems.

Acknowledgements This work was supported by a grant from the National Natural Science Foundation of China (No. 51676133).

References: [1] M.A. Siddiqi, B. Atakan, Alkanes as fluids in Rankine cycles in comparison to water, benzene and toluene, ENERGY, 45 (2012) 256-263. [2] G. Shu, X. Li, H. Tian, X. Liang, H. Wei, X. Wang, Alkanes as working fluids for high-temperature exhaust heat recovery of diesel engine using organic Rankine cycle, APPL ENERG, 119 (2014) 204-217. [3] G. Shu, Y. Gao, H. Tian, H. Wei, X. Liang, Study of mixtures based on hydrocarbons used in ORC (Organic Rankine Cycle) for engine waste heat recovery, ENERGY, 74 (2014) 428-438. [4] M. Chys, M. van den Broek, B. Vanslambrouck, M. De Paepe, Potential of zeotropic mixtures as working fluids in organic Rankine cycles, ENERGY, 44 (2012) 623-632. [5] B. Dai, M. Li, Y. Ma, Thermodynamic analysis of carbon dioxide blends with low GWP (global warming potential) working fluids-based transcritical Rankine cycles for low-grade heat energy recovery, ENERGY, 64 (2014) 942-952. [6] M. S. Kim, & D. A. Didion, Simulation of isothermal and adiabatic leak processes of zeotropic refrigerant mixtures, HVAC & R RESEARCH, 1 (1995) 3-20. [7] A. Johansson, P. Lundqvist, A method to estimate the circulated composition in refrigeration and heat pump systems using zeotropic refrigerant mixtures, International Journal of Refrigeration, 24 (2001) 798-808. [8] N. Ren, G. Yan, J. Yan, Y. Wu. Study on leakage of mixture refrigerant in auto-cascade refrigeration cycle. Refrigeration, Cryo&Supercond, 25(2007):61-64. [9] C. Chen, H. Liaw, T. Wang, C. Lin, Carbon dioxide dilution effect on flammability limits for hydrocarbons, J HAZARD MATER, 163 (2009) 795-803. [10] C. Chen, T. Wang, H. Liaw, H. Chen, Nitrogen dilution effect on the flammability limits for hydrocarbons, J HAZARD MATER, 166 (2009) 880-890. [11] G. Shu, B. Long, H. Tian, H. Wei, X. Liang, Flame temperature theory-based model for evaluation of the flammable zones of hydrocarbon-air-CO2 mixtures, J HAZARD MATER, 294 (2015) 137-144. [12] J.L. Fuentes-Bargues, C. González-Gaya, M.C. González-Cruz, V. Cabrelles-Ramírez, Risk assessment of a compound feed process based on HAZOP analysis and linguistic terms, J LOSS PREVENT PROC, 44 (2016) 44-52. [13] S. Fu, X. Yan, D. Zhang, C. Li, E. Zio, Framework for the quantitative assessment of the risk of leakage from LNG-fueled vessels by an event tree-CFD, J LOSS PREVENT PROC, 43 (2016) 42-52. [14] Z.Y. Han, W.G. Weng, An integrated quantitative risk analysis method for natural gas pipeline network, J LOSS PREVENT PROC, 23 (2010) 428-436. [15] Stefana E, Marciano F, Alberti M. Qualitative risk assessment of a Dual Fuel (LNG-Diesel) system for heavy-duty trucks J LOSS PREVENT PROC, 39 (2015) 39-58. [16] Zhang Q, Wang Y, Lian Z. Explosion hazards of LPG-air mixtures in vented enclosure with obstacles, J HAZARD MATER, 334 (2017) 59-67. [17] Ehrhardt, L., et al., Evaluation of overpressure prediction models for air blast

above the triple point, J HAZARD MATER, 311 (2016) 176-185. [18] M. Pontiggia, M. Derudi, M. Alba, M. Scaioni, R. Rota, Hazardous gas releases in urban areas: Assessment of consequences through CFD modelling, J HAZARD MATER, 176 (2010) 589-596. [19] J. Xing, Z. Liu, P. Huang, C. Feng, Y. Zhou, D. Zhang, F. Wang, Experimental and numerical study of the dispersion of carbon dioxide plume, J HAZARD MATER, 256-257 (2013) 40-48. [20] X. Liu, A. Godbole, C. Lu, G. Michal, P. Venton, Source strength and dispersion of CO2 releases from high-pressure pipelines: CFD model using real gas equation of state, APPL ENERG, 126 (2014) 56-68. [21] R.S. Nagaosa, A new numerical formulation of gas leakage and spread into a residential space in terms of hazard analysis, J HAZARD MATER, 271 (2014) 266-274. [22] M. Siddiqui, S. Jayanti, T. Swaminathan, CFD analysis of dense gas dispersion in indoor environment for risk assessment and risk mitigation, J HAZARD MATER, 209-210 (2012) 177-185. [23] He G, Liang Y, Li Y, et al. A method for simulating the entire leaking process and calculating the liquid leakage volume of a damaged pressurized pipeline[J]. Journal of Hazardous Materials, 2017, 332:19-32 [24] Jianwen Z, Xinxin Y, Yanan X, et al. Numerical investigation on three-dimensional dispersion and conversion behaviors of silicon tetrachloride release in the atmosphere.[J]. Journal of Hazardous Materials, 2015, 288:1-16. [25] S. DOROFEEV, Evaluation of safety distances related to unconfined hydrogen explosions, INT J HYDROGEN ENERG, 32 (2007) 2118-2124. [26] H. Montiel, J. A. Vılchez, J. Casal, J. Arnaldos, Mathematical modelling of accidental gas releases, J HAZARD MATER, 59(1998) 211-233 [27] Q. Liu, Y. Liu, H. Zhang, R. Yang, Z. Li, Effects of Ventilation Intensity on Leakage of Flammable Refrigerants of an Air- conditioning System. J Basic Science and Engineering,21 (2013) 1157-1166. [28] Y.D. Jo, B.J. Ahn, A simple model for the release rate of hazardous gas from a hole on high-pressure pipelines, J HAZARD MATER, 97 (2003) 31-46. [29] Y. Guo, L. He, D. Wang, S. Liu, Numerical investigation of surface conduit parallel gas pipeline explosive based on the TNT equivalent weight method, J LOSS PREVENT PROC, 44 (2016) 360-368. [30] Liu, B., et al., Computational fluid dynamics simulation of carbon dioxide dispersion in a complex environment. J LOSS PREVENT PROC, 2016. 40: p. 419-432. [31] Scargiali F, Grisafi F, Busciglio A, et al. Modeling and simulation of dense cloud dispersion in urban areas by means of computational fluid dynamics. J HAZARD MATER, 2011, 197(24):285-293. [32] J. Zhang, D. Lei, W. Feng, An approach for estimating toxic release of H2S-containing natural gas, J HAZARD MATER. 264 (2014) 350-362. [33]. Luketa-Hanlin, A., R.P. Koopman and D.L. Ermak, On the application of computational fluid dynamics codes for liquefied natural gas dispersion. J HAZARD

MATER, 2007. 140(3): 504-517. [34] Jianwen, Z., et al., Numerical investigation on three-dimensional dispersion and conversion behaviors of silicon tetrachloride release in the atmosphere. J HAZARD MATER, 2015. 288: 1-16. [35] Zhang J, Lei D, Feng W. An approach for estimating toxic releases of H 2 S-containing natural gas. J HAZARD MATER, 2014, 264(2):350. [36] Z. Wang, Y. Hu, J. Jiang, Numerical investigation of leaking and dispersion of carbon dioxide indoor under ventilation condition, ENERG BUILDINGS, 66 (2013) 461-466. [37] B.Y. Zhang, H.H. Li, W. Wang, Numerical study of dynamic response and failure analysis of spherical storage tanks under external blast loading, J LOSS PREVENT PROC, 34 (2015) 209-217. [38] H. Si, H. Ji, X. Zeng, Quantitative risk assessment model of hazardous chemicals leakage and application, SAFETY SCI, 50 (2012) 1452-1461. [39] Y. Du, F Zhang, A Zhang, et al, Consequence assessment of explosion in pipes using coupled FEM-SPH method, J LOSS PREVENT PROC, 43 (2016) 549-558. [40] F.P. Lees, Lees’ Loss Prevention in the Process Industries, second ed, Butterworth-Heinemann, Oxford, UK, 1996. [41] Pérez, J.F.S., et al., New Probit equations for the calculation of thermal effects on humans. PROCESS SAF ENVIRON, 88 (2010) 109-113. [42] Finney DJ. Probit analysis. Cambridge University Press, 1974.

Fig.1.Schematic diagram of the ORC Quantitative Risk-Assessment System.

Fig.2. QR-AS modeling flow diagram

Fig.3. Comparison of simulated results with reference [36]

Fig.4.Photograph of ORC system in the lab

Fig.5.Schematic diagram of the physical mode

Fig.6 .Grids in the computational domain

10s

30s

50s

100s Fig.7. Concentration distribution of propane on the plane of y=1.23

50s

150s Fig.8. Concentration distribution of propane on the plane of y=2.05m

Fig.9 Propane concentration on the plane of x=3.0m at different heights

Fig.10. Propane concentration on the line y=1.1m, z=2m at different time.

10s

50s

100s

150s Fig.11. Explosive zone at different time when the wind speed is 0.5m/s

Fig.12.Dangerous radius at different time

Fig.13.Relationship between overpressure and distances from explosion site at different time

(a) Injury probability with respect to distance from the explosion center.

(b).Relationship between overpressure and distances from explosion site Fig.14 Explosive damage degree to individuals and buildings in Scenario B

Plane y=2.05m

Plane y=2.55m Fig.15.Concentration distribution of propane on plane y=2.05m and y=2.55m when the wind speed is 2m/s

Fig.16. Explosive zone under different wind speed

Fig.17. Explosive zone under different ratio when the wind speed is 0.5m/s

Table 1 Boundary conditions used in the simulation Parameters Air inlet Gas inlet Out Wall Air temperature Leakage gas temperature

Setup Velocity inlet, 1m/s Velocity inlet, 245m/s Pressure outlet No-slip, adiabatic 293K 307K

Table 2 Overpressure intervals relating to expected damage [38] Overpressure (kPa) Expected damages 1.0-1.5 Window glass cracks 3.5-7.6 Minor damage in some buildings 7.6-12.4 Metal panels deform 12.4-20 Concrete walls damage >35 Wooden construction building demolition 27.5-48 Major damage on steel construction objects 40-60 Heavy damage on reinforced concrete buildings 70-80 Probable demolition of most buildings

Table 3 Transformation of Probit to Percentage [42] % 0 10 20 30 40 50 60 70 80 90 — 99

0 — 3.72 4.16 4.48 4.75 5 5.25 5.52 5.84 6.28 0 7.33

1 2.67 3.77 4.19 4.5 4.77 5.03 5.28 5.55 5.88 6.34 0.1 7.37

2 2.95 3.82 4.23 4.53 4.8 5.05 5.31 5.58 5.92 6.41 0.2 7.41

3 3.12 3.87 4.26 4.56 4.82 5.08 5.33 5.61 5.95 6.48 0.3 7.46

4 3.25 3.92 4.29 4.59 4.85 5.1 5.36 5.64 5.99 6.55 0.4 7.51

5 3.36 3.96 4.33 4.61 4.87 5.13 5.39 5.67 6.04 6.64 0.5 7.58

6 3.45 4.01 4.36 4.64 4.9 5.15 5.41 5.71 6.08 6.75 0.6 7.65

7 3.52 4.05 4.39 4.67 4.92 5.18 5.44 5.74 6.13 6.88 0.7 7.75

8 3.59 4.08 4.42 4.69 4.95 5.2 5.47 5.77 6.18 7.05 0.8 7.88

9 3.66 4.12 4.45 4.72 4.97 5.23 5.5 5.81 6.23 7.33 0.9 8.09

t (s) 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200

Appendix Evaluation results of the case in Section 3 R0.5 Rd,0.5 Rd,0.01 △p-top △p-bottom △p-left △p-right △p-front △p-back (m) 0.44 0.57 0.66 0.73 0.80 0.85 0.90 0.95 0.99 1.03 1.07 1.10 1.14 1.17 1.20 1.23 1.25 1.28 1.31 1.33

(m) 1.75 2.21 2.53 2.79 3.00 3.19 3.36 3.51 3.65 3.78 3.90 4.02 4.13 4.23 4.33 4.42 4.51 4.60 4.68 4.76

(m) 3.15 3.97 4.55 5.01 5.39 5.73 6.03 6.31 6.56 6.79 7.01 7.22 7.41 7.60 7.78 7.94 8.11 8.26 8.41 8.56

(kPa) 39.01 58.66 75.22 90.15 104.02 117.12 129.64 141.68 153.34 164.67 175.71 186.51 197.09 207.47 217.68 227.73 237.63 247.40 257.05 266.59

(kPa) 628.71 966.91 1211.24 1436.44 1650.12 1855.97 2056.11 2251.87 2444.14 2633.56 2820.60 3005.62 3188.90 3370.66 3551.08 3730.32 3908.49 4085.70 4262.05 4437.60

(kPa) 50.27 76.62 99.06 119.42 138.42 156.43 173.69 190.35 206.50 222.23 237.59 252.64 267.40 281.91 296.19 310.27 324.16 337.87 351.43 364.84

(kPa) 22.75 33.21 41.84 49.50 56.55 63.15 69.41 75.40 81.16 86.74 92.16 97.43 102.58 107.62 112.56 117.41 122.18 126.87 131.49 136.05

(kPa) 94.40 148.70 196.01 239.56 280.65 319.94 357.84 394.64 430.50 465.59 500.00 533.82 567.12 599.95 632.35 664.37 696.04 727.39 758.44 787.42

(kPa) 33.95 50.67 64.68 77.26 88.92 99.92 110.40 120.47 130.21 139.65 148.86 157.85 166.65 175.28 183.76 192.10 200.32 208.42 216.42 224.31

R0.5 is the radius of death zone Rd, 0.5, is the radius of serious injury zone Rd, 0.01 is the radius of slight injury zone respectively △p-top is the overpressure at the top wall, similarly to the meaning of △p-bottom, △p-left,△p-right,△p-front,△p-back

37