Molecular investigation on the anomalous phenomenon at liquid desiccant surfaces for air conditioning

It is necessary to disclose the two-phase interphase behavior in the liquid desiccant dehumidifier/regenerator applicable for air conditioning, but the present investigation is far from enough. In this paper, the surface structure of liquid desiccant solution is analyzed by molecular dynamics simulations. LiBr-H2O is chosen as the working solution with a concentration of 1 M and the system model is built with Gromacs. System temperatures vary from 300 to 350 K covering the temperature range of liquid desiccant dehumidification and regeneration. Density profiles of ions and water molecules are plotted along the vertical directions, and their distribution preferences on the solution surface are discussed. With the molecular simulation, it is found there is an ions-vapor layer with a thickness of 6–9 Å between the saturated vapor and bulk solution, which is not shown in the traditional macroscopic models. The results show that the density of water remains stable in the bulk while decreases sharply on the solution surface. However, the salt ions, i.e. Li+ and Br-, have a peak density on the surface. This ions-vapor layer behaves like a buffer to transfer water molecules from/to the bulk solution. More research will be required to investigate how to control the ions-vapor layer, so that air dehumidification and solution regeneration can be easily operated, which provides significant energy savings for the liquid desiccant air conditioning.


Introduction
To address the energy crisis and climate changes in the path of continuous development, China has proposed ten key energy-saving projects to encourage energy saving, such as initiatives in public transport, high-performance appliances, energy-saving buildings, and so on. As the main contribution of total energy consumption, energy-consuming in buildings has drawn great attention all over the world (Sun et al. 2018;Kim et al. 2019), in which the energy consumed by air conditioning occupies ~30% in China for its low efficiency. To enhance the efficiency of air conditioning, researchers have been exploring new ways, among which the temperature and humidity independent control (THIC) has great potential of energy saving (Liu et al. 2017). In the THIC system, the liquid desiccant is widely used to achieve the humidity control of air, such as LiBr aqueous solution. The strong solution with a relatively lower temperature can absorb the moisture from the processed air (i.e. liquid desiccant dehumidification) and becomes weak solution. Then the weak solution is heated to a higher temperature, which desorb the moisture into the ambient air (i.e. liquid desiccant regeneration), and changes back to strong solution. The performance of the dehumidifier directly determines the efficiency of humidity control of THIC system. And, the liquid desiccant regeneration is an energy-intensive process, which requires external energy sources to heat the weak solution and it is the most challenging part of the THIC system. Therefore, to achieve an energy-efficient regeneration process, it is necessary to BUILD SIMUL (2020) 13: 599-608 https://doi.org/10.1007/s12273-020-0611-8 Luo et al. / Building Simulation / Vol. 13,No. 3 600 have a comprehensive knowledge of the heat and mass transfer process in the dehumidifier/regenerator. Numerous studies have been conducted experimentally and theoretically on the heat and mass transfer process at air-liquid desiccant interface in the dehumidifier/regenerator from the macroscopic scale (Luo et al. 2015; Mohaisen and Ma 2015;Ge et al. 2012;Wen et al. 2018;Wen and Lu 2019). The theoretical models could be classified to the finite difference model (Luo et al. 2011(Luo et al. , 2012Wang et al. 2019), Effectiveness NTU (ε-NTU) model (Zhang et al. , 2018, and the simplified model (Peng and Zhang 2011;Dong et al. 2017). All of the above theoretical models are based on the double layer theory, which assumes that there is a saturated vapor layer above the solution surface and the mass transfer occurs between the saturated vapor layer and ambient air (Luo et al. 2014). In the double layer theory, it is supposed: (1) at the interface, there is a state of equilibrium for the two phase, (2) two laminar sub layers exist at both side of the interface and the main transport resistance lies in the two thin and stagnant films, (3) the mass transfer process in the two films can be described by diffusion, (4) the concentrations are uniform for the bulk phases. The specific transport mechanism is shown in Fig. 1. In the figure, P B is the partial pressure of component B in the gas phase and X B is the mole fraction of component B in the liquid phase.
Even though it is convenient to use the double layer theory for investigation, it has its own defects, including inappropriate assumption of the relationship between the mass transfer rate and the molecular diffusivity, inaccurate estimation of the layer thickness, and applicable only for the steady mass transfer process (Luo et al. 2014). Ultimately, all of the above drawbacks could be attributed to the lack of the knowledge of the gas-liquid interface behavior.  (Luo et al. 2014) The fast development of computer technology makes it more available to apply the molecular dynamics simulation technology to investigate the complicated behaviors at the two-phase interface and the corresponding microscopic characteristic of mass transfer in the devices. Till date, molecular dynamics simulation technology has been successfully used in the field of aqueous solution. However, microscopic simulations and experiments on the air-solution interfaces are very limited, and most of the research is focused on sea water rather than liquid desiccant solution Knipping et al. 2000;Garrett 2004;Ghosal et al. 2005;Jungwirth and Tobias 2006;Caleman et al. 2011). The research on seawater had indicated that the double layer model could not show the real surface structure of seawater. Knipping et al. (2000) investigated a 6 M aqueous NaCl solution. It was found that 10%-15% of the accessible surface area was taken up by chloride ions, while sodium ions were excluded from the topmost layer of liquid effectively. Ghosal et al. (2005) analyzed the composition at the liquid/vapor interface of KBr and KI with an x-ray photoelectron spectroscopy at near ambient pressure. In both of the above two cases, the halide anion had higher concentration at the surface of the saturated solution than within the bulk solution. Caleman et al. (2011) presented a quantitative analysis of the energetics of ion solvation on the basis of molecular simulations of all stable alkali and halide ions in water droplets. Employing a polarizable potential model, the concentration dependence of Na + and Cl − ions solvated in a water slab had been investigated . For all conditions, the sodium ions were located in the interior exclusively, but the chloride anions took up a great portion at the slab surface. Recently, to investigate the absorption phenomenon in absorption refrigeration systems, Weng et al. (2016) conducted molecular simulation to a traditional liquid desiccant (LiBr + H 2 O) and a novel mixture (NaCOOH + LiBr + H 2 O). However, the real surface structure was not obtained since the nonpolarizable water model (TIP3P) was used for simulation.
As shown above, most of the previous microscopic research on air-solution interface behavior with dynamics simulation technology focused on the seawater at room temperature (~20 °C) for aerosol reactions. However, there was rare microscopic investigation from the molecular level on the two-phase interphase behavior in the liquid desiccant dehumidifier/regenerator applicable for air conditioning. It is important to understand the microscopic surface structure of liquid desiccants to achieve energy-efficient air conditioning. Therefore, this paper analyzes the surface structure of a popular liquid desiccant (LiBr-H 2 O) with molecular dynamics simulation that is a powerful technology to disclose fundamental mechanisms (She et al. 2016).
The system model is built in Gromacs with a solution concentration of 1 M under various system temperatures (300 K, 325 K, and 350 K). Density profiles of ions and water molecules are plotted along the vertical direction of the solution surface, and their distribution preferences on the solution surface analyzed. Besides, saturated vapor pressure at the solution surface is predicted with the model and then compared with the empirical data. The advantage of present model lay in that it provided a brand-new view to observe the mass transfer process in the liquid desiccant dehumidifier. It could be used to investigate the dehumidification/ regeneration behavior from the microscopic molecular level and uncover the underlying mechanism of the mass transfer in the dehumidifier/regenerator, which might be good pointcuts for enhancing the performance of the devices.  system is located at the center of the rectangular box, and the initial surface position of the aqueous solution is at z = ±2 nm. Periodic boundary conditions are adopted for three dimensions. The velocity rescaling with a stochastic term is used to control the system temperature (Bussi et al. 2007), which has been proved to sample canonical ensemble. The Leap-Frog integration scheme is applied to solve Newton's equation of motion with a time step of 1 fs. The smooth Particle-Mesh Ewald method is adopted to calculate the long-range electrostatics (Essmann et al. 1995), and the short-range electrostatic potential is truncated and shifted at 12 Å. The Van der Waals forces are smoothly switched to zero between 10 and 12 Å, and long-range corrections for energy and pressure are performed. The SETTLE algorithm is utilized to restrict the water molecules (Miyamoto and Kollman 1992). The slab of aqueous solution is simulated for 500 ps with NVT ensemble to reach a steady state before it is placed into the rectangular box. For the solution surface system, each simulation consists of 1000 ps equilibration and 1000 ps data production. To analyze the system property profiles, the system is divided into 100 bins in the z-direction, with a thickness of 2 Å for each bin.
The selection for molecular models of aqueous solution (i.e., water and ions) is crucial for the molecular simulation. The most commonly used water models are nonpolarizable, which has a permanent dipole and assumes that water molecules are uniformly polarized by the surroundings. Hence, they are only applicable for reproducing bulk properties of isotropic liquid. However, for the solution surface system, water molecules are non-isotropic at the surface. Therefore, in this paper, the SWM4-NDP polarizable model (Lamoureux et al. 2006) is selected for water molecules. The main characteristic is that it has an induced dipole besides a permanent dipole, which has been presented to be necessary for describing the solution surface. Correspondingly, a polarizable model (Yu et al. 2010) is employed for ions.
The SWM4-NDP water model includes 5 sites: an oxygen site (O), two hydrogen sites (H1, H 2 ), a Drude particle (D) and a massless site (M), as shown in Fig. 3. The molecular structure is a ball-and-stick model with a separate M site. The Drude particle is so close to the oxygen site that it cannot The Drude particle has a negative charge and the oxygen atom has the same but positive charge; they are linked through elastic force and form an induced dipole, which is not available in other water models. The Drude particles are minimized at each time step until the root-mean-square force on these particles is lower than 0.001 kJ/mol (van Maaren and van der Spoel 2001). The molecular polarizability (oxygen-Drude particle) is 0.97825 Å 3 . The hydrogen atom has a positive charge, which is half of the M site (negative charge); they form a permanent dipole. The M site is in the bisector line of the angle H-O-H with a fixed distance from the oxygen atom. Parameters of the SWM4-NDP water model are shown in Table 1.
For the polarizable model of ions, each ion is made up of one Drude particle (d) and one ionic core (I). The Drude particle has a negative charge, and the ionic core has the same but positive charge. They are linked through elastic force and form an induced dipole. The parameters of the polarizable ions are shown in Table 2.
The potential energy between water molecules includes elastic potential, coulombic potential and Lennard-Jones potential. It is computed as follows.
where, s and sʹ represent the five sites of the water molecule (O, H 1 , H 2 , M, D). N w is the number of the water molecules. q s and q sʹ are the charge of the five sites, and k D is the spring constant between O and D sites. ε and  OO are water parameters of the Lennard-Jones potential. r is the coordinate vector of the particles. For the interaction between ions and water molecules, it consists of coulombic potential and Lennard-Jones potential, and is calculated by, where, q I and q d are the charges of the ionic core and Drude particle, respectively. N I represents the number of ions. The Lennard-Jones parameters between ions and water are combined according to the Lorentz-Berthelot rule.
The potential energy between ions is composed of elastic potential, coulombic potential and Lennard-Jones potential, and is computed in the following.
( ) I  I   I   2  2  I-I  d  ,I  ,d  1  ,  ,  ,   6  12  II  II  II  ,I  ,I  ,I  ,I where, t and tʹ represent the two sites of the ion (I, d). k d is the spring constant between I and d sites. ε and  II are the ion parameters of the Lennard-Jones potential. The global scalar pressure P is shown below (the first term is the kinetic energy and the second term is the inner virial): where N is the number of atoms in the system, K B is the Boltzmann constant, T is the temperature, 3 is the dimensionality of the system (2 or 3 for 2d/3d), V is the system volume (or area in 2d), r i denotes the vector from the origin to the atom i, and F i denotes the vector of total force acting on atom i.

