Multicomponent Shale Oil Flow in Real Kerogen Structures via Molecular Dynamic Simulation

: As an unconventional energy source, the development of shale oil plays a positive role in global energy, while shale oil is widespread in organic nanopores. Kerogen is the main organic matter component in shale and a ﬀ ects the ﬂow behaviour in nanoscale-conﬁned spaces. In this work, a molecular dynamic simulation was conducted to study the transport behaviour of shale oil within kerogen nanoslits. The segment ﬁtting method was used to characterise the velocity and ﬂow rate. The heterogeneous density distributions of shale oil and its di ﬀ erent components were assessed, and the e ﬀ ects of di ﬀ erent driving forces and temperatures on its ﬂow behaviours were examined. Due to the scattering e ﬀ ect of the kerogen wall on high-speed ﬂuid, the heavy components (asphaltene) increased in bulk phase regions, and the light components, such as methane, were concentrated in boundary layers. As the driving force increased, the velocity proﬁle demonstrated plug ﬂow in the bulk regions and a half-parabolic distribution in the boundary layers. Increasing the driving force facilitated the desorption of asphaltene on kerogen walls, but increasing the temperature had a negative impact on the ﬂow velocity.


Introduction
The development of shale oil has ameliorated global energy shortages, and many countries have launched programmes to investigate various development approaches [1][2][3]. However, the complex structure of nanopores in shale and the different components in shale oil hinder the further study of shale oil flow behaviours [4,5].
The sizes of pores in shale range from nanometres to micrometres, but experimental nanoscale studies are restricted to laboratory assessments [6][7][8][9][10][11][12]. Molecular dynamic (MD) simulations are often used to study the fluid behaviours within nano-confined spaces [13]. MD simulations are based on Newtonian mechanics and can be used to calculate a system's thermodynamic quantities and other macroscopic properties. Shale is mainly composed of organic and inorganic matter [14][15][16]. The inorganic matter contains quartz and clay minerals. As a source of oil, kerogen is the main component of organic matter [17]. Thus, the flow behaviour of shale oil in realistic kerogen channels needs to be investigated.
Bousige et al. [18] produced a realistic molecular kerogen model using the generalised phonon densities of states. Ungerer et al. [19] classified kerogen molecules into four types of maturity and constructed corresponding structures. It is difficult for shale oil to flow in a kerogen matrix as a result of kerogen's dense matrix structure [20], so the reasonable construction of kerogen channels is necessary. During the kerogen matrix generation process, the insertion of repulsive dummy particles can produce porous models [21,22]. Perez et al. [8] blended a fluid mixture with kerogen monomers and simulated the annealing process, and the fluid molecules gathered as clusters. Kerogen monomers form organic frameworks. However, it is difficult to obtain clear transport characteristics in irregular channels. Carbon nanotubes (CNTs) and graphene slits are common flow channel models [7,23], but the flow rate can be overestimated on ultra-smooth surfaces. In multicomponent fluids, heterogeneous density distribution characteristics cannot be produced on smooth surfaces. Due to the driving force's effect, the velocity profile tends to plug in smooth channels rather than parabolas on rough walls [24]. Thus, constructing flow channels with kerogen is a more reasonable choice.
The flow channel's boundary condition affects the slip length [25]. Due to the presence of slippage, Darcy's law is invalid in nanoscale fluid flows [20,26]. The slippage effect can be characterised using the slip length concept, which is defined as the distance between the solid surface and the point where the velocity extrapolation is zero [27]. Martini et al. [27] found, at a low driving force, the movement of only a few molecules along a wall, called a "defect slip" in a nonlinear mode. At a high driving force, all of the fluid molecules adjacent to the liquid-solid interface contribute to "global slip." Conversely, on a rough kerogen wall, the boundary's cavity can trap some molecules. In addition, fluid molecules impinge protrusions on the interface, decelerating the velocity in the boundary layer [28][29][30]. However, to the best of our knowledge, very few reports in the literature discuss multicomponent shale oil slippage.
In this study, we constructed a kerogen slit [31][32][33] and used a nonequilibrium molecular dynamic (NEMD) simulation to analyse the flow behaviour of shale oil. We also examined the effects of the driving force and temperature, and the simulation results were fitted using the Poiseuille formula by dividing the channel into two sections. We calculated the flow rate using boundary layers under multicomponent shale oil conditions.

