Hongbo Tan,Boshi Shao,Na Wen

Department of Refrigeration and Cryogenic Engineering, School of Energy and Power Engineering, Xi ’an Jiaotong University, Xi’an 710049, China

Keywords:Numerical simulation Greenhouse gas Condensation De-sublimation Vapor deep removal Non-condensable gas

ABSTRACT Nowadays,the limits on greenhouse gas emissions are becoming increasingly stringent.In present research,a two-dimensional numerical model was established to simulate the deep removal of 1,1,1,2-tetrafluoroethane (R134a) from the non-condensable gas (NCG) mixture by cryogenic condensation and de-sublimation.The wall condensation method was compiled into the Fluent software to calculate the condensation of R134a from the gas mixture.Besides,the saturated thermodynamic properties of R134a under its triple point were extrapolated by the equation of state.The simulation of the steam condensation with NCG was conducted to verify the validity of the model,the results matched well with the experimental data.Subsequently,the condensation characteristics of R134a with NCG and the thermodynamic parameters affecting condensation were studied.The results show that the section with relatively higher removal efficiency is usually near the inlet.The cold wall temperature has a great influence on the R134a removal performance,e.g.,a 15 K reduction of the wall temperature brings a reduction in the outlet R134a molar fraction by 85.43%.The effect of changing mass flow rate on R134a removal is mainly reflected at the outlet,where an increase in mass flow rate of 12.6%can aggravate the outlet molar fraction to 210.3%of the original.The research can provide a valuable reference for the simulation of the deep removal of various low-concentration gas using condensation and de-sublimation methods.

1.Introduction

The HFC-134a (1,1,1,2-tetrafluoroethane,R134a),is still one of the most widely used volatile organic compounds (VOCs) in the field of the domestic refrigerator,automobile air conditioning and industrial chillers [1].The restrictions on the emission of R134a have become increasingly stringent since the global warming potential (GWP) of R134a is up to 1301.Moreover,the direct emissions resulting in waste of resource gas simultaneously.Thus,the removal and recovery of R134a from the air has become a particularly significant engineering topic.

Many approaches,such as the condensation and de-sublimation method(CDM),can be employed to process the exhausted gas mixture.The CDM is highly universal and environmentally friendly without chemical reaction and desorption.Regarding a typical R134a/nitrogen binary mixture gas,the R134a condenses first while the high concentration of nitrogen plays as noncondensable gas (NCG).Thus,the study of the condensation and de-sublimation characteristics of R134a in the presence of NCG is very important for the condensation and de-sublimation removal and recovery of VOCs.

Many researchers had analyzed the condensation of vapor in the presence of NCG by theoretical and experimental methods.The two main methods for theoretical calculation are the diffusion layer model [2,3] and the boundary layer model [4-6].However,the application of the theoretical model is limited by the difficulty of obtaining the exact solution of some complex mathematical equations [7].

The computational fluid dynamics(CFD)simulation has become an effective method to solve the vapor/NCG condensation process.Although only pure vapor condensation models are built in the commercial software such as Fluent and CFX,vapor/NCG gas mixture condensation calculation can be implemented by specifying the source item terms in the control equations by the userdefined function (UDF).The Lee model is a common and classical method to simulate the vapor/NCG condensation.Yinet al.[8]conducted the numerical simulation of the steam laminar film condensation with NCG in a horizontal mini-tube based on the Lee model and volume of fluid(VOF)method.They found that the heat transfer coefficient of NCG has a great influence on the heat flux.Gao[9]obtained the distribution of temperature and velocity fields in the vertical tube and the variation of heat transfer coefficient with the NCG concentration by CFD simulation based on the Lee model and VOF method.Genget al.[10] carried out a numerical study by Lee model and concluded that the heat transfer suppression of NCG on steam condensation is ameliorated under high pressure.

