Get Complete Project Material File(s) Now! »

## Coarse-grained-Monte Carlo simulation (CG-MC)

Monte Carlo methods are based on a random sampling of a set of configurations, noted 𝛤𝑖, of a molecular system. A large number of representative configurations are generated in order to calculate the average properties of the system at equilibrium. In the algorithm introduced by Metropolis et al. [65], sampling is not totally random but bias are introduced to favor the most likely configurations. The probability of a configuration is calculated from its Boltzmann probability, 𝑃𝐵𝑜𝑙𝑡𝑧(𝛤), and depends on the statistical ensemble chosen for the simulations. Statistical ensembles used in this work for CG-MC simulations are described in greater details in the next sections.

In the Metropolis algorithm, configurations are generated using a Markov chain. A Markov chain must respect the condition of reversibility of a system. It means that the appearance rate of a configuration 𝛤 must be equal to the disappearance rate of this configuration. Thus, the transition rate, Π(𝛤→𝛤′), from a configuration 𝛤 to 𝛤′ is equal to the transition rate, Π(𝛤′→𝛤), from a configuration 𝛤′ to 𝛤, as shown in the following expression:

Π(𝛤→𝛤′)=Π(𝛤′→𝛤) ∀ {𝛤,𝛤′}.

### Simulations boxes and boundary conditions

DPD and CG-MC simulations are carried out in orthorhombic simulation box for which boundary conditions are added. Two types of boundary conditions can be imposed:

– The walls of the simulation box are impassable and beads are forced to remain in the space initially defined.

– Periodic boundary conditions are applied allowing simulating an infinite space in which beads are free to move.

In the second case, simulation box is considered as the primitive cell and it is reproduced in the three dimensions of space. Thus, when a bead passes through one side on the primitive cell, it reappears on the opposite side with the same velocity.

#### Coarse-grained units

Simulations based on a coarse-grained model are performed using reduced units. The reference variables are the mass which is fixed at 𝑚̅=1 for all beads, the cutoff radius which is equal to 𝑟̅ 𝑐=1 and the energy which is fixed such that 𝑘𝐵𝑇=1. In this document, all variables expressed in reduced units are indicated by an upper bar (𝑎̅). Then, the other variables can be converted in the international system (SI) from the references. A table that summarizes the conversion between the reduced units and the units in the international system for thermodynamic and physicochemical properties is given in Table 12.

**State of the art of parameterization of coarse-grained simulations**

The most widely used method to determine interaction parameters between two identical beads (like beads) has been developed in 1997 by Groot and Warren [63]. They proposed a relationship between the isothermal compressibility and the interaction parameter between like beads as shown in equations (52) and (53). 𝑎𝑖𝑖=𝜅−1−12𝛼𝜌̅𝑘𝐵𝑇, (52) with 𝜅−1=1𝑛𝑘𝐵𝑇𝜅𝑇, (53).

where 𝜅−1 is the dimensionless isothermal compressibility, 𝜌̅ the DPD number density, 𝑛 the number density of molecules and 𝛼 is a constant equal to 0.101 as proposed by Groot and Warren. Note that, in the DPD model, interaction parameters are temperature dependent. The dimensionless isothermal compressibility of water 𝜅𝑤𝑎𝑡𝑒𝑟−1 under ambient conditions with 𝜌̅=3 is approximately equal to 16 which leads to the value 𝑎𝑖𝑖=25𝑘𝐵𝑇 commonly used in the literature. However, alternative formulations introduced the degree of coarse-graining, 𝑁𝑚, which is the number of water molecules grouped within a DPD bead, to express the interaction parameter between like beads: 𝑎𝑖𝑖=𝜅−1𝑁𝑚−12𝛼𝜌̅𝑘𝐵𝑇. (54).

For example using 𝑁𝑚=3, the interaction parameter between two water beads is 𝑎𝑖𝑖=78.3𝑘𝐵𝑇. In this chapter, the two approaches will be compared in order to highlight the effects of parameters 𝑁𝑚 on the reproduction of the miscibility of compounds and the variation of IFT on liquid-liquid equilibria.

Interaction parameters between different particles (unlike beads) have been derived from several macroscopic properties in the literature. Groot and Warren [63] have related interaction parameters with the Flory-Huggins (FH) parameters (𝜒) from the Flory-Huggins theory of polymer solutions [96, 97]. However, this approach is based on an important approximation, it requires that all liquids have approximatively the same isothermal compressibility (𝑎𝑖𝑖=𝑎𝑗𝑗). In this way, interaction between like beads can be considered as the reference energy, so that interactions between unlike beads 𝑎𝑖𝑗 can be expressed in terms of excess energy Δ𝑎 compared to the reference: Δ𝑎=𝑎𝑖𝑗−𝑎𝑖𝑖.

**Statistical ensembles – methodology to compute IFT**

