DETERMINATION OF ELASTIC CONSTANTS OF ANISOTROPIC HEAVY PETROLEUM PRODUCT USING MOLECULAR DYNAMICS SIMULATION

M . S t e t s e n k o PhD candidate Fleet Technical Operation Department Odesa National Maritime Academy Didrikhson Str., 8, Odesa, Ukraine, 65029 E-mail: ms.stetsenko@gmail.com Досліджені пружні властивості деяких граничних вуглеводнів та асоціатів смолисто-асфальтенових сполук шляхом моделювання молекулярної динаміки. Отримані числові результати модуля Юнга, модуля зсуву та коефіцієнта Пуассона для цих речовин у широкому діапазоні температур. Автором запропоновано оцінювати реологічну поведінку важких нафтопродуктів матрицею пружності трансверсально-ізотропної речовини. Показано, що числові значення матриці прямо пропорційні хімічному складу нафтопродукту Ключові слова: нормальний парафін, смола, асфальтен, молекулярна динаміка, силове поле, пружні константи


Introduction
Crude oil is a key natural resource in the world, because its processing products are used in almost all spheres of human activity.Currently the major application is the use of oil products as a fuel in internal combustion engines and the working fluid in a variety of devices and mechanisms.This determines the high interest of domestic and foreign researchers in improving of understanding of petroleum and petroleum products physical and chemical properties.Nevertheless, the fundamental properties of the oil as an organized complex substance with a large free energy reserved still remain poorly understood.This is explained by the complex composition of petroleum and its boiling fractions, which can be summarized as a mixture of alkanes, cycloalkanes, aromatics and heteroatomic compounds known as resin-solvated asphaltene substances.The ability of components of heavy oils to association leads to the formation of the complex structural units.Therefore, petroleum and its products are referred to disperse oil mediums.Each fundamental property of such mediums is a function of their chemical composition.Recently, along with the study of volumetric and calorimetric characteristics, characteristics of cohesion, temperature transitions and solubility scientists are paying special attention to properties, which appear in the force fields: mechanical, electrical, electromagnetic.Data continue to emerge about the permittivity, electrical conductivity, optical refraction index, shear viscosity and other properties of various petroleum products.
In many cases, the way of understanding properties of a complex oil substance extends through the numerous researches in properties of its separate components.

Literature review
According to modern concepts in regard to the structural properties of petroleum and petroleum products, it can be stated that the high viscous oil products (such as marine heavy fuel) belong to periodic dispersion mediums [1,2].Solid phase in those is formed by the normal and dehydrated paraffin molecules with a high number of carbon atoms and asphaltene micelles which are stabilized by adsorbed resins.
Let us look at the situation with normal paraffins.There are three low-temperature crystalline modifications of n-paraffins: orthorhombic, monoclinic and triclinic.Their diversity is defined by a number of carbon atoms (length of the molecule) and the parity of the number (the symmetry of molecular chains).The authors in [3] have presented generalized information about the dimensions and type of symmetry of the unit cells of saturated hydrocarbons in the wide range of carbon number.
Rheological properties of n-paraffins have been a subject for research as well.The authors in [4,5] have reported results of yield stress study for some n-paraffins.
At the same time, physical properties of resins and asphaltenes are studied less.However, we can note a few important works.For example, first observation of liquid crystals formed by the asphaltene aggregates in petroleum fraction reported in [6]; the authors claim that viscosity of heavy petroleum is a function of asphaltene concentration due to the micellization effect.
The number of research papers devoted to paraffins and asphaltenes, where experimental part is done by means of technical means of cognition, is much less for the last.Asphaltenes are mostly studied using computer simulation methods due to difficulties in preparing samples from petroleum fractions.Complexity of the structure is another reason for that.Some of the hypothetical models of asphaltenes and resins were reported in literature [7][8][9].

The aim and objectives of the research
The aim of this study is experimental confirmation of the assumption that the elastic properties of heavy petroleum can be characterized by the elasticity matrix of for a transversely isotropic medium.In order to achieve this aim, it is necessary to solve the following problems: 1) identify the major components of the heavy oils molecular composition; 2) obtain numerical values of the elastic constants of the heavy petroleum molecular structures according to their symmetry through a computer based molecular dynamics simulations; 3) define the elasticity matrix of a particular petroleum product according to it molecular composition.

