Molecular dynamics study on viscosity coefficient of working fluid in supercritical CO2 Brayton cycle: Effect of trace gas

Molecular dynamics study on viscosity coefficient of working fluid in supercritical CO2 Brayton cycle: Effect of trace gas

Journal of CO₂ Utilization 38 (2020) 177–186 Contents lists available at ScienceDirect Journal of CO2 Utilization journal homepage: www.elsevier.com...

2MB Sizes 0 Downloads 21 Views

Journal of CO₂ Utilization 38 (2020) 177–186

Contents lists available at ScienceDirect

Journal of CO2 Utilization journal homepage: www.elsevier.com/locate/jcou

Molecular dynamics study on viscosity coefficient of working fluid in supercritical CO2 Brayton cycle: Effect of trace gas

T

Zhenyu Dua,b, Shuai Denga,b,*, Li Zhaoa,*, Xianhua Niea, Shuangjun Lia,b, Yue Zhanga, Jie Zhaoa,b, Nan Zhengc a

Key Laboratory of Efficient Utilization of Low and Medium Grade Energy (Tianjin University), MOE, Tianjin University, Tianjin, 300072, China International Cooperation Research Centre of Carbon Capture in Ultra-low Energy-consumption, China c School of Chemical Engineering and Technology, Xi'an Jiaotong University, Xi'an, 710049, China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Viscosity coefficient Molecular dynamics simulation Supercritical brayton cycle CO2-based Mixture working fluid

Supercritical CO2 Brayton cycle has been highly recommended as a promising power cycle. As a method to improve the cycle performance, some trace gases could be introduced in the system. Due to rare viscosity coefficient data of CO2-based mixtures over the supercritical regions of CO2, the effect of trace gas on the characteristics of pressure drop has not been clarified yet. Thus, it leads to a challenge in the design and optimization of the supercritical CO2 Brayton cycle when mixing working fluid is considered. In this work, microamount Xe, the mol fraction of which ranges from 1 % to 5 %, is considered as an additive gas. The viscosity coefficient of CO2-Xe was calculated using molecular dynamics simulation. A reliable simulation strategy for the viscosity coefficient calculation was formed, which was proved available since the average relative error of the simulated viscosity coefficient of pure CO2 is around 1 % compared to the REFPROP database. Then, the influence of the addition of Xe on the viscosity coefficient and the pressure drop of the recuperator was discussed. The results demonstrate that the viscosity coefficient would change a lot when CO2 is mixed with Xe. A relative change of viscosity coefficient of 66 % can be observed. Consequently, the pressure drop of the recuperator could enlarge a lot reaching up to 12 %. This work suggests that attention must be paid to the changed viscosity during the design for a reasonable prediction when the trace gas is introduced into the system.

1. Introduction With compact structure, high cycle efficiency and wide range of application, supercritical CO2 (S−CO2) Brayton cycle has been emerging as a promising power cycle for energy conversion systems [1,2]. In an S−CO2 Brayton cycle, the operating conditions of the compressor are just located above the critical point of the working fluid, which brings an extremely low compressor work, and hence a relatively high cycle efficiency could be obtained. It has been reported that the cycle efficiency could reach 47 % when the maximum and minimum temperature of cycle are 873 K and 293 K, respectively [3]. To improve the performance of the cycle, rather than pure CO2, the working fluid used in the supercritical Brayton cycle (SBC) is often mixed with trace gases in theoretical studies [4,5]. The trace gases could be admixtures of substance or impurities presented in the cycle.

The properties of working fluid, including thermodynamics properties and transport properties, would be altered by these trace gases. Intensive researches of effect of trace gases on the performance of SBC have been carried out over the past years. Most of them focus on the effect of shifted critical points on the cycle performance [6–8]. The specific properties of CO2 near critical point are one of key reasons why the S-CO2 Brayton cycle is so charming. Near the critical point, the density of S-CO2 is extremely high, which leads to a significant reduction of compression work. Consequently, the critical point of working fluid generally determines the minimum operating parameters of the cycle, by which the lower minimum temperature would lead to a higher theoretical cycle efficiency. When the trace gas presented in the system, the critical point could be shifted away from the inherent characteristics of pure CO2. From this idea, due to lower critical temperature compared to pure CO2, Xenon (Xe) is recommended as an additive in a

Abbreviations: V, volume, Å 3; S−CO2, Supercritical CO2; x, mol fraction; EoS, Equation of state ⁎ Corresponding authors at: Key Laboratory of Efficient Utilization of Low and Medium Grade Energy (Tianjin University), MOE, Tianjin University, Tianjin 300072, China E-mail addresses: [email protected] (S. Deng), [email protected] (L. Zhao). https://doi.org/10.1016/j.jcou.2020.01.023 Received 15 November 2019; Received in revised form 21 January 2020; Accepted 26 January 2020 2212-9820/ © 2020 Published by Elsevier Ltd.

Journal of CO₂ Utilization 38 (2020) 177–186

D. Zhenyu, et al.

