Molecular Dynamics Simulation of the Soret Effect on Two Binary Liquid Solutions with Equimolar n-Alkane Mixtures.

Molecular dynamics is employed to simulate the Soret effect on two binary liquid solutions with equimolar mixtures: normal pentane (n-pentane, nC-5) and normal heptane (n-heptane, nC-7) molecules plus normal decane (n-decane, nC-10) and normal pentane molecules. Moreover, two coarse-grained force field (the CG-FF) potentials, which may depict inter-/intramolecular interactions fairly well among n-alkane molecules, are developed to fulfill such investigations. In addition, thermal diffusion for the mass fraction of each of these n-alkane molecules is simulated under an effect of a weak thermal gradient (temperature difference) exerting on solution systems from their hot to cold boundary sides. Finally, quantities of the Soret coefficient (SC) for two binary solutions are calculated by means of the developed CG-FF potentials, so as to improve the calculation rationality. As a result, first, it is found that molecules with light molar masses will migrate toward the hot boundary side, while those with heavy molar masses will migrate toward the cold boundary one ; second, the SC quantities indicate that they match relevant experimental determinations fairly well, i.e., trends of these SC quantities show inverse proportionality to the thermal gradient on the systems.

Introduction

Explorations of petroleum resources are of great significance to the development of global economy in human society. In practice, how to effectively delimit a crude oil reservoir underneath still mainly depends upon accurate predictions of distributions for mass fractions of various chemical components in crude oils. According to this, several research institutions in the European Union (EU) have focused on an application of the so-called Soret effect on the oil explorations.[1−15] The Soret effect (the SE), known as one of the important phenomena in thermal diffusion, has been investigated for tens of years. It refers to that, in an isolated and isotropic liquid reservoir, due to the sole effect of a thermal gradient exerting on it (without an effect of gravitation), all the components in the reservoir will evolve their mass fraction distributions accordingly.[1−3] Early lab works have revealed that the Soret coefficient (the SC) quantities used for testing mass fractions of chemical components in the SE should play the main role in exploring crude oil resources.[4−7] However, in practice, it is doubted that distributions for mass fractions of various chemical components in crude oils were influenced by both a thermal gradient and a gravitation.[8−15] Therefore in 2016, “the Sino-EU Cooperative Soret Coefficients in Crude Oil Project” synthesized six sample cells mixing several normal alkane components (n-alkane molecules, just having a carbon bough in each molecular conformation) extracted from crude oils. Then, they were sent into space by launching a Chinese SJ-10 retrievable satellite (the SJ-10 Project), so as to make clear if the SC did playthe main role in evolving distributions for mass fractions of n-alkane components in crude oils under a very weak gravity field.[16,17]

On the other hand, during the past two decades, with rapid developments of computational algorithms, some modeling methodologies such as molecular dynamics (MD), the Monte Carlo (MC), and the finite element methods (FEM) have shown to be very competent to investigate the SE mechanism in several liquid solutions mixing various chemical components and worked out some amazing outputs. For example, Montel et al. investigated the isothermal and nonisothermal properties of some organic liquids by using the MD simulations with the Lennard-Jones (L-J) potentials.[18] They found out that, in the SE, all the components in liquids may have their mass fractions distributed diversely. Touzet et al. applied the MD simulations with some simple models (rigid sphere L-J potentials) to a ternary liquid solution mixing methane (nC-1), n-butane (nC-4), and n-dodecane (nC-12) molecules. They also observed that, qualitatively, nC-1 and nC-4 molecules would migrate toward the hot region in solution, while nC-12 molecules would migrate toward the cold one.[19] Still, from the SJ-10 Project, Galliéro et al. investigated the thermal diffusion for methane (nC-1), n-butane (nC-4), and n-dodecane (nC-12) molecules in those six sample cells by means of the gas chromatography (GC) and the MD simulations with several simple isotropic plus multipolar L-J potentials. They also found that, qualitatively, nC-1 and nC-4 molecules would gather at the hot-temperature area, but nC-12 molecules would gather at the cold one.[17,20] In addition, Mozaffari et al. applied the MD simulations with some simple L-J potentials to the Soret effect on a ternary liquid solution mixing equimolarly n-dodecane (nC-12), tetra-hydro-naphthalene, and iso-butylbenzene. They calculated the SC quantities and compared them with experimental data from the SJ-10 Project.[21] Unfortunately, they lost some important sample cells from the SJ-10 Project. Thus, their SC quantities could hardly be validated. Overall, all the molecular potentials in the above studies looked quite simple, and relevant simulation objects were too small, so their MD outputs could hardly reveal the essentials of the SC mechanism.