Molecular Models
In this study, type Ⅱ-C kerogen monomers were used to construct realistic organic slits because the Ⅱ-C kerogen tended to be present in organic-rich shales [19]. We used a simplified shale oil component model including different molecular weight models (C1, C4, C8, C12, methylbenzene, asphaltene) [8,33,34]. The details and structures of the kerogen molecules and shale oil models are provided in the Supporting Information, and the model of the slit system is shown in Figure 1. There were 26,899 atoms in the entire system. The dimensions of the simulation box were 40 Å 3 × 50 Å 3 × 200 Å 3 , and each kerogen matrix contained 15 kerogen monomers. Figure 1. Molecular model of the kerogen slits constructed with a kerogen matrix and multicomponent shale oil. Molecular colours: red, kerogen; black, asphaltene; grey, methylbenzene; purple, n-dodecane; mauve, n-octane; pink, n-butane; green, propane; yellow, ethane; blue, methane. slits constructed with a kerogen matrix and multicomponent shale oil. Molecular colours: red, kerogen; black, asphaltene; grey, methylbenzene; purple, n-dodecane; mauve, n-octane; pink, n-butane; green, propane; yellow, ethane; blue, methane. We used a polymer consistent force field (PCFF) for all of the molecules [19,35,36], and described the Van der Waals interactions using Lennard-Jones (LJ) potentials. The parameters between different molecules were calculated using Waldman-Hagler combining rules [37]. The Ewald method was adopted to compute the electrostatic potential [38]. The boundary conditions in three directions were periodic, and the cut-off distance was 12 Å.

Simulation Methods
The large-scale atomic/molecular massively parallel simulator (LAMMPS) package distributed by Sandia National Laboratories was used for the molecular simulations [39]. The MD simulation in the NVE (The number of atoms, volume and energy of system remain constant.) ensemble was conducted for 50 ps to relax the system's configuration. The NPT (The number of atoms, pressure and temperature of system remain constant.) ensemble was then performed (P = 30.0 MPa and T = 300 K) for 100 ps. We put the simulation system into the NVT (The number of atoms, volume and temperature of system remain constant.) ensemble at T = 300 K for 200 ps to achieve an equilibrium system. The timesteps of the NVE, NPT, and NVT ensemble simulations were 0.5 fs, 0.5 fs, and 1.0 fs, respectively.
We conducted NEMD simulations with the driving force along the streaming direction. A constant force was applied to each atom rather than a pressure gradient because longitudinally homogeneous pressure can be maintained with a constant field [40]. The NEMD simulations were conducted in the NVT ensemble at 300 K for 18 ns, and we collected the last 6 ns of trajectories for analysis. The temperature was maintained by a Nosé-Hoover thermostat perpendicular to the flow direction [41,42], which reduced the effect of thermal motion on the velocity in the flow direction. During the simulation, as shown in Figure 2, we froze the kerogen matrix layer (KML) and controlled the temperature (300 K) of the rough adsorption layer (RAL) [29,43], which made the kerogen interface flexible [44]. We discussed these layers in the Section 3.1. The free slit layer (FSL) was divided into a bulk phase region and boundary layer, which are discussed in Section 3.3. By averaging the velocity and number of atoms in the bins that were parallel to the flow direction, we obtained the velocity and density profiles [27]. To reduce errors, the simulations were calculated three times independently with the same initial condition.
Energies 2020, 13, x FOR PEER REVIEW 3 of 12 We used a polymer consistent force field (PCFF) for all of the molecules [19,35,36], and described the Van der Waals interactions using Lennard-Jones (LJ) potentials. The parameters between different molecules were calculated using Waldman-Hagler combining rules [37]. The Ewald method was adopted to compute the electrostatic potential [38]. The boundary conditions in three directions were periodic, and the cut-off distance was 12 Å.