Model validation
It is necessary to validate the system model before doing further analyses. Here, pure water is chosen as the solution and molecular dynamic simulation is conducted following the above procedures. Time evolutions of the energy and temperature of the water system are shown in Fig. 4. Both of the energy and temperature reach equilibrium quickly in 1000 ps, indicating that it is ready for data collection. Figure 5 shows the snapshots of water systems at steady states with the system temperature at 300 K, 325 K, and 350 K, respectively. It can be seen that water molecules at the surface tend to move away from the bulk water. What's more, with the increase of the system temperature, more water molecules can be found above the surface. This is because a higher system temperature leads to higher kinetic  Radial distribution function or pair correlation function, g(r), shows the probability of finding a particle at a distance of r away from a given reference particle, relative to that for an ideal gas. The first peak in g(r) indicates the bond length. As the r is long enough, g(r) should approach 1.0, indicating that the two particles are not interacted and behave like an ideal gas. The radial distribution functions of water at 300 K are shown in Fig. 6. From the first peaks of O-H and H-H curves, the O-H bond length is 0.95 Å and the closest distance between H atoms is 1.5 Å, which matches well with the water model (l OH = 0.9572 Å, θ HOH = 104.52°). Besides, the closest distance between O atoms is 0.28 Å based on the first peak of O-O curve. Finally, it is reasonable to say that the system model is valid for further analyses.  Figure 7 shows the snapshot of the surface condition of LiBr-H 2 O when the system temperature is 300 K. It is found that more bromide ions (yellow) occur on the surface, compared to the lithium ions (purple). The above phenomena will be analyzed in detail in the following. Fig. 7 Snapshots of surface conditions of LiBr-H2O with system temperature at 300 K. Coloring scheme: water oxygen, red; water hydrogen, white; lithium ions, purple; chloride ions, blue; bromide ions, yellow 4.1 Density distribution of particles in and above the solution Figure 8 demonstrates the density distribution of various particles (water, lithium, and bromide) in the LiBr-H 2 O system along the vertical direction (i.e., z-direction) with different system temperatures (300 K, 325 K, and 350 K). Figure 8(a) shows that at the system temperature of 300 K, water density remains unchanged in the bulk, while it decreases significantly on the surface (z is around ±2 nm). Besides, it is found that the water density of the bulk has a slight decrease when the system temperature increases. This is because a higher system temperature contributes to higher kinetic energy of molecules, which drives more water molecules into the surface. Correspondingly, there will be a slight increase of the water density on the surface with the increase of the system temperature.  Fig. 8(b), it is found that at 300 K, the density of lithium ions in the bulk has small fluctuation, and a peak value appears on the surface. The result means that the lithium ions tend to occupy the surface instead of the bulk, which is different from water molecules which prefer to stay in the bulk. When the system temperature is lower, the distribution of lithium ions is less sensitive than that at a higher system temperature. The density of lithium ions has a small decrease in the bulk when the system temperature increases from 300 K to 325 K. Nevertheless, as the system temperature increases from 325 K to 350 K, the density of lithium ions has an obvious decrease in the bulk while an obvious increase is observed at the surface. Figure 8(c) presents the density distribution of bromide ions. As the system temperature increases from 300 K to 325 K, the density of bromide ions on the surface increases significantly, and the increase rate is more obvious than lithium ions. It seems that bromide ions on the surface are more sensitive to the system temperature than lithium ions, particularly when the system temperature is relatively low. Besides, it should be noted that ions are more sensitive to the system temperature than the water molecules, which may be due to the existence of hydrogen bonds between water molecules.