ε EMD σ PACF σαα RD τ AARD η RDF

Nomenclature kB c N f p input q i, j r mix s output T

Boltzmann constant form pressure losses number of state points friction losses pressure, MPa input of channel charge of partical, e particle types distance between atoms, Å mixture standard deviation output of channel temperature, K

Depth of the potential well equilibrium molecular dynamics Collision diameter Pressure autocorrelation function press tensor Relative deviation correlation time Average absolute relative deviation viscosity coefficient Radial distribution function

Subscripts PCHE a Re

Printed circuit heat exchanger acceleration losses Reynolds number

Greek symbols MS

Molecular simulation

macro thermophysical properties, has been emerging as an important tool to make up for the lack in thermophysical property data of supercritical working fluid [17–19]. Differing from EoSs, MS has a solid physical foundation and can be extended to a wider range. As one of the branches of molecular simulation, molecular dynamics (MD) simulation has been wildly utilized to calculate the viscosity of the supercritical fluid. Using seven different potential models, Aimoli et al. [20] calculate the viscosity of pure CO2 for temperatures between 273.15 and 573.15 K and pressures up to 800 MPa. The result indicates that, with EPM2 model, the predicted viscosity could perfectly agree with the reference (overall average absolute relative deviation is below 5 %). However, when it comes to the calculation of viscosity of S-CO2 mixtures, a few pieces of research are performed using molecular dynamics simulation over the supercritical condition of CO2. To address the effect of the Xe with a diverse amount on the performance of the S-CO2 Brayton cycle, the viscosity coefficient is calculated through a serial of MD simulations. A proper MD simulation strategy is concluded in the first step. Then, the viscosity of CO2-Xe mixtures was determined. The calculation was conducted at pressure between 8 and 20 MPa and temperature up to 873 K. Based on the simulated viscosity data, the influence of Xe on the pressure drop of recuperator in SBC was discussed further.

number of publications [6,9,10]. However, limited attention was paid to the effect of transport properties like viscosity on the cycle performance. In fact, the trace gas introduced in CO2 not only affect the critical parameter but also the transport properties like viscosity coefficient. The recuperator, which is in great request for SBC, is inevitably influenced by the changed viscosity coefficient [11]. It has been studied that though the mixture of CO2-cyclohexane, CO2-butane, and CO2-isobutane would lead to a higher cycle performance compared to the pure CO2 for a recompression cycle, however, the pressure drop of these mixtures is higher than that of CO2 in the recuperator [12]. The changed performance of the recuperator largely determines the design and optimization of cycle performance further. And for pure Xe, the viscosity coefficient and density are much higher than those of CO2. As a result, the pressure drop would be definetly affected with the mixing of Xe. The accurate viscosity coefficient of the supercritical CO2-based mixture is the foundation of the above research. But few experimental data is available for most mixture over such an extreme condition [6,13,14]. An alternative is to use the equation of state (EoS) [15]. To author’s knowledge, nearly all the related analyses of exiting researches about the effect of additive on the performance of SBC are performed using an in-house cycle analysis code based on the thermophysical properties from EoSs in REFPROP program [16]. The accuracy of REFPROP properties database plays a decisive role in the feasibility of the research on the additives. The question is that the data from REFPROP program are not always reliable. The EoSs, which own lame physical significance, may not work well in the prediction of thermophysical properties beyond the scope of fitting conditions. As pointed in our previous work [17], completely opposite results may be derived from different versions of REPROP in the prediction on viscosity coefficient of CO2-based mixture. In Fig. 1, the viscosity of CO2-air obtained from different versions of REFPROP is compared. It can be seen that the relative deviations and even their trends obtained by version 9.1 and 10 are not the same. With the increase of the pressure, the relative deviation decreases monotonously for version 10, while behaving as a quadratic curve for version 9.1. Furthermore, the viscosity of the mixture can be lower than the pure for version 10, which is not observed in version 9.1. What’s worse, using REFPROP on the calculation of thermophysical properties of CO2-Xe binary mixture generated severe warning. Binary interaction parameters are not presently available for this mixture. It will be a common case when conducting the calculation of the multi-component mixture. All these calls for a reliable physical properties calculation tool. Molecular simulation (MS), bridging the nanostructure and the

Fig. 1. Comparison of viscosity of impure supercritical CO2 obtained from different versions of REFPROP. 178

Journal of CO₂ Utilization 38 (2020) 177–186

D. Zhenyu, et al.

2. The layout of S-CO2 Brayton cycle

Table 1 The working conditions studied in this work.

