- Email: [email protected]
A non-equilibrium molecular dynamics study of subcritical, supercritical and transcritical mixing of liquid-gas systems
Journal Pre-proofs A Non-Equilibrium Molecular Dynamics Study of Subcritical, Supercritical and Transcritical Mixing of Liquid-Gas Systems Farzin Rahmani, Timothy Weathers, Ashvin Hosangadi, Yee C. Chiew PII: DOI: Reference:
S0009-2509(19)30914-5 https://doi.org/10.1016/j.ces.2019.115424 CES 115424
To appear in:
Chemical Engineering Science
Received Date: Revised Date: Accepted Date:
12 September 2019 22 November 2019 9 December 2019
Please cite this article as: F. Rahmani, T. Weathers, A. Hosangadi, Y.C. Chiew, A Non-Equilibrium Molecular Dynamics Study of Subcritical, Supercritical and Transcritical Mixing of Liquid-Gas Systems, Chemical Engineering Science (2019), doi: https://doi.org/10.1016/j.ces.2019.115424
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2019 Published by Elsevier Ltd.
A Non-Equilibrium Molecular Dynamics Study of Subcritical, Supercritical and Transcritical Mixing of Liquid-Gas Systems
Farzin Rahmani(a), Timothy Weathers(b), Ashvin Hosangadi(b), Yee C. Chiew*(a), (a)
Department of Chemical & Biochemical Engineering, Rutgers University – New Brunswick, Piscataway NJ 08854 (b) Combustion Research and Flow Technology, Inc., Pipersville, PA 18947
* Corresponding author. Email: [email protected]
1
Abstract A non-equilibrium molecular dynamics (NEMD) simulation method has been developed to simulate mixing of a liquid with a vapor and identify the characteristics of the liquid-vapor interface in the mixing process. Our results show that for the case of subcritical mixing, local phase equilibrium is established, and saturated liquid and saturated vapor are formed on either side of the interface at the prevailing saturation temperature that is fixed by the pressure of the system independent of the temperature difference across the interface. For transcritical mixing, significant clustering of molecules is found at the “transitional diffuse” interface. We found that the local density and local temperature of the transcritical interface can be directly mapped to the near critical region of the Ar/Ar and Kr/Ar vapor-liquid phase diagram. Finally, for the case of supercritical mixing, our simulations show a gradual change of density consistent with diffusioncontrolled fluid mixing. Keywords: Non-equilibrium molecular dynamics simulation, Subcritical Mixing, Transcritical Mixing, Supercritical Mixing, Liquid-vapor Interface
2
1. Introduction Understanding mixing and evaporation of an injected liquid into a high-pressure gaseous environment is of relevance and importance to a wide range of technological applications such as propulsion systems including liquid rocket engines, gas turbine combustors, and diesel engines. Depending on the ambient pressure relative to the critical pressure of the liquid, mixing of the injected liquid can exhibit different mixing processes: subcritical, supercritical, and transcritical mixing mechanisms(Givler, Abraham et al. 1996, Bellan and science 2000, Yang 2000). Highpressure experimental studies(H. Mayer, A. Schik et al. 1998, Segal and Polikhov 2008, Chehroudi 2012) have provided insights on the behavior of mixing of multi-phase systems. When liquid fuel is injected into a chamber with pressure below the critical pressure of the fuel, classical jet breakup and droplet formation are observed. A well-defined molecular interface separates the injected liquid from ambient gases due to surface tension, and interactions between shear and surface tension forces promote droplet formation and secondary jet breakup. In contrast, at pressures exceeding the liquid critical pressure, the mixing process can become quite different. Under these conditions, the liquid-gas interface diminishes and is replaced by one-phase dense mixing layer which is affected by thermodynamics and transport processes.
Transcritical mixing that transitions from two-phase mixing regime to one-phase supercritical mixing has been observed in several studies.(Huai, Koyama et al. 2005, Chen, Zhang et al. 2013) Segal and Polikhov(Segal and Polikhov 2008) reported that transcritical mixing is characterized by droplet formation and irregularly shaped material break-up from the injectant. Detailed analysis suggests that the transitional changes are controlled by the multi-component
3
nature of the thermodynamics, transport processes and stability of gas-liquid interface(Oschwald*, Smith et al. 2006, Chehroudi 2012). Computational fluid dynamics and continuum-based models have been developed to model supercritical mixing. Recent theoretical and modeling studies (Mo and Qiao 2017) have employed real gas equations of state to quantify the thermodynamic properties and critical point loci of binary and multi-component fluid mixture within an equilibrium framework. Dahms and Oefelein (Dahms and Oefelein 2013) proposed a criterion, based on the thickness of equilibrium liquid-gas interface, to characterize an interface as subcritical or supercritical in behavior. Although these studies have provided insights on transcritical mixing, fundamental understanding of this complex process is still lacking due to considerable challenges involved in characterizing the thermodynamics properties and vapor-liquid equilibrium of multi-component mixtures especially in the vicinity of the mixture critical point. Moreover, numerical macroscopic models are not able to capture the fine scales needed to resolve the microscopic structure of the liquid-gas interface that in part governs the underlying processes subject to large density and temperature gradients.
Nonequilibrium molecular dynamics (NEMD) simulation(Meland 2002, Meland 2003, Bedeaux and Kjelstrup 2004, Meland, Frezzotti et al. 2004, Simon, Bedeaux et al. 2006, Xu, Kjelstrup et al. 2006) is an invaluable tool that complements macroscopic continuum methods for probing the small length scale and short time scale phenomena associated with the physical changes of mixing interface in response to variation of temperature and pressure. In contrast to continuum-based methods where a priori assumptions about the process are necessary, molecular dynamics simulations do not require assumed models for equations of state, equilibrium liquid/vapor interface, thermodynamic and transport properties. It requires only inter-molecular
4
potential energy function as input and generates the evolutionary process of mixing during the simulation. Such a study can elucidate the physics of the mixing mechanisms and test the assumptions used in different modeling approaches. Additionally, it provides information at the microscopic molecular scale, e.g., clustering of molecules and particle correlation, that are not possible from macroscopic models. Molecular dynamics simulation is computationally expensive, and the system studied is limited to several hundred to thousand nanometers and over a time duration of several hundred nanoseconds.
A number of researchers have used molecular dynamics to model the interface of onecomponent subcritical liquid-gas systems at equilibrium and nonequilibrium situations. They showed that a key parameter that controls the interface characteristic is the surface tension which is governed by the temperature at the interface. The aim of these NEMD studies was to verify the validity of proposed theories for predicting the temperature, transport properties, condensation and evaporation rates in the vicinity of interface(Kjelstrup and Hafskjold 1996, Røsjorde, Kjelstrup et al. 2001, Simon, Kjelstrup et al. 2004, Tsuruta and Nagayama 2005, Xu, Kjelstrup et al. 2006, Wilhelmsen, Trinh et al. 2015, Savin, Schweizer et al. 2016, Schweizer, Öttinger et al. 2016, Schweizer, Öttinger et al. 2016). For example, Rosjorde et al.(Røsjorde, Fossmo et al. 2000, Røsjorde, Kjelstrup et al. 2001) determined the coefficients for transfer of heat and mass across the liquid-vapor interface of one-component fluid by NEMD simulation. They used the SoaveRedlich-Kwong equation of state to test the assumption of local equilibrium and showed that local equilibrium is established not only in the bulk phases but also at the liquid-vapor interface. In recent years, NEMD simulations have been carried out to study the liquid-vapor interface of twocomponent mixtures(Olivier, Rollier et al. 2002, Bedeaux and Kjelstrup 2018). Inzoli et al.(Inzoli,
5
Kjelstrup et al. 2011) used NEMD simulation to study the transfer coefficients for liquid-vapor interface of binary mixtures of Lennard-Jones particles with the same diameter and mass, but with different interaction energy. Based on their simulation results for a two-component mixture, they found that the coupling coefficients for heat and mass transfer are negligible because of phase transition or reaction at the dividing surface.
Most of the NEMD simulations have focused on improving understanding the structure of and processes occurring near the liquid-vapor interface of subcritical systems, the study of mixing of a low temperature liquid with a high temperature gas at pressures ranging from subcritical to supercritical values have been largely overlooked. Very recently Mo and Qiao(Mo and Qiao 2017) performed molecular dynamics simulations to explore the transient behavior of a thin film of liquid n-alkanes vaporizing into gaseous nitrogen at subcritical and supercritical pressures. In these simulations, the mass and thickness of the liquid alkane decrease with time as the simulation progresses until the liquid is completely vaporized. These authors obtained the evolutionary profile of a liquid-vapor interface for a thin film of vaporizing liquid alkane exposed to nitrogen at subcritical and supercritical pressures. The authors define the interface characteristics in terms of a dimensionless transition time where the time taken for the liquid to reach supercritical temperatures is normalized by the liquid lifetime. However, it is not entirely clear how dependent these results are on the initial length scale of the liquid slab. In this study, we are interested in investigating the evaporation and mixing of a liquid with a vapor whose thermodynamic properties are held fixed. Specifically, we have developed a nonequilibrium molecular dynamics (NEMD) method to study the mixing process and the interface formed between a liquid phase at (𝑇1, 𝑃1) and vapor phase at (𝑇2, 𝑃2) whose
6
thermodynamic states are held constant through the course of the simulation. Schematic representation of the NEMD method is illustrated in Figure 1. The simulation system consists of a central mixing region, a liquid and a vapor control volume situated at each end of the box. The thermodynamic state of each control volume is held constant during the simulation. Since each box contains a pure liquid or pure vapor, its thermodynamic states can be specified by any two intensive variables, such as [𝑇,𝑃], [𝑇, 𝜌(𝑇,𝑃)] or [𝑇, 𝜇(𝑇,𝑃)]. Here 𝜌(𝑇,𝑃) and 𝜇(𝑇,𝑃) respectively represent the fluid’s density and chemical potential which depend on temperature and pressure.
Figure 1. Schematic representation of the NEMD box used to simulate mixing of a liquid at 𝑇1, 𝑃1 and vapor at 𝑇2, 𝑃2. The thermodynamic state or properties of the liquid and vapor are maintained at the respective targeted temperature and pressure in the control volumes at either end of the simulation cell. The fluids are mixed in the central mixing region.
7
The boundary between each control volume and the central region is permeable to the particles; hence, fluid particles can diffuse from the control volumes into the central region and vice versa. The fluids mix in the central region and can establish a liquid-vapor interface depending on whether the mixing is subcritical or supercritical. Since fluid particles can diffusive out of and into the control volumes, the density of the fluid within a control volume changes with time. To maintain the density of the control volume to its target density (consistent with the specified temperature and pressure), we applied the Grand Canonical Monte Carlo (GCMC) to each control volume which exchanged particles with an implicit bulk fluid. As an example, during the course of simulation, the density 𝜌1(𝑇1,𝑃1) of the liquid phase (control volume 1 in Figure 1) is maintained constant by exchanging (adding or deleting) particles with an implicit reservoir at temperature 𝑇1 and chemical potential 𝜇1. The temperature 𝑇1 of control volume 1 during the MD simulation is set constant by using a thermostat. Similarly, the thermodynamic state [𝑇2, 𝜌2(𝑇2,𝑃2) ] of the vapor phase (control volume 2) was held constant by applying the same procedures. Mixing in the central region was carried out using NVE-MD simulation.
In this study we used hybrid NEMD/GCMC method to simulate mixing of monatomic fluids whose inter-molecular interactions can be represented by the simple Lennard-Jones (LJ) potential energy. Our objectives were to test the efficacy and efficiency of the NEMD/GCMC method on simple fluid systems and gain fundamental insights and understanding of mixing processes under different conditions as a stepping stone to model propulsion fuel systems consisting of complex molecules such as aliphatic and aromatic hydrocarbons. To this end, we have investigated two cases: (i) mixing of liquid Argon (Ar) with its own vapor under subcritical, supercritical and transcritical conditions. This system can also be regarded as an ideal mixture
8
consisting of two chemical species with identical molecular characteristics (size and energy) whose thermodynamic volume and enthalpy change of mixing are zero; (ii) mixing of liquid Krypton (Kr) with vapor Argon (Ar) at subcritical, supercritical and transcritical conditions. This system differs from case (i) in that their molecular sizes and energies, and critical properties are different. The vapor-liquid equilibrium phase diagram of Kr/Ar binary mixture is more complex than that of a pure component Argon in that the former fluid depends on Argon mole fraction, an additional degree of free that does not apply to pure Argon. We would like to note that the results we obtained explicitly for Ar/Ar and Kr/Ar LJ systems could be generalized to other LJ fluids through the dimensionless scale LJ variables 𝑇 ∗ =
𝑘𝐵𝑇 𝜀
,𝑃 ∗ =
𝑃𝜎3 𝜀
𝐸
,𝜌 ∗ = 𝜌𝜎3, 𝐸 ∗ = 𝜀 .
This paper is organized as follows. Section 2 describes the simulation method and details. Sections 3 presents the results for Ar/Ar and Kr/Ar systems. Specifically, we discussed the microscopic structure of the liquid/vapor interface, mixing processes, and the relationship between local density / local temperature variations from the liquid phase to the vapor phase and the thermodynamic vapor-liquid equilibrium of the fluid. These results provide significant insight to the mixing processes under subcritical to supercritical conditions. Finally, Section 4 presents a summary of the results obtained in this study.
2.
Simulation Details
2.1.
Systems Studied The NEMD/GCMC simulations were carried out using the Large-scale Atomic/Molecular
Massively Parallel Simulator (LAMMPS) package(Plimpton 1995). We have performed simulations on six Ar/Ar systems and three Kr/Ar systems. Details on these nine systems are 9
shown in Table 1. The critical temperature (𝑇𝑐) and critical pressure (𝑃𝑐) of Argon are 150.86K and 48.55atm, respectively, and the critical temperature and critical pressure of Kr are 209.4K and 54atm, respectively. Lennard-Jones inter-molecular potential energy (LJ) was used to model the interactions between Krypton/Argon particles. The LJ parameters used are taken from Hirschfelder et al (Talu, Myers et al. 2001) with ɛ/𝑘𝐵 and 𝜎 for Ar given by 119.8K and 3.405 (Å) and ɛ/𝑘𝐵 and 𝜎 for Kr are 166.2K and 3.636 (Å), respectively. The Argon/Argon systems listed in Table 1 were selected to model subcritical (system 1 and 4), transcritical (system 2 and 5), and supercritical (system 3 and 6) mixing. The Krypton/Argon systems listed in this table were selected to model subcritical (system 7), transcritical (system 8), and supercritical (system 9) mixing. Each system is specified by the temperature and pressure of the liquid and vapor phases. Table 1 also displays the density and fugacity coefficient of Argon (=Krypton) that correspond to the specified temperature and pressure. These properties for Argon and Krypton are required in the NEMD/GCMC simulation and must be obtained in advance to set up the simulations. The calculations of these properties (density, number of molecules, and fugacity coefficient) are presented in the Supplementary Information section.
10
Table 1. Details of nine liquid (Argon or Krypton)-vapor (Argon) systems studied Liquid Species
Tri*
T (K)
Pr*Ar
Ar Ar Ar Ar Ar Ar Kr Kr Kr
0.85 0.85 0.85 0.85 0.85 0.85 0.85 0.85 0.85
128 128 128 128 128 128 178 178 178
0.65 1.2 2.3 0.65 1.2 2.3 0.65 1.2 2.3
1 2 3 4 5 6 7 8 9
* 𝑇𝑟𝑖 =
T Tci
P (atm) 31 58 111 31 58 111 31 58 111
Vapor Density (g/cm3) 1.05 1.08 1.13 1.05 1.08 1.13 1.76 1.84 1.90
N 237,692 244,760 255,958 237,692 244,760 255,958 190,365 197,695 205,349
Fugacity Species Tr coeff 0.48 Ar 1.03 0.28 Ar 1.03 0.17 Ar 1.03 0.48 Ar 2.06 0.28 Ar 2.06 0.17 Ar 2.06 0.52 Ar 2.06 0.31 Ar 2.06 0.18 Ar 2.06
P
and 𝑃𝑟𝑖 = Pci
11
T (K)
Pr*Ar
164 164 164 310 310 310 310 310 310
0.55 1.1 2.2 0.55 1.1 2.2 0.55 1.1 2.2
P (atm) 27 54 107 27 54 107 27 54 107
Density (g/cm3) 0.10 0.29 0.72 0.04 0.08 0.17 0.04 0.08 0.17
Fugacity coeff 45,685 0.87 133,436 0.72 330,056 0.49 19,133 0.98 38,583 0.98 77,348 0.97 19,133 0.98 38,583 0.98 77,348 0.97 N
2.2. Implementation of NEMD/GCMC Method Several steps are used to set up the NEMD/GCMC simulation. In the first step (system setup and equilibration), a rectangular box (with periodicity in the x and y directions) of dimensions 10×10×450 nm3 is created (Figure 2-a); the box consists of a vapor and a liquid region separated by a Lennard-Jones wall:
vapor region (10×10×300 nm3)
liquid region (10×10×150 nm3).
(a)
(b) Figure 2. Simulation box of liquid/vapor boxes. (a) equilibration step and (b) NEMD step
To set up the simulation box (Figure 2-a), the number of molecules (𝑁) that corresponds to its density are randomly inserted into the vapor and liquid region. Since the density 𝜌(𝑇,𝑃) for each system at temperature 𝑇 and pressure 𝑃 has been previously determined, the number of atoms can be calculated with following equation: (3)
𝑁 = 𝑉𝜌𝑁𝐴/𝑀𝑤
12
where V (m3) is the volume of liquid-vapor region, ρ (g/m3) is the density, and, Mw (g/mol) represents the molecular weight. The calculated number of argon and krypton molecules for each thermodynamics state are shown in Table 1. Next, energy minimization is performed using the conjugate gradient method to relax the system, followed by equilibration of the system at the target temperature in a 𝑁𝑉𝑇 ensemble for 10 ns (Deaven and Ho 1995). The Langevin thermostat is used to control the temperature(Schneider and Stoll 1978, Dünweg and Paul 1991). Figure 3-a shows the temperature profile along the z-axis of simulation box for system-1 before and after temperature equilibration. By examining the temperature profile for system-1 in figure 3-a it is evident that the system has adequately equilibrated after 10 ns of NVT simulation.
(b)
(a)
Figure 3. (a) Temperature profile along the 𝑧-axis of simulation box for system 1 before and after temperature equilibration; (b) Density profile along the z-axis of simulation box for system 1 before and after density equilibration.
In the next step we equilibrated the density of the vapor and liquid regions by performing NVT/GCMC simulations for the vapor and liquid regions. In the grand canonical Monte Carlo (GCMC) procedure, exchange of molecules with an implicit bulk fluid at the specified temperature
13
and chemical potential (𝜇) or fugacity (𝑇,𝑃) was performed. This is possible using the fugacity coefficient 𝜙(𝑇,𝑃) = 𝑓(𝑇,𝑃)/𝑃 that has been previously calculated and displayed in Table 1. Figure 4-b shows the density profile along the 𝑧-axis of simulation box for system-1 before and after density equilibration. Examination of the density profile in Figure 3-b indicates that the density of the system has adequately equilibrated after 10 ns of NVT/GCMC simulation.
Once the system is equilibrated at the target density and temperature, the wall separating the vapor and liquid regions was removed (Figure 2-b). To implement the NEMD/GCMC procedure to simulate mixing of the liquid and vapor phases, the simulation box was divided into three regions:
vapor control volume (10×10×50 nm3)
mixing region (10×10×350 nm3)
liquid control volume (10×10×50 nm3)
In the liquid and vapor control volumes, both NVT MD (Nose Hoover thermostat) and GCMC procedure were applied to ensure that the temperature and density of the control volumes are maintained at their respective target values. In the mixing region, NVE simulation was performed to model liquid-vapor mixing. The GCMC is applied every 100 MD steps while the GCMC attempts is 1000 with 1000 of MC moves. A MD time step of 2 fs was used. The trajectory file was saved every 20 ps and analyzed to generate density profile, temperature profile, diffusion coefficient, and cluster size statistics.
14
3. Results and Discussion 3.1.
Mixing of Liquid Argon with Vapor Argon In this section, the NEMD simulation results for mixing of liquid Ar with its own vapor
(Ar/Ar) are presented. System-1 and system-4 represent subcritical mixing, system-2 and system-5 are expected to exhibit transcritical mixing, and finally system-3 and system-6 follow supercritical mixing processes. 3.1.1.
Subcritical mixing of liquid Argon with vapor Argon In figure 4-a, the temperature profile and density profile of the red molecules (initially
vapor Ar), blue molecules (initially liquid Ar), and all Ar molecules along the 𝑧-axis of the simulation box are plotted for system-1 at 200 ns. The initial reduced temperature (𝑇𝑟 = 𝑇/𝑇𝑐) and initial reduced pressure (𝑃𝑟 = 𝑃/𝑃𝑐) of liquid Ar for system-1 are 0.85 and 0.65, respectively. The initial 𝑇𝑟 and initial 𝑃𝑟 of vapor Ar are respectively 1.03 and 0.55. The temperature and density of the vapor and liquid control volumes were maintained at their respective target values by using the NVT/GCMC procedure. As seen in this figure, the density profile of all particles decreases gradually away from the liquid control volume towards the liquid/vapor interface. Consistent with subcritical system, a dramatic drop in the particle density is observed at the sharp interface. The figure shows red particles penetrating (“condensation”) into the blue liquid region and blue particles “evaporating” into the vapor phase from the species density profiles. Within a distance of approximately 30 nm from the liquid/vapor interface (figure 4-b), the density and temperature are fairly constant. The variations of density and temperature along the z-axis of the simulation box starting from the vapor control volume to the liquid control volume are plotted on densitytemperature space with the vapor-liquid phase diagram (VLE) of LJ Argon. The figure indicates
15
that the temperature of subcritical interface is approximately 135K and that it is at local phase equilibrium for which saturated liquid and saturated vapor are established on either side of the interface at the prevailing saturation temperature. This result confirms that local phase equilibrium is established across a subcritical interface even in nonequilibrium multiphase system(Røsjorde, Fossmo et al. 2000).
(b)
(a)
(c) Figure 4. (a) Temperature profile and density profile of the red molecules, blue molecules, and all Ar molecules along the z-axis of simulation box in system-1 at 200 ns, (b) Temperature profile and density profile of the red molecules, blue molecules, and all Ar molecules from z = 145 nm to z = 225 nm in system-1 at 200 ns, (c) Variation of temperature-density profile in system-1 at 200 ns.
16
We next examined the evolution of the mixing process for another subcritical system-4. The temperature profile and density profile of the red molecules, blue molecules, and all argon molecules at 200ns along the 𝑧-axis of simulation box are shown in figure 5-a. For this system, the initial Tr and Pr of liquid Ar are the same as those for system-1. The initial pressure 𝑃𝑟 of vapor Ar is 0.55 (same as system-1) while its initial 𝑇𝑟 is 2.06, two times higher than that of system-1. With a higher temperature, the initial density of vapor in system-4 is lower than that of system-1. Figure 5-a shows the temperature profile decreases with 𝑧 with different slopes for the vapor and liquid phases. Despite having a higher vapor control volume temperature, the liquid-vapor interface formed in system-4 is similar to that of system-1. Since the pressure of the simulation box 𝑃 ≈ 28 atm is the same as the pressure for system-1, the temperature of the interface (which is in vapor-liquid phase equilibrium on either side of the interface) must be ~135K (figure 4-c and figure 5-b) which corresponds to the saturation temperature of Ar at a vapor pressure of 28 atm. These results suggest that the interface temperature is governed by the value of the saturation pressure, which is identical for these two systems. This finding enables us to estimate the temperature of a subcritical liquid/vapor interface given the system pressure; this in turn allow us to predict the coefficients of transport properties using kinetic theory since they depend critically on the temperature and pressure of the interface(Røsjorde, Fossmo et al. 2000, Røsjorde, Kjelstrup et al. 2001, Segal and Polikhov 2008).
17
(b)
(a)
Figure 5. (a) Temperature profile and density profile of the red molecules, blue molecules, and all Ar molecules along the z-axis of simulation box in system-4 at 200 ns, (b) Variation of temperaturedensity profile in system-4 at 200 ns.
3.1.2. Transcritical mixing of liquid Argon with vapor Argon We proceed to consider the evolution process for transcritical mixing. The temperature and density profiles for the red particles, blue particles, and all Argon particles along the 𝑧-axis of the simulation box are plotted for system-2 and system-5 at 200 ns in figure 6-a and figure 6-b respectively. For system-2, the initial (𝑇𝑟,𝑃𝑟) for liquid Argon and vapor Argon are (0.85,1.20) and (1.03,1.10) respectively. The initial 𝑇𝑟 and 𝑃𝑟 for the liquid and vapor phases for system-5 are the same as system-2 except that the initial vapor phase temperature 𝑇𝑟 = 2.06, two times higher than that for system-2. The profile for the density of Argon particles (green curve) for both systems show the formation of an interface that is more diffused than that observed for subcritical systems. The figures also show that there exists substantial penetration of the red particles (that is, mixing or “condensation” of red on the liquid side of the interface) into the blue liquid particles compared with subcritical systems. Specifically, a relatively larger number of red particles penetrate the blue
18
particle region in system-2 than those in system-5. Concurrently, a larger number of blue particles diffuse (evaporate) into the red particles (vapor side of the interface) in system-2 than those in system-5. Therefore, it is observed that mixing of liquid-vapor is more efficient in system-2 than system-5 which can be attributed to the fact that the vapor phase in system-2 is closer to the Argon critical point than that of system-5. Furthermore, these figures show that the interface in system-2 is thicker and more diffused than that in system-5.
System-2
(a)
(b)
Figure 6. Temperature profile and density profile of the red molecules, blue molecules, and all Ar molecules along the z-axis of simulation box at 200 ns for (a) system-2 and (b) system-5. It is interesting to explore how the density (𝜌 vs. 𝑧) and temperature (𝑇 vs. 𝑧) profiles of Argon displayed in figure 6 relate to the Argon critical point and vapor-liquid phase diagram. As shown in figure 7-a and 7-b, for these transcritical systems, the density and temperature of the interface traverses the Argon critical region on the ρ-T plane. The temperature of the interface for both systems fluctuates between 145K to 155K around the critical temperature of Argon of ~150.86K indicating that the interface is “unstable” fluctuating in and out of the two-phase dome near the critical point. There exist large density fluctuations (expected in the vicinity of the critical
19
point) and formation of significant number of molecular clusters of varying sizes. Despite the higher temperature of vapor control volume in system-5, the density and temperature of the interface are similar to system-2 (with temperature ranging between 145K and 155 K near the critical point) suggesting that the establishment of the interface to be strongly dependent on the pressure of the systems at 𝑃 = 55 atm, slightly above the critical pressure of Argon. Additionally, cluster analysis was performed to investigate the behavior of interface structure. Using a cut-off distance of 3.81Å to define cluster formation, which corresponds to the location of the first peak of the Ar/Ar radial distribution function, we were able to determine the clustering behavior of the Argon particles. In figure 7-c and figure 7-d, clusters with size exceeding 10 particles in interfacial region are shown for system-2 and system-5. Note that we have omitted the largest “macroscopic” cluster for visual clarity. These figures show that large clusters are formed especially on the vapor side of the interface suggesting that, unlike subcritical evaporation, Argon particles move to the vapor as small particle clusters rather than as individual particles.
20
(a)
(b)
(c)
(d) Z = 150 nm
Z = 250 nm
Figure 7. Variation of temperature-density profile at 200 ns in (a) system-2 and (b) system-5. Clusters of Ar atoms (c) system-2 and (d) system-5 at 200 ns. The macroscopic cluster and clusters having sizes smaller than 10 particles have been omitted for visual clarity.
3.1.3. Supercritical mixing of liquid Argon with vapor Argon We now consider the results for the case of supercritical mixing (system 3 and 6). The temperature and density profiles at 200 ns are depicted in figure 8-a and figure 8-b, respectively. For system-3, the initial (𝑇𝑟,𝑃𝑟) of the vapor and liquid phases are respectively (1.03, 2.2) and (0.85, 2.3) while those for system-6 are (2.06, 2.2) and (0.85, 2.3). For system-3, figure 8-a shows linear Argon density profile and temperature profile along the 𝑧-axis dimension of the simulation
21
box. We have also plotted the density-temperature profile on the 𝜌 ― 𝑇 plane in figure 8-c and 8d. Included in the figure is the Widom line which demarcates the supercritical space of Argon into liquid-like and vapor-like regions near the critical point. It is known that near the critical point of a pure fluid, the supercritical region is not homogenous; the density of a supercritical fluid changes sharply across the Widom line. Figure 8-c indicates that the density-temperature profile for system3 are located entirely in the supercritical “liquid-like” region without a marked change in Argon density. This is not the case for system-6 in which nonlinear temperature and density profiles are observed as shown in figure 8-b resulting in an “interface-like” density profile and the formation of significant number of molecular clusters of varying sizes (figure 8-e). With the temperature and density of Argon plotted on the 𝜌 ― 𝑇 plane, the profile crosses the Widom line transitioning from supercritical liquid-like states to vapor-like states. These results show that even though fluid is in its supercritical state, it can exhibit large density inhomogeneity depending on whether the fluid crosses the Widom line.
22
System-3
System-6
(b)
(a)
(d)
(c)
(e) Z = 175 nm
Z = 275 nm
Figure 8. Temperature and density profiles of the red molecules, blue molecules, and all Ar molecules along the z-axis of simulation box at 200 ns for (a) system-3 and (b) system-6. Variation of temperature-density profile at 200 ns in (c) system-3 and (d) system-6. Clusters of Ar atoms (e) system-6 at 200 ns. The macroscopic cluster and clusters having sizes smaller than 10 particles have been omitted for visual clarity.
23
3.1.4. Voronoi Tessellation Analysis We have calculated the Voronoi tessellation volume of Argon particles for system-1 (subcritical), system-2 (transcritical) and system-3 (supercritical) at 200 ns. Voronoi tessellation cells may be regarded as a measure of the “molecular volume” of each LJ particle. For subcritical mixing, we anticipate a well-defined sharp change in the volume size of Voronoi cells across the liquid-vapor interface. A gradual change of Voronoi volumes is expected for supercritical mixing. Transcritical mixing represents a domain when both subcritical and supercritical mixing types are present. Figure 9-a shows the Voronoi cell volume or “molecular volume” of Ar molecules for system-1 (subcritical system). The magnitude of Argon volume changes dramatically at the liquidvapor interface which indicate a well-defined interface. In the case of a transcritical system-2 (figure 9-b), superposition of a sharp interface to a diffusive interface is observed. While large variation or change in Voronoi volume is observed, the interface is smoothed due to the decrease in the surface tension. Particle clusters are formed at the interface which resulting in a wavy surface. Finally, with increase in the pressure of the system (figure 19-c), the gradient in volume of Ar molecules decreases monotonically, and fluid/fluid mixing (supercritical mixing) is observed. Similar observations have been reported in previous experimental study by Segal and Polikhov(Segal and Polikhov 2008).
24
750 (Å3) 60 (Å3)
(a)
3
300 (Å ) 60 (Å3) (b)
3
100 (Å ) 60 (Å3) (c) Z = 250 nm
Z = 150 nm
Figure 9. Ar molecular volume (a) subcritical system-1, (b) transcritical system-2, and (c) supercritical system-3 at 200 ns.
3.2.
Liquid Krypton/Vapor Argon System Although the Ar/Ar system described in Section 3.1 is a one-component system, it can be
regarded as an ideal binary mixture consisting of identical molecular properties (molecular size and interaction energy). Despite its simplicity, the results obtained have provided significant insights and fundamental understanding on the evolution of fluid/fluid interface and the influence of thermodynamics and phase equilibria on subcritical, supercritical, and transcritical mixing. We have also demonstrated that the NEMD/GCMC method can generate self-establishing liquid/vapor interfaces and model mixing of a liquid with a gas at different temperature and pressure. We will apply this new NEMD procedure to study mixing of liquid Krypton (Kr) with vapor Argon (Ar) – a binary mixture of monatomic species with different molecular sizes and interaction energies.
25
This study investigates the effects of dissimilar molecular properties on the evolution of mixing processes at subcritical, transcritical, and supercritical conditions.
3.2.1. Subcritical mixing of liquid Krypton with vapor Argon In figure 10, the temperature profile and density profile of the Ar molecules (initially vapor Ar), Kr molecules (initially liquid Kr) and Ar/Kr molecules along the z-axis of simulation box are plotted for system-7 at 200 ns. The initial temperature and initial pressure of liquid Kr are 178K and 31atm, respectively. The initial temperature and initial pressure of vapor Ar are respectively 310K and 27atm. The temperature and pressure of the vapor and liquid control volumes are maintained at the target values by using NVT/GCMC simulation. Similar to system-1 and system-4, the density profile of Ar/Kr particle decreases slightly outside the liquid control volume until the liquid-vapor interface. A dramatic drop in particle density is observed at the interface. Consequently, a sharp interface is established for subcritical liquid-vapor mixing. The figure shows that the Ar particles penetrate (“condense”) into the Kr liquid particles and the Kr particles “evaporate” into the vapor phase. Within a distance of approximately 10 nm away from the liquidvapor interface, the density and temperature are fairly constant.
26
Figure 10. Temperature profile and density profile of the Ar molecules, Kr molecules, and Ar/Kr molecules along the z-axis of simulation box in system-7 at 200 ns. Figure 11-a shows the variation of density-temperature-Ar mole fraction profile from vapor control volume to the liquid control volume (at 150 ns, 175 ns and, 200 ns) and the vapor-liquid phase diagram (VLE) of LJ Kr/Ar. The figure shows when these profiles are plotted in thermodynamic space, they collapse onto the same curve indicating that they are thermodynamically identical. To examine the thermodynamic equilibrium state of mixture at the interface for system-7, the VLE plot of Kr/Ar binary system at 170 K (interface temperature, see figure 10) is shown in Figure 11b. The figure shows that the vapor-liquid compositions of interface at 150 ns, 175 ns, and 200 ns are located on equilibrium line at pressure of 28 bar. This pressure is the average pressure of simulation box which suggests that, similar to the Ar/Ar system, the subcritical interface exhibits local phase equilibrium in which saturated liquid and saturated vapor are established on either side of the interface at the prevailing saturation temperature. This verifies that an equilibrium state has been achieved and that the equilibrated interface is being convected downstream maintaining the equilibrium composition on either side of the interface without any
27
change. The ability of the NEMD procedure to achieve an equilibrium state for a binary fluid system is a fundamental validation of the NEMD framework.
(b)
(a)
Figure 11. (a) Variation of temperature-density-Ar mole fraction profile in system-7, (b) the vaporliquid phase diagram (VLE) of the LJ Krypton/Argon at 170 K.
3.2.2. Transcritical mixing of liquid Krypton with vapor Argon The temperature profile and density profile of the Ar molecules (initially vapor Ar), Kr molecules (initially liquid Kr), and Ar/Kr molecules along the z-axis of simulation box are plotted for system-8 at 200 ns in figure 12-a. The initial temperature and initial pressure of liquid Kr for system-8 were 178 K and 58 atm, respectively. The initial temperature and initial pressure of vapor Ar were respectively 310 K and 54 atm. The density profile for Ar/Kr particles in system-8 shows the existence of an interface that is more diffused than that observed for subcritical system. The figure also shows that there exists substantial penetration of Ar particles (that is, mixing or “condensation” of Ar and Kr particles on the liquid side of the interface) into the Kr liquid particles
28
compared with subcritical systems. To investigate the phase behavior of interface in this system, the variation of density-temperature-mole fraction profile from vapor control volume to the liquid control volume at 200 ns and the vapor-liquid phase diagram (VLE) of the LJ Kr/Ar are shown in figure 12-b. This figure shows the density-temperature-mole fraction profile of transcritical systems traversing the critical region on the ρ-T-x space in interfacial region as the interface evolves in time. This behavior suggests large density fluctuations (expected in the vicinity of the critical point) and formation of significant number of molecular clusters of varying sizes.
(b)
(a)
Figure 12. (a) Temperature profile and density profile of the Ar molecules, Kr molecules, and Kr/Ar molecules along the z-axis of simulation box in system-8 at 200 ns, (b) Variation of temperature-density-mole fraction profile in system-8.
3.2.3. Supercritical mixing of liquid Krypton with vapor Argon In the NEMD simulation of the supercritical system, Krypton and Argon mix to establish a liquid-vapor interface as illustrated in figure 13-a. In figure 13-a, the temperature profile and density profile of the Ar molecules (initially vapor Ar), Kr molecules (initially liquid Kr), and total density profile of Kr/Ar along the z-axis are plotted at 200 ns. The initial temperature and initial 29
pressure of liquid Kr for system-9 were 178 K and 111 atm, respectively. The initial temperature and initial pressure of vapor Ar were respectively 310 K and 107 atm. The temperature and pressure of the vapor and liquid control volumes were maintained at the target values by using NVT/GCMC.
(b)
(a)
Figure 13. (a) Temperature profile and density profile of the red molecules, blue molecules, and all molecules along the z-axis of simulation box in system-9 at 200 ns, (b) Variation of temperature-density-mole fraction profile in system-9.
The evolution of the temperature and densities of Kr and Ar are captured in the NEMD simulation. Figure 13-a shows nonlinear temperature and density profiles are observed after 200 ns. The density profiles reveal diffusion-dominated mixing with Ar particles (red curve) penetrating deep inside the liquid region and Kr particles (blue curve) into the vapor region. It is interesting to note significant variation in total density (green curve) resulting in an “interface-like” density profile. Figure 13-b shows the variation of density-temperature-mole fraction profile from vapor control volume to the liquid control volume at 200 ns and the vapor-liquid phase diagram (VLE) of the LJ 30
Kr/Ar. Our NEMD result shows that the entire curve lies outside of the Kr/Ar vapor-liquid equilibrium saturation dome. Also, the temperature and density profiles cross the Widom line which differentiates gas-like and liquid-like states of supercritical Argon/Krypton mixture. We did not compute and display the Widom line, which is known to exist, for this mixture. For these supercritical profiles we note that the interface will continue to diffuse and spread in time and eventually a “steady state” profile will be obtained across the two control volumes that are maintained at fixed conditions.
4. Conclusion We have developed a NEMD simulation method to model subcritical, transcritical, and supercritical mixing of liquid and gas at high pressures. This method ensures that the temperature, pressure, and density of the liquid and vapor phases are held constant while mixing takes place between them. It can model a self-establishing liquid/vapor interface depending on whether the mixing is subcritical, supercritical, or transcritical. The procedure generates thermodynamic temperature, density, and mole fraction profiles as a function of spatial position and time and produces information at the microscopic level including the structure of liquid/vapor interface, molecular clustering statistics, molecular volumes, etc. The simulation method was applied to mixing of (i) liquid Argon with its own vapor and (ii) liquid Krypton with vapor Argon at pressures ranging from subcritical to supercritical values. Our simulations demonstrate that subcritical interface exhibits local phase equilibrium in which saturated liquid and saturated vapor are established on either side of the interface at the prevailing saturation temperature that is fixed by the pressure of the system. For transcritical systems, the interface is more diffused than that of subcritical systems. The interface consists of large particle clusters and its density-temperature
31
profile traverses the critical region. Similar to subcritical systems, irrespective of the temperatures of the liquid and vapor phases, the interface temperature is determined by the pressure of system. Finally, the results for supercritical systems show diffusion-controlled fluid mixing without the formation of an interface. However, the density-temperature variation cuts across the Widom line, that divides the supercritical region into “liquid-like” and “vapor-like” supercritical states, resulting in large density variation the resembles an “interface-like” profile.
To summarize, key contributions of this study include: (a) development and demonstration of a hybrid NEMD/GCMC method for simulating mixing of liquid and vapor at different temperature and pressure (boundary conditions); (b) this unique procedure produces for the firsttime trajectories of the liquid-vapor mixture that can be mapped onto the thermodynamic space (temperature-density-composition); and (c) it generates microscopic and molecular level structural information on mixing at different mixing conditions. Even though the Ar/Ar and Kr/Ar systems studied are simple monatomic molecular systems, the results obtained have provided significant insights and fundamental knowledge and information on mixing at subcritical, supercritical, and transcritical conditions that are not possible from macroscopic continuum-based models. The proposed NEMD method is sufficiently robust and versatile, and can be employed to study mixing of fluid mixtures consisting of realistic hydrocarbon fuel with combustible gases. We are currently applying the NEMD method to study mixing of propulsion fuels.
32
Acknowledgements This work is being performed under Air Force Phase II STTR contract FA9550-17-C-0015 with Dr. Malissa Lightfoot as the technical monitor. We gratefully acknowledge support from AFOSR for this STTR program. We would also like to thank Drs. Chiping Li, Malissa Lightfoot, Jeff Mills, and Venke Sankaran for their technical feedback during the course of this work.
33
Reference Bedeaux, D. and Kjelstrup, S. (2004). "Irreversible thermodynamics—a tool to describe phase transitions far from global equilibrium." Chemical Engineering Science 59(1): 109-118. Bedeaux, D. and Kjelstrup, S. (2018). "Fluid-Fluid Interfaces of Multi-Component Mixtures in Local Equilibrium." Entropy 20(4): 250. Bellan, J. (2000). "Supercritical (and subcritical) fluid behavior and modeling: drops, streams, shear and mixing layers, jets and sprays." Progress in Energy and Combustion Science 26(4-6): 329-366. Chehroudi, B. (2012). "Recent experimental efforts on high-pressure supercritical injection for liquid rockets and their implications." International Journal of Aerospace Engineering 2012. Chen, L., Zhang, X. R., Okajima, J. and Maruyama, S. (2013). "Numerical investigation of nearcritical fluid convective flow mixing in microchannels." Chemical Engineering Science 97: 6780. Dahms, R. N. and Oefelein, J. C. (2013). "On the transition between two-phase and single-phase interface dynamics in multicomponent fluids at supercritical pressures." Physics of Fluids 25(9): 092103. Deaven, D. M. and Ho, K. M. (1995). "Molecular geometry optimization with a genetic algorithm." Physical Review Letters 75(2): 288. Dünweg, B. and Paul, W. (1991). "Brownian dynamics simulations without Gaussian random numbers." International Journal of Modern Physics C 2(03): 817-827. Givler, S. D., Abraham, J. (1996). "Supercritical droplet vaporization and combustion studies." Progress in Energy and Combustion Science 22(1): 1-28.
34
H. Mayer, W. O., Schik, A. H. A., Vielle, B., Chauveau, C., Gokalp, I., Talley, D. G. and Woodward, R. D. (1998). "Atomization and breakup of cryogenic propellants under highpressure subcritical and supercritical conditions." Journal of Propulsion and Power 14(5): 835842. Huai, X., S. Koyama and Zhao, T. (2005). "An experimental study of flow and heat transfer of supercritical carbon dioxide in multi-port mini channels under cooling conditions." Chemical Engineering Science 60(12): 3337-3345. Inzoli, I., Kjelstrup, S., Bedeaux, D. and Simon, J. (2011). "Transfer coefficients for the liquid– vapor interface of a two-component mixture." Chemical engineering science 66(20): 4533-4548. Kjelstrup, S. and Hafskjold, B. (1996). "Nonequilibrium molecular dynamics simulations of steady-state heat and mass transport in distillation." Industrial & engineering chemistry research 35(11): 4203-4213. Meland, R. (2002). "Molecular exchange and its influence on the condensation coefficient." The Journal of chemical physics 117(15): 7254-7258. Meland, R. (2003). "Molecular dynamics simulation of the inverted temperature gradient phenomenon." Physics of Fluids 15(10): 3244-3247. Meland, R., Frezzotti, A., Ytrehus, T. and Hafskjold, B. (2004). "Nonequilibrium moleculardynamics simulation of net evaporation and net condensation, and evaluation of the gas-kinetic boundary condition at the interphase." Physics of Fluids 16(2): 223-243. Mo, G. and Qiao, L. (2017). "A molecular dynamics investigation of n-alkanes vaporizing into nitrogen: transition from subcritical to supercritical." Combustion and Flame 176: 60-71.
35
Olivier, M. L., Rollier, J. D. and Kjelstrup, S. (2002). "Equilibrium properties and surface transfer coefficients from molecular dynamics simulations of two-component fluids." Colloids and Surfaces A: Physicochemical and Engineering Aspects 210(2-3): 199-222. Oschwald, M., Smith, J., Branam, R., Hussong, J., Schik, A., Chehroudi, B. and Talley, D. (2006). "Injection of fluids into supercritical environments." Combustion Science and Technology 178(1-3): 49-100. Plimpton, S. (1995). "Fast parallel algorithms for short-range molecular dynamics." Journal of computational physics 117(1): 1-19. Røsjorde, A., Fossmo, D., Bedeaux, D., Kjelstrup, S. and Hafskjold, B. (2000). "Nonequilibrium molecular dynamics simulations of steady-state heat and mass transport in condensation: I. Local equilibrium." Journal of colloid and interface science 232(1): 178-185. Røsjorde, A., Kjelstrup, S., Bedeaux, D. and Hafskjold, B. (2001). "Nonequilibrium molecular dynamics simulations of steady-state heat and mass transport in condensation. II. Transfer coefficients." Journal of colloid and interface science 240(1): 355-364. Savin, T., Schweizer, M. and Öttinger, H. C. (2016). Nonequilibrium thermodynamics of an interface. APS Meeting Abstracts. Schneider, T. and Stoll, E. (1978). "Molecular-dynamics study of a three-dimensional onecomponent model for distortive phase transitions." Physical Review B 17(3): 1302. Schweizer, M., Öttinger, H. C. and Savin, T. (2016). "Nonequilibrium thermodynamics of an interface." Physical Review E 93(5): 052803. Schweizer, M., Öttinger, H. C. and Savin, T. (2016). "Nonequilibrium thermodynamics of an interface." Physical Review E 93(5): 052803.
36
Segal, C. and Polikhov, S. (2008). "Subcritical to supercritical mixing." Physics of Fluids 20(5): 052101. Simon, J. M., Bedeaux, D., Kjelstrup, S., Xu, J. and Johannessen, E. (2006). "Interface film resistivities for heat and mass transfers integral relations verified by non-equilibrium molecular dynamics." The Journal of Physical Chemistry B 110(37): 18528-18536. Simon, J. M., Kjelstrup, S., Bedeaux, D. and Hafskjold, B. (2004). "Thermal flux through a surface of n-octane. A non-equilibrium molecular dynamics study." The Journal of Physical Chemistry B 108(22): 7186-7195. Talu, O. and Myers, A. L. (2001). "Reference potentials for adsorption of helium, argon, methane, and krypton in high-silica zeolites." Colloids and Surfaces A: Physicochemical and Engineering Aspects 187: 83-93. Tsuruta, T. and Nagayama, G. (2005). "A microscopic formulation of condensation coefficient and interface transport phenomena." Energy 30(6): 795-805. Wilhelmsen, Ø., Trinh, T. T., Kjelstrup, S. and D. Bedeaux (2015). "Influence of curvature on the transfer coefficients for evaporation and condensation of Lennard-Jones fluid from squaregradient theory and nonequilibrium molecular dynamics." The Journal of Physical Chemistry C 119(15): 8160-8173. Xu, J., Kjelstrup, S., Bedeaux, D., Røsjorde, A. and L. Rekvig (2006). "Verification of Onsager's reciprocal relations for evaporation and condensation using non-equilibrium molecular dynamics." Journal of colloid and interface science 299(1): 452-463. Yang, V. (2000). "Modeling of supercritical vaporization, mixing, and combustion processes in liquid-fueled propulsion systems." Proceedings of the Combustion Institute 28(1): 925-942.
37
Farzin Rahmani: Methodology, Software, Writing - Original Draft, Investigation Timothy Weathers: Writing - Review & Editing, Investigation Ashvin Hosangadi: Conceptualization, Writing - Review & Editing Yee C. Chiew: Conceptualization, Methodology, Writing - Original Draft, Supervision,
Declaration of interests ☒ 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. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
38
Highlights: A non-equilibrium molecular dynamics (NEMD)/ Grand Canonical Monte Carlo (GCMC) method is employed to simulate subcritical, transcritical and, supercritical mixing. Local phase equilibrium is established at subcritical interface. Local density and local temperature of the transcritical interface can be directly mapped to the near critical region of the Ar/Ar and Kr/Ar vapor-liquid phase diagram. For transcritical mixing, superposition of a sharp interface to a diffusive interface is
39
observed. For transcritical mixing, significant clustering of molecules is found at the “transitional diffuse” interface.
40
S0009-2509(19)30914-5 https://doi.org/10.1016/j.ces.2019.115424 CES 115424
To appear in:
Chemical Engineering Science
Received Date: Revised Date: Accepted Date:
12 September 2019 22 November 2019 9 December 2019
Please cite this article as: F. Rahmani, T. Weathers, A. Hosangadi, Y.C. Chiew, A Non-Equilibrium Molecular Dynamics Study of Subcritical, Supercritical and Transcritical Mixing of Liquid-Gas Systems, Chemical Engineering Science (2019), doi: https://doi.org/10.1016/j.ces.2019.115424
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2019 Published by Elsevier Ltd.
A Non-Equilibrium Molecular Dynamics Study of Subcritical, Supercritical and Transcritical Mixing of Liquid-Gas Systems
Farzin Rahmani(a), Timothy Weathers(b), Ashvin Hosangadi(b), Yee C. Chiew*(a), (a)
Department of Chemical & Biochemical Engineering, Rutgers University – New Brunswick, Piscataway NJ 08854 (b) Combustion Research and Flow Technology, Inc., Pipersville, PA 18947
* Corresponding author. Email: [email protected]
1
Abstract A non-equilibrium molecular dynamics (NEMD) simulation method has been developed to simulate mixing of a liquid with a vapor and identify the characteristics of the liquid-vapor interface in the mixing process. Our results show that for the case of subcritical mixing, local phase equilibrium is established, and saturated liquid and saturated vapor are formed on either side of the interface at the prevailing saturation temperature that is fixed by the pressure of the system independent of the temperature difference across the interface. For transcritical mixing, significant clustering of molecules is found at the “transitional diffuse” interface. We found that the local density and local temperature of the transcritical interface can be directly mapped to the near critical region of the Ar/Ar and Kr/Ar vapor-liquid phase diagram. Finally, for the case of supercritical mixing, our simulations show a gradual change of density consistent with diffusioncontrolled fluid mixing. Keywords: Non-equilibrium molecular dynamics simulation, Subcritical Mixing, Transcritical Mixing, Supercritical Mixing, Liquid-vapor Interface
2
1. Introduction Understanding mixing and evaporation of an injected liquid into a high-pressure gaseous environment is of relevance and importance to a wide range of technological applications such as propulsion systems including liquid rocket engines, gas turbine combustors, and diesel engines. Depending on the ambient pressure relative to the critical pressure of the liquid, mixing of the injected liquid can exhibit different mixing processes: subcritical, supercritical, and transcritical mixing mechanisms(Givler, Abraham et al. 1996, Bellan and science 2000, Yang 2000). Highpressure experimental studies(H. Mayer, A. Schik et al. 1998, Segal and Polikhov 2008, Chehroudi 2012) have provided insights on the behavior of mixing of multi-phase systems. When liquid fuel is injected into a chamber with pressure below the critical pressure of the fuel, classical jet breakup and droplet formation are observed. A well-defined molecular interface separates the injected liquid from ambient gases due to surface tension, and interactions between shear and surface tension forces promote droplet formation and secondary jet breakup. In contrast, at pressures exceeding the liquid critical pressure, the mixing process can become quite different. Under these conditions, the liquid-gas interface diminishes and is replaced by one-phase dense mixing layer which is affected by thermodynamics and transport processes.
Transcritical mixing that transitions from two-phase mixing regime to one-phase supercritical mixing has been observed in several studies.(Huai, Koyama et al. 2005, Chen, Zhang et al. 2013) Segal and Polikhov(Segal and Polikhov 2008) reported that transcritical mixing is characterized by droplet formation and irregularly shaped material break-up from the injectant. Detailed analysis suggests that the transitional changes are controlled by the multi-component
3
nature of the thermodynamics, transport processes and stability of gas-liquid interface(Oschwald*, Smith et al. 2006, Chehroudi 2012). Computational fluid dynamics and continuum-based models have been developed to model supercritical mixing. Recent theoretical and modeling studies (Mo and Qiao 2017) have employed real gas equations of state to quantify the thermodynamic properties and critical point loci of binary and multi-component fluid mixture within an equilibrium framework. Dahms and Oefelein (Dahms and Oefelein 2013) proposed a criterion, based on the thickness of equilibrium liquid-gas interface, to characterize an interface as subcritical or supercritical in behavior. Although these studies have provided insights on transcritical mixing, fundamental understanding of this complex process is still lacking due to considerable challenges involved in characterizing the thermodynamics properties and vapor-liquid equilibrium of multi-component mixtures especially in the vicinity of the mixture critical point. Moreover, numerical macroscopic models are not able to capture the fine scales needed to resolve the microscopic structure of the liquid-gas interface that in part governs the underlying processes subject to large density and temperature gradients.
Nonequilibrium molecular dynamics (NEMD) simulation(Meland 2002, Meland 2003, Bedeaux and Kjelstrup 2004, Meland, Frezzotti et al. 2004, Simon, Bedeaux et al. 2006, Xu, Kjelstrup et al. 2006) is an invaluable tool that complements macroscopic continuum methods for probing the small length scale and short time scale phenomena associated with the physical changes of mixing interface in response to variation of temperature and pressure. In contrast to continuum-based methods where a priori assumptions about the process are necessary, molecular dynamics simulations do not require assumed models for equations of state, equilibrium liquid/vapor interface, thermodynamic and transport properties. It requires only inter-molecular
4
potential energy function as input and generates the evolutionary process of mixing during the simulation. Such a study can elucidate the physics of the mixing mechanisms and test the assumptions used in different modeling approaches. Additionally, it provides information at the microscopic molecular scale, e.g., clustering of molecules and particle correlation, that are not possible from macroscopic models. Molecular dynamics simulation is computationally expensive, and the system studied is limited to several hundred to thousand nanometers and over a time duration of several hundred nanoseconds.
A number of researchers have used molecular dynamics to model the interface of onecomponent subcritical liquid-gas systems at equilibrium and nonequilibrium situations. They showed that a key parameter that controls the interface characteristic is the surface tension which is governed by the temperature at the interface. The aim of these NEMD studies was to verify the validity of proposed theories for predicting the temperature, transport properties, condensation and evaporation rates in the vicinity of interface(Kjelstrup and Hafskjold 1996, Røsjorde, Kjelstrup et al. 2001, Simon, Kjelstrup et al. 2004, Tsuruta and Nagayama 2005, Xu, Kjelstrup et al. 2006, Wilhelmsen, Trinh et al. 2015, Savin, Schweizer et al. 2016, Schweizer, Öttinger et al. 2016, Schweizer, Öttinger et al. 2016). For example, Rosjorde et al.(Røsjorde, Fossmo et al. 2000, Røsjorde, Kjelstrup et al. 2001) determined the coefficients for transfer of heat and mass across the liquid-vapor interface of one-component fluid by NEMD simulation. They used the SoaveRedlich-Kwong equation of state to test the assumption of local equilibrium and showed that local equilibrium is established not only in the bulk phases but also at the liquid-vapor interface. In recent years, NEMD simulations have been carried out to study the liquid-vapor interface of twocomponent mixtures(Olivier, Rollier et al. 2002, Bedeaux and Kjelstrup 2018). Inzoli et al.(Inzoli,
5
Kjelstrup et al. 2011) used NEMD simulation to study the transfer coefficients for liquid-vapor interface of binary mixtures of Lennard-Jones particles with the same diameter and mass, but with different interaction energy. Based on their simulation results for a two-component mixture, they found that the coupling coefficients for heat and mass transfer are negligible because of phase transition or reaction at the dividing surface.
Most of the NEMD simulations have focused on improving understanding the structure of and processes occurring near the liquid-vapor interface of subcritical systems, the study of mixing of a low temperature liquid with a high temperature gas at pressures ranging from subcritical to supercritical values have been largely overlooked. Very recently Mo and Qiao(Mo and Qiao 2017) performed molecular dynamics simulations to explore the transient behavior of a thin film of liquid n-alkanes vaporizing into gaseous nitrogen at subcritical and supercritical pressures. In these simulations, the mass and thickness of the liquid alkane decrease with time as the simulation progresses until the liquid is completely vaporized. These authors obtained the evolutionary profile of a liquid-vapor interface for a thin film of vaporizing liquid alkane exposed to nitrogen at subcritical and supercritical pressures. The authors define the interface characteristics in terms of a dimensionless transition time where the time taken for the liquid to reach supercritical temperatures is normalized by the liquid lifetime. However, it is not entirely clear how dependent these results are on the initial length scale of the liquid slab. In this study, we are interested in investigating the evaporation and mixing of a liquid with a vapor whose thermodynamic properties are held fixed. Specifically, we have developed a nonequilibrium molecular dynamics (NEMD) method to study the mixing process and the interface formed between a liquid phase at (𝑇1, 𝑃1) and vapor phase at (𝑇2, 𝑃2) whose
6
thermodynamic states are held constant through the course of the simulation. Schematic representation of the NEMD method is illustrated in Figure 1. The simulation system consists of a central mixing region, a liquid and a vapor control volume situated at each end of the box. The thermodynamic state of each control volume is held constant during the simulation. Since each box contains a pure liquid or pure vapor, its thermodynamic states can be specified by any two intensive variables, such as [𝑇,𝑃], [𝑇, 𝜌(𝑇,𝑃)] or [𝑇, 𝜇(𝑇,𝑃)]. Here 𝜌(𝑇,𝑃) and 𝜇(𝑇,𝑃) respectively represent the fluid’s density and chemical potential which depend on temperature and pressure.
Figure 1. Schematic representation of the NEMD box used to simulate mixing of a liquid at 𝑇1, 𝑃1 and vapor at 𝑇2, 𝑃2. The thermodynamic state or properties of the liquid and vapor are maintained at the respective targeted temperature and pressure in the control volumes at either end of the simulation cell. The fluids are mixed in the central mixing region.
7
The boundary between each control volume and the central region is permeable to the particles; hence, fluid particles can diffuse from the control volumes into the central region and vice versa. The fluids mix in the central region and can establish a liquid-vapor interface depending on whether the mixing is subcritical or supercritical. Since fluid particles can diffusive out of and into the control volumes, the density of the fluid within a control volume changes with time. To maintain the density of the control volume to its target density (consistent with the specified temperature and pressure), we applied the Grand Canonical Monte Carlo (GCMC) to each control volume which exchanged particles with an implicit bulk fluid. As an example, during the course of simulation, the density 𝜌1(𝑇1,𝑃1) of the liquid phase (control volume 1 in Figure 1) is maintained constant by exchanging (adding or deleting) particles with an implicit reservoir at temperature 𝑇1 and chemical potential 𝜇1. The temperature 𝑇1 of control volume 1 during the MD simulation is set constant by using a thermostat. Similarly, the thermodynamic state [𝑇2, 𝜌2(𝑇2,𝑃2) ] of the vapor phase (control volume 2) was held constant by applying the same procedures. Mixing in the central region was carried out using NVE-MD simulation.
In this study we used hybrid NEMD/GCMC method to simulate mixing of monatomic fluids whose inter-molecular interactions can be represented by the simple Lennard-Jones (LJ) potential energy. Our objectives were to test the efficacy and efficiency of the NEMD/GCMC method on simple fluid systems and gain fundamental insights and understanding of mixing processes under different conditions as a stepping stone to model propulsion fuel systems consisting of complex molecules such as aliphatic and aromatic hydrocarbons. To this end, we have investigated two cases: (i) mixing of liquid Argon (Ar) with its own vapor under subcritical, supercritical and transcritical conditions. This system can also be regarded as an ideal mixture
8
consisting of two chemical species with identical molecular characteristics (size and energy) whose thermodynamic volume and enthalpy change of mixing are zero; (ii) mixing of liquid Krypton (Kr) with vapor Argon (Ar) at subcritical, supercritical and transcritical conditions. This system differs from case (i) in that their molecular sizes and energies, and critical properties are different. The vapor-liquid equilibrium phase diagram of Kr/Ar binary mixture is more complex than that of a pure component Argon in that the former fluid depends on Argon mole fraction, an additional degree of free that does not apply to pure Argon. We would like to note that the results we obtained explicitly for Ar/Ar and Kr/Ar LJ systems could be generalized to other LJ fluids through the dimensionless scale LJ variables 𝑇 ∗ =
𝑘𝐵𝑇 𝜀
,𝑃 ∗ =
𝑃𝜎3 𝜀
𝐸
,𝜌 ∗ = 𝜌𝜎3, 𝐸 ∗ = 𝜀 .
This paper is organized as follows. Section 2 describes the simulation method and details. Sections 3 presents the results for Ar/Ar and Kr/Ar systems. Specifically, we discussed the microscopic structure of the liquid/vapor interface, mixing processes, and the relationship between local density / local temperature variations from the liquid phase to the vapor phase and the thermodynamic vapor-liquid equilibrium of the fluid. These results provide significant insight to the mixing processes under subcritical to supercritical conditions. Finally, Section 4 presents a summary of the results obtained in this study.
2.
Simulation Details
2.1.
Systems Studied The NEMD/GCMC simulations were carried out using the Large-scale Atomic/Molecular
Massively Parallel Simulator (LAMMPS) package(Plimpton 1995). We have performed simulations on six Ar/Ar systems and three Kr/Ar systems. Details on these nine systems are 9
shown in Table 1. The critical temperature (𝑇𝑐) and critical pressure (𝑃𝑐) of Argon are 150.86K and 48.55atm, respectively, and the critical temperature and critical pressure of Kr are 209.4K and 54atm, respectively. Lennard-Jones inter-molecular potential energy (LJ) was used to model the interactions between Krypton/Argon particles. The LJ parameters used are taken from Hirschfelder et al (Talu, Myers et al. 2001) with ɛ/𝑘𝐵 and 𝜎 for Ar given by 119.8K and 3.405 (Å) and ɛ/𝑘𝐵 and 𝜎 for Kr are 166.2K and 3.636 (Å), respectively. The Argon/Argon systems listed in Table 1 were selected to model subcritical (system 1 and 4), transcritical (system 2 and 5), and supercritical (system 3 and 6) mixing. The Krypton/Argon systems listed in this table were selected to model subcritical (system 7), transcritical (system 8), and supercritical (system 9) mixing. Each system is specified by the temperature and pressure of the liquid and vapor phases. Table 1 also displays the density and fugacity coefficient of Argon (=Krypton) that correspond to the specified temperature and pressure. These properties for Argon and Krypton are required in the NEMD/GCMC simulation and must be obtained in advance to set up the simulations. The calculations of these properties (density, number of molecules, and fugacity coefficient) are presented in the Supplementary Information section.
10
Table 1. Details of nine liquid (Argon or Krypton)-vapor (Argon) systems studied Liquid Species
Tri*
T (K)
Pr*Ar
Ar Ar Ar Ar Ar Ar Kr Kr Kr
0.85 0.85 0.85 0.85 0.85 0.85 0.85 0.85 0.85
128 128 128 128 128 128 178 178 178
0.65 1.2 2.3 0.65 1.2 2.3 0.65 1.2 2.3
1 2 3 4 5 6 7 8 9
* 𝑇𝑟𝑖 =
T Tci
P (atm) 31 58 111 31 58 111 31 58 111
Vapor Density (g/cm3) 1.05 1.08 1.13 1.05 1.08 1.13 1.76 1.84 1.90
N 237,692 244,760 255,958 237,692 244,760 255,958 190,365 197,695 205,349
Fugacity Species Tr coeff 0.48 Ar 1.03 0.28 Ar 1.03 0.17 Ar 1.03 0.48 Ar 2.06 0.28 Ar 2.06 0.17 Ar 2.06 0.52 Ar 2.06 0.31 Ar 2.06 0.18 Ar 2.06
P
and 𝑃𝑟𝑖 = Pci
11
T (K)
Pr*Ar
164 164 164 310 310 310 310 310 310
0.55 1.1 2.2 0.55 1.1 2.2 0.55 1.1 2.2
P (atm) 27 54 107 27 54 107 27 54 107
Density (g/cm3) 0.10 0.29 0.72 0.04 0.08 0.17 0.04 0.08 0.17
Fugacity coeff 45,685 0.87 133,436 0.72 330,056 0.49 19,133 0.98 38,583 0.98 77,348 0.97 19,133 0.98 38,583 0.98 77,348 0.97 N
2.2. Implementation of NEMD/GCMC Method Several steps are used to set up the NEMD/GCMC simulation. In the first step (system setup and equilibration), a rectangular box (with periodicity in the x and y directions) of dimensions 10×10×450 nm3 is created (Figure 2-a); the box consists of a vapor and a liquid region separated by a Lennard-Jones wall:
vapor region (10×10×300 nm3)
liquid region (10×10×150 nm3).
(a)
(b) Figure 2. Simulation box of liquid/vapor boxes. (a) equilibration step and (b) NEMD step
To set up the simulation box (Figure 2-a), the number of molecules (𝑁) that corresponds to its density are randomly inserted into the vapor and liquid region. Since the density 𝜌(𝑇,𝑃) for each system at temperature 𝑇 and pressure 𝑃 has been previously determined, the number of atoms can be calculated with following equation: (3)
𝑁 = 𝑉𝜌𝑁𝐴/𝑀𝑤
12
where V (m3) is the volume of liquid-vapor region, ρ (g/m3) is the density, and, Mw (g/mol) represents the molecular weight. The calculated number of argon and krypton molecules for each thermodynamics state are shown in Table 1. Next, energy minimization is performed using the conjugate gradient method to relax the system, followed by equilibration of the system at the target temperature in a 𝑁𝑉𝑇 ensemble for 10 ns (Deaven and Ho 1995). The Langevin thermostat is used to control the temperature(Schneider and Stoll 1978, Dünweg and Paul 1991). Figure 3-a shows the temperature profile along the z-axis of simulation box for system-1 before and after temperature equilibration. By examining the temperature profile for system-1 in figure 3-a it is evident that the system has adequately equilibrated after 10 ns of NVT simulation.
(b)
(a)
Figure 3. (a) Temperature profile along the 𝑧-axis of simulation box for system 1 before and after temperature equilibration; (b) Density profile along the z-axis of simulation box for system 1 before and after density equilibration.
In the next step we equilibrated the density of the vapor and liquid regions by performing NVT/GCMC simulations for the vapor and liquid regions. In the grand canonical Monte Carlo (GCMC) procedure, exchange of molecules with an implicit bulk fluid at the specified temperature
13
and chemical potential (𝜇) or fugacity (𝑇,𝑃) was performed. This is possible using the fugacity coefficient 𝜙(𝑇,𝑃) = 𝑓(𝑇,𝑃)/𝑃 that has been previously calculated and displayed in Table 1. Figure 4-b shows the density profile along the 𝑧-axis of simulation box for system-1 before and after density equilibration. Examination of the density profile in Figure 3-b indicates that the density of the system has adequately equilibrated after 10 ns of NVT/GCMC simulation.
Once the system is equilibrated at the target density and temperature, the wall separating the vapor and liquid regions was removed (Figure 2-b). To implement the NEMD/GCMC procedure to simulate mixing of the liquid and vapor phases, the simulation box was divided into three regions:
vapor control volume (10×10×50 nm3)
mixing region (10×10×350 nm3)
liquid control volume (10×10×50 nm3)
In the liquid and vapor control volumes, both NVT MD (Nose Hoover thermostat) and GCMC procedure were applied to ensure that the temperature and density of the control volumes are maintained at their respective target values. In the mixing region, NVE simulation was performed to model liquid-vapor mixing. The GCMC is applied every 100 MD steps while the GCMC attempts is 1000 with 1000 of MC moves. A MD time step of 2 fs was used. The trajectory file was saved every 20 ps and analyzed to generate density profile, temperature profile, diffusion coefficient, and cluster size statistics.
14
3. Results and Discussion 3.1.
Mixing of Liquid Argon with Vapor Argon In this section, the NEMD simulation results for mixing of liquid Ar with its own vapor
(Ar/Ar) are presented. System-1 and system-4 represent subcritical mixing, system-2 and system-5 are expected to exhibit transcritical mixing, and finally system-3 and system-6 follow supercritical mixing processes. 3.1.1.
Subcritical mixing of liquid Argon with vapor Argon In figure 4-a, the temperature profile and density profile of the red molecules (initially
vapor Ar), blue molecules (initially liquid Ar), and all Ar molecules along the 𝑧-axis of the simulation box are plotted for system-1 at 200 ns. The initial reduced temperature (𝑇𝑟 = 𝑇/𝑇𝑐) and initial reduced pressure (𝑃𝑟 = 𝑃/𝑃𝑐) of liquid Ar for system-1 are 0.85 and 0.65, respectively. The initial 𝑇𝑟 and initial 𝑃𝑟 of vapor Ar are respectively 1.03 and 0.55. The temperature and density of the vapor and liquid control volumes were maintained at their respective target values by using the NVT/GCMC procedure. As seen in this figure, the density profile of all particles decreases gradually away from the liquid control volume towards the liquid/vapor interface. Consistent with subcritical system, a dramatic drop in the particle density is observed at the sharp interface. The figure shows red particles penetrating (“condensation”) into the blue liquid region and blue particles “evaporating” into the vapor phase from the species density profiles. Within a distance of approximately 30 nm from the liquid/vapor interface (figure 4-b), the density and temperature are fairly constant. The variations of density and temperature along the z-axis of the simulation box starting from the vapor control volume to the liquid control volume are plotted on densitytemperature space with the vapor-liquid phase diagram (VLE) of LJ Argon. The figure indicates
15
that the temperature of subcritical interface is approximately 135K and that it is at local phase equilibrium for which saturated liquid and saturated vapor are established on either side of the interface at the prevailing saturation temperature. This result confirms that local phase equilibrium is established across a subcritical interface even in nonequilibrium multiphase system(Røsjorde, Fossmo et al. 2000).
(b)
(a)
(c) Figure 4. (a) Temperature profile and density profile of the red molecules, blue molecules, and all Ar molecules along the z-axis of simulation box in system-1 at 200 ns, (b) Temperature profile and density profile of the red molecules, blue molecules, and all Ar molecules from z = 145 nm to z = 225 nm in system-1 at 200 ns, (c) Variation of temperature-density profile in system-1 at 200 ns.
16
We next examined the evolution of the mixing process for another subcritical system-4. The temperature profile and density profile of the red molecules, blue molecules, and all argon molecules at 200ns along the 𝑧-axis of simulation box are shown in figure 5-a. For this system, the initial Tr and Pr of liquid Ar are the same as those for system-1. The initial pressure 𝑃𝑟 of vapor Ar is 0.55 (same as system-1) while its initial 𝑇𝑟 is 2.06, two times higher than that of system-1. With a higher temperature, the initial density of vapor in system-4 is lower than that of system-1. Figure 5-a shows the temperature profile decreases with 𝑧 with different slopes for the vapor and liquid phases. Despite having a higher vapor control volume temperature, the liquid-vapor interface formed in system-4 is similar to that of system-1. Since the pressure of the simulation box 𝑃 ≈ 28 atm is the same as the pressure for system-1, the temperature of the interface (which is in vapor-liquid phase equilibrium on either side of the interface) must be ~135K (figure 4-c and figure 5-b) which corresponds to the saturation temperature of Ar at a vapor pressure of 28 atm. These results suggest that the interface temperature is governed by the value of the saturation pressure, which is identical for these two systems. This finding enables us to estimate the temperature of a subcritical liquid/vapor interface given the system pressure; this in turn allow us to predict the coefficients of transport properties using kinetic theory since they depend critically on the temperature and pressure of the interface(Røsjorde, Fossmo et al. 2000, Røsjorde, Kjelstrup et al. 2001, Segal and Polikhov 2008).
17
(b)
(a)
Figure 5. (a) Temperature profile and density profile of the red molecules, blue molecules, and all Ar molecules along the z-axis of simulation box in system-4 at 200 ns, (b) Variation of temperaturedensity profile in system-4 at 200 ns.
3.1.2. Transcritical mixing of liquid Argon with vapor Argon We proceed to consider the evolution process for transcritical mixing. The temperature and density profiles for the red particles, blue particles, and all Argon particles along the 𝑧-axis of the simulation box are plotted for system-2 and system-5 at 200 ns in figure 6-a and figure 6-b respectively. For system-2, the initial (𝑇𝑟,𝑃𝑟) for liquid Argon and vapor Argon are (0.85,1.20) and (1.03,1.10) respectively. The initial 𝑇𝑟 and 𝑃𝑟 for the liquid and vapor phases for system-5 are the same as system-2 except that the initial vapor phase temperature 𝑇𝑟 = 2.06, two times higher than that for system-2. The profile for the density of Argon particles (green curve) for both systems show the formation of an interface that is more diffused than that observed for subcritical systems. The figures also show that there exists substantial penetration of the red particles (that is, mixing or “condensation” of red on the liquid side of the interface) into the blue liquid particles compared with subcritical systems. Specifically, a relatively larger number of red particles penetrate the blue
18
particle region in system-2 than those in system-5. Concurrently, a larger number of blue particles diffuse (evaporate) into the red particles (vapor side of the interface) in system-2 than those in system-5. Therefore, it is observed that mixing of liquid-vapor is more efficient in system-2 than system-5 which can be attributed to the fact that the vapor phase in system-2 is closer to the Argon critical point than that of system-5. Furthermore, these figures show that the interface in system-2 is thicker and more diffused than that in system-5.
System-2
(a)
(b)
Figure 6. Temperature profile and density profile of the red molecules, blue molecules, and all Ar molecules along the z-axis of simulation box at 200 ns for (a) system-2 and (b) system-5. It is interesting to explore how the density (𝜌 vs. 𝑧) and temperature (𝑇 vs. 𝑧) profiles of Argon displayed in figure 6 relate to the Argon critical point and vapor-liquid phase diagram. As shown in figure 7-a and 7-b, for these transcritical systems, the density and temperature of the interface traverses the Argon critical region on the ρ-T plane. The temperature of the interface for both systems fluctuates between 145K to 155K around the critical temperature of Argon of ~150.86K indicating that the interface is “unstable” fluctuating in and out of the two-phase dome near the critical point. There exist large density fluctuations (expected in the vicinity of the critical
19
point) and formation of significant number of molecular clusters of varying sizes. Despite the higher temperature of vapor control volume in system-5, the density and temperature of the interface are similar to system-2 (with temperature ranging between 145K and 155 K near the critical point) suggesting that the establishment of the interface to be strongly dependent on the pressure of the systems at 𝑃 = 55 atm, slightly above the critical pressure of Argon. Additionally, cluster analysis was performed to investigate the behavior of interface structure. Using a cut-off distance of 3.81Å to define cluster formation, which corresponds to the location of the first peak of the Ar/Ar radial distribution function, we were able to determine the clustering behavior of the Argon particles. In figure 7-c and figure 7-d, clusters with size exceeding 10 particles in interfacial region are shown for system-2 and system-5. Note that we have omitted the largest “macroscopic” cluster for visual clarity. These figures show that large clusters are formed especially on the vapor side of the interface suggesting that, unlike subcritical evaporation, Argon particles move to the vapor as small particle clusters rather than as individual particles.
20
(a)
(b)
(c)
(d) Z = 150 nm
Z = 250 nm
Figure 7. Variation of temperature-density profile at 200 ns in (a) system-2 and (b) system-5. Clusters of Ar atoms (c) system-2 and (d) system-5 at 200 ns. The macroscopic cluster and clusters having sizes smaller than 10 particles have been omitted for visual clarity.
3.1.3. Supercritical mixing of liquid Argon with vapor Argon We now consider the results for the case of supercritical mixing (system 3 and 6). The temperature and density profiles at 200 ns are depicted in figure 8-a and figure 8-b, respectively. For system-3, the initial (𝑇𝑟,𝑃𝑟) of the vapor and liquid phases are respectively (1.03, 2.2) and (0.85, 2.3) while those for system-6 are (2.06, 2.2) and (0.85, 2.3). For system-3, figure 8-a shows linear Argon density profile and temperature profile along the 𝑧-axis dimension of the simulation
21
box. We have also plotted the density-temperature profile on the 𝜌 ― 𝑇 plane in figure 8-c and 8d. Included in the figure is the Widom line which demarcates the supercritical space of Argon into liquid-like and vapor-like regions near the critical point. It is known that near the critical point of a pure fluid, the supercritical region is not homogenous; the density of a supercritical fluid changes sharply across the Widom line. Figure 8-c indicates that the density-temperature profile for system3 are located entirely in the supercritical “liquid-like” region without a marked change in Argon density. This is not the case for system-6 in which nonlinear temperature and density profiles are observed as shown in figure 8-b resulting in an “interface-like” density profile and the formation of significant number of molecular clusters of varying sizes (figure 8-e). With the temperature and density of Argon plotted on the 𝜌 ― 𝑇 plane, the profile crosses the Widom line transitioning from supercritical liquid-like states to vapor-like states. These results show that even though fluid is in its supercritical state, it can exhibit large density inhomogeneity depending on whether the fluid crosses the Widom line.
22
System-3
System-6
(b)
(a)
(d)
(c)
(e) Z = 175 nm
Z = 275 nm
Figure 8. Temperature and density profiles of the red molecules, blue molecules, and all Ar molecules along the z-axis of simulation box at 200 ns for (a) system-3 and (b) system-6. Variation of temperature-density profile at 200 ns in (c) system-3 and (d) system-6. Clusters of Ar atoms (e) system-6 at 200 ns. The macroscopic cluster and clusters having sizes smaller than 10 particles have been omitted for visual clarity.
23
3.1.4. Voronoi Tessellation Analysis We have calculated the Voronoi tessellation volume of Argon particles for system-1 (subcritical), system-2 (transcritical) and system-3 (supercritical) at 200 ns. Voronoi tessellation cells may be regarded as a measure of the “molecular volume” of each LJ particle. For subcritical mixing, we anticipate a well-defined sharp change in the volume size of Voronoi cells across the liquid-vapor interface. A gradual change of Voronoi volumes is expected for supercritical mixing. Transcritical mixing represents a domain when both subcritical and supercritical mixing types are present. Figure 9-a shows the Voronoi cell volume or “molecular volume” of Ar molecules for system-1 (subcritical system). The magnitude of Argon volume changes dramatically at the liquidvapor interface which indicate a well-defined interface. In the case of a transcritical system-2 (figure 9-b), superposition of a sharp interface to a diffusive interface is observed. While large variation or change in Voronoi volume is observed, the interface is smoothed due to the decrease in the surface tension. Particle clusters are formed at the interface which resulting in a wavy surface. Finally, with increase in the pressure of the system (figure 19-c), the gradient in volume of Ar molecules decreases monotonically, and fluid/fluid mixing (supercritical mixing) is observed. Similar observations have been reported in previous experimental study by Segal and Polikhov(Segal and Polikhov 2008).
24
750 (Å3) 60 (Å3)
(a)
3
300 (Å ) 60 (Å3) (b)
3
100 (Å ) 60 (Å3) (c) Z = 250 nm
Z = 150 nm
Figure 9. Ar molecular volume (a) subcritical system-1, (b) transcritical system-2, and (c) supercritical system-3 at 200 ns.
3.2.
Liquid Krypton/Vapor Argon System Although the Ar/Ar system described in Section 3.1 is a one-component system, it can be
regarded as an ideal binary mixture consisting of identical molecular properties (molecular size and interaction energy). Despite its simplicity, the results obtained have provided significant insights and fundamental understanding on the evolution of fluid/fluid interface and the influence of thermodynamics and phase equilibria on subcritical, supercritical, and transcritical mixing. We have also demonstrated that the NEMD/GCMC method can generate self-establishing liquid/vapor interfaces and model mixing of a liquid with a gas at different temperature and pressure. We will apply this new NEMD procedure to study mixing of liquid Krypton (Kr) with vapor Argon (Ar) – a binary mixture of monatomic species with different molecular sizes and interaction energies.
25
This study investigates the effects of dissimilar molecular properties on the evolution of mixing processes at subcritical, transcritical, and supercritical conditions.
3.2.1. Subcritical mixing of liquid Krypton with vapor Argon In figure 10, the temperature profile and density profile of the Ar molecules (initially vapor Ar), Kr molecules (initially liquid Kr) and Ar/Kr molecules along the z-axis of simulation box are plotted for system-7 at 200 ns. The initial temperature and initial pressure of liquid Kr are 178K and 31atm, respectively. The initial temperature and initial pressure of vapor Ar are respectively 310K and 27atm. The temperature and pressure of the vapor and liquid control volumes are maintained at the target values by using NVT/GCMC simulation. Similar to system-1 and system-4, the density profile of Ar/Kr particle decreases slightly outside the liquid control volume until the liquid-vapor interface. A dramatic drop in particle density is observed at the interface. Consequently, a sharp interface is established for subcritical liquid-vapor mixing. The figure shows that the Ar particles penetrate (“condense”) into the Kr liquid particles and the Kr particles “evaporate” into the vapor phase. Within a distance of approximately 10 nm away from the liquidvapor interface, the density and temperature are fairly constant.
26
Figure 10. Temperature profile and density profile of the Ar molecules, Kr molecules, and Ar/Kr molecules along the z-axis of simulation box in system-7 at 200 ns. Figure 11-a shows the variation of density-temperature-Ar mole fraction profile from vapor control volume to the liquid control volume (at 150 ns, 175 ns and, 200 ns) and the vapor-liquid phase diagram (VLE) of LJ Kr/Ar. The figure shows when these profiles are plotted in thermodynamic space, they collapse onto the same curve indicating that they are thermodynamically identical. To examine the thermodynamic equilibrium state of mixture at the interface for system-7, the VLE plot of Kr/Ar binary system at 170 K (interface temperature, see figure 10) is shown in Figure 11b. The figure shows that the vapor-liquid compositions of interface at 150 ns, 175 ns, and 200 ns are located on equilibrium line at pressure of 28 bar. This pressure is the average pressure of simulation box which suggests that, similar to the Ar/Ar system, the subcritical interface exhibits local phase equilibrium in which saturated liquid and saturated vapor are established on either side of the interface at the prevailing saturation temperature. This verifies that an equilibrium state has been achieved and that the equilibrated interface is being convected downstream maintaining the equilibrium composition on either side of the interface without any
27
change. The ability of the NEMD procedure to achieve an equilibrium state for a binary fluid system is a fundamental validation of the NEMD framework.
(b)
(a)
Figure 11. (a) Variation of temperature-density-Ar mole fraction profile in system-7, (b) the vaporliquid phase diagram (VLE) of the LJ Krypton/Argon at 170 K.
3.2.2. Transcritical mixing of liquid Krypton with vapor Argon The temperature profile and density profile of the Ar molecules (initially vapor Ar), Kr molecules (initially liquid Kr), and Ar/Kr molecules along the z-axis of simulation box are plotted for system-8 at 200 ns in figure 12-a. The initial temperature and initial pressure of liquid Kr for system-8 were 178 K and 58 atm, respectively. The initial temperature and initial pressure of vapor Ar were respectively 310 K and 54 atm. The density profile for Ar/Kr particles in system-8 shows the existence of an interface that is more diffused than that observed for subcritical system. The figure also shows that there exists substantial penetration of Ar particles (that is, mixing or “condensation” of Ar and Kr particles on the liquid side of the interface) into the Kr liquid particles
28
compared with subcritical systems. To investigate the phase behavior of interface in this system, the variation of density-temperature-mole fraction profile from vapor control volume to the liquid control volume at 200 ns and the vapor-liquid phase diagram (VLE) of the LJ Kr/Ar are shown in figure 12-b. This figure shows the density-temperature-mole fraction profile of transcritical systems traversing the critical region on the ρ-T-x space in interfacial region as the interface evolves in time. This behavior suggests large density fluctuations (expected in the vicinity of the critical point) and formation of significant number of molecular clusters of varying sizes.
(b)
(a)
Figure 12. (a) Temperature profile and density profile of the Ar molecules, Kr molecules, and Kr/Ar molecules along the z-axis of simulation box in system-8 at 200 ns, (b) Variation of temperature-density-mole fraction profile in system-8.
3.2.3. Supercritical mixing of liquid Krypton with vapor Argon In the NEMD simulation of the supercritical system, Krypton and Argon mix to establish a liquid-vapor interface as illustrated in figure 13-a. In figure 13-a, the temperature profile and density profile of the Ar molecules (initially vapor Ar), Kr molecules (initially liquid Kr), and total density profile of Kr/Ar along the z-axis are plotted at 200 ns. The initial temperature and initial 29
pressure of liquid Kr for system-9 were 178 K and 111 atm, respectively. The initial temperature and initial pressure of vapor Ar were respectively 310 K and 107 atm. The temperature and pressure of the vapor and liquid control volumes were maintained at the target values by using NVT/GCMC.
(b)
(a)
Figure 13. (a) Temperature profile and density profile of the red molecules, blue molecules, and all molecules along the z-axis of simulation box in system-9 at 200 ns, (b) Variation of temperature-density-mole fraction profile in system-9.
The evolution of the temperature and densities of Kr and Ar are captured in the NEMD simulation. Figure 13-a shows nonlinear temperature and density profiles are observed after 200 ns. The density profiles reveal diffusion-dominated mixing with Ar particles (red curve) penetrating deep inside the liquid region and Kr particles (blue curve) into the vapor region. It is interesting to note significant variation in total density (green curve) resulting in an “interface-like” density profile. Figure 13-b shows the variation of density-temperature-mole fraction profile from vapor control volume to the liquid control volume at 200 ns and the vapor-liquid phase diagram (VLE) of the LJ 30
Kr/Ar. Our NEMD result shows that the entire curve lies outside of the Kr/Ar vapor-liquid equilibrium saturation dome. Also, the temperature and density profiles cross the Widom line which differentiates gas-like and liquid-like states of supercritical Argon/Krypton mixture. We did not compute and display the Widom line, which is known to exist, for this mixture. For these supercritical profiles we note that the interface will continue to diffuse and spread in time and eventually a “steady state” profile will be obtained across the two control volumes that are maintained at fixed conditions.
4. Conclusion We have developed a NEMD simulation method to model subcritical, transcritical, and supercritical mixing of liquid and gas at high pressures. This method ensures that the temperature, pressure, and density of the liquid and vapor phases are held constant while mixing takes place between them. It can model a self-establishing liquid/vapor interface depending on whether the mixing is subcritical, supercritical, or transcritical. The procedure generates thermodynamic temperature, density, and mole fraction profiles as a function of spatial position and time and produces information at the microscopic level including the structure of liquid/vapor interface, molecular clustering statistics, molecular volumes, etc. The simulation method was applied to mixing of (i) liquid Argon with its own vapor and (ii) liquid Krypton with vapor Argon at pressures ranging from subcritical to supercritical values. Our simulations demonstrate that subcritical interface exhibits local phase equilibrium in which saturated liquid and saturated vapor are established on either side of the interface at the prevailing saturation temperature that is fixed by the pressure of the system. For transcritical systems, the interface is more diffused than that of subcritical systems. The interface consists of large particle clusters and its density-temperature
31
profile traverses the critical region. Similar to subcritical systems, irrespective of the temperatures of the liquid and vapor phases, the interface temperature is determined by the pressure of system. Finally, the results for supercritical systems show diffusion-controlled fluid mixing without the formation of an interface. However, the density-temperature variation cuts across the Widom line, that divides the supercritical region into “liquid-like” and “vapor-like” supercritical states, resulting in large density variation the resembles an “interface-like” profile.
To summarize, key contributions of this study include: (a) development and demonstration of a hybrid NEMD/GCMC method for simulating mixing of liquid and vapor at different temperature and pressure (boundary conditions); (b) this unique procedure produces for the firsttime trajectories of the liquid-vapor mixture that can be mapped onto the thermodynamic space (temperature-density-composition); and (c) it generates microscopic and molecular level structural information on mixing at different mixing conditions. Even though the Ar/Ar and Kr/Ar systems studied are simple monatomic molecular systems, the results obtained have provided significant insights and fundamental knowledge and information on mixing at subcritical, supercritical, and transcritical conditions that are not possible from macroscopic continuum-based models. The proposed NEMD method is sufficiently robust and versatile, and can be employed to study mixing of fluid mixtures consisting of realistic hydrocarbon fuel with combustible gases. We are currently applying the NEMD method to study mixing of propulsion fuels.
32
Acknowledgements This work is being performed under Air Force Phase II STTR contract FA9550-17-C-0015 with Dr. Malissa Lightfoot as the technical monitor. We gratefully acknowledge support from AFOSR for this STTR program. We would also like to thank Drs. Chiping Li, Malissa Lightfoot, Jeff Mills, and Venke Sankaran for their technical feedback during the course of this work.
33
Reference Bedeaux, D. and Kjelstrup, S. (2004). "Irreversible thermodynamics—a tool to describe phase transitions far from global equilibrium." Chemical Engineering Science 59(1): 109-118. Bedeaux, D. and Kjelstrup, S. (2018). "Fluid-Fluid Interfaces of Multi-Component Mixtures in Local Equilibrium." Entropy 20(4): 250. Bellan, J. (2000). "Supercritical (and subcritical) fluid behavior and modeling: drops, streams, shear and mixing layers, jets and sprays." Progress in Energy and Combustion Science 26(4-6): 329-366. Chehroudi, B. (2012). "Recent experimental efforts on high-pressure supercritical injection for liquid rockets and their implications." International Journal of Aerospace Engineering 2012. Chen, L., Zhang, X. R., Okajima, J. and Maruyama, S. (2013). "Numerical investigation of nearcritical fluid convective flow mixing in microchannels." Chemical Engineering Science 97: 6780. Dahms, R. N. and Oefelein, J. C. (2013). "On the transition between two-phase and single-phase interface dynamics in multicomponent fluids at supercritical pressures." Physics of Fluids 25(9): 092103. Deaven, D. M. and Ho, K. M. (1995). "Molecular geometry optimization with a genetic algorithm." Physical Review Letters 75(2): 288. Dünweg, B. and Paul, W. (1991). "Brownian dynamics simulations without Gaussian random numbers." International Journal of Modern Physics C 2(03): 817-827. Givler, S. D., Abraham, J. (1996). "Supercritical droplet vaporization and combustion studies." Progress in Energy and Combustion Science 22(1): 1-28.
34
H. Mayer, W. O., Schik, A. H. A., Vielle, B., Chauveau, C., Gokalp, I., Talley, D. G. and Woodward, R. D. (1998). "Atomization and breakup of cryogenic propellants under highpressure subcritical and supercritical conditions." Journal of Propulsion and Power 14(5): 835842. Huai, X., S. Koyama and Zhao, T. (2005). "An experimental study of flow and heat transfer of supercritical carbon dioxide in multi-port mini channels under cooling conditions." Chemical Engineering Science 60(12): 3337-3345. Inzoli, I., Kjelstrup, S., Bedeaux, D. and Simon, J. (2011). "Transfer coefficients for the liquid– vapor interface of a two-component mixture." Chemical engineering science 66(20): 4533-4548. Kjelstrup, S. and Hafskjold, B. (1996). "Nonequilibrium molecular dynamics simulations of steady-state heat and mass transport in distillation." Industrial & engineering chemistry research 35(11): 4203-4213. Meland, R. (2002). "Molecular exchange and its influence on the condensation coefficient." The Journal of chemical physics 117(15): 7254-7258. Meland, R. (2003). "Molecular dynamics simulation of the inverted temperature gradient phenomenon." Physics of Fluids 15(10): 3244-3247. Meland, R., Frezzotti, A., Ytrehus, T. and Hafskjold, B. (2004). "Nonequilibrium moleculardynamics simulation of net evaporation and net condensation, and evaluation of the gas-kinetic boundary condition at the interphase." Physics of Fluids 16(2): 223-243. Mo, G. and Qiao, L. (2017). "A molecular dynamics investigation of n-alkanes vaporizing into nitrogen: transition from subcritical to supercritical." Combustion and Flame 176: 60-71.
35
Olivier, M. L., Rollier, J. D. and Kjelstrup, S. (2002). "Equilibrium properties and surface transfer coefficients from molecular dynamics simulations of two-component fluids." Colloids and Surfaces A: Physicochemical and Engineering Aspects 210(2-3): 199-222. Oschwald, M., Smith, J., Branam, R., Hussong, J., Schik, A., Chehroudi, B. and Talley, D. (2006). "Injection of fluids into supercritical environments." Combustion Science and Technology 178(1-3): 49-100. Plimpton, S. (1995). "Fast parallel algorithms for short-range molecular dynamics." Journal of computational physics 117(1): 1-19. Røsjorde, A., Fossmo, D., Bedeaux, D., Kjelstrup, S. and Hafskjold, B. (2000). "Nonequilibrium molecular dynamics simulations of steady-state heat and mass transport in condensation: I. Local equilibrium." Journal of colloid and interface science 232(1): 178-185. Røsjorde, A., Kjelstrup, S., Bedeaux, D. and Hafskjold, B. (2001). "Nonequilibrium molecular dynamics simulations of steady-state heat and mass transport in condensation. II. Transfer coefficients." Journal of colloid and interface science 240(1): 355-364. Savin, T., Schweizer, M. and Öttinger, H. C. (2016). Nonequilibrium thermodynamics of an interface. APS Meeting Abstracts. Schneider, T. and Stoll, E. (1978). "Molecular-dynamics study of a three-dimensional onecomponent model for distortive phase transitions." Physical Review B 17(3): 1302. Schweizer, M., Öttinger, H. C. and Savin, T. (2016). "Nonequilibrium thermodynamics of an interface." Physical Review E 93(5): 052803. Schweizer, M., Öttinger, H. C. and Savin, T. (2016). "Nonequilibrium thermodynamics of an interface." Physical Review E 93(5): 052803.
36
Segal, C. and Polikhov, S. (2008). "Subcritical to supercritical mixing." Physics of Fluids 20(5): 052101. Simon, J. M., Bedeaux, D., Kjelstrup, S., Xu, J. and Johannessen, E. (2006). "Interface film resistivities for heat and mass transfers integral relations verified by non-equilibrium molecular dynamics." The Journal of Physical Chemistry B 110(37): 18528-18536. Simon, J. M., Kjelstrup, S., Bedeaux, D. and Hafskjold, B. (2004). "Thermal flux through a surface of n-octane. A non-equilibrium molecular dynamics study." The Journal of Physical Chemistry B 108(22): 7186-7195. Talu, O. and Myers, A. L. (2001). "Reference potentials for adsorption of helium, argon, methane, and krypton in high-silica zeolites." Colloids and Surfaces A: Physicochemical and Engineering Aspects 187: 83-93. Tsuruta, T. and Nagayama, G. (2005). "A microscopic formulation of condensation coefficient and interface transport phenomena." Energy 30(6): 795-805. Wilhelmsen, Ø., Trinh, T. T., Kjelstrup, S. and D. Bedeaux (2015). "Influence of curvature on the transfer coefficients for evaporation and condensation of Lennard-Jones fluid from squaregradient theory and nonequilibrium molecular dynamics." The Journal of Physical Chemistry C 119(15): 8160-8173. Xu, J., Kjelstrup, S., Bedeaux, D., Røsjorde, A. and L. Rekvig (2006). "Verification of Onsager's reciprocal relations for evaporation and condensation using non-equilibrium molecular dynamics." Journal of colloid and interface science 299(1): 452-463. Yang, V. (2000). "Modeling of supercritical vaporization, mixing, and combustion processes in liquid-fueled propulsion systems." Proceedings of the Combustion Institute 28(1): 925-942.
37
Farzin Rahmani: Methodology, Software, Writing - Original Draft, Investigation Timothy Weathers: Writing - Review & Editing, Investigation Ashvin Hosangadi: Conceptualization, Writing - Review & Editing Yee C. Chiew: Conceptualization, Methodology, Writing - Original Draft, Supervision,
Declaration of interests ☒ 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. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
38
Highlights: A non-equilibrium molecular dynamics (NEMD)/ Grand Canonical Monte Carlo (GCMC) method is employed to simulate subcritical, transcritical and, supercritical mixing. Local phase equilibrium is established at subcritical interface. Local density and local temperature of the transcritical interface can be directly mapped to the near critical region of the Ar/Ar and Kr/Ar vapor-liquid phase diagram. For transcritical mixing, superposition of a sharp interface to a diffusive interface is
39
observed. For transcritical mixing, significant clustering of molecules is found at the “transitional diffuse” interface.
40