Apparently, one of most important factors to run the MD simulations successfully was how to constitute reliable models of inter-/intramolecular potentials for relevant research objects. For example, during the past two decades, the so-called anisotropic united-atom (AUA) model, which took −CH2, −CH3, and CH4 as unified interaction sites (hence united atoms) on the conformation of an alkane molecule, has shown excellent competence in many studies.[22−36] To improve the prediction accuracy of inter-/intramolecular potentials among organic molecules, in the MD simulations, Antoun et al. investigated the SE on a binary liquid solution mixing nC-5 and nC-10 molecules even by means of a model based upon the all-atom force field.[37] In general, inter-/intramolecular potentials in above studies consisted of five terms: radial stretching, angular bending, torsion (spatial geometry deforming), and van der Waals and electrostatic interactions. Moreover, the AUA model can be viewed as a simplified version of the all-atom force field. However, in practice, for electroneutral n-alkane molecules, both torsion and electrostatic terms in the AUA model can be ignored because they played obvious roles only in those alkane isomers (i-alkane molecules, having a carbon bough and twigs in each molecular conformation),[36] i.e., for a specific i-alkane molecule, the torsion term was found to behave obviously just at the joints where carbon twigs hinged.[39] However, for those macro n-alkane (larger than n-decane) molecules, adopting the AUA model would increase computational overheads when carrying out the MD simulations. According to this, if we ran the MD simulations for some super large-scale systems (about hundreds of nanometers in size), which were composed of hundreds of thousands of macro alkane molecules in a very long simulation time (hundreds of millions of MD time steps), then we must take a further coarse-graining operation on the conventional AUA model, so as to improve the whole computational efficiency.

In addition, a subsidiary key factor to sustain the MD simulations for the SE on research systems was the so-called nonequilibrium molecular dynamics (the NEMD), which usually adopted an algorithm of the heat-flow exchange between systems and their outside environments (the ΔQ, also the heat convection) in fulfilling simulation processes.[30−35] However, some of the NEMD simulations indicated that the SC quantities showed severe fluctuations at the initial stage when the ΔQ was very small. Then, they would become increasingly stable along with the ΔQ increase, meaning that the ΔQ quantity may have a positive impact on the SC calculations if it increased to quite higher than a certain level.[35] Nevertheless, such a conclusion was quite doubted because thermal fluctuations due to heat convection in experiments usually severely disturbed the accuracy of determining the SC quantity no matter whether the ΔQ quantity was obviously large or not.[38] Therefore, such the NEMD conclusion needed further corroboration.

During the school year 2011, Zhang et al. provided some concise coarse-grained beads for depicting real conformations of macro alkane molecules.[39] In their framework, three units of conventional united atoms could even be extended to three larger new ones: −(CH2)2–CH3, −(CH2)3–CH3, and −(CH2)4–CH3. The reason was that, under the standard state, most of the macro alkane molecules were at a liquid state, but conventional united atoms (−CH2, −CH3, and CH4) were usually at a gaseous state. Thus, the bonding properties of those macro alkane molecules could hardly be represented by conventional ones. However, still under the standard state, −(CH2)2–CH3, −(CH2)3–CH3, and –(CH2)4–CH3 were likely liquefied and were widely regarded as one of the basic debris among macro alkane molecules at a liquid state in petrochemical engineering. According to this, here, new united atoms could be represented by three new coarse-grained beads: A1 (−(CH2)2–CH3), A2 (−(CH2)3–CH3), and A3 (CH3–(CH2)3–CH3), which may well splice those macro alkane (both n-alkane and i-alkane) molecules at a liquid state through their effective combinations. Reference (39) indicated that not only may such a coarse-graining strategy depict the colloidal stability of petroleum engineering fairly well but save more computational resources in the MD simulations.