Brief description of the molecular dynamics method
Accurate representation of the potential energy lies at the heart of all simulations of real materials.Accurate potentials are required for molecular simulations to predict the behavior and properties of materials accurately and even qualitative conclusions drawn from simulations employing inaccurate or unvalidated potentials are problematic.
Elastic constants characterize the stiffness of a material.The formal definition is provided by the linear relation that holds between the stress and strain tensors in the limit of infinitesimal deformation.In tensor notation, this is expressed as where the repeated indices imply summation.s ij are the elements of the symmetric stress tensor, e kl are the elements of the symmetric strain tensor.C ijkl are the elements of the fourth rank tensor of elastic constants.In three dimensions, this tensor has 3 4 =81 elements.Using Voigt notation, the tensor can be written as a 6x6 matrix, where c ij is now the derivative of s i with respect to e i .Since s i is itself a derivative with respect to e i , it follows that c ij is also symmetric, and must contain 7 6 / 2 × =21 distinct elements.The number of elastic constants is farther reduced if the crystal structure of the material has a certain type of symmetry [10].
At zero temperature, it is easy to estimate derivatives c ij by deforming the cell in one of the six directions and measuring the change in the stress tensor, as this option is normally available in any molecular dynamics simulation code.
Calculating elastic constants at finite temperature is more challenging, because it is necessary to carry out a simulation that performs time averages of differential properties.One way to do this is to measure the change in average stress tensor in normal volume and temperature (NVT) simula-tion, when the cell volume undergoes a finite deformation.In order to balance the systematic and statistical errors in this method, the magnitude of the deformation must be chosen judiciously, and care must be taken to fully equilibrate the deformed cell before sampling the stress tensor.Another approach is to sample the triclinic cell fluctuations that occur in a normal pressure and temperature (NPT) simulation.This method can also be slow to converge and requires careful post-processing [11].
All simulations in this research were performed by using the NPT method in the form of commercial molecular dynamics computer program.
In the most general sense, molecular dynamics code integrates Newton's equations of motion for collections of atoms, molecules, or macroscopic particles that interact via short-or long-range forces with a variety of initial and/or boundary conditions.
Newton's equation of motion for an arbitrary atom can be written as meaning that force F i applied to the i'th atom can be found by differentiation the potential energy E with respect to the corresponding coordinates of the atom.The form of (1) has to be modified, in order to include both pressure and temperature terms.Thus, molecular dynamics code (in our case) uses extended Nose-Hover method of adopting equation of motion (1) [12,13].The potential energy in ( 1) is usually calculated by using a force field method.

1. Mathematical model of the DREIDING force field
In general, force fields can be divided into three categories: (a) force fields parameterized based upon a broad training set of molecules such as small organic molecules, peptides, or amino acids including AMBER, COMPASS, OPLS-AA and CHARMM; (b) generic potentials such as DREIDING [14] and UNIVERSAL that are not parameterized to reproduce properties of any particular set of molecules; (c) specialized force fields carefully parameterized to reproduce properties of a specific compound [11].During the past ten years, the number of research works in organic molecules structure and behavior has significantly increased.In order to provide good quality organic molecule simulations several force fields have been developed.
Various forms of classical potentials (force fields) for polymers and organic materials can be found in the literature [11].
As the DREIDING force field is characterized as useful for predicting organic crystal structures and dynamics [7,14], it was used in all simulation of the research.
The potential energy for an arbitrary geometry of a molecule is expressed as a superposition of valence (or bonded) interactions (E val ) that depend on the specific connections (bonds) of the structure and nonbonded interactions (E nb ) that depend only on the distance between the atoms [11,14] val nb In DREIDING the valence interactions consist of bond stretch (E B , two-body), bond-angle bend (E A , three-body), Технологии органических и неорганических веществ dihedral angle torsion (E T , four-body), and inversion terms (E I , four-body) while the nonbonded interactions consist of van der Waals or dispersion (E vdw ), electrostatic (E Q ) and explicit hydrogen bonds (E hb ) terms The total bond energy in ( 3) is a sum over all bonds.Number of bonds is typically very close to the number of atoms.DREIDING describes the bond stretch interaction either as a simple harmonic oscillator [14]  -α - For two bonds IJ and JK sharing a common atom the three-body angle bend terms are all taken of the harmonic cosine form ( ) where θ is the angle between bonds IJ and JK.The equilibrium angles are assumed independent of I and K.The torsion interaction for two bonds IJ and KL connected via a common bond JK is taken of the form where φ is the dihedral angle (angle between the IJK and JKL planes), n JK is the periodicity (an integer), V JK is the barrier to rotation (always positive), and 0 JK ϕ is the equilibrium angle.
In DREIDING the torsional parameters are based on hybridization and are independent of the particular atoms involved [14].
The most complex term in a typical organic force field is the inversion term, which is added to ensure that a particular atom I, which is bonded to three other atoms J, K, L, remains planar or non-planar.
It is necessary to include an energy term describing how difficult it is to force all three bonds into the same plane (inversion) or how favorable it is to keep the bonds in the same plane.
Denoting the angle between the IL bond and the JIK plane as ψ, energy expression can have a form as Both van der Waals and electrostatic interactions in (4) are calculated over pairs of atoms, so they are usually done concurrently [11,14].