Although the Lee model is simple in form and robust in convergence,it contains an empirical coefficient determined by the experiment,which is difficult for wide applications.Besides,the multi-phase model must be combined with the Lee model,which increases the complexity of the simulation.In 2013,Dehbiet al.[11] first integrated another model without empirical coefficients named wall condensation model (WCM) into the calculation of the CFD source terms for the vapor/NCG condensation simulation.They ignored the multi-phase model and compared several sets of CFD simulation using WCM with the experimental data.The results showed well consistency.Based on the above characteristics,the WCM was chosen to predict the condensation and desublimation of R134a/nitrogen mixture in this paper.Punethaet al.[12] simulated the steam condensation in the air and hydrogen mixture in the nuclear power plant containments using WCM.The simulation worked quite well compared with the experiments.Zhanget al.[13] numerically researched the condensation process of propane/methane mixture on the vertical plate by WCM and found that the heat transfer coefficient decreased by more than 90% at the NCG methane with the high mole fraction of 80%-95%compared with the pure vapor.Leeet al.[14] conjugated the heat transfer simulations with WCM in the presence of NCG whose results showed high agreement with the CONAN experimental database.Liuet al.[15] revealed the stratification characteristics of helium-air-steam mixture gas under steam condensation with wide range of steam concentration by the simulation using WCM.The details of the simulations mentioned above were summarized in Table 1.

Table 1 Summary of simulations studies of the vapor/NCG condensation

From the above researches,it can be clearly found that few studies have focused on the deep removal of R134a under low concentration by CDM.For the deep removal of R134a from the gas mixture (the molar fraction of R134a in the mixture less than 0.01%),the temperature of the mixture gas should be decreased low enough for the CDM.Consequently,this study established a 2-D numerical model for the deep removal of R134a with low inlet molar fraction (2%) from nitrogen based on the CDM.The WCM was adopted and compiled into ANSYS Fluent by UDF to simulate the phase change process.The multiphase flow process was ignored because of the low concentration of R134a.The saturated thermodynamic properties of R134a above the triple point were obtain from the National Institute of Standards and Technology(NIST) database,while those below the triple point were extrapolated and estimated by the Peng-Robinson equation [16] to simulate the R134a de-sublimation process.In addition,the effect of different parameters on the removal process was discussed.The research can provide a valuable reference for the simulation of the deep removal of various low-concentration VOCs using CDM.

2.CFD Modeling

2.1.Geometry

The heat transfer channel in the parallel plate heat exchanger is composed of multiple groups of hot and cold fluid heat transfer units.Fig.1 shows a 2-D geometric model of a heat transfer unit of the parallel plate heat exchanger used in this study.The temperature reduction of the mixture gas was achieved by countercurrent heat exchange with the coolant flows in the adjacent channel of a parallel plate heat exchanger.An axisymmetric region which was enclosed by the red dotted line in Fig.1 was selected as the calculation section for the simulation of the deep removal of the low-concentration R134a from R134a/nitrogen mixture gas to simplify the calculations.Therefore,the cooling effect of the coolant was reflected in the cold wall temperature of the channel in the simulation.It is assumed that the condensation produces dispersed droplets rather than a continuous liquid film due to the low initial molar fraction of R134a and the slow condensation rate determines the difficulty of fully wetting the wall.The single channel size was estimated as 20 m × 1 m × 0.01 m(L×W×H) and the thickness of the partition between the channels was chosen as 1 mm.

Fig.1.The sketch of the geometry model and the fluid domain.

The calculation conditions of the simulation can be seen in Table 2.The outlet temperature of the coolant was assumed and the outlet temperature of the mixture gas was calculated by the law of conservation of energy.

Table 2 Calculation conditions of the mixture gas and the coolant

2.2.Numerical methods

2.2.1.Governing equations

The continuity equation,momentum conservation equation,energy conservation equation and species conservation equation involved in the study for the two-dimension steady state model are as follow.

The continuity equation can be written as:

where ρ is the density of the mixture gas,kg·m-3;uand v are velocity in the abscissa and ordinate respectively,m·s-1;Smis the mass source term,kg·m-3·s-1.

The continuity equation can be written as:

where η is the dynamic viscosity of the mixture gas,Pa·s;SuandSvare the momentum source term in the abscissa and ordinate respectively,N·m-3.The pressure gradient term ?pin the equations is taken as zero since the pressure loss is ignored during the research.