In this article, at first, we constitute two coarse-grained force field (the CG-FF) potentials by means of new united atoms in ref (39). Then, as one of the working extensions in ref (39), the Soret effect on two binary liquid systems with n-alkane molecules (nC-5 and nC-7 molecules plus nC-5 and nC-10 molecules) is investigated by using the MD simulations with the developed CG-FF potentials. Our goal is that, if the developed CG-FF potentials here can succeed in carrying out all the MD simulations, then in the future, such techniques for constituting intra-/intermolecular potentials can also be competent to depict other complex molecules (e.g., those super macro alkane isomers) and may bridge very large-scale molecular systems and continuum media.

Figure shows a schematic of real structures of three selected n-alkane (nC-5, nC-7, and nC-10) molecules and their geometric modifications spliced by three selected coarse-grained beads (A1, A2, and A3) as mentioned above.

Figure 1

Coarse-grained modifications for real geometrics of three selected n-alkane molecules: nC-5, nC-7, and nC-10.

Table lists some physical parameters of the three selected n-alkane molecules in nature. In the past, since conventional united atoms can merely provide few pieces of physical and chemical information (e.g., vapor pressures, vaporization enthalpies, and liquid densities) for constituting some relevant inter-/intramolecular potentials in studies,[22−24,26−28] here, we adopt the first-principles calculations (e.g., calculations by using density functional theory, DFT) in optimizing those coarse-grained conformations in Figure , which would enrich the fitting information for constituting the CG-FF potentials, but ref (39) and other works did not adopt them before. That is one novel point in this article.

Table 1

Some Physical Parameters of n-Pentane (nC-5), n-Heptane (nC-7), and n-Decane (nC-10)a

All relevant data came from refs (39) and (40).

In our MD framework, to overcome the shortcomings in the NEMD, i.e., to eliminate thermal noise during the NEMD simulations at low temperature levels, we exerted two different temperature quantities (T1 and T2, T1 > T2) on a super cuboid reservoir containing hundreds of thousands of n-alkane molecules, from its left through right boundary sides. After running very long MD time steps, both hot- and cold-edge zones could be constructed at two boundary sides so that a quite steady thermal gradient (a temperature difference, DT = T1 – T2) can be generated to distribute through the edge zones,[41] which may play an equal role to the heat convection in the NEMD. That is another novel point in this article. Please note, here, we do not add an effect of gravitation into all the MD simulations. Moreover, all the MD simulations are under the standard state (e.g., initial pressure = 1.00 atm and initial T0 = room temperature (299 K)).

This article is organized as follows: Section narrates the MD simulation procedures, Section discusses the MD simulation outputs, Section summarizes the MD simulation results, and Section introduces the methodologies.

Simulation Procedures

In this article, all the MD simulations were carried out by using the LAMMPS, which is competent to run very large-scale molecular systems.[52] In details, ① initially, each of the MD systems with equimolar mixtures was constructed in a cuboid reservoir, e.g., A3 (nC-5):A1–A2 (nC-7) = 1:1 in a SET-1 system and A3 (nC-5):A1–A2–A1 (nC-10) = 1:1 in a SET-2 system. For the SET-1 system with 36,864 molecules (A3⊕A1–A2), its volume of 600 × 140 × 140 Å3 was stacked by a dozen of cuboid cells (each cell length in the x-direction: 40 Å). Meanwhile, for the SET-2 system with 73,728 molecules (A3⊕A1–A2–A1), its volume of 792 × 176 × 176 Å3 was stacked by a dozen of cuboid cells (each cell length in the x-direction: 44 Å); ② under a periodic boundary condition in the x-, y-, and z-directions, each of the MD systems was equilibrated for approximately 100,000,000 time steps at a desired temperature of 299 K (1 time step = 0.001 ps), plus a hydrostatic pressure (1.00 atm) exerting on the system. This was called the NPT ensemble simulation, in which all the relevant coupling parameters are as follows: τ = 100 and τ = 1000 for the SET-1 system and τ = 200 and τ = 2000 for the SET-2 system.[41,52] After this step, the resulting volume of the SET-1 system would become 586 × 131 × 131 Å3, and that of the SET-2 system would become 782 × 174 × 174 Å3; ③ all the intra-/intermolecular potentials in each of the MD systems were described by the developed CG-FF models in Sections . The effectiveness of these potentials is shown in Table . ④ After thermal equilibrium at a desired temperature (299 K), for each of the MD systems, a hydrostatic pressure (1.00 atm) was still exerting on the system in the y- and z-directions; meanwhile, in the x-direction, there was no pressure to control the system, but two different temperatures (T1 and T2, T1 > T2, DT = T1 – T2) were applied onto its left and right boundary sides so that hot- and cold-edge zones can be constructed near these two boundary sides to generate a thermal gradient along the x-direction.