Fig. 2 gives a simple recuperated S-CO2 Brayton cycle studied in this paper. Five components are included in the cycle: a turbine, a recuperator, a heat source, a compressor, and a precooler. The operating conditions of each components were also shown in the picture. The system in Fig. 2 was treated as a reference to deduce the possible working conditions of the working fluid. In recuperator, the heat of working fluid exiting the turbine was recovered, while the high-pressure stream exiting the compressor was preheated. To detect the effect of changing viscosity on the pressure drop characteristics in recuperator, the molecular dynamics simulations were performed over the working conditions of the recuperator. The simulated conditions were divided into two ranges including the high-pressure side and the low-pressure side. For the high-pressure side, which refers to the stream exiting the compressor and being heated, two possible pressure points, 15 MPa and 20 MPa, were considered. And 8 MPa and 9 MPa were taken into consideration to mimic the condition of low-pressure side. In both high-pressure and low-pressure sides, the temperature of streams changes a lot through recuperator. Thus, there are two ranges of temperatures for both two streams. High temperature is represented by 673 K, 773 K, and 873 K, while low temperature includes 373 K only. Since the molecular simulation near the critical point of working fluid is not stable, the low temperature was not set lower. The simulated conditions were concluded as Table. 1.

Streams

p (MPa)

Tinput (K)

Toutput (K)

High-pressure stream Low-pressure stream

15, 20 8, 9

673, 773, 873 373

373 673, 773, 873

Table 2 Potential parameters for the CO2, and Xe used in this work.

ε/kB (K) σ (Å) charge (e) length of the bond (Å) degree of the angle (deg)

U = 4εij [(

σij rij

)12 − (

σij rij

)6] +

CO2

Xe

28.129 (C-C) 80.507 (O-O) 2.757 (C-C) 3.033 (O-O) +0.6512 (C) −0.3256 (O) 1.149 (C-O) 180 (C-O-C)

221.0 – 4.1 – – – – –

qi qj 4πε0 rij

(1)

where εij, σij and rij are LJ well depth, LJ size, and separation distance respectively; ε0, qi, and qj are the dielectric constant of vacuum and partial charges for atom i and j, respectively. All potential parameters for the CO2, and Xe are described in Table. 2. Besides, the LJ parameters between the unlike atoms were obtained from the Lorentz-Berthlot combing rule:

3. Method 3.1. Force field

εij =

The molecular models of CO2, and Xe were taken from Harris et al. [21], and Rappe et al. [22], respectively. All models owning the same number of LJ centers as corresponding molecule nucleus have been considered in this study. And they are rigid with fixed bond length and fixed angle for there is limited even no improvement of accuracy with flexible models for the calculation of viscosity coefficient [20]. As a result, only the intermolecular potentials including Lennard-Jones 12-6 potential and electrostatic interaction potential are considered here to describe the system, defined as

σij =

εi εj

(2)

σi + σj 2

(3)

3.2. Viscosity coefficient calculation The transport properties in this work were obtained from equilibrium molecular dynamics (EMD) simulation. Viscosity values are determined by the Green–Kubo formulation [23], which based on the time

Fig. 2. The layout and the operation parameters of the reference recuperated S-CO2 Brayton cycle. 179

Journal of CO₂ Utilization 38 (2020) 177–186

D. Zhenyu, et al.

integral of the pressure autocorrelation function (PACF) of the press tensor as follows [24]:

η=

V kB T



∫ 〈pαβ (t0 + τ ) pαβ (t0 ) 〉 dτ 0

(4)

where α and β refer to a pair of Cartesian coordinates x, y, and z. V, T, and kB stand for the volume, temperature, and Boltzmann constant, respectively. To improve the statistics accuracy of the simulation, all elements of press tensor, including the diagonal elements and off-diagonal elements, were used in Eq. 4 to calculate the viscosity. The pαβ is derived as following [24]:

1 3

pαα = σαα −

pαβ =

∑ σγγ γ

(5)

(σαβ + σβα ) (6)

2

where the σαβ represent the elements of the pressure tensor. So the viscosity expression could be written as:

η=

V 10kB T

Fig. 3. The simulated density of pure CO2 under studied working conditions. The state points was classified into two regions according to the density. The point with lowset density in regions was illustrated in the right part of figure.



∫ 〈pαβ (t0 + τ ) pαβ (t0 ) 〉 dτ 0

(7)