Distribution preference of particles at the solution surface
The distribution preference (i.e., normalized density) is defined to be the densities of ions/water oxygen atoms (ρ) normalized by their bulk densities (ρ b ). Figure 9 shows the distribution preferences of various particles (O, Li and Br) at the surface (~2 nm) of LiBr-H 2 O solution. When the system temperature is 300 K, as shown in Fig. 9(a), the normalized density of oxygen atoms at the surface decreases significantly compared to the bulk density, indicating that water molecules tend to stay in the bulk rather at the surface. However, for the ions of Li and Br, a peak is observed around the surface, showing that the ions prefer to occupy the surface area. What's more, bromide ions have a much higher normalized density on the surface, compared to lithium ions. This indicates that the bromide ions have more favorite to occur on the surface than the lithium ions. Figure 9(b) shows that the normalized density of lithium ions at the surface increases a little as the system temperature increases to 325 K. This is because a higher system temperature results in higher kinetic energy of particles, which has a more probability to move to the surface. However, for the bromide ions, the normalized density at the surface decreases a little compared to that at 300 K. This is mainly affected by the bulk density. The differences of normalized density on the surface between bromide and lithium ions become smaller with the increase of the system temperature, as shown in Fig. 9(c). It should be noted that even with a higher system temperature (350 K), bromide ions still have a larger normalized density than the lithium ions on the surface. In the previous research of LiCl-H 2 O, it was reported that at the same system temperature (350 K), little difference of distribution preference was found between the chloride ions and lithium ions (She et al. 2019), which indicates that chloride ions and lithium ions have similar preference to take up the surface.

