Numerical simulation for evaporation of chemical agent droplet

Numerical simulation for evaporation of chemical agent droplet V International Conference on Computational Methods for Coupled Problems in Science and Engineering COUPLED PROBLEMS 2013 S. Idelsohn, M. Papadrakakis and B. Schrefler (Eds) NUMERICAL SIMULATION FOR EVAPORATION OF CHEMICAL AGENT DROPLET H. HAYAKAWA*, M. YAMAMOTO†, M. KAWAI‡ AND K. SAKAMOTO§ * Graduate School of Mechanical Engineering Kudan-Campus, Tokyo University of Science 1-14-6 Kudankita, Chiyoda-ku, Tokyo 102-0073, Japan e-mail: [email protected] † Department of Mechanical Engineering Kudan-Campus, Tokyo University of Science 1-14-6 Kudankita, Chiyoda-ku, Tokyo 102-0073, Japan email: [email protected] ‡ Numerical Engineering Dept, IHI 1, Shin-Nakahara-Cho, Isogo-ku, Yokohama 235-8501, Japan email: [email protected] § Ordnance Engineering Department Defense Systems Division Aero-Engine & Space Operations, IHI 229, Tonogaya, Mizuho-machi, Nishitama-gun, Tokyo 190-1297, Japan email: [email protected] Key words: Computational Fluid Dynamics, Chemical agent droplet, Evaporation. Abstract. Recently, terrorism is one of most daunted dangers in the world. Various terrorisms can come about. Biological and chemical terrors have particularly high mortality rate, and tend to wreak the secondary disaster. For these reasons, research and development of measures for decontaminations are emergent and of key issue on a world scale. Sulfur mustard (HD) is one of chemical agent used to commit terrorism. HD gas has severe toxicity and long-period damage because of persistent agent. In addition, as HD gas achieves the toxicity even if it only contacts on a skin, it is important that HD-contaminated objects are decontaminated quickly. In the present study, HD gas emission from a surface of HD droplet is numerically researched. HD gas emission is usually investigated with experiments, but the experiment is so dangerous because of the toxicity. Therefore, in this study, CFD (computational fluid dynamics) is used to reproduce the emission. We propose a volatilization model which can estimate the volatilization of a chemical agent with the saturated vapor pressure of chemical agent and physical quantities of the air. We apply our three dimensional volatilization prediction code to a HD droplet placed on a wind tunnel wall. Comparing the numerical results with the experimental data, it is confirmed that the numerical data are in good agreement with the experimental data. In addition, volatilization volume is coupled with droplet geometry to simulate the temporal change of droplet shape and residual mass. It is reasonably predicted that the droplet volume gradually decreases as volatilization proceeds 1191 H. Hayakawa, M. Yamamoto, M. Kawai and K. Sakamoto. 2 until the droplet disappears. 1 INTRODUCTION There have been various terror attacks around the world since 9.11 in 2001 [1]. Terrorists hold extremely wide-ranging means of terrors such as bombs and toxic gases. Above all, nuclear, biological and chemical weapons have particularly high mortality rate so that they are recognized as being the most dangerous in all other ways of terror. Among these weapons, biological and chemical weapons are relatively-easy to be produced so that they are at high risk of being used in terrorisms. In the past, there have been some cases by the use of chemical weapons. In 1978, a refugee from Bulgaria was killed by the pellet containing ricin. In Japan, sarin was used to commit indiscriminate murders in Matsumoto city in 1994 and in the Tokyo subway in1995. Sarin is extremely dangerous nerve agent so that 5,510 people were injured and 12 people were killed in the Tokyo subway attack. In terror attacks, many civilians can be victimized. Furthermore, it is hard to predict when and where it happens so that research and development against unanticipated terror attacks are emergent. However, research and development of chemical weapons are of great difficulty since the experiments are highly risky and are legally-prohibited in many cases. Therefore, numerical simulation is very useful to research it. HD (Sulfur mustard) is a type of blister agent. It affects eyes, skins and respiratory organs. In addition, HD is a persistent agent so that the damage continues over a long period. In the present study, computational fluid dynamics is used to simulate the evaporation of HD droplets. Fick’s law is used to calculate the droplet’s volatilization rate and volatilization volume is coupled with droplet geometry to simulate the temporal change of droplet shape and residual mass. In addition, the versatility of our simulation code is confirmed by applying it to various conditions of temperatures and wind speeds. It is reasonably predicted that the droplet volume gradually decreases as volatilization proceeds until the droplet disappears. 2 NUMERICAL PROCEDURE In this study, volatilization of chemical agent droplet placed in a steady flow is investigated. The flow and concentration fields with chemical agent vapor in a wind tunnel shown in Figure 1 is computed. Since the flow is assumed to be multiphase flow consisting of air and chemical agent, we use a homogenous fluid model. As the treatment of multiphase flow, we apply the two-way coupling method, that is, the concentration and the flow strongly affects each other. The numerical procedures for the flow field, the concentration field, the volatilization model and the coupling method of volatilization volume and droplet geometry will be described below. 2.1 Flow field It is assumed that the flow field is three-dimensional and incompressible. We use the conservation equations of mass and momentum expressed by Eqs.(1) and (2) as the governing equations. For pressure-velocity coupling, the Low-Mach-Number approximated MAC method (marker and cell method) is used. Since we treat two-phase flow consisting of air and chemical agent, the mixture density calculated by Eqs.(3) and (4) at each grid points are considered in the flow computations. No-slip condition between air and chemical agent is 1192 H. Hayakawa, M. Yamamoto, M. Kawai and K. Sakamoto. 3 imposed to calculate the mixture gas velocity, mix u  .   0   mixmix mix u t    (1)   gupuu t u mixsmixairmixmixmix mixmix   )( 2      (2) HDHDairairmix   (3) 1 HDair  (4) In the above equations,  is the volume fraction,  is the density, u is the flow velocity, t is the time, p is the pressure,  is the viscosity coefficient, g is the gravity acceleration, the subscript air is for the air, HD is for the chemical agent (HD), mix is for the mixture gas and s is for the average value of overall computed domain. Since the volume fraction of chemical agent is low enough, typically 0.03 percent, the viscosity of air is used in Eq.(2). The gravity term is calculated by the difference between the average density and the local density. A finite difference method is applied to discretize the governing equations. The Kawamura-Kuwahara [2] third-order upwind scheme is used for the convection terms, the Euler explicit scheme for the time marching terms, and the second order central difference scheme for the other terms. 2.2 Concentration field For the concentration field, the governing equation is given as HDHDHDHD HD CDCu t C 2    (5) where C is mol concentration and D is the diffusion coefficient. D is calculated from the estimate equation by Fuller et al. [3]. The numerical difference schemes are employed in similar to the flow field. 2.3 Volatilization model In the volatilization volume calculation, at first, the diffusion flux is calculated by Fick’s law expressed by Eq.(6). The volatilization volume is calculated by diffusion flux integrated by time and surface area. It is expressed by Eq.(7). dx dC DJ HD  (6) dxdtJm   (7) In these equations, J is the diffusion flux, x is the coordinate normal to the droplet surface, and m is the volatilization volume. In the volatilization model, whether chemical agent volatilize or not is judged in the following way. At first, the partial pressure of chemical agent is calculated by the state 1193 H. Hayakawa, M. Yamamoto, M. Kawai and K. Sakamoto. 4 equation. Next, if the partial pressure is lower than the saturated vapor pressure, the chemical agent is volatilized. On the other hand, if the partial pressure is higher than the saturated vapor pressure, the chemical agent is not volatilized. This judgement is made for each grid points on the droplet surface. Finally, volatilization volume is calculated. 2.4 Coupling method In this study, temporal change of droplet volume with the progress of evaporation is reproduced by coupling of volatilization volume and droplet geometry. In response to the deformation of droplet geometry, the flow and concentration fields change. Furthermore, the change in concentration field affects the volatilization rate. Thus, in this evaporation system, concentration field, volatilization rate and droplet geometry are affected among others. HD evaporation is long-term phenomena due to its persistent property. On the other hand, the time scale of flow and concentration fields is short. Therefore, the time scales of the two are remarkably different. For this reason, it is difficult that long-term time integration is performed in accordance with the short time scale of flow and concentration fields. Therefore, in this study, droplet geometry is estimated by using the steady solution of flow and concentration fields as long as the change of droplet geometry is not large enough to remarkably affects the flow and concentration. In other words, their mutual effects are considered by the iteration of the cycle consisting of “the flow and concentration fields” and “the change of droplet geometry by evaporation”. This method enables us to simulate the evaporation phenomena without enormous computational cost. Figure 1: Computational domain, grid and examples of flow and concentration field Chemical agent droplet z x z 25[mm] 50[mm] 120[mm] 60[mm] Outlet Chemical agent droplet Inlet 6[mm] y x (a) Computational domain and grid (b) Velocity vector (c) Distribution of HD concentration 0.00 0.12 [m/s] 0.0 12.7 [mmol/m3] 1194 H. Hayakawa, M. Yamamoto, M. Kawai and K. Sakamoto. 5 3 COMPUTATIONAL CONDITION Computations are performed for the flow in the half geometry of 5cm × 5cm wind tunnel as shown in Figure 1. We employed this geometry so as to match the condition of the experiment performed by Danberg [4]. To reduce the computational costs, only the half area of the wind tunnel is computed. The chemical agent is considered to be HD, and the droplet is placed on the wind tunnel wall as shown in Figure 1. The number of grid points is 240 × 90 × 38, clustered near the droplet, walls and inlet boundaries, where the gradients of physical quantities are expected to be steep. On the wall boundary, no-slip condition is imposed. On the upper surface, slip boundary condition is applied. On the inlet boundary, velocity is fixed and others are extrapolated from the computational domain. In order to investigate the effect of temperature and flow speed, the temperature is varied within the range of 15 degrees to 50 degrees, and the inflow velocity is changed from 0.12 to 2.57 m/s. Validity of this numerical code was carried out with the combination of three nondimensional numbers defined by Eqs.(8) to (10).   Dcc LN Sh w    (8)  uL Re (9) D Sc   (10) In these definitions, N is the volatilization rate, L is the characteristic length,  cc w is the concentration difference of the evaporating vapor between the surface and the edge of the evaporating vapor dispersion layer, u is the characteristic velocity and  is the kinetic viscosity. Sh, Re and Sc are Sherwood number, Reynolds number and Schmidt number, respectively. The computational conditions are listed in Table 1. We set Case7 as the reference case. 4 NUMERICAL RESULT AND DISCUSSION 4.1 Validation of developed code The temporal change of residual and volatilization rate of the reference case (i.e. Case7) is exhibited in Figure 2. It indicates that volatilization rate become constant after about one second. 1195 H. Hayakawa, M. Yamamoto, M. Kawai and K. Sakamoto. 6 Table 1: Computational conditions Case Temperature [degrees] Mainstream velocity [m/s] Diffusivity [m 2 /s] Saturated vapor pressure [N/m 2 ] 1 15 0.12 6.56×10-6 5.01 2 15 0.26 6.56×10-6 5.01 3 15 0.64 6.56×10-6 5.01 4 15 1.28 6.56×10-6 5.01 5 15 1.92 6.56×10-6 5.01 6 15 2.57 6.56×10-6 5.01 7 35 0.12 7.38×10-6 32.4 8 35 0.26 7.38×10-6 32.4 9 35 0.64 7.38×10-6 32.4 10 35 1.28 7.38×10-6 32.4 11 35 1.92 7.38×10-6 32.4 12 35 2.57 7.38×10-6 32.4 13 50 0.12 8.02×10-6 105.4 14 50 0.26 8.02×10-6 105.4 15 50 0.64 8.02×10-6 105.4 16 50 1.28 8.02×10-6 105.4 17 50 1.92 8.02×10-6 105.4 18 50 2.57 8.02×10-6 105.4 Figure 2: Temporal change of residual and volatilization rate of reference case 0 2 4 10 −8 10 −6 10 −4 10 −2 0 2 4 0 5 10 Time [sec] Residual Volatilization Rate [mg/m 2 s] Residual Volatilization Rate 1196 H. Hayakawa, M. Yamamoto, M. Kawai and K. Sakamoto. 7 Figure 3 shows the relation among Sh, Re and Sc. The line indicates the experiment and the symbols do the present computation. 3 1 3 2 Re ScGSh  (11) Based on the experimental data, the mean gradient G is 0.828. Our computation gives G is 0.836. Hence, it is confirmed that our numerical code is reasonably validated. Figure 4(a) shows the steady distribution of HD concentration in the reference case (Case7). The molecular mass of HD is much heavier than that of air so that HD settles near the bottom of the tunnel soon after the volatilization. Furthermore, HD concentration near the droplet surface increases remarkably. This trend affects the volatilization rate of the downstream. Figure 4(b) exhibits the steady distribution of volatilization rate over the droplet surface in the reference case. The volatilization rate of the downstream side is relatively lower than that of the upstream side. This is because that the concentration of HD near the downstream side is increased by volatilized HD vapor at the upstream, and thus the volatilization rate of the downstream side decreases. Figure 3: Relation among non-dimensional parameters Figure 4: Numerical Results of the reference case (Case 7) 010 20 30 0 10 20 30 Sh Re2/3 Sc1/3 Experiment (G=0.828) Numerical Results ・15 degrees ・35 degrees ・50 degrees z x 0.00 6.33 12.7 [mmol/m3] (a) Distribution of HD concentration (b) Distribution of Volatilization rate x y 29.8 0.00 59.6 [mg/m2s] 1197 H. Hayakawa, M. Yamamoto, M. Kawai and K. Sakamoto. 8 4.2 Evaporation simulation Figure 5 shows the distributions of HD concentration and volatilization rate. In these figures, the droplet becomes progressively smaller along with time. Computations were performed until the droplet volume becomes zero. The computational result suggests that it takes over 13 hours to dry off 9 micro liters droplet. It comes from the strong stability of HD. Figure 6 plots the temporal change of droplet volume. It indicates that the change of droplet volume is non-linear. The smaller droplet volume is, the smaller the surface area becomes. As a result, volatilization rate for the whole surface area decreases. 4.3 Effect of temperature and main stream speed Figure 7(a) shows the temporal changes of droplet volume for the cases that the inflow velocity is fixed at 0.12 m/s and the temperature is varied from 15 to 50 degrees. It is clear that evaporation speed increases with increasing the temperature. This trend is caused by the dependence of saturated vapor pressure on temperature. The value of saturated vapor pressure affects the suppression of volatilization. Figure 7(b) shows the temporal changes of droplet volume for the cases that temperature is fixed at 35 degrees and the inflow velocity is varied from 0.12 to 0.64 m/s. From this graph, evaporation speed increases with increasing the inflow velocity. This is because that HD vapor is effectively blown at high speed and it results in the low concentration of HD near the droplet surface. Through these computations, we confirmed that our numerical code can simulate the evaporation of chemical agent droplet in line with the principle of evaporation within the applied range of temperature and wind speed. Figure 5: Numerical Results of evaporation simulation z x 12.7 x y 2.0 [sec] 1,2000 [sec] 2,4000 [sec] 3,6000 [sec] 4,9500 [sec] [mmol/m3] 6.33 0.00 59.6 [mg/m2s] 29.8 0.00 (a) Distribution of HD concentration (b) Distribution of volatilization rate 1198 H. Hayakawa, M. Yamamoto, M. Kawai and K. Sakamoto. 9 Figure 6: Temporal change of droplet volume Figure 7: Temporal change of droplet volume under different conditions U inflow = 0.12[m/s] (constant) 15 degrees 35 degrees 50 degrees 35 degrees (constant) Uinflow = 0.12 [m/s] Uinflow = 0.26 [m/s] Uinflow = 0.64 [m/s] (b) Effect of inflow velocity (a) Effect of temperature 020000 40000 0 5 Droplet volume [μL] Time [sec] 020000 40000 0 5 Droplet volume [μL] Time [sec] 020000 40000 0 5 Droplet volume [μL] Time [sec] 1199