Copyright© 1998, American Institute of Aeronautics and Astronautics, Inc.
MOLECULAR DYNAMICS STUDIES OF THERMOPHYSICAL
PROPERTIES OF SUPERCRITICAL ETHYLENE
Obika C. Nwobi? Lyle N. Long! Michael M. Miccif
Department of Aerospace Engineering
The Pennsylvania State University
University Park, PA 16802
Abstract
This work involves the determination of transport coefficients and equation of state of supercritical fluids by
molecular dynamics (MD) simulations on parallel computers using the GreenKubo formulae and the virial
equation of state, respectively. The MD program uses
an effective LennardJones potential, linked cell lists for
efficient sorting of molecules, periodic boundary conditions, and a modified velocity Verlet algorithm for particle displacement. Previously, simulations had been carried out on pure argon, nitrogen and oxygen, and this
has now been extended to ethylene, C^H^, at various
supercritical conditions, with shear viscosity and thermal conductivity coefficients, and pressures computed
for most of the conditions. The results compare well
with experimental and National Institute of Standards
and Technology (NIST) values.
Nomenclature
a = acceleration vector
aza = acceleration vector of a given atom
d = molecular bond length
f = intermolecular force
g = constraint force
i = unit tensor
J p = stress tensor
J7 = heat tensor
ke = Boltzmann constant
m = atomic mass
n = number of atoms per molecule
N = number of molecules in the system
0 = time origin
P = pressure
r = position vector
r^ = position vector of a given atom
"Research Assistant, Student Member, AIAA
t Professor. Associate Fellow, AIAA
°© 1998 by Nwobi, Long and Micci. Published by the
American Institute of Aeronautics and Astronautics, Inc. with
permission.
r
a/3 = position vector between two atoms i, j
RQ = c.m. velocity of a given molecule
Ra/3 = c.m. position vector between two molecules
t = time delay
T — temperature
V = volume
v = velocity vector
v^ = velocity vector of a given atom
a,/3 = given molecules
St = timestep size
e = LennardJones energy parameter
7Q = constraint between atoms
A = thermal conductivity
H = shear viscosity
$ = intermolecular potential
p = density
Introduction
Current rocket motors, gas turbines, and many projected advanced combustor designs operate supercritically or near critical conditions.1 Computational Fluid
Dynamics (CFD) has been used to model droplet evaporation and combustion in supercritical and near critical
environments by several researchers including Yang et
al.2 and Delplanque and Sirignano.3 Quantitative information on the coefficients of mass, momentum and
energy transport, and the equation of state are very important in these studies. However, in many complex systems such as rarefied flows, supercritical environments,
shock regions, and highlyenergized plasmas, these coefficients are not known accurately and the experiments
to measure them are difficult and costly. As a result
the necessary information, is lacking in the use of the
CFD methods for simulating supercritical evaporation
and combustion. Moreover, these techniques make use of
limiting approximations such as spherical droplet shape
and constant density which are not necessarily valid under these supercritical conditions. Hence, alternative
methods for modeling the evaporation of supercritical
droplets have to be considered.
Copyright© 1998, American Institute of Aeronautics and Astronautics, Inc.
Molecular dynamics (MD), the technique used in the
present study, involves the solution of the equations of
motion of a system of molecules that interact with each
other through an intermolecular potential.4 Provided
computer program to compute the relevant properties as
the droplet vaporizes. Hence, this method is not utilized
in the current study.
MD has been used to compute the equation of state
that an accurate potential can be found for the system
of interest, MD can be used regardless of the phase and
thermodynamic conditions of the substances involved.
Although computationally intensive, this method requires no a priori assumptions regarding geometrical
symmetry, transport properties, or thermodynamic behavior. All calculations are performed from first principles based on intermolecular potentials, and thermodynamic and transport properties are results instead of
assumptions. MD has been used to support research in
such fields as chemistry, biology, physics, and material
science.5"10 It is also currently being used to model the
evaporation of a submicron droplet in supercritical environments. 11  12
MD simulations are computationally intensive since
the interaction between each particle and all its neighbors must be taken into account. The speed of these
simulations can be greatly increased by using parallel
computers. The natural parallelism in MD comes from
the fact that the force calculations and position/velocity
updates can be done simultaneously for all the molecules.
The method of parallelization currently being used is
the particle decomposition method.13 This technique is
implemented by distributing the particles evenly among
the different processors. Each processor calculates the
forces and displacements for the particles for which it is
responsible. Then the new positions and velocities from
each processor are broadcast to all other processors with
which the next time step begins. Using the particle decomposition method, almost perfect load balancing has
been achieved.11
of some fluids by several researchers including Barker
et al.,21 Singer et al.,22 and Vogelsang and Hoheisel,23
for only a few state points. No comprehensive study
has been done to calculate the pressures of ethylene and
other fluids over a wide range of states.
An attempt is made herein to compute the transport
properties and pressures of supercritical fluids which are
relevant for droplet evaporation studies using the GreenKubo formulae, and the virial equation of state, respectively. Previous studies had been done on argon, nitrogen and oxygen, and work has now been extended to
ethylene in this paper. 24  25 Successful implementation of
the code will impact the understanding and modeling of
the processes occurring in current and future aerospace
propulsion systems.
Molecular Dynamics
Molecular dynamics (MD) has become a widely used
direct simulation technique. The advantage of this technique for the solution of complex reacting flow problems
is that all attendant physical phenomena  shear viscosity, thermal conductivities, surface tension, evaporation,
and reaction rates  come from the force potential between the particles. MD does not require the use of
empirical submodels for each physical process or temperature dependent physical properties. It is an extremely
powerful and simple method of simulation. The accuracy
of a given MD model depends mainly on the chosen intermolecular potential, and the finite difference scheme.
Potentials
Transport coefficients can be computed either by
the use of GreenKubo formulae or Einstein relations
during equilibrium molecular dynamics (EMD) simulations, 4 ' 14 ' 15 or by conducting suitable nonequilibrium
molecular dynamics (NEMD) simulations.16'17 The
GreenKubo formulae have been used successfully to
compute the transport coefficients of atomic fluids, while
only limited success has been achieved for molecular
fluids which are essential for supercritical droplet studies.18"20 At various times this might be due to the inadequate modeling of the longtime tail behavior of the relevant autocorrelation functions, the use of small system
sizes because of limited computer resources available, or
due to the incorrect use of potential model parameters.
The NEMD emerged as a suitable simulation technique
for computing the transport coefficients of fluids due to
the above drawbacks. Normally to compute shear viscosities or thermal conductivities, separate computer experiments have to be carried out and the codes cannot
be incorporated easily into an MD droplet evaporation
Molecular dynamics generally adopts a classical approach in the modeling of matter. Typically, atoms or
molecules are represented as point masses interacting
through forces (potentials) that depend on the separation of the particles. The cornerstone of all molecular
dynamics simulations is the intermolecular potential. In
order to calculate the positions, velocities, and accelerations of each molecule in the system, the forces acting
on each of the molecules are required. These forces are
derived from the intermolecular potentials. All results
and properties (thermodynamic, structural and transport) determined with MD simulations, also follow directly from the potential.
An appropriate potential for MD simulations must
capture the important features of intermolecular interactions: a strong repulsive force at short distances due
to orbital overlap, a negative potential well causing regular structure in liquids at low temperatures, and an attractive tail at large distances due to dispersion forces.
Copyright© 1998, American Institute of Aeronautics and Astronautics, Inc.
The most commonly used potential which captures the
important features of intermolecular interactions is the
LennardJones (LJ) 126 potential.4
(1)
where <£ is the energy of interaction between the centers
of two atoms i and j which are a distance r^ apart. The
two fundamental constants in the LennardJones potential are e, which accounts for energy and a, which accounts for size. The values of these parameters depend
on the particular chemical specie that is being modeled.
a represents the zero energy separation distance and has
units of length, while 6 represents the depth of the potential well and has units of energy. The LJ potential
can be used for modeling solid, liquid or gaseous states.
The attraction part of the potential, 1/r6, is used to
bind the system for the solid and liquid states, while the
repulsive part. —1/r 12 prevents collapse.
For molecular fluids additional measures have to be
taken. Molecules can be modeled by one of two
means: treating each atom in the molecule as a separate
LennardJones site (sitesite), or by defining a central
site with corrections for angular dependency. The central site method starts with the NewtonEuler equations
of motion. Because singularities appear when written in
terms of Euler angles, the equations are written in terms
of quaternions which are then solved in matrix form.4
The sitesite method consists of the solution of the ordinary equations of motion for each atomic site subject
to a constraint which is the fixed bond length.4'26 The
sites are located at the positions of the nuclei of the constituent atoms. Diatomics such as oxygen and nitrogen
as well as nondiatomics such as ethylene (63 #4) can
be modeled as 2site molecules. Because of the fact that
each carbon atom is much heavier than a hydrogen atom,
each group of CH^ atoms in ethylene can be treated as
a single site (composed of just the carbon atom) with
the effects of the hydrogen atoms accounted for by the
use of an effective LennardJones potential. For modeling systems at relatively low temperatures, fixing the
bond length is an excellent approximation. At such low
temperatures, very little vibration occurs and the bonds
are essentially rigid.27
The molecular interaction generally consists of pairwise additive atomatom interactions. For two diatomic
molecules, there are four interactions between the four
sites and only interactions between atoms of different
molecules are considered.
The intermolecular force, f , is obtained from the spatial gradient of the potential.
^
.
.
,
f = _V$ = — i— J— k
dx
dy
dz
2
Figure 1 shows the LJ potential and LJ force as functions of the separation distance between two particles.
The quantities have been nondimensionalized using the
two parameters of the LJ potential a and e. It can
be seen that there are no significant long range forces
between two given particles with the potential at a radial distance of 2.5cr having a value of only 1.6% of the
potential well depth value.26 Hence, for most LJ potential systems the interactions are truncated at a relatively short range in order to reduce the computational
demand.
1
2
Distance Bewteen Particles (r/c)
Figure 1: LennardJones 126 potential and force as
functions of particle separation distance.
Finite Difference Scheme
As mentioned earlier, the force on one particle due to
another is the gradient of the potential function. Since
pairwise additivity is assumed, the resultant force, f£
on an atom, i, of a given molecule, a, is simply the
sum of all the pairwise forces due to the surrounding
n x N — 1 atoms. Using Newton's Second Law applied
to each atom, the classical equations of motion of a system of n x N atoms can be expressed in the form
"»X = %
(3)
where a^ denotes the acceleration of a given atom. A total of 3(n x N) secondorder differential equations in the
radial position of the atoms, rla are generated, where N
is the total number of molecules in the system, and n is
the number of atoms per molecule. Given the total force
on each particle, the acceleration can be determined and
Copyright© 1998, American Institute of Aeronautics and Astronautics, Inc.
the velocity obtained. Using these velocities and accelerations, all the particles in the system are moved a finite
amount at each time step. An appropriate finite difference scheme must be fast, require little memory, stable,
accurate, and permit long time step size. The scheme
should also satisfy the known conservation laws for energy and momentum, and must reproduce certain timeand spacedependent correlations to a sufficient degree of
accuracy. Several schemes have been used in MD calculations including Gears predictorcorrector and the Verlet
algorithms.4
The Verlet schemes are the most widely used algorithms.4 There are several versions of the schemes including the StoemerVerlet, the leapfrog Verlet and the
velocityVerlet. The velocityVerlet algorithm is utilized in the present study and it requires the storage of
only atomic positions, velocities and accelerations. The
scheme is a dualstep process. In the first step the positions are advanced a full step, while the velocities are
advanced only a half step. After the accelerations are
updated at the current time step, the velocities are then
advanced another half a step.
= rla(t) + 6t <(*) + l5t2 aj,(0
5t
vj,(t + 6t) =
(4)
(5)
6t)
(6)
where aj,, vla and rla are the acceleration, velocity,
and radial position, respectively, of atom i of a given
molecule a. £ is the simulation time and 6t is the time
step size. The velocity Verlet algorithm stores positions
and velocities at a given time (i), and the accelerations
at time (t + 6 t ) . This scheme minimizes roundoff error
and the atomic velocities can be easily scaled in order to
control the temperature of the system.
The sitesite method with constraint dynamics is used
in this study for modeling rigid diatomic molecular systems. Each atom of a given molecule is treated as a separate LennardJones site. The technique involves writing the equations of motion to include the constraint
force due to the rigid bond length. The constraints are
the fixed bond lengths which are given by the following
equation
j
1
—0
vfa —
— r
r a — rr a r — & =
(7)
dij denotes the desired separation distance between
atoms i and j of a given molecule a. r^ and r£ are the
radial positions of the two atoms of a given molecules a.
The equations of motion are now written as
(8)
where fa, ga, ma and ra denote the intermolecular force,
constraint force, mass and radial position of atom i of
a given molecule a. The constraint forces help to keep
the desired bond lengths constant and require a few subiterations at each time step.
Performance Enhancements
Because of the computationally intense nature of MD
simulations, the system sizes tend to be relatively small
compared to those of macroscopic systems. A major
obstacle to these simulations is the large fraction of
molecules which lie on the surface boundaries. Whether
or not the simulation box is surrounded by a containing
wall, molecules on the surface experience quite different
forces from molecules in the core of the system. The
surface effects could lead to error in the properties and
quantities computed during simulations and should be
minimized. This can be achieved by implementing periodic boundary conditions.4
One great limitation of MD is the computational demand. Generally, each molecule experiences a pairwise
contribution from every additional molecule and its infinite images in the system due to the use of periodic
boundaries. This equates to an infinite number of computations which quickly exceeds the capability of the
largest supercomputers. Hence, many techniques have
been used to increase the speed and efficiency of the
code, and these vary depending on both the system size
and type.
Firstly, a "minimum image" convention is used to reduce the number of force calculations, due to the incorporation of periodic boundary conditions. This is
achieved by the use of a potential cutoff radius. Since
intermolecular forces decay rapidly beyond just a few
atomic diameters, the effect of such forces may be neglected for all molecules which lie outside the interaction
sphere of a molecule defined by the cutoff radius. The
number of force calculations is reduced from infinity to
 AT xn(N x n 1).
Secondly, a cell index method is used for checking
which molecules need to be included in the force calculation of each molecule. This is achieved by dividing the
simulation box into small cells into which the molecules
are sorted. Since interactions are only possible between
particles that are either in the same cell or in immediately adjacent cell, only molecules in the surrounding
cells need to be considered for force calculations.
Thirdly, a "linked list" method is used for fast sorting
of the molecules into cells.4 In the cell index method,
linked lists are used to associate molecules with the cells
in which they reside at any given instant; a separate list
is required for each cell.26 The "linked list" refers to
the linking of a list array with the head array during
Copyright© 1998, American Institute of Aeronautics and Astronautics, Inc.
the neighbor search.4 The routine is rapid and may be
performed at every timestep. This method requires only
the additional storage of two relatively small array pairs,
and its efficient use of memory is applied to track the
neighbor locations. The use of above schemes with a
potential cutoff reduces the force calculations to roughly
250(N x n).
Property Calculations
All calculations in MD are performed from first principles based on the intermolecular potential, and various thermodynamic, structural and transport properties
can be obtained from the simulations. The thermodynamic properties include pressure, temperature, internal
energy; structural functions include the radial distribution function; while the transport properties include various coefficients of diffusion, shear viscosity, and thermal
conductivity.
Pressure Calculation
Pressure is related to the force per unit area acting
normal to an imaginary surface in a system. By Newton's second law, the force can be related to the momentum that crosses the surface per unit time. The pressure
of the system may be calculated during an MD simulation by means of the virial equation of state: 4 ' 28 For
molecular fluids pressure is given by the following equation
pkBT
1
„. ...0
a=1/3=1 i=1 j=i
This mav be rewritten as
(10)
PV = NkBT a=l/3=l i=l j=l
where T. P, V and p denote the temperature, pressure,
volume and density of the system. N, n are the number
of molecules in the system and the number of atoms
per molecule, respectively, while kg is the Boltzmann
constant. Ra/j denotes the difference vector between the
center of mass position vectors of two given molecules a
and Q. and $ is the centertocenter potential function.
The first term on the righthandside of the above
equation denotes the idealgas contributions. It arises
because of the motion of the particles, while the second term accounts for intermolecular forces, assuming a
pairwise additive potential. For molecular systems, all
the intramolecular contributions including those along
the bonds have to be included. At low densities the second term is negligible, but as the density increases, its
contributions becomes significant.
Correlation Functions
In computer simulations, transport coefficients can be
calculated from equilibrium correlation functions, by observing Einstein relations, or by conducting a suitable
nonequilibrium MD simulation.4 In the present study
the transport coefficients were computed using GreenKubo formulae in which a correlation function is integrated over time.4'14'15 The autocorrelation functions
(ACFs) give insight into the microscopic behavior of fluids and they help in the understanding of the dynamics
of liquids. An ACF can be derived from the Enskog
kinetic equation and the corresponding transport coefficient then obtained via the linear response theory using the hydrodynamic equations.15 The ACF of a real
dynamic variable A measures the relation between the
values of the variable at two different times: A(0) and
A(t). It is given by the following expression:
C(t) = (A(Q)A(t))
(11)
where 0 and t denote an arbitrary time (time origin)
and the time difference (delay time), respectively, and
(...) denotes a time average.
GreenKubo Formulae
The shear viscosity coefficient, /j,, is a measure of the
internal fluid friction which tends to oppose any dynamic
change in the fluid motion. It gives an indication of the
shear stress induced by an applied velocity. The corresponding correlation function contains the nondiagonal
elements of the stress tensor and it is a collective property of the entire fluid.15 The GreenKubo formula is
given by
1
3VkBT
(J p (0)J p (i)> dt
(12)
where Jp is the stress tensor which depends on the velocities, radial positions and potentials of the particles of the
system. T and V denote the temperature and volume of
the system, respectively, and kg denotes Boltzmann constant. Statistical precision is improved, when appropriate, by averaging over all the three terms (fJxyiVxnfJyz)
that result from the stress tensor. In MD simulations,
the shear viscosity is computed by averaging over a finite
number of time origins.
M
(13)
Copyright© 1998, American Institute of Aeronautics and Astronautics, Inc.
where M is the number of time origins. For molecular
fluids the stress tensor is given by
Jp = m
0=1/3=1 i=l j=l
0=1
(14)
where RQ denotes the center of mass velocity of
each molecule, while RQ, denotes the difference vector between the center of mass position vectors of the
molecules. $ indicates the pair interaction potential
between the atoms of two molecules a and /3, while n
and AT are the number of atoms per molecule and the
total number of molecules in the system, respectively.
The shear stress ACF is composed of a kinetic term,
which measures the correlation of momentum transport
caused by atomic motions; a potential term which measures the correlation of momentum transport caused by
interatomic forces, and a cross term, which measures the
coupling of atomic motions and forces. Generally, the kinetic term dominates the ACF at low densities (gases),
while the potential term dominates at higher densities
(liquids). The potential and cross terms are generally
computed within the main force subroutine in the MD
computer code.
The thermal conductivity coefficient, A, measures the
transport of heat in a system. The correlation function is
obtained from the heat current and it is also a collective
property of the entire fluid.15 The GreenKubo formula
is given by
A=
(15)
WkBT2
where J, is the heat current. Statistical precision is improved, when appropriate, by averaging over the three
terms (A x , Ay, Xz) that result from the heat current. In
MD simulations, the thermal conductivity is computed
by averaging over a finite number of time origins.
M
(16)
where M is the number of time origins.
For molecular fluids the heat current is given by
N
a=l
a=l t=l
0=1 j=l
where i denotes a unit tensor.
The ACF is also composed of kinetic, potential, and cross
(kineticpotential) terms. The computation of thermal
conductivity takes more time than the computation of
shear viscosity because of the additional terms in the formula. As in the shear viscosity case, the potential and
cross terms are computed within the main force subroutine in the MD computer code.
Molecular Dynamics Simulations
The MD computer code is written in Fortran77 and
it uses the message passing interface (MPI). All the runs
were made on an IBM SP2 computer using 8, 16 or 32
processors. The simulations were carried out in the conventional constant number of moleculesenergyvolume
(NEV) ensemble. The equations of motion were solved
using the modified velocity Verlet algorithm with a time
step of 2xlO~ 15 s, periodic boundary conditions, and a
potential cutoff radius of 2.5cr. The simulations were
done using cubical geometries, with the length of each
side about 12 to 45fr depending on the density and the
number of molecules in the system. These lengths are
much larger than the critical point correlation length.20
The system sizes of the simulations carried out ranged
between 1000 to 1200 molecules. The LennardJones 126 (LJ) potential were used for the pair interactions and
the LJ parameters of ethylene used in the simulations,
as well as the critical temperature and pressure of ethylene are given in Table 1.
Table 1: Relevant Parameters of Ethylene
PC
MPa
5.04
Tc
K
A
282.0
3.33
d
(7
K
137.7
A
2.46
Starting from an initial configuration with the
molecules placed on the positions of a facecentered cubic
(FCC) lattice, the first 5 x 104 time steps were used for
equilibrating the system to the chosen thermodynamic
state. During the equilibration run, the velocities of the
molecules are frequently rescaled in order to obtain the
desired temperature. Normally, at the end of this period, the final temperature of the system differs by less
than 1 % of the desired value, while the total energy of
the system would have stabilized with a fluctuation of
less than 0.2 % from the mean value. The data runs
were about 2 x 106 time steps for the transport coefficients and 2 x 105 for pressures. Statistical uncertainties
in these calculations were estimated by evaluating independent averages over blocks of 1 x 105, and 1 x 104
time steps for the transport coefficients and pressures,
respectively.9
The number of molecules used in the simulation, the
time step size and the length of the runs have been
Copyright© 1998, American Institute of Aeronautics and Astronautics, Inc.
250
5000
10000
15000
Number of Molecules
Figure 2: Total computer run times for 1000 time steps
for different system sizes for oxygen at 15 MPa and
250 K.
chosen to allow for accurate evaluation of the GreenKubo time correlation functions and the computation
of pressures. Each simulation typically takes about 36
hours and 12 hours, respectively, for the transport coefficients and pressure determinations. Figure 2 shows
the total time required to carry out 1000 time steps using the transport coefficients and pressure computation
codes for various system sizes for pure oxygen at 15 MPa
and 250 K. The total time includes the computation,
communication, linked list, waiting and miscellaneous
times. The codes can be seen to scale with the number
of molecules in the system.
The calculated transport coefficient results are compared to values obtained from the National Institute of
Standards and Technology (NIST) software which uses
equations to determine the transport coefficients and
pressures of substances.29 The parameters of these equations for the transport coefficients were determined by
experiments.30'31 The average uncertainties for these
values are 2 and 4 per cent for shear viscosity and thermal conductivity, respectively. The NIST results are
only available for temperatures up to 400 K. The pressure results are compared to experimental values determined by Sychev et al.32
Results and Discussion
The main objective of this research is to develop tools
for predicting the transport properties and pressures of
pure fluids and fluid mixtures which are relevant for
modeling supercritical droplet evaporation. Simulations
have been performed using pure ethylene at various
subcritical, supercritical and near critical conditions.24
Pressures, shear viscosity, and thermal conductivities
have been computed.
Systematic error in an MD simulation is minimized
by extrapolating results to the thermodynamic limit by
using a large number of particles N and large system volume V while keeping the density N / V constant. The
effects of the number of particles in the system on the
computed transport properties were investigated for two
different states and are presented in Figures 3 and 4.
The plots show results for system sizes of about 300 to
3200, and they demonstrate clearly the independence of
the shear viscosity and thermal conductivity values on
the number of particles for the subcritical and near critical cases studied. This indicates that these system sizes
are sufficient for accurate computation of the transport
coefficients. It is believed that these results would be
the same at supercritical conditions. Most of the simulations in this study were done with about 1000 to 1200
molecules. The effect of the number of molecules in the
system N on the computed pressures is also presented
in Figure 5. The plots show results for system sizes of
about 500 to 3200 molecules. Like in the case of transport properties, the calculated pressures are independent
of the number of molecules for the range studied. Most
of the pressure results in this study were obtained using
system sizes of 10001200 molecules.
Figures 6, 7, 8 and 9 show plots of the shear viscosity,
as functions of temperature and pressure for ethylene.
Figures 10, 11, 12 and 13 show plots of the thermal conductivity, as functions of temperature and pressure for
ethylene. The MD results compare well with the NIST
values for most of the cases. However the NIST software
does not compute values for temperatures higher than
400 K, hence no comparison can be carried out with MD
values at these high temperatures. The MD thermal conductivity results do not match the NIST values very well
for the 15 MPa case.
Figure 14 shows plots of calculated pressures against
temperature for various densities for ethylene. Simulations have been done under various conditions between densities of 50300 kg/m3 and temperatures of
100600 K. For most of the cases, there is excellent agreement between the MD and experimental results.32
Copyright© 1998, American Institute of Aeronautics and Astronautics, Inc.
60
50
03
CO
45
3MPa, 160K(NIST)
40 
3 MPa, 160K(MD)
SMPa, 160 K (MIST)
5MPa, 160K(MD)
CO
Q.
30
£30
CO
CO
03
o
% 20
200 kg/m3, 600 K (exp)
200 kg/m3, 600 K (MD)
15
10
1000
2000
4000
3000
Number of Molecules
Figure 3: Effect of the number of molecules on shear
viscosity for oxygen.
1000
2000
3000
Number of Molecules
4000
Figure 5: Effect of the number of molecules on calculated
pressure for ethylene.
Conclusions
ou
I
:
,
20
n

I
,
,
i

1
,
.
,

1
1
I
,
————— 3 MPa, 160K(NIST)
•
3 MPa, 160K(MD)
60
I4"
I
....._.._.... — . 5 MPa, 160K(NIST)
•
5 MPa, 160K(MD)
x" I T
:
1
.
I
1
1
,
1000
,
r
*
,
I
1
,
2000
.
!
!
1
~
1
3000
!
,
The transport properties and pressures of ethylene
are being computed via molecular dynamics (MD) simulations using the GreenKubo formulae and the virial
equation of state, respectively. Simulations have been
carried out for a large number of subcritical, near critical and supercritical conditions. Results show that the
numbers of molecules in the simulation systems are sufficient for obtaining accurate results for the subcritical
and nearcritical conditions studied and the same result
should hold for supercritical conditions. Computed shear
viscosity and thermal conductivity coefficients compare
well with NIST results, while the pressure results compare well with experimental values. Future studies will
include the computation of the transport properties of
other hydrocarbons typically found in propulsion systems, and the properties at the surface of evaporating
droplets in supercritical environments.33 A successful
implementation of the code should have a great impact
on the understanding of the processes occurring in current and future aerospace propulsion systems.
,
4000
Acknowledgments
Number of Molecules
Figure 4: Effect of the number of molecules on thermal
conductivity for oxygen.
This work was supported by AFOSR Grant F496209410133 and NASA Grants NAGW1356 Supplement
10, and NGT10034. The authors would like to thank
Dr. Jeff Little and Dr. Teresa Kaltz for their contributions.
Copyright© 1998, American Institute of Aeronautics and Astronautics, Inc.
80
80

•
3 MPa (NIST)
3MPa(MD)
60
60 
03
03
Q_
Q.
<D
<D
O
0
20
I I
20
300
400
600
500
300
400
500
600
Temperature (K)
Temperature (K)
Figure 6: Shear viscosity values from MD and NIST for
ethylene at a pressure of 3 MPa.
Figure 8: Shear viscosity values from MD and NIST for
ethylene at a pressure of 10 MPa.
80rr
80
5 MPa (NIST)
5MPa(MD)
60
60
co
CO
03
03
Q_
CO
Q.
40
40
o
1
20
300
400
500
20
600
ou 300
400
500
600
Temperature (K)
Temperature (K)
Figure 7: Shear viscosity values from MD and NIST for
ethylene at a pressure of 5 MPa.
Figure 9: Shear viscosity values from MD and NIST for
ethylene at a pressure of 15 MPa.
Copyright© 1998, American Institute of Aeronautics and Astronautics, Inc.
i i i i i i i i i i i i i . i i
IUU
80
1
j=60

————— 3 MPa (NIST)
•
3MPa(MD)
IUU
1
1 

J
1
'
1
i

1
y
80 
I
I
I
c

I
1
1
1

————— 10 MPa (NIST)
•*
10MPa(MD)
r
J=60
:\
§40 
T
"
"
20 ~~
0
300
400
500

§,40
~
20
0 i
60C
T T T f
i i
300
i
i
.
i i
400
i
.
,
.
, i
60C
500
Temperature (K)
Temperature (K)
Figure 10: Thermal conductivity values from MD and
NIST for ethylene at a pressure of 3 MPa.
Figure 12: Thermal conductivity values from MD and
NIST for ethylene at a pressure of 10 MPa.
100
IUU
 5 MPa (NIST)
A
80
5MPa(MD)
£
7
I
60
20
500
'
' i
'
i '
'
i
'
'
'
'
i
 <
•
0
600
V
r
.
§,40
20
400
;
'
80
^60
300
.'
T
o °

<L
o o
:
————— 15 MPa (NIST)
•
15MPa(MD)
l
300
i
,
i
i
l
400
,
i
,
i
i
500
i
i
i
l
60C
Temperature (K)
Temperature (K)
Figure 11: Thermal conductivity values from MD and
NIST for ethylene at a pressure of 5 MPa.
Figure 13: Thermal conductivity values from MD and
NIST for ethylene at a pressure of 15 MPa.
10
Copyright© 1998, American Institute of Aeronautics and Astronautics, Inc.
[5] Garrison, B. J. Molecular dynamics simulation of
80
surface chemical reaction. Chem. Soc. Review, 21,
pp. 155162, 1992.
MD (symbols) vs Experiment (lines)
[6] Thompson, S. M., Gubbins,
60
K. E.,
Walton,
J. P. R. B., and Rowlinson, J. S. A molecular dy
namics study of liquid drops. Journal of Chemical
Physics, 81, pp. 530542, 1984.
300 kg/m3..
[7] Koplick. J., and Banavar, J. R. Continuum deduc
CO
CO
0)
tions from molecular hydrodynamics. Annual Re
view of Fluid Mechanics, 27, pp. 257292, 1995.
20
300
400
500
[8] Ohara, T., and Aihara, T. Molecular dynamics study on structure of nearcritical water. The
JapanU.S. Seminar on Molecular and Microscale
Transport Phenomena (J14), 1993.
600
Temperature (K)
[9] Stassen, H., and Steele, W. A. Simulation studies of
shear viscosity timecorrelationfunctions. Journal
of Chemical Physics, 102(2), pp. 932938, 1995.
Figure 14: Pressure calculations from MD and experiment for ethylene.
[10] Rapaport, D. C.
Microscale hydrodynamics:
Discreteparticle simulation of evolving flow pat
References
tern. Physical Review A, 36, pp. 32883299, 1987.
[1] Long, L. N., Micci, M. M., Kaltz, T. L., Little, J. K.,
[11] Little, J. K. Simulation of Droplet Evaporation in
and Wong, B. C. Submicron droplet modeling using
Supercritical Environments Using Parallel Molecu
molecular dynamics. AIAA Paper 950412, January
lar Dynamics. Ph.D. thesis, The Pennsylvania State
1995.
University, May 1996.
[2] Yang, V., Lin, N., and Shuen, J. S. Vaporization
[12] Kaltz, T. L. Simulation of Transcritical Oxygen Va
of liquid oxygen droplets in supercritical hydrogen
porization using Molecular Dynamics. Ph.D. thesis,
environments. Combustion Science and Technology,
The Pennsylvania State University, August 1998.
97(4), pp. 247270, 1994.
[13] Plimpton, S. Fast parallel algorithms for short
range molecular dynamics. /. Comp. Phys., 117,
[3] Delplanque, J. P., and Sirignano, W. A. Numeri
pp. 119, 1995.
cal study of the transient vaporization of an oxygen
droplet at sub and supercritical conditions. Inter
[14] Hoheisel, C. Theoretical Treatment of Liquids and
national Journal of Heat and Mass Transfer, 36(2),
Liquid Mixtures. Elsevier, Amsterdam, 1st edition,
pp. 303, 1993.
1993.
[4] Alien, M. P., and Tildesley, D. J. Computer Sim
[15] Hoheisel, C. Transport properties of molecular liq
ulation of Liquids. Clarendon Press, Oxford, 1st
uids. Physics Reports, 245(3), pp. 111157, 1994.
edition, 1987.
11
Copyright© 1998, American Institute of Aeronautics and Astronautics, Inc.
[16] Evans, D. J., and Morriss, G. P. Nonnewtonian
critical fluids. Journal of Thermophysics and Heat
molecular dynamics. Computer Physics Reports, 1,
Transfer, 12(2), pp. 16, 1998.
pp. 297343, 1984.
[26] Rapaport. D. C. The Art of Molecular Dynamics
[17] Millat, J., Dymond, J. H., and Nieto de Castro,
Simulation. Cambridge University Press, 1st edi
C. A. Transport Properties of Fluids: Their Correlation, Prediction and Estimation. Cambridge University Press, 1st edition, 1996.
tion, 1995.
[27] Vincenti, W. G., and Kruger, C. H. Introduction to
Physical Gas Dynamics. Krieger Publishing Com,
[18] Cheung, P. S. Y., and Powles, J. G. The properties
pany, Malabar, Florida, 1986.
of liquid nitrogen, part iv. a computer simulation.
[28] Haile, J. M. Molecular Dynamics Simulation: Ele
Molecular Physics, 30(3), pp. 921949, 1975.
mentary Methods. John Wiley, New York, 1st edi
tion, 1992.
[19] Cheung, P. S. Y., and Powles, J. G. The properties
of liquid nitrogen, part v. a computer simulation
with quadrupole interaction.
[29] National Institute of Standards and Technology,
Molecular Physics,
Gaithersburg, MD. NIST Thermophysical Proper
32(5), pp. 13831405, 1976.
ties of Pure Fluids Database: Version 3.0, 1996.
[20] Hoheisel, C., and Vogelsang, R. Thermal transport
[30] Hanley, H. J. M., McCarty, R. D., and Haynes,
W. M. The viscosity and thermal conductivity coefficients for dense gaseous and liquid argon, krypton,
xenon, and oxygen. Journal of Physical Chemistry
Ref. Data, 3(4), pp. 9791018, 1974.
coefficients for onecomponent and twocomponent
liquids from time correlation functions computed
by molecular dynamics. Computer Physics Report,
18(8), pp. 169, 1988.
[21] Barker, J. A., Fisher, R. A., and Watts, R. 0. Liq
[31] Younglove, B. A., and Hanley, H. J. M. The viscos
uid argon: Monte carlo and molecular dynamics cal
ity and thermal conductivity coefficients of gaseous
culations. Molecular Physics, 21(4), pp. 657673,
and liquid argon. Journal of Physical Chemistry
1971.
Reference Data, 15(4), pp. 13231337, 1986.
[22] Singer, K., Taylor, A., and Singer, J. V. L. Thermo
[32] Sychev, V. V., Vasserman, A. A., Golovsky,
dynamic and structural properties of liquids mod
E. A. Kozlov, A. D., Spiridonov, G. A., and Tsy
elled by '2lennardjones centres' pair potentials.
marny, V. A. Thermodynamic Properties of Ethylene. National Standard Reference Data Service of
the USSR, A Series of Property Tables. Hemisphere
Publishing Company, New York, 1987.
Molecular Physics, 33(6), pp. 17571795, 1977.
[23] Vogelsang, R., and Hoheisel, C. Comparison of various potential models for simulation of the pressure
of liquid and fluid nitrogen. Physics and Chemistry
[33] Nwobi, 0. C. Molecular Dynamics Studies of Transport Properties and Equation of State of Supercritical Fluids. Ph.D. thesis, The Pennsylvania State
University, August 1998.
of Liquids, 16(3), pp. 189203, 1987.
[24] Nwobi, 0. C., Long, L. N., and Micci. M. M. Molecular dynamics studies of transport properties of su
percritical fluids.
AIAA Paper 970598, January
1997.
[25] Nwobi, 0. C.. Long, L. N., and Michael, M. M.
Molecular dynamics studies of properties of super12