Table 2

Effectiveness of the CG-FF Models for Two Binary Solutions with Equimolar Mixtures: The SET-1 and the SET-2 Systems

parameters	system	the CG-FF	other works
coefficient of volume thermal expansion (×10–3, K–1)	SET-1	0.7063	0.6926a
	SET-2	0.4931	0.5065b
mass density (g/m3, 20 °C)	SET-1	0.6760	0.6600a
	SET-2	0.7010	0.6900c

Data came from ref (40).

Data came from ref (53).

Data came from ref (58).

In details, when simulating an MD system, the x-direction was under a free boundary condition, i.e., two vacuum spaces with their own lengths larger than two cutoff distances were set outside its left and right boundary sides, while in the y- and z-directions, it was still under periodic boundary conditions. As a result, a net thermal pressure due to the DT increase may appear in one middle region between two edge zones to drive molecular migrations along the x-direction. For example, along the x-direction, both hot- and cold-edge zones were controlled by two different external temperatures (e.g., T1 = 302, 305, and 308 K for the hot-edge zone, but T2 = 299 K for the cold-edge zone only. Thus, the DT between the two edge zones could be set at 3, 6, and 9 K, respectively). However, the middle region in the system would not perform temperature control but allow heat exchange and diffusion along the x-direction, which was realized by a constant-energy MD simulation,[41] see Figure .

Figure 2

Schematic of the MD simulation geometry.

In Figure , both Zone-1 and Zone-2 were named hot- and cold-edge zones, which were controlled by T1 and T2, respectively. The middle region between them would perform the heat evolution. In details, for the SET-1 system, along the x-direction, each of the two edge zones was set to 93 Å in length, and the middle region was set to 400 Å in length comprising 10 cuboid cells (each cell length along the x-direction: 40 Å). Please note, inside the middle region, two areas (Area-1 and Area-2, each of which was set to 120 Å in length) were adjacent to two edge zones, respectively, so that molecular numbers in these two areas could be counted during thermal diffusion. Similarly, for the SET-2 system, along the x-direction, each of its two edge zones was set to 172 Å in length, and the middle region was set to 440 Å in length comprising 10 cuboid cells (each cell length along the x-direction: 44 Å). Also, inside the middle region, each of the two areas was set to 132 Å in length. Finally, after 500,000,000 time steps of the thermal constraint on the two edge zones, all the MD systems were found to be under the steady state, and the molecular number in each of the 10 cuboid cells in the middle region was counted by an average of 4000 time steps, just for the SC calculations in Section .

Results and Discussion

Distribution Trends for Mass Fractions of n-Alkane Molecules in Binary Liquid Solutions

Figure a shows three distribution (regression) curves representing temperature spatial profiles in the middle region along the x-direction in response to three DT values between hot- and cold-edge zones in three SET-1 systems (DT = 3, 6, and 9 K, and T2 = 299 K in the cold-edge zone), respectively. Also, Figure b shows three distribution (regression) curves representing temperature spatial profiles in the middle region along the x-direction in response to three DT values between hot- and cold-edge zones in three SET-2 systems (DT = 3, 6, and 9 K, and T2 = 299 K in the cold-edge zone), respectively. In Figure , every temperature point was achieved by an average of 4000 time steps after 500,000,000 time steps of the thermal constraint on the two edge zones. Moreover, the solid line was for T1 = 308 K (DT = 9 K), the dashed line was for T1 = 305 K (DT = 6 K), and the dash-dotted line was for T1 = 302 K (DT = 3 K). From Figure , it may be concluded that, under the steady state, the increase in the thermal gradient was well exerting on those SET systems to promote molecular migrations during the whole thermal diffusion.

Figure 3

Distribution (regression) curves representing temperature spatial profiles in the middle region along the x-direction corresponding to three temperature differences (DT of 3, 6, and 9 K): (a) three SET-1 systems and (b) three SET-2 systems.