simulations were performed beforehand to determine a proper correlation time and the number of independent runs. The effect of correlation time was investigated first. As demonstrated above, the correlation time was selected according to density and the density of pure CO2 under supercritical conditions as given in Fig. 3. To minimize the total length of simulation time, working conditions are divided into two regions, namely “Region 1” surrounded by the red dash line and “Region 2” surrounded by the blue dash line. For each region, the correlation time is set by the preliminary simulation of the state point with the lowest density within the region. From this principle, working condition with temperature of 873 K and pressure of 8 MPa for “Region 1”, and working condition with temperature of 873 K and pressure of 15 MPa for “Region 2” are selected to perform preliminary simulations. The density of the selected state point is also shown in the figure and surrounded by lines owning corresponding color. The PACF and the calculated viscosity of two selected state points using diverse correlation time are given in Fig. 4. For PACF, lines with different colors denote results ⋅from different elements of press tensor. While for the viscosity coefficient, lines with different colors indicate calculated values derived from different correlation time. It can be observed that the PACFs of all elements of press tensor begins decaying to zero at correlation time equal to about 3000 fs for the state point in the “Region 1”, which means that the correlation time used to determine the viscosity coefficient should at least excesses 3000 fs for the “Region 1”. The Fig. 4(b) demonstrates the same conclusion. When correlation time is relatively low, the calculated viscosity values keep increasing with correlation time. Until correlation time is high enough equal to about 3000 fs, the viscosity coefficient does not change with the correlation time. Taking potential statistics error into consideration, the correlation time equal to 3500 fs was chosen to determine the viscosity coefficient for state points in the “Region 1”. To get sufficiently reliable calculated viscosity coefficient data, 7 ns simulation time was performed during the production run with sampling the press tensor every 5 fs. In total, 400 viscosity coefficient values were generated during an individual production run. Compared with the PACFs of selected state point in the “Region 1”, the PACFs from “Region 2” decays to zero in a shorter time equal to about 2000 fs. Similarly, a final correlation time of 2500 fs was selected to calculate the viscosity coefficient for state points from “Region 2”. To ensure the quality and amount of data consistent with that in the “Region 1”, a production run of 5 ns was exerted for state points in the “Region 2”. Then, to detect a proper number of runs, 30 individual runs were performed over the condition for systems owning the highest density

The reason why a factor 10 instead of 9 is presented in the denominator is that the integral of the diagonal component yields a value equal to 4η/3. 3.3. Simulation details All simulations were performed with the LAMMPS package. The standard velocity-Verlet algorithm was applied to integrate the equations of motion. For CO2, a linear three-atom molecule, the algorithm proposed by the Kamberaj et al. [25] was utilized to keep the molecule rigid. Standard Ewald summation method was chosen to compute the long-range Coulombic interactions with a desired relative error in forces of 1 × 10−4, and the cutoff for Lennard Jones interactions was taken to be 12 Å. A time step of 1 fs, as well as usual periodic boundary conditions in three directions x, y, and z, were employed for all simulations. The simulated system consists of 1000 particles which are enough to ensure the accuracy of calculations. A run of 1 ns was performed to reach the equilibration before the production runs for viscosity coefficient calculations. The NpT ensemble was applied during all the runs. It has been shown that EMD simulations yield a value for the viscosity coefficient in the NpT ensemble similar to that in the NVT ensemble [26]. The viscosity calculation in EMD simulations is challenging since it is subject to high statistical uncertainties [27]. Theoretically, the PACF decays to zero in the long time limit and the time integral of PACF calculated by Eq. (4) reaches a constant value, which is exactly the calculated viscosity. The correlation time τ in the PACF determines the accuracy and the reliability of the calculation. It has been shown that the correlation time has a lot of things to do with the density of the studied system [17,28]. With the same correlation time, PACF of a system with high density can completely converge, while that of a system with a low density may not converge at all. But it does not mean a high correlation time is the best all the time. Fewer data are generated with a higher correlation time, which leads to larger statistical uncertainty, and hence more production runs for viscosity calculation are needed in simulations [29]. On the other hand, due to the large statistical noise at a long time, the integral of PACF does not necessarily converge to a constant value. Generally, a long simulated trajectory could be produced to improve the statistics. However, this method tends to be too time-consuming. An alternative is to run multiple independent trajectories for each point and the average of the running integrals is regarded as the final viscosity value. To obtain a reliable viscosity value, a series of preliminary 180

Journal of CO₂ Utilization 38 (2020) 177–186

D. Zhenyu, et al.

Fig. 5. The effect of the number of run on the calculation of viscosity coefficient: (a), (b): the state point with the lowest density; (c), (d): the state point with the highest density.

Fig. 4. The PACF and the calculated viscosity using diverse correlation time of two selected state points: (a), (b): the state point with the lowest density in the “Region 1”; (c), (d): the state point with the lowest density in the “Region 2”.

181

Journal of CO₂ Utilization 38 (2020) 177–186

D. Zhenyu, et al.

(673 K and 20 MPa) and lowest density (873 K and 8 MPa), respectively. What should be pressed is that the random seed of velocity is distinct in each run during the initialization of configuration. The calculated viscosity coefficient of each individual run for two working conditions is given in Fig. 5. For each individual run, the result fluctuates a lot. And all the curves deviate from one another to a certain degree. The average of all curves is also presented in Fig. 5 shown as a red line. It can be seen that the curve of average results is quite smooth compared to that from individual run. On the other hand, a plateau is somehow identified for each run within production time. Compared with the simulated system with low density, the calculated viscosity coefficient of system with high density using a smaller correlation time tend to convergence to a constant value in a shorter time. Consequently, the viscosity coefficient data during the final proper length of production run were averaged. For state points belong to two different regions presented in Fig. 3, a proper length of the final production run was intercepted to make the amount of valid data equal to 100. The viscosity value of an individual run is concluded that way. Further, taking an average to such viscosity values of all individual runs, a final viscosity can be derived. Following this principle, the standard deviation of simulation could also be determined. To minimize the total length of running time as much as possible, the effect of the individual running number on the final viscosity coefficient under two conditions was given in Fig. 5(b) & (d), where the grey point is the viscosity coefficient value of individual run and the blue point is the average value of individual runs. And the red line, equal to the average value of 30 individual runs, is shown as a reference line. The confidential interval with a 95 % confidence of the average value of 30 individual runs is calculated using the Student’s t-value and represented as the red region. Viscosity coefficient values of the individual run are different from each other. It can noticed that the final average viscosity values vary with running numbers during the initial stage. For both two working conditions, the final average viscosity value tends to enter the confidential interval when the number of run increases to 10 times, and then reach a constant once the number of run is higher than 15 times. Thus, 15 individual runs at each working conditions were performed and the final viscosity was determined by averaging the viscosity values of 15 individual runs.