The energy conservation equation can be written as:

whereTis the temperature of the mixture gas,K;λ is the thermal conductivity of the mixture gas,W·m-1·K-1;cpis the specific heat capacity of the mixture gas,J·kg-1·K-1;Seis the energy source term,W·m-3.

The species conservation equation can be written as:

where ωvis the mass fraction of the vapor;Dvis the diffusion coefficient of the vapor,m2·s-1;Ssis the species source term,kg·m-3·s-1.

2.2.2.Wall condensation method

The WCM was used to simulate the mass and energy transfer behaviors at the cold wall of the mixture gas during the condensation of vapor/NCG.The method was implemented by compiling the source terms involved in the governing equations mentioned above into the ANSYS Fluent through the UDF.

Since the liquid and solid film was neglected in this study,the cold wall surface was regarded as the equivalent phase interface.The mass diffusion flux with a suction correction factor θcat the interface when the condensation and de-sublimation occurs can be expressed as [11,17]:

where ωv,iis the saturated mass fraction of the vapor corresponding to the interface temperature;nis the normal direction of the cold wall;the differential termrepresents the mass fraction gradient within the first layer of the grid immediately above the wall;θcis the suction correction factor,which is defined as:

whereBis the Bird suction parameter [18] can be expressed as:

where ωv,bis the mass fraction of the vapor in the bulk flow.Based on the value of θcis 1.0071 in the present calculation condition,the suction effect represented by θccan be ignored reasonably.

The condensation or de-sublimation rate of the vapor should be equal to the mass diffusion flux through the interface when the phase change process keep equilibrium.Besides,the gradientis related to the difference between the vapor mass fraction in the first layer of the grid adjacent to the interface (ωv) and the saturated vapor mass fraction at the phase interface (ωv,i) when the height of the first grid (l) is low enough,which means that:

where ˙mpcindicates the phase change rate,kg·m-3·s-1;ωvis the vapor mass fraction in the first layer of the mesh adjacent to the cold wall;lis the vertical distance from the center of the first layer of the mesh to the cold wall,m.The difference between ωvand ωv,irepresents the driving force for mass transfer.

The cold wall is impermeable to the mass diffusion flux based on the no-slip boundary condition.As a result,the vapor condensation/sublimation is implemented with the source terms in the governing equations in the first cells adjacent to the cold wall.Hence the source terms in Eqs.(1)-(5) can be expressed as a universal form using the condensation or de-sublimation rate mentioned above:

whereAfaceis the projected area of the grid on the cold wall,m2;Vcellis the volume of the grid adjacent to the cold wall,m3;? is a generic variable whose values are listed in the Table 3.

Table 3 The value of ? for different source term calculations

Table 4 Values of the parameters involved in the Eq.(11)

Table 5 Values of the parameters involved in the Eq.(12)

The Eq.(10) was compiled into ANSYS Fluent by the UDF to implement the function of the phase change.In summary,the vapor phase change will occur and Eq.(10) will take a non-zero value when the temperature in the first cells adjacent to the cold wall (Tcell) is lower than or equal to the saturation temperature under the partial pressure of the vapor (Tsat),while conversely,the source term will keep zero.The computational procedure of the UDF is shown in Fig.2:

Fig.2.The computational procedure of the WCM by the UDF.

2.3.Parameter determination

2.3.1.Saturation temperature and pressure

The determination of saturation temperature (Tsat) and saturation pressure(psat)under different pressures or temperature is crucial to the condensation and de-sublimation simulation.The segmental fitting equation is usually used to compile theT-prelationship into the UDF.Nevertheless,the values are usually different at the junction points of the segmented fitting curves,resulting in the saltation in the results.In addition,the exactT-prelationship of R134a when the de-sublimation occurs is vague because of the lack of relevant experimental data.

To solve these problems,the difference of the two fitting equations at the junction points was added to one equation as a constant term.The equation added the constant term is generally larger in value to reduce the error.Besides,the different slopes at the junction point will also cause the saltation.This effect will be attenuated by increasing the number of the segment equations.And then,in order to simulate the de-sublimation of R134a,it is necessary to extrapolate theT-prelationship below its triple point.Consequently,the Peng-Robinson (P-R) equation was chosen to implement the extrapolation based on its high accuracy in the description of the gas-solid phase change of the sulfur in the sour gas [19].