In this chapter, three different statistical ensembles are used depending on phenomena under investigation. A workflow summarizing the proposed methodology is shown in Figure 17.

– The DPD model is combined with coarse-grained Monte Carlo (CG-MC) technique in order to simulate systems in the Gibbs (NVT) ensemble. This approach has already been used and validated by Wijmans et al. [98] with beads and soft potentials from DPD models. In the Gibbs ensemble, two separated simulation boxes that can exchange particles are used with a constant total volume 𝑉. Thus, it is possible to describe equilibrium between two phases without considering explicitly the interface. Gibbs (NVT) ensemble simulations are used in our work in order to compute phase diagram and, thus, to check the relevance of parameterization methods to reproduce composition in bulk phases. Three different types of Monte Carlo moves are used: (1) translation of beads, (2) transfer of beads between the two boxes and (3) concerted change of volume of each box. In addition, for hexane molecules which are represented by two beads, rigid body rotation and configurational regrowth moves are added. These movements are described in greater details in the section 3.3.3. During CG-MC simulations in the Gibbs (NVT) ensemble, chemical potential of each species are calculated using Widom insertion test method [109].

– Some CG-MC simulations are also performed in the osmotic (μsoluteNsolventPzzT) ensemble in order to describe a system with a constant number of solvent particles (Nsolvent) and a variable number of solute particles, fixing the chemical potential of the solute (μsolute). Osmotic ensemble has already been used by Rekvig et al. [84] to compute the number of surfactants necessary to reach an imposed IFT value between water and oil phases. CG-MC simulations in osmotic (μsoluteNsolventPzzT) ensemble with an explicit interface allow to predict the solute concentration at the interface from known bulk compositions. The imposed chemical potentials are obtained from previous simulations in the Gibbs (NVT) ensemble. Three different types of Monte Carlo moves are used: (1) translation of beads, (2) change of volume along z-axis (perpendicular to the interface) and (3) insertion or removal of solute beads. For hexane molecules, rigid body rotation and configurational regrowth moves are added. These movements are described in greater details in the section 3.3.3.

– Finally, DPD simulations are performed in the NVT ensemble in order to compute the interfacial tension.

In this chapter, CG-MC in the Gibbs (NVT) ensemble and DPD simulations are performed at constant density (𝜌̅=3). Therefore, the total pressure of the system varies depending on the composition. This choice was made to simplify the parameterization procedure. However, it is important to notice that an alternative parameterization of DPD interactions can be done by working at constant pressure. In this case, calculation of interaction parameters between like and unlike beads is dependent on the total density of the system (see equations (52), (56) and (57)). Consequently, additional bulk phase density data are required to obtain interaction parameters of pure components. Noting that, on this basis, any phase property should be calculated at constant pressure using a reference value.

CG-MC simulations were carried out with the molecular simulation package GIBBS [110]. For simulations in the Gibbs (NVT) ensemble, the two subsystems (“boxes”) have each an initial dimension of 𝐿𝑥=𝐿𝑦=𝐿𝑧=10 (in DPD units). The volume of each boxes can vary during the simulation but the total volume remains constant. The total number of beads is 6000. For simulations in the osmotic (μsoluteNsolventPzzT) ensemble, box dimensions were set to 𝐿𝑧=60,𝐿𝑥=𝐿𝑦=10 (in DPD units). Two planar water-organic compound interfaces are created normal to the z-axis. Box length in z-direction is six times larger than in the x and y-directions in order to avoid interactions between the two interfaces. Initial boxes containing a total of 18 000 DPD beads are built for different solute/solvent concentrations using the PACKMOL software package [111, 112]. DPD simulations in the NVT ensemble were performed using the molecular dynamics simulation package NEWTON [113]. Initial configurations are derived from simulations in the osmotic (μsoluteNsolventPzzT) ensemble. The area of the interface is kept constant (𝐿𝑥=𝐿𝑦=10, in DPD units). A modified version of the velocity-Verlet algorithm [63] governed the equation of motion, and the time step is fixed at δt = 0.001 in DPD units. Constants in the dissipative force γ and random force σ were set to 4.5 and 3, respectively, in order to keep the temperature fixed at 𝑘𝐵𝑇=1, thus satisfying the fluctuation-dissipation theorem (see section 3.2.4). In all simulations, periodic boundary conditions were imposed in all directions.

The IFT values are evaluated using two local methods (Irving and Kirkwood [69] and Harasima [70]) and one global method (Kirkwood-Buff [67]). More details on these methods are given in section 3.4.3.

**Table of contents :**

ACKNOWLEDGEMENT

TABLE OF CONTENTS

INTRODUCTION

**CHAPTER 1. SETTING UP STAGE **

1.1 INTRODUCTION

1.2 CRUDE OIL AND PETROLEUM FRACTIONS

1.2.1 FORMATION OF PETROLEUM

1.2.2 DESCRIPTION OF PETROLEUM

1.2.3 DESCRIPTION OF PETROLEUM FRACTIONS