AARD =

Ns

∑ RDi i=1

(9)

where the Ns is the total number of state points and the summation over RDs of all simulated results is performed. As can be seen in Table 3, the overall average absolute relative deviation of the viscosity coefficient from MD simulations under all temperatures is around 1 %. As a reference, the average error of the CO2 viscosity coefficient using EPM2 model under supercritical conditions by Aimoli et al. [20] is 4.69 %. It can be concluded that the simulation strategy employed in this research could work successfully in the reliable prediction of the viscosity coefficient over supercritical conditions of CO2.

4.2. Viscosity coefficient of the binary mixture of CO2-Xe The viscosity coefficient of the binary mixture of CO2 and Xe with different components, as well as the standard deviation, at different working conditions is given in Fig. 8. Except for the data under the highest pressure, over the working condition in this work, the viscosity coefficient of binary mixture tends to increase with the increase of temperature and pressure overall, which is basically consistent with the trend of pure CO2. It is worth noting that the viscosity coefficient of both two binary mixtures increases obviously when pressure raises from 15 MPa to 20 MPa at 373 K. Such phenomenon has not been observed in the cases of pure CO2. Radial distribution function (RDF) was calculated to illustrate the structural change of working fluid when pressure enhances from 15 MPa to 20 MPa at 373 K. RDF denotes the spherically averaged local order around a center atom, yielding the probability of finding another particle at a given distance apart from the center atom [24]. In Fig. 9, the RDF between carbon atoms of carbon dioxide for the system containing different components of Xe is presented. When pressure locates at 15 MPa, the RDFs between carbon atoms in different systems are quite similar. However, when pressure increased to 20 MPa, the first peak of RDF in the mixture systems is higher than that in the pure CO2. The first peak height indicates the local structural order of carbon dioxide molecules. A higher peak generally means the enhancement of local structure order, which leads to a higher viscosity consequently [30]. In a word, near the critical temperature of CO2, the existence of trace Xe is likely to induce the carbon dioxide molecules to form a more organized structure when pressure is relatively high. Such results explain why a surge in viscosity of CO2-Xe mixture is observed when pressure increases from 15 MPa to 20 MPa.

4. Results 4.1. Validation of pure CO2 The viscosity coefficient of pure CO2, as well as the standard deviation, at different working conditions, is given in Fig. 6. On the whole, a growth of the viscosity coefficient of pure CO2 with the increase of pressure and temperature could be observed except for the cases under 20 MPa. A few exceptions may be the joint effect of the limited inherent difference and relatively high simulated uncertainty due to measuring repeatability which is given in the Fig. S1. To demonstrate the quality of the viscosity coefficient calculated by MD simulations in this work, the data from REFPROP 10.0 was also shown in the Fig. 6. An excellent agreement between the simulated data and the reference could be observed in the figure. A quantitative comparison between the simulated data and the reference was made over the studied range further. The relative deviation (RD) between viscosity values from MD simulation and REFPROP is determined as Eq. (8) and shown in Fig. 7.

RD =

1 Ns

ηMD − ηREFPROP ηREFPROP

(8)

According to the figure, the viscosity values from MD simulations agree very well with the REFPROP database. The absolute relative errors range from 1.08 % to 3.36 %. The average absolute relative deviation (AARD) of simulated results when compared to REFPROP data under different temperatures were also calculated as follows:

Fig. 6. The viscosity coefficient of pure CO2 from MD simulations and REFPROP. 182

Journal of CO₂ Utilization 38 (2020) 177–186

D. Zhenyu, et al.

Fig. 7. The comparison between the viscosity coefficient from REFPROP 10.0 and simulations.

Table 3 AARDs of simulated viscosity coefficient of pure CO2.

AARD (%)

373 K

673 K

773 K

873 K

Total

1.51

1.55

1.33

0.64

1.26