In present study,the thermal properties of R134a above the triple point were obtained from the database of the NIST,while the extrapolation by the P-R equation was conducted by the fluid package in Aspen HYSYS.The universal form of the fitting equations of theT-prelationship of R134a can be written as:

The values of the parameters involved in the Eqs.(11) and (12)are shown in the Tables 4 and 5,respectively.

2.3.2.Condensation or de-sublimation rate

The condensation or de-sublimation rate among a section Δxin the calculation domain is an important parameter,which will be determined through the post-processing process.A section in the computation domain is represented by the point O and the rate can be determined through the difference between the molar concentrations of the R134a at the two adjacent points of the point O shown in Fig.3.

Fig.3.The sketch of the section Δx.

The condensation or de-sublimation rate among the section Δxcan be expressed as:

whereMR134ais the relative molecular mass of the R134a,g·mol-1;uis the flow velocity of the mixture gas,m·s-1;WWandWEare the molar concentration obtained easily from Fluent at the point W and E respectively,mol·m-3;wand δ are the wide and thickness of the channel,whose values in present study are 1.0 m and 0.01 m respectively.

2.3.3.Thermal properties

The diffusion coefficient of the R134a/nitrogen mixture gas at different temperature or pressure is estimated by a semi-rational formula given by Gilliland [20].

where A and B are the components of the binary mixture gas,which refer to R134a and nitrogen respectively in present research;Vis the liquid molar volume at normal boiling point,cm3·mol-1;Mis the relative molecular mass,g·mol-1.For the present simulation,the parameters are shown in the Table 6.

Table 6 The relevant parameters in Eq.(14)

The latent heat of condensation or de-sublimation is also a parameter that needs to be specified to calculate the energy source term.The latent heat of R134a condensation changes linearly with the saturation temperature approximately,which can be fitted by the data from NIST as:

Besides,the latent heat of the R134a de-sublimation can also be easily extrapolated using Eq.(15).The error from the estimation of the de-sublimation latent heat is limited since the molar fraction of R134a and the phase change rate are already quite low when the de-sublimation occurs.

The remaining thermal parameters of the mixture gas (such as density,thermal conductivity and viscosity) were estimated by weighting the mass fractions of the components as:

where φ is a generic variable.The properties of the constituents were obtained from the NIST and the fluid package of Aspen HYSYS.

2.4.Simulation conditions

The assumptions adopted in the simulation are as follows:

(1) the liquid and solid film was ignored and only single gas phase was considered due to the low R134a inlet molar fraction;

(2) the pressure loss along the channel was ignored;

(3) the thickness of the cold wall was ignored since the Biot number (Bi) in present study is much less than 0.1;

(4) only the phase change at the cold wall was considered;

(5) the evaporate and the sublimation was neglected;

(6) the radiation heat transfer was neglected;

(7) the thermodynamic equilibrium was always maintained at the interface.

2.4.1.Meshing

The 2-D fluid domain was meshed with structural grids using ANSYS ICEM.Since the heat transfer resistance of the liquid film when the NCG mass fraction is higher than 10% can be ignored[11],the assumption of ignoring the liquid film is quite reasonable for the conditions studied in this paper (the mass fraction of the NCG nitrogen is 93.08% which is much higher than 10%).As a result,only single-phase fluid domain was simulated in the model.Besides,the mass and heat transfer occur mainly near the cold wall,hence they+should be less than and close to 1 if the viscous sublayer must be simulated [21].As a result,y+was taken as 0.95 in the present simulation.The 2-D fluid domain with the partial grid can be seen in Fig.4.

Fig.4.Schematic of the fluid domain with the partial grid.

2.4.2.Boundary conditions