Simulation Methods
The large-scale atomic/molecular massively parallel simulator (LAMMPS) package distributed by Sandia National Laboratories was used for the molecular simulations [39]. The MD simulation in the NVE (The number of atoms, volume and energy of system remain constant.) ensemble was conducted for 50 ps to relax the system's configuration. The NPT (The number of atoms, pressure and temperature of system remain constant.) ensemble was then performed (P = 30.0 MPa and T = 300 K) for 100 ps. We put the simulation system into the NVT (The number of atoms, volume and temperature of system remain constant.) ensemble at T = 300 K for 200 ps to achieve an equilibrium system. The timesteps of the NVE, NPT, and NVT ensemble simulations were 0.5 fs, 0.5 fs, and 1.0 fs, respectively.
We conducted NEMD simulations with the driving force along the streaming direction. A constant force was applied to each atom rather than a pressure gradient because longitudinally homogeneous pressure can be maintained with a constant field [40]. The NEMD simulations were conducted in the NVT ensemble at 300 K for 18 ns, and we collected the last 6 ns of trajectories for analysis. The temperature was maintained by a Nosé-Hoover thermostat perpendicular to the flow direction [41,42], which reduced the effect of thermal motion on the velocity in the flow direction. During the simulation, as shown in Figure 2, we froze the kerogen matrix layer (KML) and controlled the temperature (300 K) of the rough adsorption layer (RAL) [29,43], which made the kerogen interface flexible [44]. We discussed these layers in the Section 3.1. The free slit layer (FSL) was divided into a bulk phase region and boundary layer, which are discussed in Section 3.3. By averaging the velocity and number of atoms in the bins that were parallel to the flow direction, we obtained the velocity and density profiles [27]. To reduce errors, the simulations were calculated three times independently with the same initial condition.

Heterogeneous Density Distribution Analysis
As shown in Figure 2, to more clearly describe the flow behaviour, we divided the system into three regions according to kerogen's density profile. The KML density fluctuated at 1.1-1.2 g·cm −3 [18,22,32]. The FSL was defined as the free space without kerogen, and the RAL was comprised of the cavities and protrusions of the kerogen between the KML and FSL. The widths of the FSL, RAL, and KML were 5 nm, 1.5 nm, and 2.5 nm, respectively. The roughness of the kerogen matrix's two walls differed as a result of its amorphous structure, so we only used half of the system to analyse the shale oil fluid's properties. Due to the random molecular thermal motion effects, there were always slight errors in the shale oil's distribution. The kerogen wall friction also played a role in the difference between the two halves, but the friction was not our focus.
Using NEMD simulations, we simulated the shale oil flow through the kerogen slit under the driving force's effect. Figure 3a shows that the shale oil density profile varied considerably at different driving forces. There were two density peaks in the boundary layer at a low driving force, consistent with our previous research [33]. The density profiles presented as a single peak in the middle of the slit at a larger driving force. The difference in the density was small below a driving force of 5 × 10 −4 kcal/(mol·Å). Figure 3b shows that the methane density distribution peaked in the RAL region, and the asphaltene's density peak was distributed in the FSL.

Heterogeneous Density Distribution Analysis
As shown in Figure 2, to more clearly describe the flow behaviour, we divided the system into three regions according to kerogen's density profile. The KML density fluctuated at 1.1-1.2 g·cm −3 [18,22,32]. The FSL was defined as the free space without kerogen, and the RAL was comprised of the cavities and protrusions of the kerogen between the KML and FSL. The widths of the FSL, RAL, and KML were 5 nm, 1.5 nm, and 2.5 nm, respectively. The roughness of the kerogen matrix's two walls differed as a result of its amorphous structure, so we only used half of the system to analyse the shale oil fluid's properties. Due to the random molecular thermal motion effects, there were always slight errors in the shale oil's distribution. The kerogen wall friction also played a role in the difference between the two halves, but the friction was not our focus.
Using NEMD simulations, we simulated the shale oil flow through the kerogen slit under the driving force's effect. Figure 3a shows that the shale oil density profile varied considerably at different driving forces. There were two density peaks in the boundary layer at a low driving force, consistent with our previous research [33]. The density profiles presented as a single peak in the middle of the slit at a larger driving force. The difference in the density was small below a driving force of 5 × 10 −4 kcal/(mol·Å). Figure 3b shows that the methane density distribution peaked in the RAL region, and the asphaltene's density peak was distributed in the FSL. Due to the high flow speed caused by the driving force, the shale oil molecules collided with the kerogen protrusion and had a diffuse scattering behaviour, forming a velocity component perpendicular to the streaming direction [43,45]. The perpendicular velocity component enriched the shale oil in the bulk phase, as shown in Figure 3b, especially in the asphaltene with its complex and large molecular structure. The diffuse scattering reduced the density in the interface [43].
Notably, the density profile remained steady under the different driving forces on the smooth walls [7], and it was difficult to obtain a heterogeneous distribution because the smooth wall could not provide a vertical force from the wall to the fluid. This characteristic was similar to blood flow in the blood vessels; haemoglobin concentrates in the central area and ionic liquid concentrates at the boundary [46,47]. By comparing the results of different driving forces, that the density profiles of asphaltene and methane were opposite, which was important for calculating the flow rate. Due to the high flow speed caused by the driving force, the shale oil molecules collided with the kerogen protrusion and had a diffuse scattering behaviour, forming a velocity component perpendicular to the streaming direction [43,45]. The perpendicular velocity component enriched the shale oil in the bulk phase, as shown in Figure 3b, especially in the asphaltene with its complex and large molecular structure. The diffuse scattering reduced the density in the interface [43].
Notably, the density profile remained steady under the different driving forces on the smooth walls [7], and it was difficult to obtain a heterogeneous distribution because the smooth wall could not provide a vertical force from the wall to the fluid. This characteristic was similar to blood flow in the blood vessels; haemoglobin concentrates in the central area and ionic liquid concentrates at the boundary [46,47]. By comparing the results of different driving forces, that the density profiles of asphaltene and methane were opposite, which was important for calculating the flow rate.