5. Discussion 5.1. The effect of the addition of Xe on the viscosity of working fluid As emphasized in Section. 1, the addition of Xe in the S-CO2 Brayton cycle is likely to influence the pressure drop characteristics, and the extent of change on the viscosity coefficient of working fluid is the foundation of that research. Thus the relative change of viscosity coefficient between binary mixtures CO2-Xe and pure CO2 was summarized in Fig. 10 according to the simulated viscosity coefficient listed in Section. 4. Clearly, when the temperature is relatively low at 373 K, the change on the viscosity due to the addition of Xe is extremely obvious. The highest relative change of the viscosity coefficient could be up to 66 % when 5 % of Xe is presented in the system under 20 MPa. Even at high temperature, the viscosity of working fluid could change a lot, reaching up to 13 %. Such results highlight that the attention to the changed viscosity must be paid when the trace gas is introduced in the system, otherwise, it would be counterproductive. At all the temperatures, the relative change value on viscosity assumed an increasing tendency with a higher pressure. Interestingly, when temperature is relatively high (i.e., T > 673 K), the relative change could change from positive to negative. Particularly with higher temperature, when the pressure keeps decreasing, the condition that the relative change is negative appears earlier. And the absolute relative change values keep enlarging. Such changes once again stress that the measurement of viscosity coefficient of S-CO2 mixture should be listed in agenda.

Fig. 8. The viscosity coefficient of CO2-Xe from MD simulations. (a) mol fraction of Xe equal to 0.01; (b) mol fraction of Xe equal to 0.05.

recuperator. The straight channel is used for the PCHE and the channel diameter Dh of PCHE is equal to 2 mm in this study. The pressure drop in the channel of PCHE can be calculated as Eq. (10):

Δp = Δpf + Δpc + Δpa

(10)

where Δpf, Δpc, Δpa are the friction losses, form pressure losses and acceleration losses, respectively. Among the three components of pressure drop, the friction losses are related to the viscosity tightly. In this work, only the friction losses are considered and can be determined as follows:

5.2. The effect of trace gas on the pressure drop characteristics To detect the effect of changed viscosity coefficient of working fluid on the performance of recuperator, pressure drop characteristics were analyzed using data from MD simulations. What should be emphasized is that this work is a preliminary analysis on the pressure drop only considering the change of viscosity. Actually, the pressure dorp in the heat exchanger is also influenced by several other important factors including mass flow rate and frontal area in addition to viscosity. Generally, a printed circuit heat exchanger (PCHE) was adopted as a heat exchanger in the S-CO2 Brayton cycle to reduce the size of the

Δpf = f

L v2 ρ Dh 2

(11)

In Eq. (11), L is the length of PCHE, while f is the friction factor and can be calculated with Reynolds number:

f=

183

15.78 Re,Re <2300 ⎧ ⎪ 2 1⎛ 1 ⎞⎟ , Re≥2300 ⎨ ⎜ ⎪ 4 ⎝ 1.8 log ( Re ) − 1.5 ⎠ ⎩

(12)

Journal of CO₂ Utilization 38 (2020) 177–186

D. Zhenyu, et al.

Fig. 9. RDFs between carbon atoms of carbon dioxide with a different component of Xe at 373 K under different pressures: (a) 15 MPa; (b) 20 MPa.

The velocity of the fluid in this study is set as 1 m/s. Owing to that the state points with calculated viscosity coefficient are limited, it is impractical to perform an iteration calculation on the pressure drop. In a compromising approach, the pressure drop was put aside firstly to locate the state of output so that the thermophysical properties of both input point and output point can be estimated. According to the working conditions of two kinds of streams in recuperator shown in Table 1, all 12 possible cases were generated and summarized in Table 4. Then, the pressure drop of each case can be derived. The absoulute pressure drop and the relative change of pressure drop between mixtures and pure CO2 in the recuperator is shown in Fig. 11. Due to the increase in the viscosity coefficient and density when Xe is introduced, the pressure drops of mixture working fluids are higher than those of pure CO2. The influence is limited with 1 % mol fraction of Xe mixed in the system. But when the working fluid contains 5 % of Xe, the pressure drop could be 12 % larger. It means that, even though the addition of Xe could level up the total cycle efficiency according to the existing publications, such additions may be not conducive to the S-CO2 Brayton cycle. Obviously, when the same component of Xe is mixed, the pressure drop tends to increase with higher absolute pressure of channel. However, the pressure drop changes a

Fig. 10. The effect of the addition of Xe on the viscosity coefficient of working fluid at different temperatures: (a) 373 K; (b) 673 K; (c) 773 K; (d) 873 K.

184

Journal of CO₂ Utilization 38 (2020) 177–186

D. Zhenyu, et al.

Table 4 The summary of cases for the research on the pressure drop.

Low-pressure stream

Cases

Sites

T(K)

p(MPa)

Case1

input output input output input output input output input output input output input output input output input output input output input output input output

373 673 373 773 373 873 373 673 373 773 373 873 673 373 673 373 773 373 773 373 873 373 873 373

8 8 8 8 8 8 9 9 9 9 9 9 15 15 20 20 15 15 20 20 15 15 20 20

Case2 Case 3 Case 4 Case 5 Case 6 Case 7 High-pressure stream Case 8 Case 9 Case 10 Case 11 Case 12