The inlet and outlet thermal boundary conditions of the mixture gas fluid domain can be seen in Table 2.And the symmetry boundary condition was set to the upper wall which was shown in Fig.4.The Dirichlet boundary condition was used at the cold wall shown in Fig.4.An engineering method for the design of the countercurrent heat exchanger was used to estimate the temperature along the cold wall in present research.In the method,the ratio of the heat-capacity flow rate of the hot and cold fluids,which can be seen in Eq.(17),is compared to 1 to analyze the temperature distribution in the heat exchanger.The average ratio of present research over the temperature range is 0.894 which is close to 1,hence the linear temperature distribution of the two fluids can be assumed approximately [22].

whereqmis the mass flow rate;cpis the constant pressure specific heat capacity.

The heat balance equation of the fluids in the flow channel are expressed as Eq.(18) as:

wherehmixandhcare the convective heat transfer coefficient of the mixture gas and the coolant respectively which was calculated by the Dittus-Boelter equation,W·m-2·K-1;Tmix,TcandTware the mixture gas temperature,the coolant temperature and the cold wall temperature respectively,K;hlis the specific latent heat term,W·m-2.

Since all the parameters in the Eq.(18) are of the first-power,the linear distribution of the two fluids will lead to the linear distribution in the cold wall temperatureTwalong the flow path.Therefore,the temperature distribution at the entire cold wall can be simply fitted by determining the wall temperature at the inlet and outlet using Eq.(19).It can be expressed as:

wherexis the abscissa shown in Fig.4 whose value range is from 0 to 20 m.

2.4.3.Modeling options

Turbulence was modeled by thek-ω shear stress transport (kω SST) model on account of its better near-wall accuracy than thek-ε model [11].The convective terms in the governing equations were discretized using the second-order-upwind scheme.The velocity-pressure coupling was implemented by the SIMPLE method.The species transport model in the ANSYS Fluent was activated to consider the mixture gas.

2.4.4.Grid independence

The grid independence test was carried out to ensure the computing reliability and save the computing resources.The outlet temperatures under different grid quantity are shown in Fig.5.The temperature remains stable when the grid quantity exceeds 400000,hence the grid quantity of 400000 was adopted in the building of mesh.

Fig.5.Mixture outlet temperature with different grid quantity.

3.Results and Discussion

3.1.Verification of the method

The CFD numerical and the experimental result of steam-air condensation heat flux at a cold wall is compared to verify the numerical model as shown in Fig.6.The experimental conditions of the study performed by Chenget al.[23]are displayed in Table 7.

Fig.6.Comparisons of the heat fluxes obtained by the simulation and experiment.

Table 7 Experimental conditions of Cheng’s study [23]

The simulation results agree well with the experiment for different inlet velocity.Heat fluxes are predicted with a maximum deviation of 17.32%as shown in Fig.6.Consequently,the rationality of ignoring the liquid film can be justified,and the results obtained in present paper can be regarded as a reasonable prediction for the condensation and de-sublimation of R134a-nitrogen mixture.

3.2.Condensation and de-sublimation characteristics

The simulated average outlet temperature was 170.26 K,which is close to the assumed result shown in Table 2 (170.6 K).In addition,the average R134a molar fraction at the outlet was reduced to 7.483 × 10-5,which has met the emission requirement (less than 0.01%).Consequently,the validation of the temperatures assumed in Table 2 and the estimation length (20 m) of the heat exchanger were verified.

The partial contours of the mole fraction of R134a simulated based on the aforementioned boundary conditions are shown in Fig.7.The NCG layer with lower R134a concentration expands gradually from the cold wall to the bulk flow region as the proceeding of the condensation or de-sublimation.The expansion of the NCG layer mentioned above will be repeated with the phase change along the flow channel,which leads to the continuous decrease of the average molar fraction of R134a in the bulk flow.

Fig.7.Partial contours of the mass fraction of R134a near inlet and outlet:(a)near the inlet,(b) near the outlet.