∑∑
where R 0 is the van der Waals bond length (Å), D 0 is the van der Waals well depth (kcal/mol), Q i and Q j , are the charges in electron units, R is the distance in Å, ε is the dielectric constant (usually ε=1), and 332,0637 converts E to kcal/mol.DREIDING treats hydrogens in a special way [14].Hydrogens are not given charges or van der Waals parameters.This is to eliminate errors possible when using the van der Waals parameters appropriate for non-hydrogenbonded systems.
As a result, the DREIDING force field has a special hydrogen bond term for D-H-A interactions, where D is the hydrogen bond donor, H is the hydrogen bonded to it covalently, and A is the hydrogen bond acceptor, noncovalently attached as shown in Fig. 1.The DREIDING hydrogen bond uses both a radial DA R and an angular DHA θ part:

2. Building experimental models of hydrocarbons and heteroatomic compounds
The simulation models of hydrocarbons were built with use of the software called Materials and Processes Simulation (MAPS).In some cases (for complex heteroatomic compounds) application called Accelrys was used.In distinct to the n-paraffins crystal symmetry of aromatics was not considered in modeling process.This was so as aromatic molecules exact orientation in asphaltene micelles is not well understood and could be rather arbitrary.The authors report disc-like or spherical-like micelle formation [15].Aromatic hydrocarbons used in modeling properties of asphaltene aggregates are shown in Fig. 2, a, b. 1,7-dimethylnaphthalene and methylcyclohexane were chosen for experiments, as their properties are well studied, and there are simulation experiments in asphaltenes mixtures available [9].This makes easier for us to compare and evaluate experimental results.Saturated hydrocarbons were simulated separately.There were three n-paraffins selected for simulation: C 19 H 40 , C 24 H 50 and C 36 H 74 .The reason is that those are responsible for composition of heavy petroleum fractions.Structural properties of those mentioned aromatics and n-paraffins depend on their composition and H/C ratio.As can be seen from Table 1, mass of carbon takes approximately 85 % in molecules of n-paraffin, and 92 % in aromatic molecules.