little with temperature for both low-pressure stream and high-pressure stream. It indicates that higher maximum temperature of the cycle could enhance the theoretical cycle efficiency but not enlarges the pressure drop significantly.

6. Conclusion The viscosity coefficient of CO2 mixtures is a necessity for the design and optimization of the S-CO2 Brayton cycle when trace gas is introduced in the system. In this work, the viscosity coefficient of CO2-Xe in the supercritical region of CO2 is calculated using MD simulations. A reliable simulation strategy is constructed. CO2 mixtures containing different mol fractions of Xe are considered. Based on the simulated thermophysical properties, the effect of the addition of Xe on the viscosity coefficient and the pressure drop of the recuperator is discussed in detail. Some conclusions can be drawn as follows: (1) Over the working conditions in this work, with fifteen diverse trajectories for one state point, a reliable viscosity coefficient can be determined when the correlation time for the high-density region and low-density region are equal to 2500 fs and 3500 fs, respectively. Compared to the REFPROP data, the AARDs of the viscosity coefficient of pure CO2 from MD simulations is around 1 %. (2) The viscosity coefficient could change a lot under high pressure at 373 K when Xe is mixed. The highest relative change of the viscosity coefficient could be up to 66 % when 5 % of Xe is introduced. Even at a high temperature, the relative change could reach up to 13 %. (3) Due to the change in the viscosity coefficient and density, when the working fluid containing 5 % of Xe, the pressure drop could be 12 % larger. The relative change of pressure drop in recuperator tends to enlarge with higher pressure, but it is not sensitive to the temperature. The addition of Xe may be not conducive to the S-CO2 Brayton cycle due to its impact on the characteristic of pressure drop. Such results revealed by MD simulations indicate that a further experimental measurement on the viscosity coefficient as well as other thermophysical properties of CO2-based mixtures is urgently needed if trace gas is introduced in the system. Simulated data in this work could be powerful information for the design and optimization of S-CO2 Brayton cycle.

Fig. 11. The pressure drop of recuperator: (a) The absolute pressure drops for different working fluids; (b) The effect of the addition of Xe on the pressure drop.

CRediT authorship contribution statement Zhenyu Du: Methodology, Writing - original draft, Formal analysis. Shuai Deng: Supervision, Project administration. Li Zhao: Supervision, Funding acquisition. Xianhua Nie: Methodology, Software. Shuangjun Li: Methodology, Validation. Yue Zhang: Investigation. Jie Zhao: Investigation. Nan Zheng: Funding acquisition. 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. Acknowledgment The authors are grateful for the support provided by the National Key Research and Development Program of China under Grant No. 2018YFB1501004, China's Post-doctoral Science Fund under Grant No. 185

Journal of CO₂ Utilization 38 (2020) 177–186

D. Zhenyu, et al.

2018M631155.