The variation of R134a molar fraction along the length of the channel was illustrated in Fig.8.The decrease rate of the molar fraction of R134a keeps slowing down along the length of the channel.This trend can be identified in the Fig.9,which shows the variation of the mass transfer rate near the inlet calculated by the Eq.(13).The mass transfer rate decreases gradually except for the segment near the inlet,whose increasing then decreasing trend isshown clearly in Fig.9(b).The variation in mass transfer rate is related to the subcooling between the interface and the vapor in the grid adjacent to it.The effect of the subcooling in the WCM is reflected by the differential term which can be approximated as/lin Eq.(9).A large subcooling corresponds to a large mass fraction gradient near the interface,which leads to the increasement in the mass transfer rate caused by the high diffusion rate.The difference in mass fraction(ωv-ωv,i)in the range of 0 to 1 m can be seen in Fig.10.The similarity in trend between Figs.10 and 9(b) demonstrates the high correlation between the subcooling(mass fraction gradient term in the WCM)and the mass transfer rate.

Fig.8.Variation of the molar fraction of R134a along the channel.

Fig.9.Mass transfer rates of R134a along the channel: (a) Δx in Fig.3 taken as 0.5 m,(b) Δx in Fig.3 taken as 0.1 m.

Fig.10.Mass fraction difference in the grid adjacent to the cold wall.

The de-sublimation of R134a occurs when the wall temperature is lower than the triple point temperature of R134a (169.85 K).Solid R134a accumulated by de-sublimation at the cold wall will block the channel.Therefore,the estimation of the desublimation volume rate of R134a is the focus to assess the freezing risk of the heat exchanger.The volume de-sublimation rate can be obtained by dividingqm,Δxin Eq.(13)by the density of the solid R134a.The minimum density was assumed as 318.14 kg·m-3,which is 20%of the liquid density at the triple point.Since the wall temperature reaches the triple point at about 19 m,the desublimation accumulation volume rate from 19 m can be calculated and shown in Fig.11.The blockage will happen when the accumulated solid R134a reaches the height of the channel.Hence the effective heat transfer periodtin present simulation conditions was calculated as 36.74 h by Eq.(20) as:

Fig.11.Frozen R134a accumulation rate (Δx in Fig.3 was taken as 0.1 m).

whereqV,Δx,maxis the maximum volume accumulation rate of R134a which is 7.56 mm3·s-1in Fig.11;wand δ are the same as the ones in Eq.(13);Δxis the section length shown in Fig.3 which was taken as 0.1 m in present calculation.

3.3.Effect of the cold wall temperature

The effect of the coolant will be reflected in the change of the cold wall temperature.The change in cold wall temperature was conducted by adding a constant term ΔTinto Eq.(19),which can be written as:

The variation of the molar fraction of R134a along the channel under different ΔTare shown in Fig.12.The outlet molar fraction of R134a decreases gradually and the change in the molar fraction near the inlet decreases more sharply as the wall temperature decreases.Fig.13 shows the corresponding mass transfer rates under different wall temperatures.The trend of mass transfer rate along the channel changes to monotonical decrease when ΔTis greater than zero.Besides the maximum mass transfer rate increases gradually as ΔTincreasing.Subsequently,the larger mass transfer rate dropped more rapidly due to the lower subcooling causing by intensive condensation of the R134a.For example,the mass transfer rate within one meter in Fig.13(b) when ΔTwas taken as 15 K decreases by 94.95%compared to the rate at the inlet,while the corresponding R134a molar fraction in Fig.12 was reduced to 38.75% of the initial molar fraction simultaneously.Therefore,the key factor affecting the outlet R134a molar fraction lies more in the mass transfer rate near the inlet,which also causes the distribution of the molar fraction as shown in Fig.12.

Fig.12.Molar fractions of R134a along the channel under different wall temperature.

Fig.13.Mass transfer rates of R134a along the channel under different wall temperature: (a) Δx in Fig.3 taken as 0.5 m,(b) Δx in Fig.3 taken as 0.1 m.

It can be seen in Figs.12 and 13 that increasing ΔTcan significantly increase the removal performance at various locations in the heat exchanger.When ΔTis increased by 15 K,the outlet R134a molar fraction decreases to 14.57% of the original.Therefore,the increase in the subcooling(ΔT)is proven to be beneficial for shortening the heat exchanger length or reducing the outlet molar fraction of R134a.