Figure shows three distributions for nC-5 mass fractions (Mf , also the relevant C) in the middle region along the x-direction corresponding to three DT values (DT = 3, 6, and 9 K) between hot- and cold-edge zones in three SET-1 systems, respectively. In this figure, every point was achieved from calculations in the Appendix. Three slopes representing distribution trends for the Mf were determined by the linear regression. From this figure, it can be seen that the larger the temperature difference (DT) on the system (being an effect of the thermal gradient), the steeper the Mf slope and the smaller the Mf fluctuations became, and thus, the bigger the driving effect on nC-5 molecular migrations during thermal diffusion. Overall, three Mf slopes corresponding to three DT values were shown to be declining gradually, meaning that more nC-5 molecules (with light molar mass) would migrate into the hot area (Area-1) than the cold one (Area-2) with the DT increase.

Figure 4

Distributions for nC-5 mass fractions (Mf) in the middle region in three SET-1 systems (u0 = 10–6): (a) DT = 3 K, (b) DT = 6 K, and (c) DT = 9 K.

Figure shows three distributions for nC-5 mass fractions (Mf , also the relevant C) in the middle region along the x-direction corresponding to three DT values (DT = 3, 6, and 9 K) between hot- and cold-edge zones in three SET-2 systems, respectively. In this figure, every point was achieved from calculations in the Appendix. Three slopes representing distribution trends for the Mf were determined by the linear regression. From this figure, it can also be seen that the larger the temperature difference (DT) on the system (being an effect of the thermal gradient), the steeper the Mf slope and the smaller the Mf fluctuations became, and thus, the bigger the driving effect on nC-5 molecular migrations during thermal diffusion. Overall, three Mf slopes corresponding to three DT values were shown to be declining gradually, meaning that more nC-5 molecules (with light molar mass) would migrate into the hot area (Area-1) than the cold one (Area-2) with the DT increase.

Figure 5

Distributions for nC-5 mass fractions (Mf) in the middle region in three SET-2 systems (u0 = 10–6): (a) DT = 3 K, (b) DT = 6 K, and (c) DT = 9 K.

Apparently, Figures and 5 may corroborate what the NEMD did before,[31−35] i.e., a thermal gradient (equivalent to the heat-flow exchange) had a positive impact on the Soret effect if it increased to higher than a certain level. Thus, increasing the thermal gradient could be fairly effective to the SC calculations because a linear response and constant phases in the MD simulations can be well ensured under this higher constraint, while they could hardly be well determined in relevant experimental works.

Mechanisms of thermal diffusion in Figures and 5 could be explained as follows: In each of the two binary liquid solutions here, from data in Table , it may receive the following: ε [nC-5⊕nC-7] = −0.06328 eV > [ε (nC-5) + ε (nC-7)] = −0.06740 eV and ε [nC-5⊕nC-10] = −0.06319 eV > [ε (nC-5) + ε (nC-10)] = −0.06733 eV, meaning that there were very weak hydrophobic associations between nC-5 and nC-7 molecules and between nC-5 and nC-10 molecules, even if these two binary liquid solutions looked uniform. In addition, geometric differences between nC-5 and nC-7 molecules and between nC-5 and nC-10 molecules were both more than 15%. Therefore, when the thermal gradient was increasingly applied onto each of the SET systems, all the molecules in the systems would perform their thermal migrations by means of severe molecular collisions.[25] For example, in Area-1 (near the hot-edge zone), it was quite easier for more nC-5 molecules (with light molar mass) than nC-7 and nC-10 molecules (both with heavy molar masses) to have their own thermal segregation. As a result, some of nC-7 or nC-10 molecules in Area-1 were squeezed out of this area but migrated toward inner sites of the middle region. Consequently, vacancies in Area-1 due to heavy molecular leaving would likely attract more nC-5 molecules from the middle region to fill in. Then, new vacancies would appear in the middle region due to nC-5 molecular leaving. Simultaneously, those heavy nC-7 and nC-10 molecules squeezed out of Area-1 but flowing into the middle region could partially fill in those new vacancies because of their obvious geometric differences with nC-5 molecules. Accordingly, the remaining new vacancies in the middle region would attract some other nC-5 molecules as well from Area-2 (near the cold-edge zone) to fill in. Following this migration path, some nC-7 or nC-10 molecules in the middle region were also squeezed out of this region but migrated toward Area-2. At the final steady state (after 500,000,000 time steps of the thermal constraint on the two edge zones), nC-5 molecules would be abundant in Area-1 (near hot Zone-1), while nC-7 and nC-10 molecules would be affluent in Area-2 (near cold Zone-2). Obviously, more nC-7 and more nC-10 molecules (both with heavy molar masses) in the SET systems were found to have their thermal migrations opposite to the nC-5 molecules with the DT increase, which were also observed in other works.[17−19]