[14] H. Li, Ø. Wilhelmsen, Y. Lv, W. Wang, J. Yan, Viscosities, thermal conductivities and diffusion coefficients of CO2 mixtures: review of experimental data and theoretical models, Int. J. Greenh. Gas Con. 5 (2011) 1119–1139. [15] S. Li, S. Deng, L. Zhao, R. Zhao, M. Lin, Y. Du, Y. Lian, Mathematical modeling and numerical investigation of carbon capture by adsorption: literature review and case study, Acs Appl. Energy Mater. 221 (2018) 437–449. [16] L. EW, H. ML, M. MO. REFPROP, reference fluid thermodynamic and transport properties. NIST: Standard Reference Database. [17] X. Nie, Z. Du, L. Zhao, S. Deng, Y. Zhang, Molecular dynamics study on transport properties of supercritical working fluids: literature review and case study, Acs Appl. Energy Mater. 250 (2019) 63–80. [18] P. Ungerer, C. Nieto-Draghi, B. Rousseau, G. Ahunbay, V. Lachet, Molecular simulation of the thermophysical properties of fluids: From understanding toward quantitative predictions, J. Mol. Liq. 134 (2007) 71–89. [19] X. Nie, L. Zhao, S. Deng, W. Su, Y. Zhang, A review of molecular simulation applied in vapor-liquid equilibria (VLE) estimation of thermodynamic cycles, J. Mol. Liq. 264 (2018) 652–674. [20] C.G. Aimoli, E.J. Maginn, C.R.A. Abreu, Transport properties of carbon dioxide and methane from molecular dynamics simulations, J. Chem. Phys. 141 (2014) 134101. [21] J.G. Harris, K.H. Yung, A. Abdelrasoul, Carbon dioxide’s liquid-vapor coexistence curve and critical properties As predicted by a simple molecular model, J. Phys. Chem. A 99 (1995) 12021–12024. [22] D. Dubbeldam, S. Calero, D.E. Ellis, R.Q. Snurr, RASPA: molecular simulation software for adsorption and diffusion in flexible nanoporous materials, Mol. Simulat. 42 (2015) 81–101. [23] R. Kubo, Statistical-mechanical theory of irreversible processes. I. General theory and simple applications to magnetic and conduction problems, Manag. Financ. 25 (2007) 16–29. [24] G. Raabe, Molecular Simulation Studies on Thermophysical Properties, Springer, Singapore, 2017. [25] H. Kamberaj, R.J. Low, M.P. Neal, Time reversible and symplectic integrators for molecular dynamics simulations of rigid molecules, J. Chem. Phys. 122 (2005) 224114. [26] C. Nieto-Draghi, T. de Bruin, J. Pérez-Pellitero, J. Bonet Avalos, A.D. Mackie, Thermodynamic and transport properties of carbon dioxide from molecular simulation, J. Chem. Phys. 126 (2007) 64509. [27] Y. Zhang, A. Otani, E.J. Maginn, Reliable viscosity calculation from equilibrium molecular dynamics simulations: a time decomposition method, J. Chem. Theory Comput. 11 (2015) 3537–3546. [28] U.A. Higgoda, R. Hellmann, T.M. Koller, A.P. Fröba, Self-diffusion coefficient and viscosity of methane and carbon dioxide via molecular dynamics simulations based on new ab initio-derived force fields, Fluid Phase Equilibr. 481 (2019) 15–27. [29] X. Yang, C. Duan, J. Xu, Y. Liu, B. Cao, A numerical study on the thermal conductivity of H2O/CO2/H2 mixtures in supercritical regions of water for coal supercritical water gasification system, Int. J. Heat Mass Transf. - Theory Appl. 135 (2019) 413–424. [30] N. Sohrevardi, M.R. Bozorgmehr, M.M. Heravi, M. Khanpour, Transport properties of mixtures composed of acetone, water, and supercritical carbon dioxide by molecular dynamics simulation, J. Supercrit. Fluids 130 (2017) 321–326.

Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.jcou.2020.01.023. References [1] R.S. Norhasyima, T.M.I. Mahlia, Advances in CO₂utilization technology: a patent landscape review, J. Co2 Util. 26 (2018) 323–335. [2] L. Teng, Y. Xuan, Design of a composite receiver for solar-driven supercritical CO2 Brayton cycle, J. Co2 Util. 32 (2019) 290–298. [3] W.S. Jeong, J.I. Lee, Y.H. Jeong, Potential improvements of supercritical recompression CO2 Brayton cycle by mixing other gases for power conversion system of a SFR, Nucl. Eng. Des. 241 (2011) 2128–2137. [4] L. Vesely, K.R.V. Manikantachari, S. Vasu, J. Kapat, V. Dostal, S. Martin, Effect of impurities on compressor and cooler in supercritical CO2 cycles, J. Energy Resour. Technol. 141 (2019). [5] C.M. Invernizzi, T. van der Stelt, Supercritical and real gas Brayton cycles operating with mixtures of carbon dioxide and hydrocarbons, Proceedings of the Institution of Mechanical Engineers, Arch. Proc. Inst. Mech. Eng. Part A J. Power Energy 19901996 226 (2012) 682–693. [6] Y.H. Jeong, W.S. Jeong, Performance of supercritical Brayton cycle using CO2-based binary mixture at varying critical points for SFR applications, Nucl. Eng. Des. 262 (2013) 12–20. [7] J. Guo, M. Li, Y. He, J. Xu, A study of new method and comprehensive evaluation on the improved performance of solar power tower plant with the CO2-based mixture cycles, Acs Appl. Energy Mater. 256 (2019) 113837. [8] L. Hu, D. Chen, Y. Huang, L. Li, Y. Cao, D. Yuan, J. Wang, L. Pan, Investigation on the performance of the supercritical Brayton cycle with CO2-based binary mixture as working fluid for an energy transportation system of a nuclear reactor, Energy 89 (2015) 874–886. [9] J. Guo, M. Li, J. Xu, J. Yan, K. Wang, Thermodynamic performance analysis of different supercritical Brayton cycles using CO2-based binary mixtures in the molten salt solar power tower systems, Energy 173 (2019) 785–798. [10] L. Vesely, K.R.V. Manikantachari, S. Vasu, J. Kapat, L.J. Abbott, Effect of Mixtures on Compressor and Cooler in Supercritical Carbon Dioxide Cycles, ASME Turbo expo 2018, Turbomachinery Technical Conference and Exposition, Oslo, Norway, 2018. [11] Y. Jiang, E. Liese, S.E. Zitney, D. Bhattacharyya, Design and dynamic modeling of printed circuit heat exchangers for supercritical carbon dioxide Brayton power cycles, Acs Appl. Energy Mater. 231 (2018) 1019–1032. [12] X. Liu, Z. Xu, Y. Xie, H. Yang, CO2-based mixture working fluids used for the drycooling supercritical Brayton cycle: thermodynamic evaluation, Appl. Therm. Eng. 162 (2019) 114226. [13] C.W. White, N.T. Weiland, Evaluation of property methods for modeling directsupercritical CO2 power cycles, J. Eng. Gas Turbine. Power 140 (2018) 11701.

186