Effect of the Driving Force on the Flow Behaviour in the Kerogen Slit
To investigate the driving force's effects, we conducted the NEMD at driving forces from 15 × 10 −4 kcal/(mol·Å) to 40 × 10 −4 kcal/(mol·Å). Unlike conventional parabola and plug profiles [36], Figure 4a shows that the velocity distribution was a trapezoid with two arched waists at high driving forces. The heterogeneity of the shale oil caused an interesting flow conformation. According to the previous density analysis, the asphaltene concentrated in the bulk region, while the methane enriched the boundary layer at high driving forces. In this case, as a result of cohesion, the asphaltene formed clusters with a molecular pore network and trapped some small molecules. The linear alkanes sheared easily. The orientation property of the alkanes was observed by solation force experiments [48,49]. However, it was hard to change the configuration of the asphaltene cluster with the shear friction, so the bulk velocity was flat.

Effect of the Driving Force on the Flow Behaviour in the Kerogen Slit
To investigate the driving force's effects, we conducted the NEMD at driving forces from 15 × 10 −4 kcal/(mol·Å) to 40 × 10 −4 kcal/(mol·Å). Unlike conventional parabola and plug profiles [36], Figure 4a shows that the velocity distribution was a trapezoid with two arched waists at high driving forces. The heterogeneity of the shale oil caused an interesting flow conformation. According to the previous density analysis, the asphaltene concentrated in the bulk region, while the methane enriched the boundary layer at high driving forces. In this case, as a result of cohesion, the asphaltene formed clusters with a molecular pore network and trapped some small molecules. The linear alkanes sheared easily. The orientation property of the alkanes was observed by solation force experiments [48,49]. However, it was hard to change the configuration of the asphaltene cluster with the shear friction, so the bulk velocity was flat.  In the boundary layer, the methane enrichment made the velocity profile approximately parabolic. Based on the platform and the two halves of the parabolic velocity profile, we divided the FSL into two regions (the bulk phase and boundary layers) and fit the velocity and flow rate [50]. The platform region was the bulk phase (0-12.5 Å), and the halves of the parabola were the boundary layers (12.5-25 Å). Consequently, using the Poiseuille equation, we quantitatively analysed the flow behaviour.
The velocity in the boundary layer increased homogeneously over forces of 30 × 10 −4 kcal/(mol·Å), as shown in Figure 4b. By increasing the driving force, the fluid molecules had more kinetic energy, leading to stronger scattering effects. Figure 5 demonstrates that the driving force had a positive effect on the asphaltene density in the bulk phase, and the density profile stabilised beyond 30 × 10 −4 kcal/(mol·Å). Due to the strong scattering, the asphaltene density peaks in the RAL were diminished. However, the asphaltene still adsorbed in the RAL and formed density peaks at driving forces below 30 × 10 −4 kcal/(mol·Å). Therefore, the constant distribution of the asphaltene density facilitated a steady increase in velocity. In the boundary layer, the methane enrichment made the velocity profile approximately parabolic. Based on the platform and the two halves of the parabolic velocity profile, we divided the FSL into two regions (the bulk phase and boundary layers) and fit the velocity and flow rate [50]. The platform region was the bulk phase (0-12.5 Å), and the halves of the parabola were the boundary layers (12.5-25 Å). Consequently, using the Poiseuille equation, we quantitatively analysed the flow behaviour.
The velocity in the boundary layer increased homogeneously over forces of 30 × 10 −4 kcal/(mol·Å), as shown in Figure 4b. By increasing the driving force, the fluid molecules had more kinetic energy, leading to stronger scattering effects. Figure 5 demonstrates that the driving force had a positive effect on the asphaltene density in the bulk phase, and the density profile stabilised beyond 30 × 10 −4 kcal/(mol·Å). Due to the strong scattering, the asphaltene density peaks in the RAL were diminished. However, the asphaltene still adsorbed in the RAL and formed density peaks at driving forces below 30 × 10 −4 kcal/(mol·Å). Therefore, the constant distribution of the asphaltene density facilitated a steady increase in velocity.

