Molecular simulation of crystal alkanes biology essay

A dissertation submitted toThe University of Manchesterfor the degree ofMaster of Engineeringin Chemical EngineeringbyAbhit GhelaniSupervisor: Professor Andrew Masters

School of Chemical Engineering and Analytical Science

The University of Manchester

May 2013

Acknowledgements

ABSTRACT

Contents

1Introduction

The diesel engine is becoming more popular in Europe due to their efficiency and high fuel economy. This is due to the rising gasoline fuel prices. This trend shows the importance for the development of diesel fuel efficiency. Figure – The trends of diesel and gasoline fuelThe behaviour of diesel in cold temperatures is problem faced both in the upstream and downstream. Diesel contains a number of different compounds such as naphthalene, alkyl benzene and sulphur containing compounds. When diesel is exposed to low temperatures it causes the fuel to gel and produce insoluble wax crystals (Henry S Ashbaugh 2004). It is generally accepted that it is the n-alkanes which are the primary component of the wax crystals (Hansen 1991). N-alkanes are the simplest hydrocarbons with the molecular formulate, where the n is the carbon number of the alkane. The solid wax crystals in diesel cause blockage of filters and narrow lines in the fuel system which cases operational failure of diesel and hence starvage of the engine (Filter Manufacturers Council 2012). Wax crystal form as the solubility of n-alkanes decrease with temperature, where diesel is the solvent. There are several different methods already developed to tackle wax crystals such as fuel heaters. Alternatively the most economical and commonly used method is to add cold flow improvers to diesel reduce the amount of n-alkanes crystals in the diesel below the cloud point.

Crystallisation of n-alkanes

The classic nucleation theory is prevailing theory in the homogeneous crystallisation of alkanes which is based on the work by Turnbull and Fisher (Turnbull D 1949). The classic nucleation theory states that the crystal nucleus which consists of a thermodynamically stable phase is separated from the liquid by a sharp infinitely think interface. Temperatures that are below melting point, the competition between the free energy gain of the interior of the nucleus and the free energy cost of the interface creates a free energy barrier. The time required to surpass the free energy results in the observed induction period before a nucleation event. Depending on the length and carbon number, diesel fuel usually crystallise into the following four crystal shapes which are triclinic, monoclinic, orthorhombic and hexagonal (X. G. Henry S Ashbaugh 2004). Triclinic crystals are formed when an alkane with an even carbon number between 6 and 26 crystallises, medium sized molecules with carbon number between 28 and 36 crystallise into the monoclinic structure and even carbon alkanes with a carbon number greater than 36 crystallise into orthorhombic sttructure and odd carbon number molecules all crystallise hexagonally (X. G. Henry S Ashbaugh 2004). It is observed that the waxes pack in single layers the crystallisation of alkanes has been studied thoroughly. It has been found that the composition of paraffin has a large influence on cold flow and pour point, rather than their concentration and or the the nature of the solvent used, impling that the carbon numbers from 21 to 27 have the greatest influence on cold flow problems (Khan U. H 1995).

Homogeneous crystal nucleation is hard to detect using experimental equipment as the nucleation takes place in the interior of the under-cooled liquid. Vonnegut used the droplet technique to address this problem and successfully implemented by Turnbull. The liquid is dispersed into a large number of droplets, only micron in size, thereby getting rid of the impurities in the liquid, and allowing homogenous nucleation. However the problem with this method was the due to the droplet size being so small, the crystal growth is much faster than the nucleation rate causing the droplet to immediately to crystallise. The crystallisation method can be monitor using X-ray scattering or differential scanning calorimetey allowing the estimation of the nucleation rate.

Molecular simulation of crystal alkanes

Cross polarized microscopy has led to an understanding of wax crystals at a crystal level view of the system and how the composition of a compound will affect the overall morphology of the crystal. But, it lacks the molecular level view of the molecular aggregation and the formation of single crystals by the mechanism known as co-crystallisation. Therefore molecular simulations enable us to get an insight at molecular level. Molecular simulation must take into account a number of parameters that affect crystals, such as solvent properties, crystal growth kinetics and transport properties. Conducting a molecular simulation on all the components found in diesel or crude oil would be computational intensive. Therefore, it would be wise to conduct on simpler organic systems which have already been conducted. The first computational work that was successful was work conducted by Liu and Bennema. CHINESE REFERENCE!!!!!. Liu and Bennema used stastiscal mechanics to investigate the interfacial interactions between crystal and the solute, the mechanics was also combined with kinetic growth models and the solvent-crystal physical properties. These properties were obtained from previous experimentation or molecular dynamics to predict the crystal shapes for systems. To reduce the intensity of the computation power required, in 1998 Winn and Doherty used the properites of a pure component and an attachment energy model to determine the interfacial interactions between the crystal and solute. The attachment energy model summed together the intermolecular forces such as van der waals. However this is only for one molecule as diesel in reality is a multicomponent mixture further study on other independent n-alkane lengths is required. CHINESE REFERENCE!!!!!. Nucleation is a non-equilibrium dynamic process therefore molecular dynamics simulation method is the most reliable method. The molecular dynamic simulation was used by Madel to study of nucleation using a Lennard jones system.