Proposed heat and mass transfer model at the air-solution interface
In the traditional air-solution interface model, as shown in Fig. 10(a), it is usually assumed that there is a saturated vapor layer directly above the bulk solution. However, based on the above analyses, it can be seen that there is an ions-vapor layer between the vapor and bulk solution, as shown in Fig. 9(b). In this layer, anions have a much higher quantity than cations. What's more, anions prefer to stay at the top of the layer. The thickness of the layer is 6-9 Å. This layer may play a very important role during the heat and mass transfer process, which is usually neglected in the previous model.
The ions-vapor layer behaves like a buffer to transfer water molecules from/to the bulk solution. In the process of air dehumidification, it is expected that more ions can move to the ions-vapor layer from the bulk solution, so that this layer can easily absorb water molecules from humid air even at room temperatures. In the process of solution regeneration, ions at the ions-vapor layer are expected to move to the bulk solution, so that this layer is mainly occupied by water molecules and has a large vapor pressure, which allows the operation of the regeneration process more easily even at a lower heat source temperature. Therefore, if researchers can find a mechanism to control the ions-vapor layer in the future, it be a revolution for the liquid desiccant air conditioning.

Predication of vapor pressure at the solution surface
For aqueous solutions, vapor pressure is one of the crucial parameters to evaluate their feasibility of dehumidification. Figure 11 shows the prediction of vapor pressure for pure H 2 O and LiBr-H 2 O under different system temperatures (300 K, 325 K and 350 K). The vapor pressure calculation is based on GROMACS-LS, a custom GROMACS version to compute the local stress tensor in 3D from a molecular dynamics simulation (Vanegas et al. 2014). For LiBr-H 2 O, the saturated vapor region is between 7.6 and 10.0 nm. And it is considered as the saturated vapor region. Small fluctuations are observed for LiBr-H 2 O. The result is different from that of LiCl-H 2 O, whose vapor pressure fluctuates significantly at first and then remains almost unchanged between 6.0 and 10.0 nm (She et al. 2019). Figure 12 shows the predication of saturated vapor pressure for pure H 2 O and LiBr-H 2 O based on the saturated vapor region in Fig. 11. By observation, it is found that simulation data has a good agreement with the empirical data when the system temperature is lower, such as 300 K.
As the system temperature increases, the deviations become more obvious, while the change trends of the saturated vapor pressure are identical. The deviations are mainly because the current water models are developed for biomolecules at room temperatures. Therefore, it may be not that accurate to predict the properties of aqueous solution especially at a higher temperature.