Effect of the Temperature on the Flow Behaviour in the Kerogen Slit
We compared the velocity profiles of the shale oil subjected to different temperatures, as shown in Figure 6a. In general, the fluid had a higher velocity as the temperature increased because the stronger thermal motion facilitated the reduction in viscosity. However, the velocity decreased and the asphaltene's density profiles were gentler as the temperature increased, as shown in Figure 6. The higher temperature of the single component fluid reduced its viscosity. In the confined nanoslits, the temperature affected the shale oil's heterogeneous density distribution. The asphaltene cluster behaved as a complex molecular network. The stronger thermal motion increased the distance between the asphaltene atoms, expanding the asphaltene cluster and increasing its concentration in the boundary layer, as shown in Figure 6b. Thus, the higher temperature weakened the density distribution heterogeneity, and the fluid in the boundary layer was mixed with more asphaltene molecules, increasing the friction between the kerogen and asphaltene. Therefore, the high temperature had a negative effect on the shale oil flow, with asphaltene clusters in the nanochannels.

Effect of the Temperature on the Flow Behaviour in the Kerogen Slit
We compared the velocity profiles of the shale oil subjected to different temperatures, as shown in Figure 6a. In general, the fluid had a higher velocity as the temperature increased because the stronger thermal motion facilitated the reduction in viscosity. However, the velocity decreased and the asphaltene's density profiles were gentler as the temperature increased, as shown in Figure 6.

Effect of the Temperature on the Flow Behaviour in the Kerogen Slit
We compared the velocity profiles of the shale oil subjected to different temperatures, as shown in Figure 6a. In general, the fluid had a higher velocity as the temperature increased because the stronger thermal motion facilitated the reduction in viscosity. However, the velocity decreased and the asphaltene's density profiles were gentler as the temperature increased, as shown in Figure 6. The higher temperature of the single component fluid reduced its viscosity. In the confined nanoslits, the temperature affected the shale oil's heterogeneous density distribution. The asphaltene cluster behaved as a complex molecular network. The stronger thermal motion increased the distance between the asphaltene atoms, expanding the asphaltene cluster and increasing its concentration in the boundary layer, as shown in Figure 6b. Thus, the higher temperature weakened the density distribution heterogeneity, and the fluid in the boundary layer was mixed with more asphaltene molecules, increasing the friction between the kerogen and asphaltene. Therefore, the high temperature had a negative effect on the shale oil flow, with asphaltene clusters in the nanochannels. The higher temperature of the single component fluid reduced its viscosity. In the confined nanoslits, the temperature affected the shale oil's heterogeneous density distribution. The asphaltene cluster behaved as a complex molecular network. The stronger thermal motion increased the distance between the asphaltene atoms, expanding the asphaltene cluster and increasing its concentration in the boundary layer, as shown in Figure 6b. Thus, the higher temperature weakened the density distribution heterogeneity, and the fluid in the boundary layer was mixed with more asphaltene molecules, increasing the friction between the kerogen and asphaltene. Therefore, the high temperature had a negative effect on the shale oil flow, with asphaltene clusters in the nanochannels.