Cold flow improvers

There has been a great effort to develop cold flow improver. It has been found that the inhibition by cold low improver of the crystal formation in diesel can be categorised into two different mechanisms, nucleation and growth inhibition. The cold flow improvers are polymeric materials which are able to improve the low temperature properties of diesel. The cold flow properties of diesel is characterised by three parameters cloud point, pour point and cold filter plugging point. The cloud point is the temperature at which the crystals first appear, this is dependent on the alkane’s molecular weight and the chemical nature of the alkane crystals. The pour point is the temperature at which diesel will no longer flow, this is because the diesel gets trapped by a wax crystal matrix structure. Pour point depressants improve the flow properties of the diesel fuel. The cold filter plugging point corresponds to the plugging of a 45µm filter under standardized cooling conditions due to wax crystal growth. The presence of the wax crystals in the fuel is dealt with by using cold flow Improvers. There are three different types of cold flow improvers; random copolymers, Comb type polymer and Nucleators. Random copolymers are ethylene vinyl ester co-polymers acetate, the most commonly used is ethylene vinyl acetate as it has the most effective wax crystal growth inhibition properties; it affects both the cold filter plugging point and the pour point. The comb type polymers have long n-alkyl side groups along the backbone, these co-crystallise with n-alkanes. Finally Nucleators are additives with low molecular weight polymers, these additives work by reducing the number of large crystals are formed by increasing the number of small crystals present in the diesel.

Previous work on cold flow improvers

There has been a large amount of work carried out on wax additives. There are several different types of inhibitors that have been studied but there is still a lack of deep understanding of the inhibition properties of cold flow improvers. It is known that the cold flow improvers interact with the surface of wax crystals by co-crystallising with the wax crystals. This minimises the wax cohesive energy density, the force responsible for the ordering transformation. It has been found that inhibitors can alter the solubility of n-alkanes in n-alkane solvents. The mechanism has been further investigated by Jang whom managed to show that inhibitors alter solubility by associating with the chained paraffin (Jang YH 2007). The interaction between the inhibitor poly(octadecyl acrylate), PA-18, and the n-octacosane was investigated by using two different method. The first was to adding layers of the crystal to PA-18 and the second involved rapid crystal growth by adding alkane molecules one by one to the surface of an adsoberd PA-18. The united atom representation was used for the PA-18. The simulations of PA-18 were simulated on several different surfaces of n-octacsane crysal and it was that that the wax crystal adsorber along the inter-lamellar gap, occupying the sites at which the alkanes would have been place even there was no inhibition. It was also found that when the crystal wax was grown onto the inhibited surface that these caused the system shifted so the alkanes that were added adsorbed onto the inhibitor rather than forming a larger crystal. This shows that the PA-18 modified the wax surface and created derfects and at low concentration the defects prevent further growth of the crystal. Another cold flow improver that has been investigated is the inhibitor Maleic Anhydride based copolymers. The copolymers of the inhibitor were found to decrease the wax crystal size and prevented the formation of a crystal lattice. Another variation that was used was to add both Maleic anhydride-vinyl acetate and Steramide derivatives, which showed good improvement in the cold flow improvements in reducing the CFPP. This was found to be due to the carbonyl oxygen group found in the maleic anhydride.

Ethylene Vinyl Acetate

Characteristics of EVA

There are several different co-polymers of Ethylene Vinyl acetate and there is a small amount of preciece information regarding its molecular structure. The linear EVA structure is found to an effective cold flow improver.

Previous work on EVA