Conclusions
In this paper, molecular dynamics simulations have been conducted to analyze the surface structure of liquid desiccant for air conditioning. LiBr-H 2 O is chosen as the working solution. The system model is built in Gromacs at different system temperatures (300 K, 325 K, and 350 K) with a solution concentration of 1 M. The density profiles and distribution preferences of ions and water molecules in and above the solution are analyzed. Finally, as one of the essential parameters, vapor pressure profiles are predicted, and the result is compared with the empirical data. All these studies will be helpful to understand liquid desiccant dehumidification/regeneration processes from molecular scales. The main conclusions are highlighted as follows: 1) Molecular simulation shows that there is an ions-vapor layer between the saturated vapor and bulk solution, which is not shown in the traditional macroscopic models. This ions-vapor layer behaves like a buffer to transfer water molecules from/to the bulk solution. More research will be required to investigate how to control the ions-vapor layer, so that the air dehumidification can be easily operated at room temperatures and weak solution can be easily regenerated even at a low heat source temperature. 2) Water density remains unchanged in the bulk solution, while decreases greatly on the solution surface, which indicates that water molecules prefer to stay in the bulk solution. However, for the ions (Li + and Br -), there are density peaks on the surface, showing that ions prefer to occur on the solution surface instead of the bulk solution.
Besides, it should be noted that ions are more sensitive to the system temperature than the water molecules, which may be due to the existence of hydrogen bonds between water molecules. 3) For LiBr-H 2 O, the ions-vapor layer has a thickness of 6-9 Å. In this layer, bromide ions have a more preference to occur on the surface than lithium ions for all temperatures. Therefore, it seems that the distribution layers on the solution surface from the top to the bottom should be bromide ions, lithium ions and water molecules in sequence.