The solid R134a accumulation rates under different subcooling(ΔT) are shown in Fig.14.Although the R134a starts desublimating at different locations due to diverse wall temperatures,the maximum accumulation rates keep the approximate constant.It means the effective heat transfer periodtcalculated by Eq.(20) has less correlation with the subcooling (ΔT).On the other hand,the increasement in ΔTwill result in the increasement in the total frozen rate which is shown in Fig.15.The increasement trend is more obvious when the subcooling (ΔT) is low.The total frozen rate increase by 101.3% compared to the original value when ΔTwas taken as 5 K.This phenomenon means that a larger ΔTwill make the cost for recovering the solid R134a higher on account of the larger latent heat.Nevertheless,the reduction of the cold wall temperature is beneficial to the removal of R134a on account of the better removal efficiency.

Fig.14.Solid R134a accumulation rates under different ΔT (Δx in Fig.3 taken as 0.2 m).

Fig.15.Total accumulation rates of solid R134a under different ΔT.

Fig.16.Molar fractions of R134a along the channel under different mass flow rate:(a)distribution along entire channel,(b)distribution near the inlet,(c)distribution near the outlet.

Fig.17.Mass transfer rates of R134a along the channel under different mass flow rate: (a) Δx in Fig.3 taken as 0.5 m,(b) Δx in Fig.3 taken as 0.1 m.

3.4.Effect of changing mass flow rate

The mass flow rate allocated to each channel will change when the number of parallel plate heat exchanger channels varies.Not only the inlet velocity will change,but also the outlet temperature of the mixture gas will change with the mass flow rate,thus the temperature distribution of the cold wall should be recalculated.The cold wall temperatures at different mass flow rate are summarized in Table 8.

Table 8 Wall temperature distribution at different mass flow rate

The molar fraction of R134a and the mass transfer rate at different location are shown in Figs.16 and 17 respectively (ΔTin Eq.(21) was fixed as 10 K).Although the overall molar fraction and mass transfer rate change little with the flow rate,the R134a molar fraction near the outlet rises significantly as the increase in the mass flow rate.An increase in mass flow rate of 12.6% can aggravate the outlet molar fraction to 210.3% of the original value.Therefore,the increase of flow rate is obviously unfavorable for the removal of the R134a at very low molar fraction range.Nevertheless,when the total mass flow rate and the size of the single channels remain constant,a lower mass flow rate distributed to the single channels means the requirement of a larger number of the channels,which will increase the cost of the heat exchanger.

The solid R134a accumulation rate and the total frozen rate under different mass flow rates are shown in Figs.18 and 19 respectively.As shown in Fig.18,it can be found that the location where de-sublimation occurs moves back and the maximum freezing rate increases as the mass flow rate increases.Therefore,the effective heat transfer periodtcalculated by Eq.(20) decreases slightly as the flow rate increases,which changes to 91.59% of the original value when the flow rate increases by 12.6%.The total frozen rate in Fig.19 remains the same value when the mass flow rate increases.While the de-sublimation start point will move back as the mass flow rate increases.It causes the R134a de-sublimating area becomes smaller and the total amount of freezing starts to decrease.It can be predicted that the de-sublimation will disappear and the frozen rate will decrease to zero when the calculated cold wall temperature is always above the triple point of R134a on account of the excessive mass flow rate.

Fig.18.Solid R134a accumulation rates under different mass flow rate(Δx in Fig.3 taken as 0.2 m).

Fig.19.Total accumulation rates of solid R134a under different mass flow rates.

4.Conclusions

Based on the condensation and de-sublimation method,the deep removal of a small amount of R134a from non-condensable gas (nitrogen) was simulated by the wall condensation method(WCM)compiled into Fluent.Besides,the relationship of the saturation temperature and pressure under the triple point of R134a was extrapolated by the Peng-Robinson equation.The influence of key factors on the R134a removal performance was investigated at the basis of revealing the characteristics of R134a condensation and de-sublimation in the gas mixture.The main conclusions are as follows:

(1) A non-condensable gas (NCG) layer develops near the cold wall to the bulk flow region as the temperature decreasing and the mass transfer process proceeding.The continuous decrease of the R134a molar fraction in the mixture gas is conducted by the repeated generation and expansion of this layer.The molar concentration of R134a decreases remarkably near the inlet,which determines relatively higher removal efficiency at this section.