The ethylene vinyl acetate has also shown to display nucleation and growth inhibition mechanisms. During the nucleation mechanism, EVA copolymers are found to crystallise from the solution before the wax crystals form at the melting temperature. As with the PA-18 the phase separated EVA would attract the n-alkane and adsorb then on the surface of the EVA. As the growth inhibitor, the EVA inhibits crystal formation and the temperature at which wax crystal appear. The factors that decide whether the EVA will act as a nucleating agent or crystal inhibitor were found to be dependent on the vinyl acetate content and the nature of the diesel matrix. Therefore it shows that the EVA will interact differently depending on the composition of the diesel. Another study has investigated ethylene vinyl acetate in a mixture of cyclohexane/dichlorethane (Jin Wen Qian 1997). Ethylene vinyl acetate was added to the mixture containing 10 weight% of wax, as the concentration of the cyclohexane was increased the EVA no longer had an effect on the pour point. However it was found the larger the EVA copolymer the large decrease of the pour point. Therefore it would be beneficial to further study this using molecular dynamics to understand the interaction. Molecular dynamic simulations have found that the EVA thermodynamically prefers to adsorb parallel with the wax surface. The polar VA on the surface prevent paraffin deposition as the CH3 and CH2 repel the VA. The non-polar ethelyene backbone assembled itself closely to the crystal lattice, while the polar groups extended towards the surface of the growing plane. Suggesting that the EVA modifies wax crystallisation by hindering further growth along the plane.

Project aim

The project will look into addressing the concerns with the formation of plate like n-alkane wax crystals in diesel fuel during winter. The wax crystals cause the fuel to gel and block the fuel filers in the engine. Therefore cold flow improvers are added to the diesel fuel to alter the wax crystal shapes occurring. Diesel fuel is a mixture of alkanes, however in this project will focus on Undecane. The project is designed to get insights into: The structure of odd alkane crystals by oThe behaviour of cold flow improvers with short odd alkanesThe interactions of cold flow improvers and n-alkane crystalsInteractions of cold flow improvers with the n-alkane crystalsTo achieve these objectives it is planned to: Build moleculer moedels of undecane fuel using the untied atom modelModel the ethylene vinyl acetate and add it to the C-11 alkane crystals

Simulation method

Force field models

The force field method requires that’s several different assumptions, firstly it is assumed that the nuclei and electron are taken as a single atom particle and that these atoms are spherical and a single charge. The interactions for the molecules are based on springs and classical potentials.

Description of molecular mechanics for Force field model

The approximation that is used in molecular dynamic simulation is that the energy of the molecule can be described as a function of the force field and it is dependant only on the atomic positions. This method computers the molecular potential energy as the sum of the variation of energy as a function of variations in bond lengths, bond angles and torsional angles not at equilibrium as well as the intramolecular and intermolecular forces such as Van der Waals and electrostatic interactions. The feature that makes the force field is the fact they are transferable this means that a particular atom will have the same parameters in different molecules. Such as the stretching of aliphatic carbon-hydrogen bonds are independent of the molecular environment.

Bond stretching

The is the energy function for stretching the bond between two atom A and B. It uses Taylor expansion around the equilibrium bond length. The expression that describes this isThe derivatives are evaluate at and is normally set to zero. Since the Taylor expansion is going to be considered around the equilibrium value it will also be set to zero reducing the above equation toThe simples form is the harmonic form and sufficient to determine of most equilibrium geometrics. The refence bond length and is known as the equilibrium bond lengh, however this is when all other force fiel terms are set to zero. The actual equilibrium bond lengh is when the bond length for the minimum energy of a molecule. The force of the bond between atoms is much stronger than compared to the other forces that may be present, justifying the use of the harmonic approximation. However this approximation is limited to the fact that there is no great amount of deviations from the equilibrium bond length and if there is a large deviation the harmonic approximation can no longer be used. To calculate the molecular structure and vibrational frequencies more accurately, the Morse potential can be used. The morse potential can also be used if the bond length deviates from . The morse potential is given byWhere is he dissociation energy and http://rimg. geoscienceworld. org/content/42/1/1/F3. large. jpgFigure comparison of the simple harmonic with the Morse curve (Geoscience worl n. d.)It can be seen from Figure comparison of the simple harmonic with the Morse curve (Geoscience worl n. d.) the harmonic approximation is reasonable to the at the bottom of the potential curve which correspond to bonding in ground-sate molecules. However it is less accurate away from the equilibrium as the distance of the atoms get far enough where there are no interactions between the molecules. The Morse potential models the dissociation energy between the molecules.

Bond angle bending

The energy required for bening an angle between atoms where there is a bond between 1 and 2 and between 2 and 3 is. As with the taylor series expansion is expanded around the an equilbirum bond angle and terminated at the second order. Giving the harmonic approximation as: To increase the accuracy of the bond stretching, it is best to incoporae higher order terms. Higher order terms greater than 3 are usually used for pathological simulations which are usually highly strained molecules.

Torsions