a b
Fig. 2. Three-dimensional models of some aromatics: a -1,7-dimethylnaphthalene; b -methylcyclohexane Crystal lattice of each n-paraffin was built in MAPS by multiplication of original unit cell, which was formed according to the data from literature [3] (Fig. 3).For n-paraffins C 24 H 50 and C 36 H 74 number of molecules was 216, and for C 19 H 40 : 512.These numbers should indicate an evidence of good experimental setup.The total number of atoms in simulation models (Table 2) was made at an average of two thousand.
Although the molecular structure of the resins and asphaltenes is not well known today, there are many hypothetical models that can be found in numerous research papers.For this research molecule of asphaltene C 84 H 98 N 2 S 2 O 3 was chosen (Fig. 4, a, b).Its chemical formula and molecular structure were taken from literature [16].
As can be seen in Fig. 5, a, b, structure of resin in petroleum is similar to asphaltene and may also contain nitrogen and sulfur atoms.Hipothetical model of the resin molecule C 26 H 40 NS 2 was taken from literature [8].Comparison between molar masses of the resin and asphaltene molecules, is given in Table 3.We can note that for this type of molecules the mass of carbon is about 70-80 % of the total mass.Since authors in [15] claim that asphaltenes can undergo the process of both flocculation and self-association depending on the fraction content, it was decided to create two experimental models.As it can be seen from Table 4, the first model contains molecules of resins and asphaltenes that represents the process of flocculation in the presence of n-paraffins.N-paraffins were not included in the model, as they were studied separately.The second model reproduces the process of micelle formation, or self-association of asphaltenes in an environment with high aromaticity.Model 1 was obtained by multiplication of the unit cell containing arbitrarily positioned one molecule of the asphaltene and five resin molecules.Model 2 was obtained in a similar manner from the unit cell, where one molecule of asphaltene was surrounded by 10 molecules of 1,7-dimethylnaphthalene and 20 molecules of methylcyclohexane.The final numbers of molecules in simulated volumes are given in Table 4.
Peculiarities of the atomic composition of experimental models are shown in Table 5.It may be noted that the first    Linear dimensions of the models obtained are shown in Table 6.Isotropic behavior of the models is expected as liner dimensions of the unit cell as well as its angles are the same for each model.

3. Simulation
Simulation begins with operational parameters assignment.These include the following: -cutoff radius (usually takes 10 Å); -initial pressure (usually 0 MPa); -the final pressure (in most experiments taken -101.325MPa corresponds to -1000 bar); -temperature at the beginning, during and at the end of the experiment (all three are normally the same; simulation temperature were the next: 298.15 о K, 333.15 о K and 373.15 о K); -duration of the experiment (60-120 ps); -time step (1 fs); -interval for writing to a file (usually 10 or 20 fs); -force direction (coordinates X, Y, Z).
The default boundary conditions in MAPS are the periodic boundary conditions.The idea of periodic boundary conditions (PBC) is to embed the simulation volume or simulation cell into an infinite, periodic array of replicas or images [17].
Time spent on the molecular dynamics simulation depends on how powerful computer is used.It was noted that a computer with a dual-core processor, clocked at 2.2GHz, and having 4GB of RAM has to run approximately 8.5 hours in order to simulate 100 ps of molecular model response.Thus, we can estimate the total time conducting simulation experiments.Excluding failed experiments, this time will be about 250 hours.
The example of three-dimensional view of the relaxed structure after running the molecular dynamics simulation and its original condition is shown in Fig. 6.Those simulations were performed for asphaltene aggregates with temperature set point of 100 о С.
The final stage of the experimental study is the post-processing.Data array produced by the molecular dynamics code of the MAPS program contains values of required variables in the matrix form.The number of values corresponds to a predetermined time interval.In order to calculate the elasticity, there following parameters are needed: pressure and the linear dimensions of the cells at each time point.Data array analysis is automated and performed by the firmware code, but operator is still required to set up the coordinate system and control the accuracy of the calculation.The program calculates the strain value at each point of time and returns the function of pressure versus strain of the simulated volume.Since this dependence in most cases is non-linear, a linear function is found by means of interpolation in the form y kx b.
= + Then the coefficient k is numerically equal to the value of Young's modulus.Shear modulus is calculated in the sa me way.Poisson's ratio is calculated using values of t h e longitudinal and transverse strain.
N-paraffins with triclinic, monoclinic and orthorhombic symmetry require simulation in three dimensions: x, y and z; however in this research simulation was conducted in the directions x and z only.The validity of such simplification was justified in the literature [3], where authors report that n-paraffins raise their symmetry up to hexagonal with increase of temperature.