1.3 MOLECULAR COMPOSITION OF CRUDE OIL

1.3.1 HYDROCARBONS

1.3.2 OTHER CRUDE OIL COMPONENTS

1.4 CRUDE OIL EXTRACTION

1.4.1 CRUDE OIL EXTRACTION METHODS

1.4.2 RECOVERY OF RESIDUAL OIL USING EOR TECHNIQUES

1.5 CONCLUSIONS

**CHAPTER 2. CHARACTERIZATION OF CRUDE OIL **

2.1 INTRODUCTION

2.2 PROPERTIES OF CRUDE OIL: EXPERIMENTAL METHODS

2.2.1 ELEMENTAL ANALYSIS

2.2.2 DENSITY AND SPECIFIC GRAVITY

2.2.3 MOLECULAR WEIGHT

2.2.4 BOILING TEMPERATURE AND DISTILLATION CURVES

2.2.5 CHEMICAL FAMILIES

2.2.6 INTERFACIAL TENSION

2.3 STATE OF THE ART FOR THE REPRESENTATION OF CRUDE OIL

2.3.1 FRACTIONATION APPROACHES

2.3.2 LUMPING METHODS

2.3.3 MOLECULAR RECONSTRUCTION

2.4 CONCLUSIONS

**CHAPTER 3. COARSE-GRAINED SIMULATION METHODS **

3.1 INTRODUCTION

3.2 DISSIPATIVE PARTICLES DYNAMICS (DPD)

3.2.1 CONSERVATIVE FORCE

3.2.2 INTRAMOLECULAR FORCE

3.2.3 DISSIPATIVE AND RANDOM FORCES

3.2.4 FLUCTUATION-DISSIPATION THEOREM

3.2.5 TIME INTEGRATION SCHEMES

3.3 COARSE-GRAINED-MONTE CARLO SIMULATION (CG-MC)

3.3.1 PRINCIPLE AND METROPOLIS ALGORITHM

3.3.2 STATISTICAL ENSEMBLES

3.3.3 MONTE CARLO MOVES

3.4 COARSE-GRAINED SIMULATIONS

3.4.1 COARSE-GRAINED MODEL

3.4.2 SIMULATIONS BOXES AND BOUNDARY CONDITIONS

3.4.3 INTERFACIAL TENSION CALCULATION

3.4.4 COARSE-GRAINED UNITS

3.5 CONCLUSIONS

**CHAPTER 4. PARAMETERIZATION OF LIQUID-LIQUID TERNARY MIXTURES **

4.1 INTRODUCTION

4.2 STATE OF THE ART OF PARAMETERIZATION OF COARSE-GRAINED SIMULATIONS

4.3 NEW PARAMETERIZATION APPROACH FOR LIQUID-LIQUID EQUILIBRIUM OF TERNARY SYSTEMS: METHODOLOGY AND COMPOSITIONAL DETAILS

4.3.1 PARAMETERIZATION OF INTERACTIONS

4.3.2 STATISTICAL ENSEMBLES – METHODOLOGY TO COMPUTE IFT

4.3.3 SYSTEMS STUDIED

4.4 RESULTS

4.4.1 COMPOSITION DEPENDENCE OF THE FLORY-HUGGINS INTERACTION PARAMETERS

4.4.2 LIQUID-LIQUID EQUILIBRIUM (LLE)

4.4.3 INTERFACE COMPOSITIONS

4.4.4 INTERFACIAL TENSION

4.5 THERMODYNAMIC MODELS FOR PREDICTION OF LIQUID-LIQUID EQUILIBRIA

4.5.1 PREDICTION OF THE COMPOSITION OF LIQUID-LIQUID EQUILIBRIA

4.5.2 PREDICTION OF INTERACTION PARAMETERS USING A THERMODYNAMIC MODEL

4.5.3 TRANSFERABILITY OF INTERACTION PARAMETERS FOR ALKANES

4.6 CONCLUSIONS

**CHAPTER 5. REPRESENTATION AND COARSE-GRAINED SIMULATIONS OF CRUDE OILS **

5.1 INTRODUCTION

5.2 THEORETICAL BACKGROUND FOR CRUDE OIL MOLECULAR REPRESENTATION

5.2.1 REPRESENTATION OF THE LIGHT FRACTION (C20-): THE LUMPED ALGORITHM

5.2.2 REPRESENTATION OF HEAVY FRACTIONS (C20+)

5.3 RESULTS

5.3.1 MOLECULAR MODEL OF CRUDE OIL

5.3.2 COARSE-GRAINED MODEL OF CRUDE OIL

5.3.3 PARAMETERIZATION OF INTERACTIONS FOR DPD SIMULATIONS

5.4 RESULTS OF DPD SIMULATIONS OF CRUDE OIL

5.5 CONCLUSIONS

GENERAL CONCLUSION

**BIBLIOGRAPHY**