The bond stretching and angle bending terms are known as hard degrees of freedom, as these terms require substantial energies to cause deformation from the reference values. In fact the variation in the structures is due to the interplay between the torsional and non-bonded contributions. The energy change, , occurs due to rotation around a 2-3 bond in an atom sequence where are bonded. There are two differences between and, the first is that the energy function must be periodic with torsion angle, , therefore if the bond is rotated the energy must return to the original value. Secondly the energy required to distort a molecule by rotation around a bond is low. cannot be expressed using taylor series expression, instead it must be expressed as a Fourier series so that the periodicity is still present. The constant determines the size of the barrier for the rotation around the bond. When the term n is 1 the rotation is periodic by, if n is 2 the term is periodic by. For example in ethane the most stable conformation is one in which the hydrogen atoms stagger relative to each other, so the eclipsed conformation is the energy maximum. This is due to the fact ethane’s structure contains three identical hydrogen, three energetically equivalent staggered, and three equal eclipsed conformations. Therefore the energy profile must have three maxima and three maxima. All organic molecule have a maximum number of four valences, there are a few exceptions which only have rotational profiles which only have three minima. The above EQUATION!!!!! Can qualitively produce this profile and when it is combined with non-bonded interactions can produce quantitative data. To calculate the torsional energy for a for a molecule which has the structure the equation!!!!! is commonly used: The equation contains both + and – so for one fold rotational term has a minima at the angle , the two rotational term has minima for angles and , and the three fold term minima for angles of , , and .

Non-bonded interactions

Force fields have both bonded interactions and non-bonded interactions. The non-boned interactions are between atoms in one molecule with those in other molecules. The two types of non-bonded interactions that are used in force fields are the Van der Waal interactions and electrostatic interactions.

Electrostatic interactions

Electrostatic interactions occur when there is an unequal distribution of charge in a molecule. The uneven distribution of charge can be modelled by placing point charges at each individual atomic site. The coulomb potential can model the interactions between these point chargesWhere, is the permittivity of the free space, and are atomic charges, and is the distance between atom and

Van Der Waal interactions

The Van der Waal interactions describe the repulsion and attraction between atoms which are not bonded. As the distance between molecules increases the Van der Waal forces decrease and eventually reaches zero and for small distances the force become more repulsive. This is due to the overlapping of the negatively charged electron clouds of the two atoms with each other, however there is a slight attraction at intermediate distance between the electron clouds due to the induced dipole-induced dipole interactions. The attraction is created by a transient dipole moment being created which induces a charge polarisation in the neighbouring molecule. The Van der waal force is positive at small distance, it has a minimum which is slightly negative at a distance corresponding to two atoms touching each other and tends to zero fast as the distance between the atoms gets larger. The general functional form which takes these into consideration isThe Lennard Jones potential is widely used as it usually obeys the general requirements (Jones 1924), where the replusive part is given by a dependance. This is given byThere are no are strong theoretical arguments that favour the repulsive a. The twelfth term is found to reason for rare gases but it is too steep for systems that contain hydrocarbons. However the 6-12 potential is widely used as the can be found by squaring the term. The can be calculated from the square without performing a computationally expensive square root calculation.

Types of force fields

Force fields can be classified into two different groups the classical force fields and second generation force fields. Classical force fields include the AMBER which is used for proteins and DNA (Wang ZX 2006), GROMOS which is a general purpose moleclular dynamics computer simulation package for the study of biomolecular systems (Karl-Heinz Ott 1996). Second-generation force fields include the MMFF which was developed at Merck for a broad range of chemicals.

Molecular Dynamics Simulation

The energy of a state is given byIf a 2 state system with energies and . The relative populations of molecules in state 1 over state 0 is given byIf two of more states are available then the Boltzmann distribution can be used to obtain the probability of each state. The probability of the state is given bywhere q is the partition function. The partition function can be found usingthe value is the lowest energy of the system and is the number of different energy states. This therefore causes a problem for larger systems; if more than one molecule is available there will be an exponential increase in the number of possibilities as the number of molecules increase. Calculating large number of energy states would be computational expensive and time consuming. The solution to this is to calculate the average energies without calculating the partition function. Molecular dynamics ignores the partition function by generating different states with the correct Boltzmann distribution, producing states with a probability. The mean energy is given by an average over number of states,

Molecular dynamics method