Table 5

The Fitting Parameters to Represent van der Waals Interactions among A1, A2, and A3 Beadsb

bead type	potential model	bond strength ε (eV)	bond strength σ (Å)	bead molar mass (g/mol)
A1–A1a	L-J	–0.04381	5.1276	44.107–44.107
A1–A2	L-J	–0.14420	5.5289	44.107–58.138
A1–A3	L-J	–0.06300	5.5835	44.107–72.169
A2–A2a	L-J	–0.04459	5.9302	58.138–58.138
A2–A3	L-J	–0.06356	5.9848	58.138–72.169
A3–A3a	L-J	–0.09059	6.0395	72.169–72.169

Data came from ref (39).

Please note, there is no L-J interaction between A1 and A2 beads inside each single molecule.

Calculations of the Soret Coefficients

By using eq , three SC quantities for nC-5 molecules in three SET-1 systems corresponding to three DT values in Figure were calculated. In addition, three SC quantities for nC-5 molecules in three SET-2 systems corresponding to three DT values in Figure were calculated. Please note, here, T2 = 299 K was taken as a reference point, and the Mf quantities in Area-1 and Area-2 were selected for the SC calculations by means of the optimized interpolation.[35,39]Table lists relevant results for these SC quantities. From this table, it may be concluded that the SC quantities for nC-5 molecules in both the SET-1 and the SET-2 systems would decline gradually along with the DT increase.

Table 3

Three SC Quantities for SET-1 and SET-2 Systems Corresponding to Three Temperature Differencesa

temperature difference between hot- and cold-edge zones, DT (K)		3 (T1 = 302 K)	6 (T1 = 305 K)	9 (T1 = 308 K)
the Soret coefficient (×10–3, K–1)	SET-1	–3.17 ± 0.14	–2.36 ± 0.17	–2.08 ± 0.20
	SET-2	–3.45 ± 0.19	–2.61 ± 0.20	–2.33 ± 0.13

T2 = 299 K in cold Zone-2, while T1 in hot Zone-1 are listed in the table.

Still coming back to eq , along the x-direction, since the C was the T-dependent, an indefinite integral solution for this equation: may represent a function of the variable. To approve of this point, Figure shows two trends of the SC quantities for nC-5 molecules versus the temperature difference (DT = T – T0) on the SET-1 and the SET-2 systems. Relevant data in this figure came from Table , and T0 = 299 K was taken as the reference point. From Figure , it may be observed that two trends of the SC quantities for nC-5 molecules in the two SET systems were both inversely proportional to the DT value (corresponding to the thermal gradient). Also, from this figure, the SC on the SET-2 system was found to perform more obviously than on the SET-1 system.

Figure 6

Two trends of the SC quantities for nC-5 molecules corresponding to three temperature differences in the SET-1 and SET-2 systems.

In experimental works, the SC quantity for a binary liquid solution was usually represented by its absolute value. Thus, qualitatively, Dominguez et al. reported that an expression of the Soret coefficient (absolute value) for an isotopic fractional element (the melting silicate) under a thermal gradient showed a function of the inverse of the absolute temperature if choosing the zero point as the temperature reference point.[54] Essentials of such curve distributions in Figure could be explained as follows: According to the SC definition, , in a uniform liquid solution, the D usually became blunt if the thermal gradient did not severely impact the system (e.g., just ±1 K deviation at each temperature point in Figure may represent a very low phonon excitation). However, the D was fairly sensitive to the thermal gradient even if the latter just behaved a small alteration (see obvious Mf fluctuations in Figures and 5). On the whole, the ST quantity would be declining if the thermal gradient was increasing gradually. In addition, three other works also reported that:[55−57] if a multicomponent system was under an effect of a thermal gradient, then lighter fractional elements would concentrate in the hot region, while heavier elements would concentrate in the cold one. According to them, qualitatively, trends of the SC quantity versus the temperature difference, and trends of molecular migrations under influences of the thermal gradient on the SET systems here, would be in good agreement with those of the experimental observations.