Experimental results
Known values of Young's modulus, shear modulus and Poisson's ratio are enough to characterize the elasticity of a material completely.However, we will show later how elastic constants c ij can be derived, as these are used in some calculations (like waves propagation and interaction).
Table 7 contains values obtained of those elasticity parameters for n-paraffins, at three different temperatures.Among these values particular interest is drawn by the Poisson's ratio, which possesses the negative value almost in all cases.
Organic auxetic are known in modern science.Thus, the authors of [18,19] showed that in nematic liquid crystals with parallel orientation of molecular axes such effect is possible.However, there are no other experimental proofs of n-paraffins negative Poisson's ratio available.
For the case of asphaltene associates only Young's modulus and Poisson's ratio were obtained due to isotropic symmetry.Values of these parameters are listed in Table 8.
Change in elasticity depending on the number of carbon atom in n-paraffin molecules illustrates Fig. 7. Graph in Fig. 7 is linear, and it includes value for C 8 H 18 , specially obtained for the reference purpose.
The role of temperature in normal paraffins elasticity was found as expected.As it can be seen in Fig. 8, the rise of temperature leads to the phase transition, which is indicated by the non-liner behavior of crystal cell elasticity.The liner behavior begins, when temperature reaches approximately 100 о C. In case of asphaltene aggregates elastic behavior is more liner, as it is shown in Fig. 9. Now let us search the elasticity of heavy petroleum products in the form of elasticity matrix.Since we are dealing with we expect n-paraffins and other chain molecules to have hexagonal symmetry, elasticity matrix of transversely isotropic medium can be well used in our case.
Transversely isotropic medium is characterized by the five independent elastic constants.In this medium, stiffness is unique along one of the axis.Thus, hexagonal crystal symmetry can be an example.The elasticity matrix of transversely isotropic medium can be written as where 11 c , 12 c , 13 c , 33 c , 44 c and 66 c are the elastic constants.Elastic constants in Eq. ( 5) can be found using the following well known formulas: Finally, it can be shown how elastic properties of a particular petroleum product can be estimated.

Discussion of numerical results for a particular heavy petroleum product
The author is much interested in marine heavy fuels elasticity.The literature survey has shown that heavy fuel composition in terms of elastic behavior described in this paper can look like ( ) ( ) where summands are multiplied on their predictive concentration.Then the elasticity matrix for a heavy marine fuel oil, at the temperature of 100 o C, will have the following numerical values: As it can be seen from ( 7), the transverse elasticity of marine heavy oil is three times smaller as longitudinal one.Shear elasticity is rather small and do not exceed 1 GPa.Nevertheless, shear modulus will rise with decrease of temperature, meaning that transverse waves in such a medium may become perceptible.The negative values of tangential components in (7) would look quite unusual for a liquid, but they are much expected for a liquid-crystalline substance, which heavy petroleum products, obviously, are.The stability of the materials with negative elastic constants is proven to be possible in a condition of both full and partial constraint [20].

Conclusions
The numerical results obtained for the elastic parameters of some normal paraffins and resin-solvated asphaltene aggregates, using molecular dyna m ics simulation method, suggest that their elastic proper ties are closer in nature to the properties of solids, rather than amorphous entities.
Elastic behavior of saturated an d cyclic hydrocarbons and heteroatom compounds was bot h expected and unexpected.Particular attention is drawn by the negative Poisson's ratio, which gives reason t o rank these substances to auxetics.
For the first time, numerical values obtained of marine heavy fuel elasticity matrix at the temperature of 100 о C.
The author finds proposed in this paper deductive approach to assess the elasticity of complex chemical substances, according to which new scientific knowledge about a complex matter can be obtained by examining its individual components, to be quite successful.Thus, by analyzing the molecular composition of a complex compound and defining the main structural elements of it, properties of each of them can be estimated separately, in a first place, and then a general result can be derived for a substance as a whole.
=and K I is the force constant.

Fig. 1 .
Fig. 1.Hydrogen bonding interpretation Both the radial and angular parts are set to zero beyond certain cutoff values and switching functions are used to make the transition to hb E 0 = smooth.The values of hb R and hb D depend on the convention for assigning charges.

Fig. 6 .
Fig. 6.Volumetric image of simulation models before and after relaxation, at the temperature of 100 о С

Table 1
Composition properties of some hydrocarbons

Table 2
Lattice parameters of some n-paraffins model contains a significantly higher amount of inorganic atoms.The higher ratio of H/C in the second model is explained by a significant number of aromatic rings in the molecular ensemble.

Table 3
Atomic properties of the resin and asphaltene

Table 5
Atomic mass distribution in the simulation models

Table 4
Molecular composition of simulation models

Table 6
Parameters of the simulation models

Table 7
Young's modulus, shear modulus and Poisson's ratio of some n-paraffins

Table 8
Young's modulus and Poisson's ratio of some resin-solvated and aromatics associated asphaltenes