Slip Length Analysis of the Transport Condition
The slip occurred in the nanochannels, which was confirmed by many studies [25,[51][52][53]. The slip length was affected by the slip wall location, so we selected the slip boundary at the interface between the FSL and RAL [54]. Because the kerogen protrusion created extra flow resistance in the RAL, destabilising the flow, the flow in the FSL was only affected by the viscosity and driving force.
We calculated the slip length using L s = u surf /(du/dr), where L s is the slip length, u surf is the velocity at the slip wall, and du/dr is the velocity gradient at the interface [55]. As shown in Figure 7, the slip length increased steadily over a driving force of 30 × 10 −4 kcal/(mol·Å). Below a driving force of 25 × 10 −4 kcal/(mol·Å), the slip length increased unstably. Because the asphaltene was still adsorbed at the boundary at low driving forces, the adsorption affected the slip at the boundary.

Slip Length Analysis of the Transport Condition
The slip occurred in the nanochannels, which was confirmed by many studies [25,[51][52][53]. The slip length was affected by the slip wall location, so we selected the slip boundary at the interface between the FSL and RAL [54]. Because the kerogen protrusion created extra flow resistance in the RAL, destabilising the flow, the flow in the FSL was only affected by the viscosity and driving force.
We calculated the slip length using Ls = usurf/(du/dr), where Ls is the slip length, usurf is the velocity at the slip wall, and du/dr is the velocity gradient at the interface [55]. As shown in Figure 7, the slip length increased steadily over a driving force of 30 × 10 −4 kcal/(mol·Å). Below a driving force of 25 × 10 −4 kcal/(mol·Å), the slip length increased unstably. Because the asphaltene was still adsorbed at the boundary at low driving forces, the adsorption affected the slip at the boundary.

Segment Fitting on the Flow Behaviour in the Kerogen Slit
Due to the extremely heterogeneous density and velocity, it was difficult to obtain a constant boundary layer viscosity. We attained a linear viscosity distribution by separating the regions and calculating the shear stress and shear velocity, respectively, in the boundary layer using Equations (1)-(4) [56,57], where n is the fluid density, z is the distance from the centre of the slit, P is the shear stress, F is the driving force, γ is the shear velocity, u is the velocity, and η is the viscosity. We calculated a series of

Segment Fitting on the Flow Behaviour in the Kerogen Slit
Due to the extremely heterogeneous density and velocity, it was difficult to obtain a constant boundary layer viscosity. We attained a linear viscosity distribution by separating the regions and calculating the shear stress and shear velocity, respectively, in the boundary layer using Equations (1)-(4) [56,57], where n is the fluid density, z is the distance from the centre of the slit, P is the shear stress, F is the driving force, γ is the shear velocity, u is the velocity, and η is the viscosity. We calculated a series of shale oil density values in parallel bins of 4 × 5 × w bin nm 3 using n(z i ) = N i /(4 × 5 × w bin ), where w bin is the width of the bins and N i is the number of atoms in the bins. The viscosity distribution is expressed as where a = w/2 − b, as shown in Figure 8. a is half of the length of the half-bulk phase layer, b is half of the length of the boundary layer, and w is the slit width. N bin is the number of bins, η(z a ) is the viscosity on the border between the bulk phase layer and boundary layer, and η(z w ) is the viscosity at the interface between the FSL and RAL. We obtained the viscosity distribution in the boundary layer by averaging the difference between η(z a ) and η(z w ) as expressed in Equation (4).  The velocity of the Poiseuille flow modified by the slip is expressed in Equation (5) where the w is the slit width. The velocity in the boundary layer is expressed in Equation (6).
Equation (7) expresses the velocity of the bulk phase layer where z is equal to a.
Although the linear distribution of the viscosity cannot perfectly describe the real viscosity, the velocity profile presented in Figure 8 calculated using the NEMD simulations fully fit with Equation (6) and Equation (7).
The FSL flow rate is integrated using Equation (6) and Equation (7), as expressed in Equation (8): Substituting Equation (4) into Equation (8) obtains The velocity of the Poiseuille flow modified by the slip is expressed in Equation (5) [7], where the w is the slit width. The velocity in the boundary layer is expressed in Equation (6).
Equation (7) expresses the velocity of the bulk phase layer where z is equal to a.
Although the linear distribution of the viscosity cannot perfectly describe the real viscosity, the velocity profile presented in Figure 8 calculated using the NEMD simulations fully fit with Equation (6) and Equation (7).
The FSL flow rate is integrated using Equation (6) and Equation (7), as expressed in Equation (8): Substituting Equation (4) into Equation (8) obtains where Q is the flow rate. Compared with the flow rate using the NEMD simulations, as demonstrated in Equation (10), Figure 9 shows that the separate fitting worked well for the heterogeneous shale oil flow.
Energies 2020, 13, 3815 9 of 12 in Equation (10), Figure 9 shows that the separate fitting worked well for the heterogeneous shale oil flow. For the different interfaces in the kerogen matrix, we selected half of the system to fit the flow rate. Figure 9 demonstrates that the monotonically increasing flow rate was a result of the increased driving force. However, the different slopes in the flow rate vs. the driving force at 30 × 10 −4 kcal/(mol·Å) were an unexpected finding.
As shown in Figure 5, the density peak of the asphaltene in the RAL slowly decreased, suggesting that the asphaltene molecules desorbed from the RAL and flowed into the bulk phase below 30 × 10 −4 kcal/(mol·Å). Figure 9 demonstrates that the flow rate increased at a higher rate below the driving force of 30 × 10 −4 kcal/(mol·Å), so the stronger driving force facilitated the recovery of asphaltene in the boundary layer.

Conclusions
In this work, using MD simulations, we studied the flow behaviour of multicomponent shale oil in realistic kerogen nanoslits. Compared with previous static adsorption results, we found opposite density distributions for asphaltene and methane. The asphaltene increased in the bulk phase layer, while the methane concentrated in the boundary layer. Due to kerogen's effects, high-velocity asphaltene more easily scattered in the bulk phase. The fast transportation of shale oil demonstrated a different flow behaviour. The asphaltene group was constrained in the bulk phase region and behaved as a plug flow, and the velocity profile in the boundary layer presented as a half of a parabola.
The driving force had a positive effect on the shale oil velocity. The density distribution of asphaltene remained steady beyond the driving force of 30 × 10 −4 kcal/(mol·Å), while it showed a slightly denser peak at a smaller driving force. As the temperature increased, due to the stronger thermal motion, the asphaltene's density distribution was relatively gentle. Thus, the boundary layer viscosity increased and then the velocity decreased in the slit. This characteristic is unusual in macro-flow, but the asphaltene cluster cohesion increased the internal friction of the shale oil in the nano-confined channel.
According to the different velocity profiles in the bulk phase and boundary layer, we employed the Poiseuille equation to fit the velocity and flow rate, respectively. As a result of the asphaltene desorption, the slope of the flow rate maintained a higher value below the driving force of 30×10 -4 For the different interfaces in the kerogen matrix, we selected half of the system to fit the flow rate. Figure 9 demonstrates that the monotonically increasing flow rate was a result of the increased driving force. However, the different slopes in the flow rate vs. the driving force at 30 × 10 −4 kcal/(mol·Å) were an unexpected finding.
As shown in Figure 5, the density peak of the asphaltene in the RAL slowly decreased, suggesting that the asphaltene molecules desorbed from the RAL and flowed into the bulk phase below 30 × 10 −4 kcal/(mol·Å). Figure 9 demonstrates that the flow rate increased at a higher rate below the driving force of 30 × 10 −4 kcal/(mol·Å), so the stronger driving force facilitated the recovery of asphaltene in the boundary layer.

Conclusions
In this work, using MD simulations, we studied the flow behaviour of multicomponent shale oil in realistic kerogen nanoslits. Compared with previous static adsorption results, we found opposite density distributions for asphaltene and methane. The asphaltene increased in the bulk phase layer, while the methane concentrated in the boundary layer. Due to kerogen's effects, high-velocity asphaltene more easily scattered in the bulk phase. The fast transportation of shale oil demonstrated a different flow behaviour. The asphaltene group was constrained in the bulk phase region and behaved as a plug flow, and the velocity profile in the boundary layer presented as a half of a parabola.
The driving force had a positive effect on the shale oil velocity. The density distribution of asphaltene remained steady beyond the driving force of 30 × 10 −4 kcal/(mol·Å), while it showed a slightly denser peak at a smaller driving force. As the temperature increased, due to the stronger thermal motion, the asphaltene's density distribution was relatively gentle. Thus, the boundary layer viscosity increased and then the velocity decreased in the slit. This characteristic is unusual in macro-flow, but the asphaltene cluster cohesion increased the internal friction of the shale oil in the nano-confined channel.
According to the different velocity profiles in the bulk phase and boundary layer, we employed the Poiseuille equation to fit the velocity and flow rate, respectively. As a result of the asphaltene desorption, the slope of the flow rate maintained a higher value below the driving force of 30 × 10 −4 kcal/(mol·Å). Calculating the viscosity of multicomponent shale oil remains challenging, especially in nano channels. This study proposed a new perspective on shale oil transport.