Still, in experimental works, quantitatively, Leahy-Dios and Firoozabadi reported that an SC quantity (absolute value) for a binary liquid solution with equal-weight (50 wt %) nC5–nC10 mixtures would lay in [3.66, 4.06] × 10–3 K–1, mainly dependent upon different molecular shapes and sizes (normal and isomer molecular conformations) in the system.[53] In Perronace et al.’s work,[27] starting from 300 K, one SC quantity (absolute value) for a multicomponent C5–C10 system with equimolar mixtures would lay in [2.76, 4.42] × 10–3 K–1. In addition, starting from 298 K, Mounir Bou-Ali et al. determined the SC quantities for several binary liquid solutions mixing various n-alkane molecules.[58] In their works, the SC quantity for a binary solution with equimolar nC-5 and nC-10 mixtures was 2.96 × 10–3 K–1 if an obvious thermal flow passed through the whole system. Compared with the above results, for the SET-2 system, its SC quantities obtained from the MD simulations here, (2.33–3.45) × 10–3 K–1, would match those works fairly well.

Unfortunately, there were no relevant reports about the SC quantity for the SET-1 system. However, in ref (58), the SC quantities of four equimolar mixtures, nC-7 and nC-10, nC-10 and nC-14, nC-12 and nC-17, and nC-13 and nC-18, were determined. Their values were 1.87 × 10–3, 2.17 × 10–3, 2.40 × 10–3, and 2.35 × 10–3 K–1, respectively. Since the ratio of carbon atoms between two components in each of the four mixtures above was within [1.38, 1.42], all of which were near that in the SET-1 system (nC-5⊕nC-7), the SC quantities for the SET-1 system obtained from the MD works here, (2.08 ∼ 3.17) × 10–3 K–1, would also match those findings.

Conclusions

At the nanoscale, coarse-grained molecular dynamics (MD) was fulfilled to simulate the Soret effect (SE) on two binary liquid solutions with equimolar mixtures: n-pentane (nC-5) and n-heptane (nC-7) molecules plus n-pentane (nC-5) and n-decane (nC-10) molecules. Some results can be drawn as follows:

First, before the MD simulations, two coarse-grained force field (the CG-FF) potentials, which were based upon a model of new united atoms, were developed to depict intra-/intermolecular interactions in the above liquid solutions. As a result, they were quite reliable throughout the whole MD simulations.

Second, under an effect of a weak thermal gradient on the two systems, the MD simulations indicated that nC-5 molecules with light molar mass would migrate toward the high-temperature region, while nC-7 and nC-10 molecules with heavy molar masses would migrate toward the low-temperature region, both of which were also observed in other relevant works. In addition, the higher the thermal gradient on two systems, the more obvious the SE would be. Mechanisms of such phenomena were mainly due to very weak hydrophobic associations between two different types of n-alkane molecules in each of the binary liquid solutions.

Third, still under the effect of a weak thermal gradient, quantities of the Soret coefficient (SC) for the above binary liquid solutions were calculated by means of the developed CG-FF potentials. Calculations were found to match other relevant works fairly well. In addition, trends of the SC quantity versus the temperature difference on the systems looked inversely proportional.

Finally, the fact that the SE outputs corroborated fairly well with other relevant works here may support that our developed CG-FF potentials have opened a good perspective: in the future, they may play a very promising role in simulating other relevant properties of petroleum engineering.

Methodologies

Equations of Thermal Diffusion

In a cuboid reservoir containing a uniform liquid solution mixing various components, under an effect of a thermal gradient (temperature difference) from hot to cold boundary sides of the reservoir, mass fractions of these components would distribute accordingly.[2−4] Assuming that a uniform liquid system just consists of two components, thus, a comprehensive flux for component i should obey the following diffusion law[2]where T is a thermal temperature field exerting on the system, C is a mass fraction for component i in response to the T variable, ρ is a mass density for component i, and D and D represent two rates of thermal diffusion and mass diffusion for component i, respectively. When the system reaches the steady state (i.e., the J = 0), the T and the C would distribute steadily. In this case, eq could be converted as followswhere , denoting a definition of the Soret coefficient (SC) for component i, and its unit is K–1. In practice, if considering the SE along the x-direction only, then eq will have its projection on the x-axis as an ordinary differential equation, i.e.,