(2) The mass transfer rate of the gas mixture is strongly correlated with the subcooling of the cold wall,which is expressed by the concentration gradient in the first grid adjacent to the cold wall in the WCM.And the maximum mass transfer rate tends to appear near the inlet.When the wall temperature is lower than the triple point temperature,R134a will be de-sublimated and deposited on the cold wall,which may hinder the flow after a long-term operation.The effective heat transfer period before the channel is blocked for present research was calculated as 36.74 h.

(3) The cold wall temperature has a significant effect on R134a removal effect.The reduction in the molar fraction of R134a at various locations in the channel will be promoted by dropping the wall temperature.The outlet molar fraction of R134a decreased to 14.57% of the original value when the temperature of the cold wall is dropped by 15 K,and the effective heat transfer period remained almost unchanged.Therefore,lowering the wall temperature as much as possible is obviously beneficial to the removal of R134a.

(4) The effect of mass flow rate of the mixture gas on R134a removal is mainly reflected in the outlet molar fraction of R134a.A reduction in mass flow rate of 11.6% (from 500 to 444.2 kg·h-1) can reduce the outlet molar fraction to 47.6%of the original.Therefore,reducing the mass flow rate can greatly facilitate the removal of R134a in the low molar fraction range and extend the effective heat transfer period.

In summary,the cold wall temperature and the mass flow rate distributed to each heat transfer channel are important factors to be considered when designing the heat exchangers for the deep removal of VOCs from exhaust gas by the condensation and desublimation method.The lower cold wall temperature and less mass flow rate in each channel under the cost and safety permit mean better removal performance.

Data Availability

Data will be made available on request.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

The work was funded by the National Natural Science Foundation of China (52076159).

Nomenclature

Afaceprojected area of the grid on the cold wall,m2

aparameter used to determine wall temperature,K

BBird suction parameter

BiBiot number

cpconstant pressure specific heat capacity of the mixture gas,J·kg-1·K-1

Dvdiffusion coefficient of the vapor,m2·s-1

hheat transfer coefficient,W·m-2·K-1

hllatent heat of the vapor phase change,J·kg-1

iiteration step

kslope used to determine the wall surface temperature

lvertical distance from the center of the first layer of the mesh to the cold wall,m

Mrelative molecular mass,g·mol-1

nnormal direction of the cold wall

ppressure of the mixture gas,Pa

qm,Δxphase change rate among the section ΔT,kg·s-1

Seenergy source term,W·m-3

Smmass source term,kg·m-3·s-1

Ssspecies source term,kg·m-3·s-1

Sumomentum source term in the abscissa,kg·m-2·s-2

Svmomentum source term in the ordinate,kg·m-2·s-2

Ttemperature of the mixture,K

ΔTchange in cold wall temperature,K

Tcelltemperature in the first cell adjacent to the cold wall,K

Tsatsaturation temperature under the partial pressure of the vapor,K

teffective heat transfer period,h

uabscissa velocity,m·s-1

Vliquid molar volume at normal boiling point,cm3·mol-1

Vcellvolume of the grid adjacent to the cold wall,m3

v ordinate velocity,m·s-1

Wvolumetric molar concentration,mol·m-3

wwide of the channel,m

xlength of abscissa coordinate,m

Δxsection length used in the calculation of the phase change rate,m

ρ mass density,kg·m-3

η dynamic viscosity of the mixture gas,Pa·s

ω mass fraction

θcsuction correction factor

Φ generic variable

φ generic variable

δ thickness of the channel,m

ε ratio of the heat capacity flow rate of the hot and cold fluids

Subscripts

b mixture bulk flow

c coolant (cold nitrogen in present study)

co cold fluid

E eastern point

h hot fluid

i position at gas-liquid interface

inlet flow inlet

m mass

mix mixture gas

sat saturated

ΔTchange in cold wall temperature

uni universal

V volume

v vapor

W western point

w wall

Δxsection length used in the calculation of the phase change rate