The molecular dynamics method applies the Newton’s equations of motion. Using the information of the force applied on each atom, it’s possible to calculate the acceleration of each atom. By integrating the equation of motion the trajectory that average values of the position, velocities and acceleration of each atom can be found. Molecular dynamics is a deterministic which means that the state of the system in the future can be predicted from the current state. The force can be expressed as the gradient of the force fieldCombing this equation with Newton’s second lawGivesWhere is the force applied on the particle, is the mass and is the acceleration and is the potential energy of the system. This yields the derivate potential energy as a function of time. These equations can then be solved for a system of N atoms to find the position vectors as a function of time.

Integration algorithms

The potential energy is a function of the atomic position of the atoms in the system. There is no analytical method to solve the equations motions, they are solve numerically. There are several numerical algorithms that have been developed for integrating the equations of motion by using finite difference algorithms. The velocity verlet algorithm will be used.

Statistical mechanics

Molecular dynamic simulation is used macroscopic properties of a system through microscopic simulation. Statistical mechanics are used to connect microscopic simulations with macroscopic properties by providing mathematical expressions that relate macroscopic properties to the distribution and the motion the atoms and molecules of the N-body system (Chandler 1987). The pressure, , number of particles, and the numbers are some parameters that define the thermodynamic state of a system. Any other thermodynamic properties may be found from equations of state and other thermodynamic equations. The multi-dimensional phase space coordinates the atomic postions, r and momenta mv define the microscopic state. A system containing N particles has dimension of 6N and a single point in point in the phase space, defines the state of the system. The ensemble are a set of points which statisfy the conditions for a particular thermodynamic state. In a molecular dynamics simulation a number of sequence of points are generated in phase space as a function of time. There are a number of different ensembles with different characteristics as described in the table below. EnsembleCharacteristicsMicrocanonical ensemble(NVE)Fixed number of atoms, NFixed volume, VFixed energy, ECanonical ensemble(NVT)Fixed number of atoms, NFixed volume, VFixed Temperature, TIsobaric-Isothermal ensemble(NPT)Fixed number of atoms, NFixed pressure, PFixed Temperature, TThe simplest molecular dynamic simulation that can be conducted is the NVE ensemble. Rhis ensemble is usually used for an isolate system. The NVE ensemble is beneficial as it allows for the calculation of average values for thermodynamic properties. However it is of greater interest that the simulation are carried out using the NVT and NPT ensembles. The NVT ensemble is used to calculate the thermodynamic values of a system under a fixed tempeterature. The NpT ensemble is used for systems that require constant pressure regime. It allows for the changing in simulation box size during the course of the simulation, either isotropically or anisotropically depending on the method and system being studied. Periodic boundary conditions allows a simulation to be performed using a small number of particles by making the particles experience forces as though the particles were in a bulk solution. The primary paricle forces are calculated from particles in the same box as well as the image box. The cutoff is chosen so that the particles in the primary box do not see their images in the surrounding boxes.

Force fields

The study being conducted is for n-linear alkanes and branched alkanes. These are both nonpolar and flexible in character. Force field parameters need to be transferable between different molecules for any given interaction site. This mean that the methyl group in n-undecane and undacene should have identical parameters and the force field should be transferable to different state points and different properties. There have been several representation simulations conducted using molecular dynamic models. The explicit atom (EA) representation is the most detailed representation, where each carbon and hydrogen atom on the chain is individually represented. These require sources for the bond length, bond angles, torsional angles and, non-bonded interactions. However an explicit atom simulation is computationally costly. Another force field which is simpler is the commonly used united atom representation, where CH2 and CH3 groups are represented as a single bead representation. Therefore the chains of alkanes can be regulated by a similar set of bond lengths, bond angles, torsional and non-bonded angles for interactions between CH2 and CH3 beads. Also comparisons made between UA and EA force fields found that there is only a slight conformation difference but the UA over predicts diffusion, due there being no interatomic packing between the hydrogen. The Trappe-UA force field was developed by research group of Siepmaan of the University of Minnesota Chemistry department. (Siepmaan 1998). The parameters have for n-alkanes have been optimised over many years. It was found that the TRAPP-UA is the force field that gives results that are the closest to the experiments.

TraPPE-UA

This section, Trappe-UA is discussed further in greater detail. Trappe-UA uses the carbon beads in the alkyl gorup. Both the bonded and non-bonded potentials are used in the sum of the total potential energy. The intermolecular are used only for interactions between different psudo-atoms in different molecules or in the same molecules if it is not accounted for by the intramolecular bonded potentionals. The intramolecular bonded potentials are the bond angle bending and dihedral potentialsThe non-bonded interactions between two atoms in a system scan be described by the pairwise-additive potential consisting of Lennar-Jones 12-6 potentialInteractions between unlike Lenard jones interactions can be computer using the Lorentz-Berthelot combing rules.