Generally, eq may well depict the characteristic of thermal diffusion in most of the macro systems with uniform liquids. However, if the size of a system is around the nanoscale and the thermal gradient exerting on the system is quite weak, then both the D and the D would perform gently.[42] In this case, the S could be viewed as a constant throughout the whole system. Therefore, along the x-direction, if setting up hot (T1)- and cold (T2)-edge zones at the left and right boundary sides of the system (T1 > T2) and naming C the mass fraction for component i in an area adjacent to the hot-edge zone and C the one adjacent to the cold-edge zone, by solving eq with the method of variable separation, then the S for component i can be derived as follows

Please note, at the final steady state, both the C and the C must have accumulated all information of the C evolution throughout the whole thermal diffusion.

Models Selected for the MD Simulations

In the MD simulations here, two uniform binary liquid solutions with equimolar mixtures: the SET-1 system with nC-5 and nC-7 molecules and the SET-2 system with nC-5 and nC-10 molecules, were selected for the studied objects.

As discussed in Section , both the AUA and the all-atom force field models usually consisted of two main termswhere Ebond represented strong intramolecular chemical bonds (atomic bonds in one molecular conformation) and Enonbond represented weak intermolecular interactions among all the molecular groups. In details, (a) the Ebond term consisted of three subtermswhere Eradial = kr(r – r0)2, a harmonic model denoting an effect of radial stretching (kr was a radial strength, and r0 was an initial radial quantity), Eangle = kθ(θ – θ0)2, a harmonic model denoting an effect of angular bending (kθ was an angular strength, and θ0 was an initial angular quantity), and Etorsion = Edihedral + Eoff-plane, denoting a deforming effect of spatial geometry. Here, kr, r0, kθ, and θ0 are empirical parameters.

(b) The Enonbond term consisted of two subtermswhere the EvdW denoted van der Waals interactions and the Eelectric denoted electrostatic interactions.

(c) For each of the electroneutral n-alkane molecules, both the Etorsion and the Eelectric contributions would become secular. Thus, eq can be simplified as follows

Equation may represent the CG-FF model in this article. Overall, from this model, one nC-7 molecule was spliced by an A1–A2 group, one nC-10 molecule was spliced by an A1–A2–A1 group, and one nC-5 molecule was depicted by an A3 bead only.[39,42] By means of the DFT calculations,[43−45] the empirical parameters kr, r0, kθ, and θ0 in the Eradial and the Eangle terms were determined by using the Birch–Murnaghan equation of state (BMES) at its harmonic stage.[46,47]Table lists all the fitting parameters for these two terms in A1–A2 and A1–A2–A1 groups as shown in Figure .

Table 4

The Fitting Parameters for the Bonding Status in the Coarse-Grained Beads

		bond strength		bond geometry
bead type	potential model	radial kr [eV/Å2]	angular kθ [eV/(°)2]	length r0 [eV/Å2]	angular θ0 [°]	bead molar mass [g/mol]
A1–A2	the CG-FF	28.998		4.193		44.107–58.138
A1–A2–A1	the CG-FF	28.998	0.4129	4.193	179.6	44.107–58.138–44.107
another calculationa		25.911	0.5183	4.180	180.0

Data came from ref (39).

Still, in eq , the EvdW term was represented by an L-J potential,, where ε and σ were empirical parameters. Table lists all the fitting parameters in those EvdW terms for A1, A2, and A3 bead interactions (intermolecular potentials) as shown in Figure . The cutoff distance of the potential tail in each of the EvdW terms was chosen as 12.00 Å. In Table , values of ε and σ for A1–A1, A2–A2, and A3–A3 beads were chosen from refs (39) and (47−49). However, ε values for A1–A2, A1–A3, and A2–A3 beads were obtained from those for A1–A1, A2–A2, and A3–A3 beads by using the geometric mean method. Also, σ values for A1–A2, A1–A3, and A2–A3 beads were obtained from those for A1–A1, A2–A2, and A3–A3 beads by using the arithmetic mean method. Regarding validations of these mean methods, see refs (50) and (51) .