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 Green-Kubo formulae and the virial equation of state, respectively. The MD program uses an effective Lennard-Jones 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 = time-step size e = Lennard-Jones 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 highly-energized 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 Green-Kubo formulae or Einstein relations during equilibrium molecular dynamics (EMD) simulations, 4 ' 14 ' 15 or by conducting suitable non-equilibrium molecular dynamics (NEMD) simulations.16'17 The Green-Kubo 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 long-time tail behavior of the relevant auto-correlation 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 Lennard-Jones (L-J) 12-6 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 Lennard-Jones 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 L-J 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 Lennard-Jones site (site-site), or by defining a central site with corrections for angular dependency. The central site method starts with the Newton-Euler 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 site-site 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 non-diatomics such as ethylene (63 #4) can be modeled as 2-site 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 Lennard-Jones 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 atom-atom 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 L-J potential and L-J force as functions of the separation distance between two particles. The quantities have been non-dimensionalized using the two parameters of the L-J 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 L-J 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: Lennard-Jones 12-6 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) second-order 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 space-dependent correlations to a sufficient degree of accuracy. Several schemes have been used in MD calculations including Gears predictor-corrector and the Verlet algorithms.4 The Verlet schemes are the most widely used algorithms.4 There are several versions of the schemes including the Stoemer-Verlet, the leap-frog Verlet and the velocity-Verlet. The velocity-Verlet algorithm is utilized in the present study and it requires the storage of only atomic positions, velocities and accelerations. The scheme is a dual-step 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 <(*) + l-5t2 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 round-off error and the atomic velocities can be easily scaled in order to control the temperature of the system. The site-site 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 Lennard-Jones 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 cut-off 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 cut-off 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 time-step. 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 cut-off 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 center-to-center potential function. The first term on the right-hand-side of the above equation denotes the ideal-gas 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 non-equilibrium 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. Green-Kubo 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 non-diagonal elements of the stress tensor and it is a collective property of the entire fluid.15 The Green-Kubo 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 (fJ-xyiVxnfJ-yz) 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 Green-Kubo 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 (kinetic-potential) 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 Fortran-77 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 molecules-energy-volume (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 cut-off 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 Lennard-Jones 126 (L-J) potential were used for the pair interactions and the L-J 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 face-centered 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 1000-1200 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 50-300 kg/m3 and temperatures of 100-600 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 Green-Kubo 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 near-critical 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 F4962094-1-0133 and NASA Grants NAGW-1356 Supplement 10, and NGT-10034. 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. 80r-r 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. 155-162, 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. 530-542, 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. 257-292, 1995. 20 300 400 500 [8] Ohara, T., and Aihara, T. Molecular dynamics study on structure of near-critical water. The Japan-U.S. Seminar on Molecular and Microscale Transport Phenomena (J-14), 1993. 600 Temperature (K) [9] Stassen, H., and Steele, W. A. Simulation studies of shear viscosity time-correlation-functions. Journal of Chemical Physics, 102(2), pp. 932-938, 1995. Figure 14: Pressure calculations from MD and experiment for ethylene. [10] Rapaport, D. C. Microscale hydrodynamics: Discrete-particle simulation of evolving flow pat- References tern. Physical Review A, 36, pp. 3288-3299, 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 95-0412, 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. 247-270, 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. 1-19, 1995. cal study of the transient vaporization of an oxygen droplet at sub- and super-critical 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. 111-157, 1994. edition, 1987. 11 Copyright© 1998, American Institute of Aeronautics and Astronautics, Inc. [16] Evans, D. J., and Morriss, G. P. Non-newtonian critical fluids. Journal of Thermophysics and Heat molecular dynamics. Computer Physics Reports, 1, Transfer, 12(2), pp. 1-6, 1998. pp. 297-343, 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. 921-949, 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. 1383-1405, 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. 979-1018, 1974. coefficients for one-component and two-component liquids from time correlation functions computed by molecular dynamics. Computer Physics Report, 18(8), pp. 1-69, 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. 657-673, and liquid argon. Journal of Physical Chemistry 1971. Reference Data, 15(4), pp. 1323-1337, 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 '2-lennard-jones 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. 1757-1795, 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. 189-203, 1987. [24] Nwobi, 0. C., Long, L. N., and Micci. M. M. Molecular dynamics studies of transport properties of su- percritical fluids. AIAA Paper 97-0598, January 1997. [25] Nwobi, 0. C.. Long, L. N., and Michael, M. M. Molecular dynamics studies of properties of super12