Transport Properties of the Rough Hard Sphere Fluid by Olga Kravchenko B.Sc., Moscow Institute of Physics and Technology, Russian Federation, 2000 M.Sc., University of California at Davis, USA, 2005 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Chemistry) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) January, 2013 c Olga Kravchenko 2013 Abstract This thesis examines the dynamics and physical properties of the rough hard sphere fluid (RHSF). The RHSF model consists of spherical particles with well-defined radii that exchange linear and angular momenta upon collision. The simplicity of this model allows for a precise theoretical description that provides a basis for studying fluid properties on the most fundamental level. Extensive molecular dynamics calculations were made of transport properties as functions of density, tracer particle size, and degree of rotational-translational coupling. Results were compared with the smooth hard sphere case and it was found that transport coefficients change significantly due to rotational-translational coupling, becoming stronger with an increase in coupling. Tracer diffusion coefficients were examined for a range of tracer sizes and at various densities. As tracer particles become larger, their diffusion coefficient moves from an Enskog (molecular) to a Stokes-Einstein (hydrodynamic) functional form; the latter depends upon the boundary condition at the surface of the tracer. These boundary conditions for the RHSF are directly proportional to the degree of rotational-translational energy exchange, and can be tuned from “slip” to “stick” values. The validity of several kinetic theory equations have been examined as functions of density and translational-rotational coupling. At very low densities, Boltzmann theory was accurate even at low order except for describing the dependence upon rotational-translational coupling, where low order expansions are less accurate. Enskog theory performed well at low and moderate densities but deviated at larger densities, as expected. For thermal conductivity as a function of translational-rotational coupling even the qualitative behaviour was incorrect. The Enskog predictions for diffusion were also found to be quite poor at low order. Finally, motivated by the results of the thesis, experimental diffusion coefficient data were analysed, especially for nanoparticles. It was shown that defining the correct radius is crucial for describing such systems. In addition, a new formula for predicting tracer diffusion was tested. ii Preface This thesis is based on works that have been published or are prepared to be published by the author of this thesis, Olga Kravchenko, under guidance and in co-authorship with research superviser Prof. Mark Thachuk. Olga Kravchenko is the first author and the main writer of the published works and works in preparation. Chapters 1-4 are introductory chapters. Chapter 1 gives the general overview of the research area, while Chapters 2 - 4 provide more comprehensive theoretical and computational background. Results are presented in Chapters 5, 6 and 7, followed by a summary in Chapter 8. The model description and derivations given in Chapter 2 are expanded versions of the material presented in the publication by Olga Kravchenko and Mark Thachuk, ”The effect of rotational and translational energy exchange on tracer diffusion in rough hard sphere fluid”, J. Chem. Phys. 134, 114310 (2011). In Chapter 3, versions of sections 3.5 and 3.6 have been published as part of the publication by Olga Kravchenko and Mark Thachuk, ”Transport properties of the rough hard sphere fluid”, J. Chem. Phys. 136, 044520 (2012). The derivation presented in section 3.6 was done by Prof. Mark Thachuk. Chapter 4 describes the methodology used in conducting this work. Versions of this description have been published as parts of both publications named above. In Chapters 5 and 6 results of this work are presented. A version of Chapter 5 is published by Olga Kravchenko and Mark Thachuk, ”The effect of rotational and translational energy exchange on tracer diffusion in the rough hard sphere fluid”, J. Chem. Phys. 134, 114310 (2011), and a version of Chapter 6 is published by Olga Kravchenko and Mark Thachuk, ”Transport properties of the rough hard sphere fluid”, J. Chem. Phys. 136, 044520 (2012). Chapter 7 is based on material that is currently being prepared for submission by Olga Kravchenko, as a first author, under the supervision of Prof. Mark Thachuk. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Preface 1.1 An Introductory Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 Rough Hard Sphere Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1 Collision Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 Rotational-Translational Coupling . . . . . . . . . . . . . . . . . . . . . . . 13 3 Theoretical Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.1 Knudsen Number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2 Hydrodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.3 Boltzmann Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.4 Enskog Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.5 Transport Coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.6 Notes on the Derivation of the Transport Coefficients . . . . . . . . . . . . 25 4 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.1 Computational Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.2 Transport Coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 iv 4.3 Finite Size Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Tracer Diffusion in the Rough Hard Sphere Fluid 34 . . . . . . . . . . . . . 35 5.1 Objectives and Background . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.2 Computational Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.3 Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 5.3.1 Data Collection and Convergence . . . . . . . . . . . . . . . . . . . 36 5.3.2 Finite Size Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 5.3.3 Rotational-Translational Coupling . . . . . . . . . . . . . . . . . . . 41 5.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 6 Transport Properties of the Rough Hard Sphere Fluid . . . . . . . . . . 56 6.1 Objectives and Background . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 6.2 Methods 57 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.1 Computational Details . . . . . . . . . . . . . . . . . . . . . . . . . 57 6.2.2 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6.2.3 Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 6.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 7 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 7.1 Motivation 7.2 Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 7.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 8.1 Computational Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 8.2 Experimental Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 107 8.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 v List of Tables 5.1 Values of the mean square displacement for various time differences collected over corresponding number of samples. . . . . . . . . . . . . . . . . . . . . . 5.2 38 For the “fixed I” class, reduced diffusion constants of the bath and tracer for differing tracer sizes Σ/σ and total number of particles N for two reduced densities ρ∗ . Values for N = ∞ are obtained by extrapolating to the infinite volume limit, as described in the text. . . . . . . . . . . . . . . . . . . . . . 5.3 Values of αtracer and αbath for the “fixed I”, and values of Itracer /Ibath for the “fixed α” calculations, as a function of tracer size Σ/σ. . . . . . . . . . . . . 5.4 42 44 For the “fixed α” class with α = 2/5, reduced diffusion constants of the tracer for different tracer sizes Σ/σ and total number of particles N for two reduced densities ρ∗ . Values for N = ∞ are obtained by extrapolating to the infinite volume limit, as described in the text. . . . . . . . . . . . . . . . . . . . . . 5.5 47 For the “fixed α” class with α = 2/3, reduced diffusion constants of the tracer for different tracer sizes Σ/σ and total number of particles N for the reduced density ρ∗ = 0.52. Values for N = ∞ are obtained by extrapolating to the infinite volume limit, as described in the text. . . . . . . . . . . . . . . . . . 6.1 48 Bulk viscosity coefficients ηb∗ for a range of reduced densities ρ∗ and corresponding packing fraction at different values of the dimensionless radius of gyration α. The reported values were calculated using the Eq. (4.10), as described in the text. The values have been extrapolated to the infinite size limit. The low density values denoted by (a) are predictions of Eq. (3.17). . 6.2 Shear viscosity coefficients responding packing fraction η∗ 78 for a range of reduced densities ρ∗ and corat different values of the dimensionless radius of gyration α. The values have been extrapolated to the infinite size limit. The low density values denoted by (a) and (b) are predictions of Eqs. (3.16) and (3.21), respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 vi 6.3 Thermal conductivity coefficients λ∗ for a range of reduced densities ρ∗ and corresponding packing fraction at different values of the dimensionless ra- dius of gyration α. The values have been extrapolated to the infinite size limit. The low density values denoted by (a) and (b) are predictions of Eqs. (3.15) and (3.20), respectively. . . . . . . . . . . . . . . . . . . . . . . . 6.4 Products of the density and diffusion coefficient ρ∗ D∗ densities ρ∗ and corresponding packing fraction 80 for a range of reduced at different values of the dimensionless radius of gyration α. The values have been extrapolated to the infinite size limit. The low density values denoted by (a) and (b) are predictions of Eqs. (3.14) and (3.19), respectively. . . . . . . . . . . . . . . . 6.5 Values of V k1B T Φ(Gηb ) from Eq. 4.6 densities ρ∗ , dimensionless radii of 81 in reduced units for a range of reduced gyration α and number of particles N spanning from 32 to 5300. In order to obtain the reduced bulk viscosity coefficients ηb∗ , 43 η ∗ should be subtracted from each number, where η ∗ is the corresponding reduced shear viscosity coefficient from Table 6.2. . . . . . . 6.6 Reduced shear viscosity coefficients η∗ for a range of reduced densities 82 ρ∗ , dimensionless radii of gyration α, and number of particles N spanning from 32 to 5300. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.7 Reduced thermal conductivity coefficients ρ∗ , λ∗ for a range of reduced densities dimensionless radii of gyration α, and number of particles N spanning from 32 to 5300. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.8 85 Reduced self diffusion coefficients D∗ for a range of reduced densities 88 ρ∗ , dimensionless radii of gyration α and number of particles N spanning from 32 to 5300. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 91 Experimental data for particle sizes measured by TEM with corresponding diffusion coefficients. Numbers in square parentheses are references to published sources. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 vii List of Figures 5.1 The time dependence of the tracer mean square displacement from a “fixed I” (i.e. α = 0) simulation with N = 55297, Σ/σ = 16 and N = 5120, Σ/σ = 5, at the reduced densities ρ∗ = 0.34 (upper panel) and ρ∗ = 0.52 (lower panel), respectively. The bars indicate the estimated errors in the simulation data points while the solid lines are the rational function fits to the data. The dashed lines are the first derivatives of the fit function. The scale for derivatives is not shown. . . . . . . . . . . . . . . . . . . . . . . . . 5.2 39 The time dependence of the tracer mean square displacement from simulations with α varying from 0 to 2/3, with Σ/σ = 5, at the reduced density ρ∗ = 0.52. The bars indicate the estimated errors in the simulation data points while the solid lines are the rational function fits to the data. The dashed lines are the first derivatives of the fit function and approach the diffusion constant value at long times. The scale for derivatives is not shown. 40 5.3 ∗ as a function of inverse simulation Plots of reduced diffusion constants DN box length for a “fixed α = 2/5” simulation for Σ/σ = 10.1 and for the reduced densities ρ∗ = 0.34 (upper panel) and ρ∗ = 0.52 (lower panel). The lines are the linear function fits to the data, whose y intercept gives the extrapolated value of the diffusion constant. . . . . . . . . . . . . . . . . . . 5.4 Extrapolated reduced tracer diffusion constants, ∗ , D∞ multiplied by 43 (R/σ)2 plotted as a function of tracer size, R/σ (R = (Σ + σ)/2) for the reduced densities ρ∗ = 0.34 (upper panel) and ρ∗ = 0.52 (lower panel). The solid and dashed lines denote the hydrodynamic limit with “slip” and “stick” boundary conditions, respectively. The circles, squares, and triangles denote extrapolated results from the “fixed I”, “fixed α = 2/5”, and “fixed α = 2/3” calculations, respectively. The solid curved lines represent the predictions of Eq. (5.2) for the “fixed α” calculations. . . . . . . . . . . . . . . . . . . . . 45 viii 5.5 Using the extrapolated reduced tracer diffusion constants, and employing Eq. (3.2), values of ζ were calculated and plotted as a function of tracer size, R/σ (R = (Σ + σ)/2) for the reduced densities ρ∗ = 0.34 (upper panel) and ρ∗ = 0.52 (lower panel). The solid and dashed lines denote the hydrodynamic limit with “slip” and “stick” boundary conditions, respectively. The circles, squares, and triangles denote extrapolated results from the “fixed I”, “fixed α = 2/5”, and “fixed α = 2/3” calculations, respectively. The crosses denote values calculated from the smooth hard sphere results of Sokolovskii et al. [23] 50 5.6 For the largest-sized tracers, estimates of ζ for α = 2/5 and α = 2/3 from Fig. 5.4 are plotted as a function of α in the upper panel. The smooth sphere result is also included for α = 0. The solid line is the best linear fit through the three points. In the lower panel, using the smooth hard spheres results of ∗ (α)/D ∗ (α = 0) of extrapolated Sokolovskii et al [23] for α = 0, the ratio D∞ ∞ reduced diffusion constants is plotted for α = 2/5 and α = 2/3 as a function of α. The open and filled symbols denote values for the reduced densities ρ∗ = 0.34 and ρ∗ = 0.52, respectively. The circles, upwards triangles, and downwards triangles denote values Σ/σ = 1, 3.17, and 10.1, respectively. The solid line is the function (1 + α)/(1 + 2α). . . . . . . . . . . . . . . . . . . . 5.7 51 Extrapolated reduced diffusion constants for the fixed α calculations as a function of the tracer size. Symbols have the same meaning as in Fig 5.3. The solid curves are the corresponding predictions of Eq. (5.4) with x = 1.4 6.1 54 Mean squared values < ∆G2 > as a function of time scaled for ρ∗ = 0.4298. Each panel corresponds to a different transport coefficient and contains data for four different values of α. . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 59 Plots of reduced constants of thermal conductivity, shear and bulk viscosity as a function of inverse particle number for ρ∗ = 0.2865 and α = 0.004. The lines are linear function fits to the data, whose y intersept gives the extrapolated value of the corresponding transport coefficient. . . . . . . . . 6.3 61 The solid line represents the mean squared values of the off–diagonal elements of the pressure tensor plotted as a function of time at number density ρ∗ = 0.4298 and α = 2/5. Dashed and dotted lines represent the four different fits used to extract the value of the transport coefficient. . . . . . . . . . . . . 62 ix 6.4 Values of the bulk viscosity coefficient ηb∗ multiplied by the radial distribution function at contact g and divided by the value of ηb∗ at ρ∗ = 0 plotted as a function of number density ρ∗ . Each panel corresponds to a different value of α. For the simulation data, plotted as open circles, the values for ρ∗ = 0.01 are taken to be the values for ρ∗ = 0. The solid lines are Enskog theory predictions of Eq. 3.25. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 65 Values of the shear viscosity coefficient η ∗ multiplied by the radial distribution function at contact g and divided by the value of η ∗ at ρ∗ = 0 plotted as a function of number density ρ∗ . Each panel corresponds to a different value of α. For the simulation data, plotted as open circles, the values for ρ∗ = 0.01 are taken to be the values for ρ∗ = 0. The solid and dashed lines are Enskog theory predictions of Eqs. (3.23) and (3.24), respectively. . . . . 6.6 Values of the thermal conductivity coefficient λ∗ 66 multiplied by the radial distribution function at contact g and divided by the value of λ∗ at ρ∗ = 0 plotted as a function of number density ρ∗ . Each panel corresponds to a different value of α. For the simulation data, plotted as open circles, the values for ρ∗ = 0.01 are taken to be the values for ρ∗ = 0. The solid and dashed lines are Enskog theory predictions of Eqs. (3.23) and (3.24), respectively. . 6.7 Squares represent the values of function of number density ρ∗ . ρ∗ D∗ g divided by ρ∗ D∗ at ρ∗ 68 = 0 plotted as a Circles represent the values of ρ∗ D∗ divided by ρ∗ D∗ at ρ∗ = 0 plotted as a function of number density ρ∗ . Each panel corresponds to a different value of α, as indicated. The values for ρ∗ = 0.01 are taken to be the values for ρ∗ = 0. . . . . . . . . . . . . . . . . . . . . . . 6.8 69 Values of the radial distribution function plotted as a function of r/σ for α = 0.2 at various reduced densities. These values were calculated from the simulation with N = 1360 particles. . . . . . . . . . . . . . . . . . . . . . . 6.9 Values of the bulk viscosity coefficient ηb∗ divided by the value of ηb∗ 70 for α = 0.004 plotted as a function of α. Symbols represent the values corresponding to different reduced densities, as shown in figure legend. Solid lines represent the Enskog predictions of Eq. (3.25). . . . . . . . . . . . . . . . . . . . . . . 6.10 Values of the shear viscosity coefficient η∗ divided by the value of η∗ 72 for α = 0 plotted as a function of α. For the simulation data, plotted as symbols, values for α = 0.004 are taken to be the values at α = 0. Symbols represent the values corresponding to different reduced densities, as shown in the figure legend. Solid and dashed lines represent the Enskog predictions of Eqs. (3.23) and (3.24), respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 x 6.11 Values of the thermal conductivity coefficient λ∗ divided by the value of λ∗ for α = 0 plotted as a function of α. For the simulation data, plotted as symbols, the values for α = 0.004 are taken to be the values at α = 0. Symbols represent the values corresponding to different reduced densities, as shown in the figure legend. Solid and dashed lines represent the Enskog predictions of Eqs. (3.26) and (3.27), respectively. . . . . . . . . . . . . . . . 74 6.12 Values of the self-diffusion coefficient D∗ divided by the value of D∗ for α = 0 plotted as a function of α. For the simulation data, plotted as symbols, the values for α = 0 are taken to be the values at α = 0. Symbols represent the values corresponding to different reduced densities, as shown in the figure legend. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 75 Experimental values for the reduced diffusion coefficient multiplied by R∗ in 2 the upper panel and by R∗ in the lower panel as a function of the distance between the centers of the particles, drawn for two different representation of the particle’s size. Circles, squares and crosses represent Ag [79], Ag [80] and Cu2 O [81], respectively. 7.2 . . . . . . . . . . . . . . . . . . . . . . . . . . 98 Millikan diffusion data and the fitted curves plotted as a function of the Knudsen number Kn. The circles represent the values derived from the raw Millikan data, the solid line represents Eq. (7.6) with parameters obtained by Allen and Raabe [77], the dashed line represents Eq. (5.4) with x = 1.1 and the dashed-dotted line represents Eq. (5.4) with x = 1.4. . . . . . . . . 7.3 Values 1/πD∗ η ∗ 99 for the rough hard sphere data given in Chapter 5 and smooth hard sphere data given by Sokolovskii at. al. [23] and Heyes [61] are plotted as a function of Knudsen number Kn. . . . . . . . . . . . . . . . . . 7.4 101 Deviation from the hydrodynamic diffusion coefficient as a function of Kn at various densities, represented by D/Dh , for the rough and smooth hard sphere particles at various densities and degree of rotational-translational coupling. In the top panel the radius of the bath particle is not taken into account, while in the bottom panel the radius of the bath particle is accounted for. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.5 102 Deviation from the hydrodynamic diffusion coefficient as a function of Kn at various densities, represented by (D/Dh − 1) /Kn , for the rough and smooth hard sphere particles at variuos densities and degree of rotational-translational coupling. In the top panel the radius of the bath particle is not taken into account, while in the bottom panel the radius of the bath particle is accounted for. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 xi List of Symbols electric mobility, in Chapter 1 packing fraction, in Chapter 6 q charge of a particle, in Chapter 1 q heat flux, in Chapter 3 kB Boltzmann constant T temperature n number density m mass M total mass of a colliding pair µ reduced mass of a colliding pair R distance between particles at their closest approach σ diameter of a bath particle Σ diameter of a tracer particle I moment of inertia r position v, c velocity ω angular velocity k unit vector p linear momentum R position of the center of mass P momentum of the center of mass J impulse g relative velocity χ parameter, in Chapters 2.1, 2.2 χ scattering angle, in Chapters 3.3, 3.4 α radius of gyration Ek kinetic energy ζ boundary conditions at the surface of a particle L edge length of a simulation cell N number of particles in a simulation cell D self diffusion coefficient xii DB diffusion coefficient in the Boltzmann regime DH diffusion coefficient in the hydrodynamic regime DE Enskog diffusion coefficient D12 binary diffusion coefficient λ mean free path, in Chapters 3.1, 7.2 λ thermal conductivity, in Chapters 3.5, 3.6, 4.2, 5.1, 6.3 ηB bulk viscosity η shear viscosity, in Chapter 3.2, 3.5, 3.6, 4.2, 5.2, 6.2.1 η packing fraction, in Chapter 5.4 Φ interaction potential, in Chapter 1.1 Φ perturbation, in Chapter 3.5 Kn Knudsen number f (r, c, t) distribution function c velocity F external force b second virial coefficient ρ∗ reduced density σ ∗ reduced diameter m∗ reduced mass E ∗ reduced kinetic energy D∗ reduced diffusion coefficient η ∗ reduced shear viscosity ∗ reduced bulk viscosity ηB F0 drag force V relative speed kD drag coefficient C Cunningham correction factor TEM transmission electronic microscopy DMA differential mobility analyzer xiii Acknowledgements First and foremost I would like to thank my research supervisor Mark Thachuk for his invaluable guidance throughout the course of my studies. It is with his encouragement and support I have developed deep understanding of the subject and completed this work. I would like to thank Ruslan Sokolovskii for getting me started and for developing the backbone of the computer code that was employed in this work. Dr. Roman Baranowski and Mr. Brent Gawryluik provided excellent technical support on a number of issues with computing facilities. I thank my committee members for reading my thesis and for their advice on various aspects of this thesis; faculty members, staff and graduate students of the Chemistry Department for creating an inspirational environment that made my stay here enjoyable and productive. I am gratefull to my friends Heidrun Spohr and Victoria Savalei for unforgettable outdoor adventures around beautiful British Columbia. My very special gratitude is to my parents, Lidiya and Valeriy Kravchenko, who made my education their utmost priority; and my brother, Ilya Kravchenko, who led me by example throughout my life. I lovingly thank my husband Sergey Frolov for his support and encouragements. I acknowledge financial support from the Natural Sciences and Engineering Research Council (NSERC) of Canada. All computations were performed using WestGrid computing resources, which are funded in part by the Ca nada Foundation for Innovation, Alberta Innovation and Science, BC Advanced Education, and the participating research institutions (www.westgrid.ca), as well as Scinet and U Shebrook. xiv This thesis is dedicated to my family xv Chapter 1 Introduction 1.1 An Introductory Overview The main objective of the research presented in this thesis is to study transport properties of the rough hard sphere fluid in different fluid regimes. The effects of various microscopic parameters were isolated and studied to reveal fundamental dependencies that lead to better understanding of the macroscopic behaviour of the fluid. A number of research areas, both in theory and in experiment, involve systems that can be described as a collection of particles in a fluid moving freely or under the influence of an external field. In experimental studies, the detailed description of the movement of the particles is often necessary for a correct interpretation of measured values [1–3]. For example, in drift tube experiments, atomic and molecular ions [4–7], aerosols [8] and biological molecules [9] are introduced into a buffer gas with the purpose of sorting them by their mass, size or charge in the presence of an electric field. Ions move through the tube with a velocity which is proportional to the magnitude of the electric field measured experimentally. In the low field limit, the coefficient of proportionality, ε, is related to diffusion via a simple formula ε = (q/kB T )D, where q is the charge of the particle, kB is the Boltzmann constant and T is temperature. Therefore ε, like D, may also be obtained analytically from equilibrium properties of the system [10]. Similarly, in capillary electrophoresis, ionic species are separated based on their hydrodynamic radius, charge and frictional forces while migrating through a capillary filled with bath fluid. One of the challenges of these experiments is to define the radius of the drifting particles. In the low-density limit, the radius of small spherical particles can be estimated from the Boltzmann equation for diffusion coefficient, i.e. 3 DB = 8nR2 kB T 2πµ 1/2 (1.1) in which R is the distance between the centres of the particles at their closest approach, µ is the reduced mass of the tracer and bath particles and n is the number density of the fluid [11]. As shown above, this expression is used to estimate the spherical radius of ionic species in gas-phase ion mobility experiments. On the other hand, for large tracers in the 1 hydrodynamic regime, the diffusion coefficient for spherical particles is given by the Stokes - Einstein expression D= kB T , ζπηR (1.2) in which η is the shear viscosity, R is the hydrodynamic radius and ζ represents the boundary conditions at the surface of the tracer and may adopt values from 4 to 6, referred to as “slip” and “stick” boundary conditions, respectively [11]. In experiments such as capillary electrophoresis, ζ = 6 is the universally adopted value for the boundary conditions. Even though Eqs. (1.1) and (1.2) are used routinely for the systems described above, it is not always clear which theory is the most suitable for a specific system in a given set of conditions. For example, nanoparticles may not be big enough to be described by hydrodynamic equations, yet not small enough for the molecular regime. Biological molecules in a gas phase, on the other hand, may be too big to be accurately described by the molecular theory, which is usually applied. In other cases, when particles are in the hydrodynamic regime, the choice of boundary conditions has to be considered more carefully, because six is not the only possible value. In real molecules, the radius of the particles is not well defined, that may lead to uncertainties in interpreting experimental results. One of the ways to obtain the radius is to deduce it analytically from the diffusion coefficient. Given that such ambiguities exist, the following questions come to mind. If the radius is to be found from the diffusion coefficient, which expression would then give the most appropriate description for diffusion of a specific system and allow one to interpret the radius most accurately? Which boundary conditions should be used in the hydrodynamic case for a given system? In other words, is it possible to identify a set of criteria, according to which an accurate analytical description of the system can be provided? A number of studies conducted by other researchers tried to address these questions. In the following paragraphs, an overview of these works will be given, followed by description of the common difficulties with these studies. The work in this thesis overcomes many of those difficulties and provides a solid foundation to address the questions given above. Schmidt and Skinner [12] studied the behaviour of very large and massive solutes in a dense fluid of spherical solvent particles and found that diffusion could be described by Eq. (1.2) with “slip” boundary conditions. In a later work the same authors studied the effect of the solute-solvent potential on the boundary conditions on the surface of different tracer models and found that large tracers approach “slip” boundary conditions. However, when the tracer is represented by a cluster of particles, “stick” boundary conditions are obtained when the hydrodynamic radius is considered [13]. Masters and Madden studied the diffusion of a large particle in a fluid of small light particles and found that the choice 2 of potential affects the boundary conditions at the surface of the tracer [14]. Bhattacharyya and Bagchi have studied diffusion of a tracer in a dense fluid and found that Stokes Einstein behaviour manifests itself when the solute is only 2 - 3 times bigger than the solvent molecules [15]. Li [16] studied the validity of Eq. (1.2) for cluster particles of different sizes in a Lennard-Jones system and showed that diffusion depends upon both particle size and the strength of the attractive tracer-bath interaction. It was found that Eq. (1.2) breaks down for small particles but eventually does predict the diffusion for large enough particles. Particle friction has also been studied in different systems. Uvarov and Fritzsche [17, 18] developed a semiphenomenological Fokker-Planck approach to the diffusion of small macromolecules in solution as a function of mass ratio between solute and solvent particles, and observed deviations from slip for certain mass ratios. In another work Bhattacharyya [19] employed a mode coupling theoretical description of microscopic friction to study the transition between molecular and hydrodynamics regimes as a function of size, mass, and solute - solvent interaction parameters and found that crossover to the Stokes - Einstein regime occurs sooner if the mass of the tracer is increased along with its size. Ould-Kaddour and Levesque [20] studied the behaviour of the velocity autocorrelation function for Lennard-Jones tracer particles in a fluid, and tested scaling relationships for both short and long time behaviour. The mass dependence of transport coefficients was studied in binary Lennard-Jones mixtures with high mass asymmetries by Fenz et al. [21]. They found that the results satisfied a generalized Stokes-Einstein relation. Lee and Kapral [22] studied the size and mass dependence of the friction and diffusion of a Brownian particle in a mesoscopic solvent using a hybrid MD-multiparticle collision model with truncated Lennard-Jones interactions. At last, a study by Sokolovskii et al. [23] showed that the diffusion coefficient starts to deviate from behaviour predicted for the molecular regime as soon as tracer particle becomes slightly bigger than bath particles. In most studies listed above there are common difficulties that leave room for further investigation. Particularly, studies based on non-hard sphere models where radius is not precisely defined are always susceptible to an error and may produce ambiguous results. The work of Cappelezzo et al. [24] studies the validity of the Stokes - Einstein relation for self-diffusion in a pure Lennard – Jones fluid and compared values of the boundary conditions ζ obtained by using different definitions of the hydrodynamic radius at various fluid densities. They showed that the value of ζ depends upon particular definition chosen for the hydrodynamic radius. For non-hard sphere potentials, the radius can be defined as the average distance of collision averaged over a typical distribution of collision energies, which may vary. For large enough particles the difference between these definitions may not be noticeable. But for small particles and for studying the transition from molecular to hydrodynamic regimes, defining the radius one way or the other may affect the interpretation of the results. 3 Another common problem with the studies mentioned above is accounting for dynamic finite size effects. When a particle moves through the fluid and collides with other particles the resulting correlated motion propagates throughout the whole fluid. Periodic boundary conditions used in molecular dynamics simulations cut off some of the correlations and may produce a significant error, even when the size of the simulation cell is large. Finite size effects can greatly affect tracer diffusion, yet in almost all studies, except a few that will be named further, finite size effects have not been removed. Finite size effects have been investigated by Ould-Kaddour et al. [20] who observed a dependence of the friction coefficient upon the size of the simulation cell used in simulations of massive particles suspended in a Lennard-Jones fluid. Fushiki [25] suggested a method of extrapolating the diffusion coefficient to the infinite volume. Sokolovskii et al. [23] investigated the effects of a finite system size upon diffusion in the smooth hard sphere fluid and showed that the extrapolation method introduced by Fushiki is crucial for correct interpretation of the computational results, particularly in determining the boundary conditions at the surface of the tracer. Finally, in order for results of the simulations to be reliable, a great number of computational jobs have to be run for sufficiently long time. This used to be a limitation until the last decade, when computational resources improved dramatically and made it possible to gather enough data for the accurate statistical analysis presented in this thesis while accounting for the difficulties mentioned above. In addition to the computational studies named above there were analytical works that tried to bridge the molecular and hydrodynamic regimes. Mehaffey and Cukier derived repeated-ring kinetic theory to account for correlations in the fluid, and calculated the diffusion coefficient for a large hard sphere tracer in a dense fluid of small hard spheres. They suggested the transition region could be described by the sum of the Enskog value and the hydrodynamic value with ζ = 5 [26]. Beijeren and Dorfman then pointed out that including more terms in the repeated-ring sum leads to “slip” boundary conditions [27]. To study the subject, a theoretical model that removes ambiguities in the description of the system is needed. Such a model should allow one to isolate and study the effect of certain parameters on the entire system. For example, the effect of increasing the size of the particle in comparison with the size of the bath fluid, changing the density of the bath fluid, etc.. The hard sphere fluid is a simplified model fluid often invoked to study certain properties of real fluids on a fundamental level. The model consists of spherical particles with well defined radii that move freely and collide. Particles do not interact with one another at any time except at the moment of collision. At that moment, instantaneous repulsive forces arise and particles change direction. The interaction potential Φ between the particles in 4 the hard sphere fluid is defined precisely by the following equation: Φ= ∞ ;r < σ 0 ;r ≥ σ where r is the relative distance between the particles and σ is the distance of closest approach (or the diameter of the particles when the particles are identical in size) [28]. In real fluids, the interactive potential consists of a long-range attractive part and shortrange repulsive part. Even though the hard sphere potential is lacking the attractive part, it is still a very effective model for the study of transport properties, such as diffusion, shear and bulk viscosity and thermal conductivity, for the following reason. There are two ways to transport energy and momentum through a fluid: at very low density it is done mainly by transporting the molecule itself, but in all other cases collisions play a major role in propagating the change in properties throughout the fluid. Since collisions are defined by the repulsive part of the interactive potential, and the simplest way to represent repulsion is by introducing the infinite repulsive wall, the hard sphere model is a suitable and convenient model for studying properties that rely on collisional transfers. The hard sphere fluid model allows one to study particles on a microscopic level and obtain reliable computational data for various dependencies, that may be useful in predicting the behaviour of the particles in a more realistic environment. To this point, most of the studies have been done on the smooth hard sphere fluid for which rigorous theoretical descriptions have long existed [11, 29, 30]. In this work we focus on the lesser studied rough hard sphere fluid. In the rough hard sphere fluid, unlike in the smooth hard sphere fluid, individual particles are allowed to rotate. Thus when particles collide they exchange not only translational but also rotational energy and momenta. Structurally the smooth hard sphere and rough hard sphere fluids are the same, and therefore, by comparing the two, it is possible to study the effect of particle rotation upon various macroscopic properties of the fluid. Rough hard spheres allow one to introduce the simplest internal degree of freedom, as well as the coupling between internal and translational states. Even though the rough hard sphere model is considered too inelastic in comparison with real molecules [31], it gives a good qualitative description of processes that take place in real fluids. The rough hard sphere model was first introduced by Bryan in 1894 [32] and further developed by Pidduck [33], who gave the precise kinetic description of the collision event. The model is described in detail in the works of Chapman and Cowling [31] who solved the Boltzmann equation using Chapman-Enskog theory and predicted transport coefficients for the rough hard sphere fluid. This work was expanded by Dahler and coworkers who also solved Enskog equation for the rough hard sphere fluid [34–36]. Berne extensively studied 5 the effect of the roughness of the particles on rotational - translational motion [37–40]. For large enough particles and high enough density hydrodynamic theory is appropriate [11]. The intermediate region between these two limits is not well described. In this thesis an attempt is made to bridge the two limiting cases and observe the transition from one limit to another. It is particularly important that particles of the hard sphere fluid model in general, and rough hard sphere model in particular have well-defined radii; the only thing we add to the smooth sphere model to transform it into the rough hard sphere model is rotationaltranslational coupling. In real systems, especially for non-spherical particles, the physical description of a radius is a matter of debate. The treatment of experimental data may be affected by the way of defining the radius and thus it may be difficult to obtain experimental data with a high degree of accuracy. This is why computational experiments are particularly useful in isolating and studying the effects of various conditions upon fluid properties; the results from such studies are highly accurate and reliable. This thesis is organized in the following way. Chapters 1-4 are introductory chapters. Chapter 1 provides general overview of the research area. In Chapters 2-4 comprehensive background information about theoretical and computational approaches used in this study is given. In Chapters 5-7 results of the studies are presented and discussed. In Chapter 8 results of this thesis are briefly summarized and possible future work is suggested. Chapter 5 is based on the paper “The effect of rotational and translational energy exchange on tracer diffusion in rough hard sphere fluids” published in 2011 [41]. In this paper a study of diffusion of a tracer particle in a rough hard sphere fluid as a function of particle size at two different densities is presented. As the tracer grows in size it is expected that the diffusion constant will follow the Stokes-Einstein hydrodynamic result of Eq. (1.2). It is known that smooth hard spheres adopt slip boundary conditions in the hydrodynamic regime [23]. Results of this thesis show that rough hard spheres adopt boundary conditions proportional to the degree of translational – rotational energy exchange. Spheres for which this exchange is the largest adopt stick boundary conditions while those with more intermediate exchange adopt values between slip and stick limits. This dependence is found to be almost linear. The changes in the diffusion constants as a function of this exchange have also been examined. It was found that the dependence is stronger than that suggested by the low-density, Boltzmann result. Compared with smooth hard spheres, real molecules undergo inelastic collisions and have attractive wells. Rough hard spheres model the effect of inelasticity and show that even without the presence of attractive forces, the boundary conditions for large particles can deviate from slip and approach stick. In Chapter 6, results based on the paper ”Transport properties of the rough hard sphere fluid” published in 2012 [42] are discussed. Systematic benchmark calculations of the transport properties of the rough hard sphere fluid have been conducted. Molecular dynamics 6 was used to study diffusion, shear and bulk viscosity and thermal conductivity coefficients as functions of density and degree of rotational-translational coupling. As well, the validity of several kinetic theory equations at various levels of approximation was examined. Particularly, expressions from Enskog theory using different numbers of the basis sets in the representation of the distribution function were tested. It was found that Enskog theory works well at low density but not as well at higher densities. Enskog theory as a function of degree of rotational translation coupling has also been examined and found that even at low density the agreement with simulation results is generally poor. In comparison with smooth sphere results, the transport coefficients change significantly due to rotational translational coupling. In Chapter 7, several sets of experimental data for the diffusion of nanoparticles were analysed and compared with computational results presented in this thesis in order to test the applicability of various theoretical equations in different fluid regimes. It was found, that in the case of Ag and Cu2 O particles in a nitrogen bath hydrodynamic behaviour doesn’t manifest itself which suggests that Boltzmann theory may be more appropriate to describe a fluid with a similar particle size distribution. An attempt was made to describe the transition from molecular to hydrodynamic regimes by fitting the data from the Millikan oil drop experiment with two different phenomenological expressions, one of which is well established, and the other one proposed for the smooth hard spheres by Sokolovskii et al. [23] that was known to be applicable at low density, and compare the outcomes. It appeared, that both expressions perform well in the hydrodynamic regime, giving expression by Sokolovskii et al. a wider range of applicability. As well, different definitions of the representative particle size have been tested. It was shown that including the size of the bath particle in the definition produces trends that are theoretically predicted for both molecular and hydrodynamic regimes, while excluding it does not lead to correct trends. 7 Chapter 2 Rough Hard Sphere Model 2.1 Collision Dynamics In this chapter the nature of collisional dynamics taking place in the rough hard sphere fluid will be explained. Since the smooth and rough hard sphere fluids share common properties important in this work, the generalized theoretical description of the collisions that occur in both of them will be provided, and the differences between the two cases will be outlined. The derivation presented in this chapter has been done in various forms by several other authors, [37–40, 43] but it was completely re-derived here in terms specific for this work. Consider a system of spherical particles set in motion by randomly assigning velocities in a way that corresponds to a Maxwell distribution. Particles move linearly and with constant velocity from one collision to the next, and exchange energy and momentum during collisions. In order to predict the linear trajectory of a particle the positions, velocities and angular velocities of the particle before and after collision must be related. These relations will be derived below. Let m, σ, and I be the mass, diameter, and moment of inertia, and r, v and ω the position, velocity and angular velocity, respectively, of a hard sphere particle. The relative position and velocity between two particles, labelled “1” and “2”, are given by r = r2 − r1 and v = v2 − v1 , respectively. Define k = −r/|r|, the unit vector pointing from the center of particle “2” to the center of particle “1”. Denote parameters before and after collision as unprimed and primed, respectively. The conservation of linear momentum can be written as m1 v1 + m2 v2 = m1 v 1 + m2 v 2 . (2.1) The total angular momentum consists of the contributions of each individual particle, I1 ω1 and I2 ω2 , as well as the orbital angular momentum arising from orbital movement of the particles around each other. A particle moving with the linear momentum p at the point r relative to an arbitrary origin has angular momentum r × p. Therefore, the conservation of total angular momentum for the two particles relative to an arbitrary origin can be written 8 as I1 ω1 + I2 ω2 + r1 × p1 + r2 × p2 = I1 ω 1 + I2 ω 2 + r 1 × p 1 + r 2 × p2 (2.2) To write this expression in the center of mass coordinate system, let M = m1 + m2 be the total mass, then M R = m1 r1 + m2 r2 , (2.3) where R is the position of the center of mass. The positions of the particles are then m2 r, M m1 r. = R+ M r1 = R − (2.4) r2 (2.5) Differentiating Eqs. (2.4) and (2.5) and multiplying by the corresponding masses yields p1 = p2 = m1 P − µv, M m2 P + µv, M (2.6) (2.7) where P is the momentum of the center of mass, and µ = m1 m2 /M is the reduced mass. Substituting Eqs. (2.6) and (2.7) into Eq. (2.1) shows that P = P , that is the center of mass momentum is conserved during the collision. The momenta of the particles after the collision can be written as p1 = p2 = m1 P − µv , M m2 P + µv . M (2.8) (2.9) Subtracting either Eq. (2.6) from (2.8) or Eq. (2.7) from (2.9) gives p2 − p2 = p1 − p1 = µ(v − v) = J, (2.10) where J is the impulse transferred from one sphere to another during the collision. Using relations (2.6) and (2.7) and considering that R and r do not change at the point of contact during the collision, the angular momentum conservation law Eq. (2.2) becomes I1 ω1 + I2 ω2 + r × p = I1 ω 1 + I2 ω 2 + r × p , (2.11) in which p = µv. 9 The conservation of total kinetic energy can be written as 1 1 1 1 m1 v12 + m2 v22 + I1 ω12 + I2 ω22 = 2 2 2 2 2 2 2 2 1 1 1 1 m1 v1 + m2 v2 + I1 ω1 + I2 ω2 . (2.12) 2 2 2 2 The system of Eqs. (2.1), (2.11) and (2.12), that represent conservation of linear and angular momenta and total kinetic energy, fully describe the hard sphere collision. The goal is to solve this system of equations and obtain parameters necessary to predict the motion of the particles in the fluid. In terms of J the velocities and angular velocities after the collision can be represented as J , m1 J , = v2 + m2 σ1 = ω1 + k × J, I1 σ2 = ω2 + k × J. I2 v1 = v1 − (2.13) v2 (2.14) ω1 ω2 (2.15) (2.16) Equations (2.13) and (2.14) follow directly from Eq. (2.10). Equation (2.15) and (2.16) can be derived by noting the collision creates torque at the points of contact that in turn generates a change in angular momentum given by R1 × p1 − p1 for particle “1” and R2 × p2 − p2 for particle “2”, where R1 and R2 are radius-vectors from the centres of the spheres to their points of contact at collision. Taking into account that σ1 R1 = − k, 2 σ2 R2 = k, 2 (2.17) (2.18) and using Eqs. (2.13) and (2.14) then gives Eqs. (2.15) and (2.16). Note that by substituting v1 and v2 from Eqs. (2.13) and (2.14) into Eq. (2.1), or ω1 and ω2 from Eqs. (2.15) and (2.16) into Eq. (2.11) and recalling that at the point 2) k shows that the form of Eqs. (2.13)-(2.16) guarantee linear and of collision r = − (v1 +v 2 angular momenta are conserved regardless of the value of J. The value of J is fixed by the conservation of total energy. Substituting Eqs. (2.13)-(2.16) into Eq. (2.12) and simplifying leads to χ k×J 2 2 + σ1 σ2 J2 ω1 + ω2 · k × J + + v · J = 0, 2 2 2µ (2.19) 10 in which σ12 σ22 + I1 I2 1 4 χ= . (2.20) Now define g=v− k × (σ1 ω1 + σ2 ω2 ) , 2 (2.21) which is composed of two contributions. The first is the relative velocity of the spheres and the second gives the difference in the velocities of the surface of the spheres at contact due to their rotational motion. Thus g is the total relative velocity of the points of contact of the two spheres at collision. Rewriting Eq. (2.19) in terms of g gives g·J + J2 χ + k×J 2µ 2 2 = 0. (2.22) The impulse J can be represented as a vector sum of the components parallel and perpendicular to the unit vector k, that is J = J · k k, J⊥ = J − J , (2.23) (2.24) 2 . The relative velocity g can be presented in a so that k × J = k × J⊥ and J 2 = J 2 + J⊥ similar way. Written in terms of parallel and perpendicular components, Eq. (2.22) becomes g J + J⊥ · g⊥ + 1 χ 2 2 J 2 + J⊥ + J⊥ = 0. 2µ 2 (2.25) Let’s separate expressions for the parallel and perpendicular components of J. The equation for the parallel component is g J + 1 2 J = 0, 2µ (2.26) which gives two possible outcomes for J , namely J = 0, (2.27) J = −2µg . (2.28) 11 The equation for the perpendicular components is J⊥ · g⊥ + 2 J⊥ χ 2 + J⊥ = 0, 2µ 2 (2.29) which leads to the following values for J⊥ , J⊥ = 0, J⊥ = − (2.30) 2µ g⊥ . µχ + 1 (2.31) First, note that even though the parallel component of J can be equal to 0, this solution implies J = µ v − v = 0, that is v = v , which means that the velocity along the k direction stays the same after the collision. Since this result is unphysical, it is not taken into account for constructing the description of the collision event. It can be referred to as the result of the case when collision between two particles does not occur. Therefore, only two possibilities are left. One is attributed to the smooth sphere case and described by J = −2µg = −2µv , J⊥ = 0. Using these expressions with J = µ v − v (2.32) (2.33) shows that v = −v and v⊥ = v⊥ . Thus, for a smooth sphere collision the component of the parallel velocity changes sign while the perpendicular stays unchanged. In this case, k × J = k × J = 0 and Eqs. (2.15) and (2.16) show that ω1 = ω1 and ω2 = ω2 , that is no change in angular momentum can be induced by the collision. In this case Eq. (2.12) simplifies to 2 2 1 1 1 1 m1 v12 + m2 v22 = m1 v 1 + m1 v 2 , 2 2 2 2 (2.34) corresponding to the case of perfectly elastic collisions where kinetic energy is conserved. The other possible solution represents the rough hard sphere case, which is the main subject of our studies, and is given by J = −2µg = −2µv , J⊥ = − (2.35) 2µ 2µ k g⊥ = − v⊥ − × (σ1 ω1 + σ2 ω2 ) . 1 + µχ 1 + µχ 2 (2.36) In this case, collisions are inelastic and translational and rotational conversion can occur. As with the smooth hard sphere case J = µ v − v = −2µv therefore v = −v . 12 Using Eq. (2.21) to determine an expression for g in terms of post-collisional variables leads to (along with Eqs. (2.13)-(2.16)) g −g = v −v− = = 1 k × σ1 ω1 + σ2 ω2 − k × (σ1 ω1 + σ2 ω2 ) 2 J − χk × k × J µ J µχ + 1 k+ J⊥ , µ µ (2.37) where the relation k × k × J = J k − J = −J⊥ has been used. Now rewrite this expression in terms of the components of g, that is g −g = − 2µg k+ µ µχ + 1 µ − 2µ 1 + µχ g⊥ = −2 g k + g⊥ = −2g. (2.38) From the last expression it follows that in the rough hard sphere case g = −g, that is the relative velocity at the point of contact of the spheres is reversed upon collision. 2µ Note that Eqs. (2.35) and (2.36) can be combined to give J = − 1+µχ g + µχ g · k k . Because J⊥ = 0 torques are generated and Eqs. (2.15) and (2.16) show that angular velocities can change when rough spheres collide. 2.2 Rotational-Translational Coupling As shown in Sec. 2.2, rough hard sphere collisions are inelastic and when two particles collide, they exchange both translational and rotational energy and momenta. In some cases, especially for glancing collisions, this exchange can be very large. Using Eqs. (2.13) and (2.14) the change in kinetic energy during the collision is given by ∆Ek = = = 2 2 1 1 1 1 m1 v 1 + m2 v 2 − m1 v1 2 − m2 v2 2 2 2 2 2 1 2 J +v·J 2µ 1 J⊥ · [J⊥ + 2µv⊥ ], 2µ (2.39) and based on Eqs. (2.15) and (2.16) the change in rotational angular momentum is given by ∆Iω = I1 ω1 + I2 ω2 − I1 ω1 − I2 ω2 = ( σ1 + σ2 )k × J⊥ , 2 (2.40) 13 where the final expressions in each case were obtained by noting J = −2µv and k × J = k × J⊥ . For smooth hard spheres, J⊥ = 0, and the equations above show that no change in kinetic energy or rotational angular momentum occurs in this case, as expected. Define α = 4I/mσ 2 , the reduced moment of inertia of the spherical particle. The value of α ranges from zero to 2/3 as the mass of the sphere is concentrated from its center to its surface, respectively. A sphere with uniform mass density corresponds to α = 2/5. In terms of α, the parameter χ can be written as χ= 1 1 + . m1 α1 m2 α2 (2.41) Consider the limit of the rough sphere model when α becomes zero. Specifically, let the two particles, “1” and “2” be referred to as tracer and bath particles, respectively, with the corresponding coefficients α1 and α2 , and let α1 → 0. In this limit, the first term of Eq. (2.41) dominates so lim χ = α1 →0 1 σ2 = 1, m1 α1 4I1 (2.42) and Eq. (2.36) gives lim J⊥ = − α1 →0 8I1 1 v⊥ − k × (σ1 ω1 + σ2 ω2 ) . 2 2 σ1 (2.43) When a tracer is infinitely dilute in a bath, it is possible to change its size without affecting the overall density of the fluid. Consider the case where the moment of inertia of the tracer I1 is held fixed and its size is allowed to increase. In this limit, α1 → 0 and as the equation above shows J⊥ → 0 also. This result in Eq. (2.39) then gives Ek → 0; the collisions become elastic. This means that by changing α1 we can change the degree to which translational and rotational energies are coupled. When α1 = 0 there is no exchange of energy between translations and rotations. This exchange grows as α1 becomes larger. However, it follows from Eqs. (2.40) and (2.43) that for fixed I1 and growing σ1 ∆Iω = 2I1 k × (k × ω1 ) = −2I1 ω1⊥ . (2.44) Although there is no exchange between translational and rotational energy when α1 → 0, this is not the case for the exchange of rotational angular momentum between the tracer and the bath upon collision. Once the same limit is applied to Eqs. (2.15) and (2.16) the following result is obtained 14 ω1 = ω1 + 2k × k × ω1 = ω1 − 2ω1⊥ = ω1 − ω1⊥ , ω2 = ω2 . (2.45) (2.46) This shows that the change in rotational angular momentum is transferred solely to the tracer, while the bath conserves its rotational angular momentum. At first sight this may seem to violate the conservation of angular momentum, until we recall that the total angular momentum in the system arises from two sources: (i) rotational angular momentum of the particles, and (ii) orbital angular momentum of the form r × p for each particle. As the size of the tracer increases, the magnitude of J⊥ decreases but this impulse is applied at the edge of the tracer. The increasing value of σ1 represents the growing lever arm upon which J⊥ acts, in such a proportion that the resulting contribution to the orbital angular momentum remains non-zero (and becomes independent of σ1 ). The net result is that this change in the orbital angular momentum is transferred wholly to the tracer particle, with the effect that the perpendicular component of its angular velocity is reversed. Now consider a pure rough sphere fluid with m1 = m2 = m, I1 = I2 = I, v1 = v2 = v and α1 = α2 = α, in the limit where α → 0. In this case lim J⊥ = − α→0 4I σ v⊥ − k × (ω1 + ω2 ) . 2 σ 2 (2.47) Another way of achieving the limit α → 0 is to keep σ fixed and let I → 0. This corresponds to condensing the mass of the particle to its center. When I → 0, one can obtain J⊥ = 0, ∆Ek = 0, ∆Iω = 0, and ω1,2 = ω1,2 + k k × (ω1 + ω2 ) − 2 v⊥ . σ (2.48) In this limit, there is no change in total angular momentum, and the total rotational energy vanishes as I → 0 but there are still changes in rotational velocities upon collision. In summary, the limit α → 0 does not extinguish exchange among the rotational degrees of freedom. Chapman and Cowling [31] pointed out this effect for the pure rough sphere fluid. For this reason, the rough sphere model does not approach precisely the smooth sphere one in the α → 0 limit. While this difference does affect some transport properties, such as the thermal conductivity, it is not expected to affect diffusion. This issue will expanded upon in Chapter 3. Finally, when performing simulations of a tracer particle in a bath, most of the collisions are between bath particles, in which case “1” and “2” label the same type of particle and the expressions for µ and χ simplify. In this case, σ1 = σ2 = σ, where σ is the diameter of 15 a bath particle. Only for the less frequent tracerbath collisions does one require the more general equations, and in this case, one of σ1 or σ2 will be equal to Σ, the diameter of the tracer particle. We will use these notations for the diameters of the bath and the tracer particles for the rest of the thesis. 16 Chapter 3 Theoretical Background 3.1 Knudsen Number A fluid can be treated as an ensemble of discrete particles that interact on a microscopic level or as a continuous medium, in which individual particles are ignored. In order to give the most suitable theoretical description of the fluid regime, the Knudsen number, the ratio of the mean free paths to a representative length of the particle size, that is Kn = λ , R (3.1) is defined [28]. The length can be, for example, the radius or the diameter of the particle. The Knudsen number Kn shows how far a particle travels, in comparison with it’s own size, before it collides with another particle. When the fluid density is low or the tracer particle is small in comparison to the mean free path, that is Kn ≥ 10, the fluid is considered to be a collection of particles and can be described by the Boltzmann Equation. For very large Kn values collisions between particles are rare the collision term in the Boltzmann Equation can be set to zero. When the ratio of the mean free path to the size of the tracer is very small, i.e. Kn < 0.1, the fluid can be treated as a continuum and described by the Navier-Stokes Equations. In the intermediate region of 0.1 ≤ Kn < 10 the fluid is said to be in the slip flow regime, where it can generally be treated as a continuum while accounting for microscopic properties may be useful in some cases [11]. The Knudsen number is a useful parameter for defining the transition region between the limiting cases; it allows one to identify where the behaviour of the fluid changes from molecular to hydrodynamic and thus to choose the most appropriate theory to describe it. One of the goals of this work is to study whether the Knudsen number can be used as a universal identifier of the behaviour of the fluid. 3.2 Hydrodynamics At small Knudsen numbers, the fluid can be described by hydrodynamic theory using Navier-Stokes equations which refer to the flow in the fluid rather than treat the individual particles that constitute the fluid. In hydrodynamic theory, for large enough spherical 17 particles, the diffusion coefficient is found to be [11] DH = kB T , ζπηR (3.2) in which R is the radius of the tracer particle, η is the shear viscosity coefficient of the bath, and the value of ζ depends upon the boundary conditions at the surface of the tracer: ζ = 6 is said to be “stick” boundary conditions and ζ = 4 corresponds to “slip” boundary conditions. Kinetic theories and computer simulations suggest that the boundary conditions at the surface of the tracer determine the value of ζ, and the precise manner in which they are determined is one of the subjects of study in this work. Previous simulation studies have observed the slip limit [12, 14, 15, 44–46], while others the stick limit [26, 47, 48], as described in Chapter 1. Simulation results by Sokolovskii et al. [23] have determined that ζ corresponds to the “slip” value for smooth hard sphere tracers in a smooth hard sphere bath. This is consistent with the earlier results of Schmidt and Skinner [12] who found slip boundary conditions for a spherical tracer in a Lennard-Jones fluid. Regardless of the boundary conditions employed at the surface, the diffusion coefficient is inversely proportional to a particle’s size, unlike in the Boltzmann limit, where this dependence is quadratic (this will be expanded upon in further chapters). As a tracer particle grows in size, the diffusion coefficient transitions from Enskog to the hydrodynamic limit and the functional dependence upon the size of the particle changes along the way. One of the goals of this thesis is to test how well the diffusion coefficient follows Eq. (3.2) while varying external parameters such as the density of the fluid and the size of the tracer particle. Results of our studies of this subject are presented in Chapter 5. 3.3 Boltzmann Equation The macroscopic properties of the fluid, such as momentum, energy, pressure, transport coefficients, etc. can be inferred from the velocity distribution function f (r, c, t), where r and c are the position and velocity of a given particle at time t, which in turn can be found from the Boltzmann equation when Kn is large enough. Consider a fluid that consists of particles subjected to an external force F . At time t the velocity of a particle located at r1 is equal to c1 , and the number of particles in a volume element dc1 dr1 is equal to f (c1 , r1 , t)dc1 dr1 . After time dt, in the absence of collisions, the velocity of the same particle will change to c1 + F m dt, the position becomes r1 + c1 dt, and the number of particles located in the volume element dc1 dr1 will be given by f c1 + F m dt, r1 + c1 dt, t + dt dc1 dr1 . The number of particles at time t is going to change 18 by the time t + dt due to collisions: some of the particles will be pushed out of the given volume, and some will enter. Therefore, the net change in the number of particles will be given by f c1 + F dt, r1 + c1 dt, t + dt m − f (c1 , r1 , t) dc1 dr1 = ∂e f dc1 dr1 dt. ∂t (3.3) Dividing this expression by dc1 dr1 dt and setting dt → 0, Boltzmann’s equation can be obtained [31] F ∂f ∂e f ∂f ∂f + = + c1 · · , ∂t ∂ r1 m ∂ c1 ∂t (3.4) where ∂e f /∂t represents the rate of change of the distribution function due to collisions. The Boltzmann equation is based upon two important assumptions: (1) the motion of the particles is uncorrelated, and (2) the distance between the centres of the particles is ignored, i.e. particles are treated essentially as point particles. In the first assumption, also called the “molecular chaos approximation”, it is presumed that the distribution function of a single particle is uncoupled from that of the other particles. In other words, for two colliding particles it is enough to know position and velocity of one particle in order to calculate it’s trajectory and the collision history of the other particle is ignored. This assumption works best at low density because in denser systems collisions occur more often and correlations develop. The second assumption addresses the size of the particles. The particles are considered to be very small in comparison with the mean free path and their size is not taken into account. That is, when particles collide, the collision is said to occur at a certain point and the distribution function is written only at that point rather than at two distinct coordinates for the centres of the particles. This approximation works well when the fluid is dilute, however, it breaks down when the density, or the size of the particles increases. In this case the two colliding particles have corresponding distribution functions that differ because their centres do not coincide, and this difference becomes significant. In Eqs. (3.3) and (3.4) the parameters have indices “1” because it is enough to know the position and velocity of one of the colliding particles to find the distribution function. Now let’s look closer at the collision term ∂e f /∂t which accounts for interactions between the particles during collision and is given by ∂e f = ∂t f c1 , r1 , t f c2 , r1 , t − f (c1 , r1 , t) f (c2 , r1 , t) gσ (g, χ) sin(χ)dχd dc2 . (3.5) Equation (3.5) describes the collision between two particles, labelled “1” and “2”. Unprimed 19 and primed quantities denote values before and after a collision between the two particles. The relative velocity of the two particles is g = c2 − c1 . The scattering angle χ indicates how much the final relative velocity is deflected off the initial one. Note that the distribution functions are evaluated at the same position r1 , that is the radius vector that defines the position of the first particle, consistent with the second assumption described above. The terms f1 and f2 represent the probabilities that both particles (1) and (2) are at position r1 with velocities c1 and c2 at time t; similarly the terms f1 and f2 are probabilities that both particles are at position r1 with velocities c1 and c2 at time t + dt. The second term in Eq. (3.5) which can be written in condensed notation as f1 f2 gσ sin χdχd dc2 is called the “loss term”. The loss term indicates the number of particles lost from the volume element at c1 due to collisions. The first term at Eq. (3.5), written as f1 f2 gσ sin χd dc2 is called the “gain term” and accounts for the number of particles which have velocity c1 after the collision. The difference between the two terms gives the net change in the number of particles due to collisions. Note that in the absence of external forces the term containing F will be equal to zero, which is the case of the hard sphere fluid of this study. 3.4 Enskog Theory In the case of a dense gas, the ratio of the particle’s size and the mean free path becomes smaller and the size of the particle can no longer be ignored. Boltzmann’s equation must be modified to account for a particle’s size, and this was first done by Enskog [31]. The Enskog Theory for the hard sphere fluid is based on the following assumptions. (1) as in Boltzmann theory, the motion of the particles is considered to be uncorrelated, (2) the density of the fluid is not necessarily uniform, and (3) unlike Boltzmann theory, the distance between the centres of the particles at their closest approach is not ignored. The position of each particle in the colliding pair is represented by r1 and r2 , and the distance between their centres is |r2 − r1 |. This difference can be neglected at low densities, but when the density increases, the effect of the difference in the positions of the particle becomes greater. In Enskog theory, the collision term in Eq. (3.5) is generalized as ∂e f E ∂t = σ12 k f c1 , r1 , t f c2 , r2 , t 2 σ12 − X r2 − k f (c1 , r1 , t) f (c2 , r2 , t) gσ (g, χ) sin χdχd dc1 , 2 X r1 + (3.6) where σ12 is the distance of the closest approach of the colliding particles, and the function X accounts for the non-uniformity of the local density and should be evaluated at the point of contact of the two particles. Albeit similar notation, the function X is not the same as 20 the scattering angle χ. Both r1 and r2 equally give positions of particles “1” and “2”. The right-hand side of Eq. (3.6) can be approximated by expanding functions X, f and f in a Taylor series, in which the third and higher order terms are typically ignored, i.e. ∂e f E = J 0 (f f ) + J I (f f ) + J II (f f ) . ∂t (3.7) Collision terms J 0 (f f ), J I (f f ) and J II (f f ) are given by the following expressions J 0 (f f ) = X 2 f2 f1 − f2 f1 σ12 g · k dkdc1 , J I (f f ) = σ12 X f1 ∂f2 1 ∂χ ∂f2 − f1 + f f + f2 f1 ∂r1 ∂r1 2 ∂r1 2 1 2 g · k dkdc1 , ×k1 σ12 J II (f f ) = 2 σ12 2 + X f1 ∂ 2 f2 ∂ 2 f2 + f1 ∂r1 ∂r2 ∂r1 ∂r2 1 ∂2X f f − f2 f1 4 ∂r1 ∂r2 2 1 2 k1 k2 σ12 g · k dkdc1 , (3.8) where dk = sinχdχd . In Eqs. (3.7) the values X, f1 , f2 , f1 and f2 are evaluated near corresponding r leading to a series of equations of increasing order to solve. The lowest term J 0 (f f ) is the Boltzmann collision term with a factor X multiplied by it. More detailed derivation can be found in the work of Chapman and Cowling [31]. 3.5 Transport Coefficients For systems near equilibrium, the solution of the Boltzmann equation can be obtained via the Chapman-Enskog method and appears in varying degrees several times in the literature [31, 33, 35, 49]. Since the full derivation is too long to be presented here, only the major assumptions and final results will be given in this chapter as an example. The idea of the method is the following. We start by assuming that the system is close to equilibrium. The distribution function can be expanded in a series about the equilibrium distribution, f (0) , which to lowest order is f = f (0) (1 + Φ) , (3.9) where f (0) is given by the local Maxwellian distribution and Φ is a small perturbation. The form of Φ must be, in general, Φ=− 1 n 2kB T m 1/2 A1 (C)C · ∇ ln T − 2 B1 (C)C C : ∇c0 . n (3.10) 21 in order to match the functional form rising from the left hand side of of the Boltzmann equation of Eq. (3.4). In this expression, the peculiar velocity is C = c − c0 , where c0 is the velocity of the fluid, n is the number density, C C is a second rank tensor, A = A1 (C)C is a vector function of C, and B = B1 (C)C C is a tensor. This expression of Φ is written for the spherical case. In the rough hard sphere case Eq. (3.10) contains more terms that depend upon angular velocity and translational-rotational coupling. The heat flux q in the system is given by 1 f 0 Φ · mC 2 Cdc, 2 q = −λ∇T = (3.11) where λ is the thermal conductivity. Substituting Φ from Eq. (3.10) into Eq. (3.11) gives an expression for the thermal conductivity λ= 3 kB T A, A . 2 m (3.12) Similarly, for shear viscosity we obtain 2 η = kB T [B, B] . 5 (3.13) The notation [B, B] denotes a bracket integral which involves the collision operator in the Boltzmann equation[31]. The bracket integrals in Eqs. (3.12) and (3.13) are evaluated by expanding A and (n) B in terms of a basis of Sonine polynomials Sm (x), where x = C 2 , that differ from the generalized Laguerre polynomials only by a normalization constant. They will be reduced to single terms in that expansion and these terms are evaluated. At the lowest order of approximation, Pidduck’s formulae for the transport coefficients are obtained, namely D = 3 8nσ 2 kB T πm 1+α , 1 + 2α (3.14) λ = 75 64σ 2 3T kB πm 12(1 + α)2 (37 + 151α + 50α2 ) , 25(12 + 75α + 101α2 + 102α3 ) (3.15) η = 5 16σ 2 mkB T π 2(1 + α)2 , 6 + 13α (3.16) ηb = 1 32σ 2 mkB T π (1 + α)2 , α (3.17) in which the quantities within square brackets for the self-diffusion, thermal conductivity 22 and shear viscosity are the corresponding expressions for the smooth sphere fluid. The latter has no bulk viscosity in the low density limit. Another expression for the diffusion coefficient at low densities for binary mixture of pure rough sphere fluids is given by [11] D12 = 3 8nR2 kB T 2πµ 1/2 1 + µχ , 2 + µχ (3.18) in which R = (σ +Σ)/2 and the quantity in square brackets is the binary diffusion coefficient for a bath of smooth spheres. The roughness modifies the smooth result by a multiplicative factor but does not change the dependence upon temperature or density. The value of α ranges from zero, when mass is concentrated at the center, to 2/3, when mass is confined to the surface of the sphere. Thus, the roughness decreases the diffusion coefficient relative to the smooth case but this attenuation is not too severe. For a pure rough sphere fluid, µχ = 3/2 when α = 2/3 so the attenuation factor is 5/7 for the largest value of α. For the more general case of tracer diffusion, a similar attenuation is expected. The most complete solution of the Boltzmann equation was given by Condiff, Lu, and Dahler [35] who compared the contributions to transport coefficients from different terms in a basis of Sonine polynomials. When Condiff, Lu, and Dahler examined the accuracy of the Pidduck formulae they found that, especially for diffusion, thermal conductivity and shear viscosity, the predictions could vary by almost 10% from those obtained by including additional basis functions in the expansion of the distribution function. A number of different basis functions can contribute and based upon an examination of this, the authors proposed new formulae which included the effects of the most important ones. Based upon this analysis, the following modified formulae were proposed, −1 ˜ = D 3 8nσ 2 kB T πm 1+α πα(1 + α) 1+ 1 + 2α 2(1 + 2α)(5 + 9α + 8α2 ) ˜ = λ 75 64σ 2 3T kB πm 4(1 + α)3 (1169 + 6667α + 7690α2 + 2000α3 ) ,(3.20) 25(116 + 969α + 2560α2 + 3973α3 + 3626α4 + 1360α5 ) η˜ = 5 16σ 2 mkB T π 2(1 + α)2 (3 + 10α) . (6 + 33α + 35α2 ) , (3.19) (3.21) For higher densities, McCoy, Sandler and Dahler[36] applied Enskog theory to a fluid of rough hard spheres. The equation of state for rough hard spheres is identical to that for smooth hard spheres and is given by p = nkB T (1 + bng) , (3.22) 23 in which g is the value of the equilibrium radial distribution function evaluated at the surface of the sphere (that is, the contact value) and b = 2πσ 3 /3 is the second virial coefficient for hard spheres. At the first-order level of the Chapman-Enskog approximation, the expressions for the shear viscosity, bulk viscosity and thermal conductivity are given by Eqs. (32), (33), and (35) of Ref. [36]. For the shear and bulk viscosities, these expressions retain the number of basis functions consistent with the expressions for η ∗ and ηb given (K) above. By determining explicit expressions for η (K) and ηb in Eqs. (32) and (33) of Ref. [36], and performing some algebraic manipulations, final expressions for the prediction of shear and bulk viscosities within Enskog theory are η 2 (2 + 5α) 2 (6 + 13α)(4 + 7α) (bng)2 1+ (2 + 5α)2 + bng + 2 g 5 (1 + α) 25(1 + α) π (1 + α) ∗ 2 η˜ 2 (2 + 5α) (bng) 6 (6 + 33α + 35α2 ) 2 1+ bng + (2 + 5α) + g 5 (1 + α) 25(1 + α)2 π (3 + 10α) ηb 32 α 1 + 2bρg + (bρg)2 1 + . g π (1 + α)2 ηE = η˜E = ηbE = (3.23) , (3.24) , (3.25) The expressions for thermal conductivity require more detail. As pointed out in Ref. [50], Eq. (27) of Ref. [36] contains an error on the left hand side. To correct this error the authors suggest that Eq. (26b) of Ref. [50] should replace the left hand side of Eq. (27) of Ref. [36] and Eq. (28) of Ref. [50] should replace Eq. (35) of Ref. [36]. Upon following the suggested framework, it appears that the proposed correction may also have some errors associated with it. When the appropriate corrections are made and propagated through the mathematical framework of Ref. [36], the expressions for the thermal conductivity coefficient, within Enskog theory, are found to be bng (136 + 603α + 425α2 + 150α3 ) 2(bng)2 + 2 2 3(1 + α) (37 + 151α + 50α ) 9(1 + α) (37 + 151α + 50α2 ) 24 × 96 + 491α + 459α2 + 240α3 + (12 + 75α + 101α2 + 102α3 ) , (3.26) π λE = λ g ˜E = λ ˜∗ λ g 1+ bng (1962 + 19247α + 40320α2 + 35211α3 + 16700α4 + 4000α5 ) (1 + α)2 (1169 + 6667α + 7690α2 + 2000α3 ) 2 (bng) + 3 (1 + α) (1169 + 6667α + 7690α2 + 2000α3 ) × (1 + α)(757 + 5680α + 10583α2 + 9510α3 + 5850α4 + 2000α5 16 + (116 + 969α + 2560α2 + 3973α3 + 3626α4 + 1360α5 ) . (3.27) π 1+ The details of the derivation are presented in the next section. Note that the resulting expressions depend upon the number of basis functions retained ˜ In principle, the expressions in the distribution function, and hence differ for λ and λ. 24 ˜ should be the more accurate. Finally, we are not aware of any derivations of involving λ expressions for the self-diffusion of rough spheres using Enskog theory. 3.6 Notes on the Derivation of the Transport Coefficients The kinetic theory determination of transport coefficients for the rough sphere fluid has been presented by Condiff, Lu and Dahler [35] in the low density limit, for which the Boltzmann equation was solved, and by McCoy, Sandler, and Dahler [36] in the dense limit, for which Enskog theory was employed. In a subsequent work, in which the effect of external fields are examined, Klein, Hoffman, and Dahler [50] indicate an error in some of the equations of Ref. [36]. Specifically, the left hand side of Eq. (27) of Ref. [36] should be replaced by Eq. (26b) of Ref. [50], and that the table of calculated thermal conductivity coefficients, Table I of Ref. [36] should be replaced by updated values in Table IV of Ref. [50]. Unfortunately, it appears that the equations and table of Ref. [50] may also be in error. In the limit of low density, Eq. (26b) reduces to Eq. (26a) in Ref. [50]. The resulting column matrix should agree with the corresponding one used in the low density theory of Ref. [35], which is the left hand side of Eq. (43) of that work. A comparison shows that the element in the third row of Eq. (26a), written as zero, differs from the low density value of -1. However, the calculated low density values appearing in Table IV of Ref. [50] for bn = 0 agree completely with the low density values appearing in Tabel I of Ref. [36] and Table II of Ref. [35]. In other words, calculated low density values presented in tables are all completely consistent, indicating that the low density equations are correct, and that Eq. (26a) of Ref. [50], being inconsistent with them, must be in error. Furthermore, if a value of -1 is added to the element of the third row of Eq. (26b) to make it correct in the low density limit or if the values in Eq. (26b) are used as written, and if the numerical calculations of Ref. [50] are repeated, the resulting calculated values do not agree with those appearing in Table IV. Given these inconsistencies, a rederivation of the Enskog results was made. The discrepancies were traced to Eqs. (24a)–(24c) of Ref. [36]. The corresponding, rederived equations were found to be, using the notation of Ref. [36], 3 − [5δi1 δk0 + 3δi0 δk1 + 2bngδi0 δk0 ] − n−1 4 m 2kB T 1/2 (i) (k) dτ f [0] HT · WS3/2 (W 2 )S1/2 (Ω2 ) (i) (k) = g[A; {WS3/2 (W 2 )S1/2 (Ω2 )}] , (A1) 25 − p − 2 δi0 δk0 + 5δi1 δk0 + 5δi0 δk1 nkB T 3 2 8 −n −1 m 2kB T 1/2 (i) (k) dτ f [0] HT · ΩΩ · WS3/2 (W 2 )S3/2 (Ω2 ) (i) (A2) (k) = g[A; {ΩΩ · WS3/2 (W 2 )S3/2 (Ω2 )}] , −n−1 m 2kT 1/2 (i) (k) dτ f [0] HT · (Ω × W)S3/2 (W 2 )S3/2 (Ω2 ) (i) (k) = g[A; {(Ω × W)S3/2 (W 2 )S3/2 (Ω2 )}] , (A3) (p) in which HT is given by Eq. (20) of Ref. [36], and Sq (x) are Sonine polynomials whose normalization is given by Eq. (34) of Ref. [35]. Substituting the expression for HT into Eqs. (A1)–(A3) and performing the resulting integrals (this requires use of the relations between pre and post collision values of velocity and angular momentum given by Eq. (A1) of Ref. [35]) gives − 3 4 5+ = − 3 4 (5α + 3) 2 bng δi1 δk0 + 3 + bng δi0 δk1 (α + 1) (α + 1) (i) (k) g[A; {WS3/2 (W 2 )S1/2 (Ω2 )}] 5 (5α + 3) + bng δi1 δk0 + 2 2(α + 1) (A4) , 5 1 + bng δi0 δk1 − δi0 δk0 2 (α + 1) (A5) (i) (k) WS3/2 (W 2 )S3/2 (Ω2 )}] = g[A; {ΩΩ · , √ 3 8 α (i) (k) √ bngδk0 Ii = g[A; {(Ω × W)S3/2 (W 2 )S3/2 (Ω2 )}] , 4 2π (κ + 1) (A6) in which 1 Ii = √ 8 2π 5/2 dW dW1 e−W 2 +W 2 1 (i) S3/2 (W 2 )[2GW 2 + G−1 (W 2 W12 − (W · W1 )2 )] , (A7) and G = W1 − W. Because the factors of G mix the W and W1 dependencies in a nonseparable manner, it is not possible to determine a simple analytic formula for Ii . However, in the present case, only the value for i = 0 is required, and for this case the integral is determined as I0 = 1. When considering the particular terms retained in the Sonine basis set expansion of the distribution function, Eqs. (A4)–(A6) predict that Eq. (26b) of Ref. [50] and the left hand 26 side of Eq. (27) of Ref. [36] should be 3 − 4g 5 + [(5α + 3)/(α + 1)]bng 3 + [2/(α + 1)]bng . −1 −[8/(2π)1/2 ][α1/2 /(α + 1)]bng (A8) After making this replacement, and retaining only the contributions from the terms A10 1 , 00 A01 1 and A2 , the procedure outlined in Ref. [36] was followed, using the modified equations of Eq. (28) of Ref. [50] where appropriate, to obtain the Enskog formula given by Eq. (3.27) above. Repeating the same procedure but retaining only the contributions from the terms 01 A10 1 and A1 produces Eq. (3.26) above. 27 Chapter 4 Methodology 4.1 Computational Approach In order to study physical properties of the rough hard sphere fluid a computer program has been developed that performs molecular dynamics simulations and collects the necessary data. The initial C++ code was created by R. Sokolovskii for measuring the diffusion coefficient in the smooth hard sphere fluid. The code was then modified and expanded for the ongoing study of the rough hard sphere model. First, the collision subroutine was edited in accordance with the theoretical model described in Chapter 2. That is, the angular velocity was added to the description of the particle, representing a third degree of freedom, and conservation law equations were modified accordingly. Second, calculations for various properties of the fluid were added, namely, calculation of the shear and bulk viscosity, thermal conductivity, pressure, and corresponding data acquisition and analysis routines were added. In the computational representation the fluid consists of N spherical particles, in a cubic cell of edge length L, and repeated multiple times, i.e. periodic boundary conditions are employed. For tracer calculations, N − 1 particles are identical and represent the bath, and the remaining one particle is the tracer. The bath particles have diameter σ and the tracer diameter Σ can be varied relative to σ. Another variable parameter is responsible for rotational-translational coupling and represented by α, the dimensionless radius of gyration, introduced in Chapter 2. The values of N , α, Σ/σ and the reduced density ρ∗ = nσ 3 and the value of the radius of gyration α are supplied as input parameters at the start of simulations. The tracer particle is placed into the simulation cell randomly, and bath particles are added thereafter. The bath particles are distributed in random configuration with no overlaps. At higher densities random configurations produce many overlaps. To overcome this issue the particles are positioned in a face-centred cubic lattice configuration. All particles are assigned random initial velocities. The system is equilibrated for a sufficient amount of time before statistics are gathered. The mass m of the tracer was chosen to be identical to the mass of the bath particle for consistent comparison with the studies for the smooth sphere case done by R. Sokolowskii [23], for which this condition was used. We expect the mass dependence to be weak for 28 the following reasons. For large tracers, diffusion does not depend upon mass, as shown in Eq. (3.2). For small tracers, the dependence upon mass appears in the reduced mass factor which changes at most by the factor of 2, as shown in Eq. (3.14). The “event-driven” mechanism described by Allen and Tildesley [51] takes place in the simulation. In this method, all collisions are modelled as binary, instantaneous events. Particles do not interact between collisions and no forces are calculated. At each collision, new velocities and angular velocities are calculated according to Eqs. (2.13) - (2.16), and the next collision is predicted based on updated parameters, considering that particles move linearly and with constant speed. When calculating transport properties the system is sampled at equal periods of time, generally chosen to be approximately 10 times the average collision time. Since the sampling times do not necessarily coincide with collision times, Newton’s Laws are used to correct the positions of particles to account for the difference between the sampling time and the last collision time. Thus the simulation proceeds from the previous collision to the next one predicted, until it reaches the required number of collisions, set at input. One of the computational challenges is dealing with so called “long tracks” particles that travel for a long time before they collide. When the distance a particle travels exceeds the dimensions of the simulation cell, it leaves the cell and it’s periodic image re-enters on the opposite side of the cell. In this case the next point of collision is recalculated using this image of the particle. Thus, in the case of long tracks, the point of the next collision is calculated twice, which adds to computer simulation time when the density is very low and collisions are especially rare. The parameters used to construct the simulation system are scaled. The basic scaled values inside the program are length, mass and energy. The size of the particle is σ ∗ = σ/L, where L = 1 is the size of the simulation cell. The size of the bath is chosen relative to the size of the simulation cell while the size of the tracer Σ is chosen relative the size of the bath. The mass is scaled as m∗ = m/mbath so that m∗ = 1 for the bath; for the tracer we chose m∗ = 1 as well. The kinetic energy E ∗ = E/kT = 1, meaning no thermostat is employed in the simulation. The kinetic energy and temperature are monitored during simulation. In order to obtain enough data for reliable statistical analysis, a large number of computational jobs for a particular set of parameters is submitted at the same time. All data collected during simulation is averaged over the number of samples collected and the number of job instances. A set of utility programs have been written to submit and monitor the jobs, and for further data averaging and analysis. The typical number of jobs submitted for a particular combination of N , ρ∗ , Σ/σ and α varies from 30 to 200. To insure the system is well converged, the mean square displacements for the tracer are plotted as a function of time. Usually, at short times the mean square displacement is 29 well converged, but is far from the limit required to extract reliable values for the transport coefficients. At longer times, the mean square displacement is not as well converged, due to poor statistics involved with sampling at large time differences. Fitting the raw data to a rational function of the type f (t) = (6Dt + a)t/(t + b) allows for the reliable estimation of the transport coefficients. This function is chosen in a way to ensure that it is equal to 0 at t = 0 and to 6D at infinitely long time. Specific details and examples of how the data is processed will be given in the following chapters. 4.2 Transport Coefficients The transport coefficients are calculated using Green-Kubo expressions which can be expressed in two equivalent ways. For any dynamic variable G the transport coefficient of interest can be obtained using the following relation [52] ∞ Φ(G) = G˙ (t) G˙ (t + s) ds. (4.1) 0 Alternatively, it can be obtained by taking a long-time limit in the expression [52] 1 d < [G(t + s) − G(t)]2 > . s→∞ 2 ds Φ(G) = lim (4.2) In both cases, simulations must be run for long times in order to reach convergence. Performing the integration in Eq. (4.1) becomes particularly challenging for diffusion coefficients when the velocity autocorrelation function is used. The behaviour of the velocity autocorrelation function at long times was explored by R. Sokolowskii in relation to estimating the diffusion coefficient in the smooth hard sphere model [23]. He found it to be difficult to obtain accurate estimates of the diffusion coefficient because the velocity autocorrelation function takes a long time to converge and even at very long times it still oscillates around the limit, producing large statistical errors. Hence Eq. (4.2) was used to obtain the transport coefficients. Transport coefficients can be determined from Eq. (4.2) in the following way. For diffusion, D, thermal conductivity, λ, shear viscosity, η, and bulk viscosity, ηb , the specific 30 expressions for G are, respectively, GD = x i , (4.3) N Gλ = x˙i [Ei − < E >] , (4.4) mi x˙i yi , (4.5) mi x˙i xi − pV t , (4.6) i=1 N Gη = i=1 N Gηb = i=1 in which the sums include all N particles in the system, xi and x˙i are position and velocity in the respective Cartesian coordinate of particle i, Ei and < E > are the total energy of particle i and average total energy of the system, and p and V are the pressure and volume, respectively. Because there are three degrees of freedom each in translation and rotation, the average energy in the thermodynamic limit in this particular case is < E >= 3kB T . Note that for each expression above there exist two others, analogous in form but with the other Cartesian coordinates. During the numerical calculations, all three expressions are used for each coefficient in order to improve the statistical averaging. Using the definitions above, the expressions for the transport coefficients become D = Φ(GD ) , 1 Φ(Gλ ) , λ = V kB T 2 1 η = Φ(Gη ) , V kB T 1 4 Φ(Gηb ) − η . ηb = V kB T 3 (4.7) (4.8) (4.9) (4.10) The dynamics of the system in this study is event-driven rather than continuous thus making it impossible to calculate G directly. The acceleration term x ¨ behaves like a delta function and G becomes discontinuous at the points of collision. To overcome this problem the method of Alder, Gass, and Wainwright is used [52]. The change in G is calculated using t+s G(t + s) − G(t) = ˙ )dτ . G(τ (4.11) t This integral is divided into two contributions: i) one arising from times between collisions, and ii) one arising from collision events. Between collisions particle velocities are constant, while they change in a discontinuous but predictable way at collisions. During the simulation, values for both contributions are collected and then combined to produce the final 31 statistically averaged quantities. To illustrate this process consider the viscosity coefficient for a rough sphere fluid. It can be calculated from the off-diagonal elements of the pressure tensor, that is 1 d < [Gη (t + s) − Gη (t)]2 > s→∞ 2 ds ∞ < G˙ η (t)G˙ η (t + s) > ds , = (V kT )−1 η = (V kT )−1 lim (4.12) 0 where differentiating Eq. (4.5) gives N G˙ η (t) = [mi x˙i y˙i + mi x¨i yi ] . (4.13) i=1 Imagine p collisions occur during the time period from t to t + s at the times tk , k = 1, ..., p. The integral in Eq. (4.11) is divided into smaller intervals between whose end points are the times tk . In this way, one knows that collisions occur only at the upper and lower limits of each interval but not within. For ease of notation, define t0 = t and tp+1 = t + s. Consider now the terms on the right hand side of Eq. (4.13). The first term involves the velocity components of the particles. These velocities are constant between collisions so their contributions to Eq. (4.11) are given by p N N tk+1 p mi (vxi )k (vyi )k (tk+1 − tk ), mi x˙i y˙i ds = k=0 i=1 tk (4.14) i=1 k=0 in which (vxi )k and (vyi )k are the x and y components of the velocity of particle i between the times tk and tk+1 . Because no collisions occur between tk and tk+1 , the velocities of all particles are constants, and their contributions can be factored out of the integrals on the left hand side rendering the integral trivially. The second term on the right hand side of Eq. (4.13) involves components of the acceleration. These components are zero for all times except at the points of collisions so the contributions of these terms to Eq. (4.11) are given by p N tk +δ mi x¨i yi ds , k=1 i=1 (4.15) tk −δ in which δ represents an infinitesimally small time. However, at each collision time tk , only two of N particles are colliding. Let them be denoted with the labels ik and jk . Only these pairs of particles have non-zero accelerations at tk and hence only they will contribute to Eq. (4.15). This collapses the sum over all N particles to a sum over only the colliding 32 pairs, that is Eq. (4.15) becomes p tk +δ tk +δ tk −δ k=1 ¨ik yik ds + mik x p ¨jk yjk ds mjk x tk +δ tk +δ m ik yik = tk −δ k=1 p tk −δ x ¨ik ds + mjk yjk tk −δ x ¨jk ds mik yik (vxik )k − (vxik )k−1 + mjk yjk (vxjk )k − (vxjk )k−1 = , (4.16) k=1 in which (vxik )k − (vxik )k−1 represents the difference in the x-component of the velocity of the particle with index ik after and before its collision at tk , for example. These velocity differences are directly related to the collision impulses through Eqs. (2.13)-(2.16). Associating ik and jk with particle labels ‘1’ and ‘2’ then gives mik vik − mik vik = −Jk , (4.17) mjk vjk − mjk vjk = Jk , (4.18) in which Jk is the impulse associated with the collision at time tk . Using these equations then gives p k=1 tk +δ tk −δ p tk +δ mik x ¨ik yik ds + tk −δ mjk x ¨jk yjk ds = ∆yk Jxk , (4.19) k=1 in which ∆yk = yik − yjk is the relative difference in the y components of the positions of particles with labels ik and jk colliding at the time tk . Combining the results from Eqs. (4.19) and (4.14) then gives N p G(t + s) − G(t) = p mi (vxi )k (vyi )k (tk+1 − tk ) + i=1 k=0 ∆yk Jxk , (4.20) k=1 and this is the working equation used, along with Eq. (4.12), in the simulation code. Note that equations analogous to Eq. (4.20) also exist for the other two component pairs (yz and zx) and in practice all three pairs are calculated and averaged together to improve statistical convergence. Note that Eq. (4.20) has the advantage of depending only upon relative distances so it is independent of the particular choice of coordinate system. In practice, the contributions towards G are stored with each particle during the simulation. When a collision occurs, the contributions from the terms in Eq. (4.20) are calculated and added to these stored variables so that at any time, the contributions from the full history of collision events for a given particle are accounted for. The calculation of bulk viscosity presents the additional challenge of removing the 33 quadratic dependence of time that manifests itself when Φ(Gηb ) is calculated without subtracting the pV s contribution, i.e. < N i=1 mi [x˙i (t + s)xi (t + s) − x˙i (t)xi (t)] >= pV s [52]. While in principle this quadratic dependence could be removed with the use of an appropriate fitting function, in practice the quadratic growth swamps the values, especially at large times, and obscures the underlying linear behaviour upon which the transport coefficient depends. The simplest way to deal with this quadratic dependence is to remove the factor pV s during the calculation of Φ(Gηb ). This requires that pV be known at the start of the simulation. In principle one could obtain this value on-the-fly as the simulation proceeds, or by some other estimation method. In the present case, the Carnahan-Starling formula [53] is used to estimate pV under the conditions of each simulation since it is a reasonably accurate equation of state for hard sphere fluids (both smooth and rough). However, because of small inaccuracies in this formula and because of finite size effects in the prediction of the equation of state, the entire contribution of the pV term is not removed by this formula. A quadratic time dependence in Φ(Gηb ) still remains although smaller in magnitude. This remaining dependence is removed in a fitting procedure by using a function that explicitly incorporates a quadratic dependence. 4.3 Finite Size Effects When a tracer is introduced into a bath, the velocity flow around the tracer leads to correlations in the motion of the particles that propagate throughout the fluid. If the simulation box is not large enough, some of the correlations will be cutoff when applying periodic boundary conditions, leading to significant error when estimating the diffusion coefficient [23]. To compensate for this effect, the values of Φ(G) are extrapolated to the infinite size limit. This is done using linear regression on the set of Φ(G) for simulations with different numbers of particles. In the case of the diffusion coefficient, the Fushiki method is used [25] where the diffusion coefficient is plotted as a function of 1/L and extrapolated to 1/L → 0. For large tracers, it is only with such extrapolations that accurate diffusion coefficients can be obtained. In the case of shear and bulk viscosity and thermal conductivity extrapolation is done over the entire number of particles, i.e. 1/N . Examples of these extrapolations will be shown in Chapters 5 and 6. 34 Chapter 5 Tracer Diffusion in the Rough Hard Sphere Fluid 5.1 Objectives and Background In this chapter a study of tracer diffusion in the rough hard sphere fluid is presented. Rough sphere particles possess rotational degrees of freedom which influence the behaviour of the fluid in comparison with that of a smooth sphere fluid where rotation of the individual particles is not taken into account. This additional degree of freedom affects the fluid on a microscopic level by changing the physics of a collision, which in turn influences the macroscopic properties of the fluid, such as diffusion. Since the structures of the smooth hard sphere and the rough hard sphere fluids are precisely the same, and particles in both types of fluid have precisely defined radii, it is possible to isolate and examine the sole effect of rotation of individual particles on the behaviour of the fluid. The work is focused on the following objectives. First, the dependence of the diffusion coefficient upon particle size will be studied across different fluid regimes, from molecular to hydrodynamic, for different densities. The results will be compared with theoretically predicted limits. Second, the boundary conditions adopted by the tracer surface in the hydrodynamic regime will be determined. Third, the dependence of the diffusion coefficient upon the degree of translational-rotational coupling will be examined and various dependencies observed. The detailed explanation of the results will follow. To accomplish this study, extensive molecular dynamics simulations were performed using the rough hard sphere model described in Chapters 2 and 3. In the current chapter, the specific computational details of the simulation will be provided, followed by the data analysis. Results and discussion follow thereafter. 5.2 Computational Details Event driven molecular dynamics simulations were performed consisting of N rough hard sphere particles enclosed into a cubic cell. The diameter of the tracer, Σ, was varied relative to that of the bath, σ. Simulations were done for Σ/σ varying from 0.125 to 16 and for N 35 varying from 128 to 55 297. The ratio Σ/σ is chosen in such a way as to maintain integer values of the input parameter n in a relation Σ = 2n/3 σ. Values of n are −9, −3, 0, 3, 5, 10, 12. The mass of the tracer was set equal to the mass of the bath particles. The diameter dependence of tracer diffusion was explored for reduced densities ρ∗ = nσ 3 of ρ∗ = 0.34 and 0.52 since the transition to hydrodynamic behaviour is expected to be most pronounced in this range. For each set of parameters, 30 to 200 simulation runs were carried out to increase statistical averaging. Each system was equilibrated for 400,000 collisions and production runs used 4 × 108 collisions. Diffusion coefficients are reported as scaled values, D∗ = Dt0 /σ 2 , where the time unit t0 is defined as t0 = σ are also scaled as η∗ m/kB T . Viscosities = ησt0 /m where m is the mass of the bath gas. The event-driven algorithm employed is described in Chapter 4. The total number of of jobs ran during this work is approximately 120, 000. While jobs with the smaller number of particles and small tracers take hours of CPU time, the amount of CPU time goes up significantly with increasing N and tracer size. The longest CPU time was approximately 16 weeks for the case of Σ/σ = 16 and N = 55297. The data occupy approximately 800GB of disk space. 5.3 5.3.1 Data Analysis Data Collection and Convergence As each simulation progresses the positions of the particles are collected at equal periods of time, as described in Chapter 4, and these, along with mean squared displacements ∆R2 , are calculated and stored in separate files labelled by corresponding ρ∗ , Σ/σ and N . For example, 100 jobs submitted for ρ∗ = 0.52, Σ/σ = 16, N = 55297 result in 100 separate data files, each containing the data sampled at equal time periods throughout the simulation for the desired number of collisions. For small time differences, millions of samples are collected, but at longer time differences the data is sampled less frequently, and therefore the number of samples reduces to hundreds of thousands, and then to tens of thousands at very long times, for each sampling point. The data is then averaged for each trajectory at each sampling point and stored in a human-readable data file that is used for further analysis. Table 5.1 is an example of such a data file. In this table, the first column represents the time difference at which the mean square displacement is calculated; the second column is the mean square displacement averaged over all number of samples collected for this time difference; and the third column is the number of samples used for the averaging. The values in the second column are then averaged over the total number of jobs submitted for a given set of parameters and the final value for a given time difference is calculated, along with the corresponding error, which is estimated as a standard deviation 36 among the different jobs. For the data to be reliable the simulation system has to be well converged. The convergence is examined by plotting the mean square displacements for the tracer as a function of time. Figures 5.1 and 5.2 are examples of the convergence plots. Figure 5.1 represents the largest tracer (Σ/σ = 16) in a simulation using 55297 particles at the reduced densities of ρ∗ = 0.34 and 0.52. At short times, the mean square displacement is well converged but it is far from the limit required to extract the diffusion coefficient. At longer times, the mean square displacement is not as well converged due to the poorer statistics involved with sampling at large time differences. While directly differentiating the mean square displacement curve allows for estimation of the diffusion coefficient, this procedure would produce poor results due to statistical uncertainties. In order to overcome this problem the mean square displacements were fit to the rational function f (t) = (6Dt + a)t/(t + b). This function was chosen in a way to ensure that it is equal to 0 at t = 0 and to 6D at infinitely long time. This function is represented by the solid lines in Fig. 5.1 where it can be seen to fit the data well. The limit of this function should therefore be reliable for estimating the diffusion coefficient. The dashed lines on the graph are the first derivatives of the fit function and approach the diffusion constant value at long times. Similarly, Fig. 5.2 is a sample plot for convergence of the data at different values of the translational-rotational coupling coefficient α. For each value of α the data is well converged and provides a reliable estimate for the diffusion coefficient. The reason to use the rational function rather than straight line is to obtain the best value at the longest time. In the case of the very good data, both straight line and rational function would produce the same result. However, slight deviations in the data may affect the linearity at long times. The rational function is more flexible and takes into account such deviations, producing a more precise fit with a very small error. 5.3.2 Finite Size Effects The finite size effects, described in Chapter 4 are removed by using Fushiki extrapolation [25]. Diffusion coefficients calculated for different N (and hence simulation cell size) are plotted as a function of 1/L and extrapolated to 1/L → 0 for each tracer size at each density to obtain the final value of the diffusion coefficient D∞ . Figure 5.3 shows some sample results of such an extrapolation for tracers with Σ/σ = 10.1 for reduced densities ρ∗ = 0.34 and 0.52. The same extrapolation is performed for every tracer size. Using extrapolation to extract the diffusion coefficient for the large tracers is crucial, as the finite size effects are even more pronounced when the tracer grows in size. When the size of the tracer increases, the volume left for the rest of the particles decreases which leads to increase in the bath density. The difference between the values of the diffusion coefficient from the 37 Table 5.1: Values of the mean square displacement for various time differences collected over corresponding number of samples. t < ∆R2 (t) > number of samples 0.112495 0.00109231 7772 0.224991 0.0024072 7771 0.337486 0.00380467 7770 0.449982 0.0052554 7769 0.562477 0.0067404 7768 0.674973 0.00822421 7767 0.787468 0.00973216 7766 0.0112602 7765 0.899964 1.01246 0.0128011 7764 1.12495 0.0143493 7763 2.24991 0.0307179 776 3.37486 0.0472742 775 0.0637194 774 4.49982 5.62477 0.0807459 773 6.74973 0.0985555 772 7.87468 0.116026 771 8.99964 0.133932 770 10.1246 0.151978 769 0.170979 768 11.2495 22.4991 0.369434 76 33.7486 0.542966 75 44.9982 0.726084 74 56.2477 0.904204 73 67.4973 1.05503 72 78.7468 1.18746 71 89.9964 1.37241 70 1.57692 69 101.246 112.495 1.7848 68 224.991 3.50937 6 337.486 4.88218 5 449.982 5.38599 4 562.477 6.86265 3 674.973 13.0796 2 787.468 11.4814 1 38 350 ∗ 50 ∗ ρ = 0.34, N = 5120, Σ/σ = 5 ρ = 0.34, N=55297, Σ/σ = 16 300 40 150 20 2 30 <∆ R (t)> 200 2 <∆R (t)> 250 100 10 50 0 ∗ ∗ ρ = 0.52, N = 5120, Σ/σ = 5 75 ρ = 0.52, N=55197, Σ/σ = 16 0 12 10 8 2 45 <∆ R (t)> 2 <∆R (t)> 60 6 30 4 15 2 0 0 0 500 1000 t/τf 1500 0 500 1000 t/τf 1500 2000 Figure 5.1: The time dependence of the tracer mean square displacement from a “fixed I” (i.e. α = 0) simulation with N = 55297, Σ/σ = 16 and N = 5120, Σ/σ = 5, at the reduced densities ρ∗ = 0.34 (upper panel) and ρ∗ = 0.52 (lower panel), respectively. The bars indicate the estimated errors in the simulation data points while the solid lines are the rational function fits to the data. The dashed lines are the first derivatives of the fit function. The scale for derivatives is not shown. 39 800 2 <∆R (t)> 600 400 200 α = 0, Σ/σ = 5 0 2 <∆R (t)> 60 40 20 α = 2/5, Σ/σ = 5 400 0 2 <∆R (t)> 300 200 100 α = 2/3, Σ/σ = 5 0 0 4000 8000 12000 16000 20000 t/τf Figure 5.2: The time dependence of the tracer mean square displacement from simulations with α varying from 0 to 2/3, with Σ/σ = 5, at the reduced density ρ∗ = 0.52. The bars indicate the estimated errors in the simulation data points while the solid lines are the rational function fits to the data. The dashed lines are the first derivatives of the fit function and approach the diffusion constant value at long times. The scale for derivatives is not shown. 40 largest simulation cell size and extrapolated values for large tracers can be up to 30%, according to Table 5.2. Iin the work of R. Sokolowskii [23] a comparison was made of the diffusion coefficients obtained by using Fushiki extrapolation against the data obtained by simple volume correction. He showed that for increasing tracer sizes the difference between the diffusion coefficients is significant. The volume correction works fairly well for static changes in fluid density, but it doesn’t take into account the dynamic finite size effects, and therefore the use of Fushiki extrapolation is necessaryl for obtaining the accurate results for large tracers. In order to obtain the values of shear viscosity, separate calculations were run and a new set of data was produced. A similar extrapolation procedure was used on this new set of data to extract the shear viscosity coefficient. The reduced viscosity values obtained were 0.445 and 0.860 for ρ∗ = 0.34 and 0.52, respectively. 5.3.3 Rotational-Translational Coupling Compared with the smooth hard sphere case, the rough hard sphere model contains one additional degree of freedom, rotational motion, whose effect needs to be examined. Several different approaches can be taken. Using the definition of α, and recalling that the masses of the tracer and bath were set equal in all calculations, gives Itracer αtracer = Ibath αbath Σ σ 2 . (5.1) Two different classes of calculations were performed which were labelled “fixed I” and “fixed α”. In the “fixed I” class, all calculations enforce the moment of inertia of the tracer to always be equal to the moment of inertia of the bath, i.e. Itracer = Ibath , regardless of the tracer size. In the rotational degree of freedom, this would be analogous to setting the masses of the tracer and bath equal in the translational degree of freedom. In other words, the tracer rotates with the same average angular speed regardless of its size. In addition, αbath was set to 2/5, which corresponds to uniform mass distribution within the bath spheres. With these constraints, Eq. (5.1) shows that as the tracer grows in size, Σ/σ increases and αtracer decreases. Thus, in the “fixed I” class, αtracer gets closer to zero as the tracer grows in size. Conversely, for tracer sizes smaller than the bath, αtracer must grow. However, the value of α has an upper limit of 2/3, corresponding to the mass in the sphere being concentrated solely at its surface. According to Eq. (5.1), once Σ/σ < 3/5 it would be necessary for αtracer to exceed 2/3 in order to maintain Itracer = Ibath so this tracer size represents the limit on the lower end. In order to move beyond this, the value of αbath is changed instead when Σ < σ. More specifically, αtracer is fixed at 2/5 and αbath is decreased 41 Table 5.2: For the “fixed I” class, reduced diffusion constants of the bath and tracer for differing tracer sizes Σ/σ and total number of particles N for two reduced densities ρ∗ . Values for N = ∞ are obtained by extrapolating to the infinite volume limit, as described in the text. ρ∗ = 0.34 ρ∗ = 0.34 ρ∗ = 0.52 ρ∗ = 0.52 3 ∗ ∗ ∗ ∗ N tracerDN bathDN tracerDN bathDN (Σ/σ) 1/512 128 1.57 0.4256 0.885 0.2125 1/512 256 1.56 0.4348 0.888 0.2186 1/512 512 1.58 0.4420 0.894 0.2239 1/512 1024 1.60 0.4479 0.860 0.2280 1/512 5120 1.58 0.4571 0.907 0.2349 1/512 ∞ 1.62 0.4700 0.910 0.2439 1/8 128 0.690 0.3797 0.357 0.1857 1/8 256 0.701 0.3891 0.362 0.1918 512 0.707 0.3964 0.368 0.1968 1/8 1/8 1024 0.703 0.4022 0.373 0.2006 1/8 5120 0.714 0.4114 0.366 0.2072 1/8 ∞ 0.738 0.4246 0.390 0.2154 1 128 0.305 0.3053 0.144 0.1432 1 256 0.315 0.3146 0.163 0.1489 1 512 0.322 0.3216 0.153 0.1532 1 1024 0.325 0.3269 0.156 0.1567 1 5120 0.324 0.3351 0.158 0.1621 1 ∞ 0.348 0.3481 0.171 0.1699 32 128 0.0714 0.2701 0.0254 0.1105 512 0.0892 0.3123 0.0377 0.1443 32 32 1024 0.0950 0.3222 0.0416 0.1521 32 5120 0.107 0.3342 0.0480 0.1611 32 ∞ 0.119 N/A 0.0535 N/A 128 512 0.0449 0.2883 0.0173 0.1214 128 1024 0.0502 0.3099 0.0209 0.1401 128 5120 0.0598 0.3317 0.0264 0.1587 128 13824 0.0623 0.3370 0.0294 0.1629 128 ∞ 0.0710 N/A 0.0345 N/A 1024 1024 0.0130 0.2128 N/A N/A 1024 13824 0.0261 0.3291 0.0118 0.1551 1024 32768 0.0279 0.3365 N/A 0.1616 1024 55297 0.0284 0.3389 0.0145 0.1638 1024 ∞ 0.0355 0 N/A 0.0191 N/A 4096 13824 0.0126 0.3030 0.00510 0.1287 4096 32768 0.0151 0.3256 0.00690 0.1511 55297 0.0146 0.3324 0.00755 0.1574 4096 4096 ∞ 0.0225 N/A 0.0109 N/A 42 0.03 0.025 D * N 0.02 0.015 0.01 ρ = 0.34 0.005 0 0.014 0.012 D * N 0.01 0.008 0.006 0.004 ρ = 0.52 0.002 0 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 1/L ∗ as a function of inverse simulation Figure 5.3: Plots of reduced diffusion constants DN box length for a “fixed α = 2/5” simulation for Σ/σ = 10.1 and for the reduced densities ρ∗ = 0.34 (upper panel) and ρ∗ = 0.52 (lower panel). The lines are the linear function fits to the data, whose y intercept gives the extrapolated value of the diffusion constant. 43 in such a way as to keep Itracer = Ibath . Thus, as the tracer gets smaller, the value of αbath approaches zero. Table 5.3 lists the values of α for both the tracer and bath as a function of tracer size for the “fixed I” class of calculations. In the “fixed α” class, all calculations enforce αtracer = αbath . In this case, Eq. (5.1) shows that the ratio Itracer /Ibath = (Σ/σ)2 so that the moment of inertia of the tracer (relative to the bath) varies with its size. In other words, the angular speed of the tracer decreases as its size increases, and vice versa. Calculations were performed in this class for α = 2/5 and 2/3. Table 5.3 lists the values of Itracer /Ibath for different tracer sizes for the “fixed α” class of calculations. Table 5.3: Values of αtracer and αbath for the “fixed I”, and values of Itracer /Ibath for the “fixed α” calculations, as a function of tracer size Σ/σ. fixed I fixed I fixed α 3 Σ/σ αtracer αbath Itracer /Ibath (Σ/σ) 1/512 1/8 0.40 0.0063 0.016 1/8 1/2 0.40 0.10 0.25 1 1 0.40 0.40 1.00 32 3.17 0.040 0.40 4.00 128 5.04 0.016 0.40 10.1 1024 10.1 0.0040 0.40 102 16 0.0016 0.40 256 4096 5.4 Results and Discussion For the “fixed I” class, calculations were performed at reduced densities of ρ∗ = 0.34 and 0.52 for a range of tracer sizes. Reduced diffusion constants for both the tracer and bath using different numbers of particles are shown in Table 5.2. Also included are tracer values extrapolated to the infinite volume limit. These data are also plotted in Fig. 5.4 where the trend against tracer size can be compared with the hydrodynamic “slip” and “stick” values. The hydrodynamic lines are estimated according to Stokes-Einstein equation mentioned in the introduction with the shear viscosity calculated according to Eq. (4.12) using the procedure described in Chapter 4. A number of observations can be made. Consider first the values of the bath diffusion constants in Table 5.2. The pure bath (Σ = σ) reduced diffusion constants are 0.348 and 0.170 for ρ∗ = 0.34 and 0.52, respectively. For fixed N , as the tracer becomes larger, the bath diffusion constants begin to deviate ever more from these values. This deviation occurs because the tracer removes a significant amount of free volume from the cell. The bath particles exist in an effective smaller volume, leading to a higher effective density and hence 44 2 * 2 D R /! 2 * 2 D R /! 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 # " = 0.34 # " = 0.52 1 2 3 4 5 6 7 8 9 R/! ∗ , multiplied by (R/σ)2 Figure 5.4: Extrapolated reduced tracer diffusion constants, D∞ plotted as a function of tracer size, R/σ (R = (Σ + σ)/2) for the reduced densities ρ∗ = 0.34 (upper panel) and ρ∗ = 0.52 (lower panel). The solid and dashed lines denote the hydrodynamic limit with “slip” and “stick” boundary conditions, respectively. The circles, squares, and triangles denote extrapolated results from the “fixed I”, “fixed α = 2/5”, and “fixed α = 2/3” calculations, respectively. The solid curved lines represent the predictions of Eq. (5.2) for the “fixed α” calculations. 45 a lower diffusion constant. For the largest tracer, this deviation is quite significant even when 55297 particles are in the simulation cell. This is just another manifestation of finite size effects and shows the importance of their removal by extrapolation. Note that attempting to correct for the effective densities of the bath particles by scaling the calculated diffusion coefficients destroys the linear behaviour necessary for the extrapolation as a function of ‘1/L’. As the tracer becomes smaller than the bath, the bath diffusion constants increase as a function of N , as expected from the reasoning given above, but their values begin to exceed the pure bath values. This is a manifestation of the changing value for αbath , as seen in Table II, as the tracer size decreases. In this limit, αbath moves closer to zero, that is as discussed above, the bath approaches the smooth limit. The pure bath gas reduced diffusion constants for the smooth hard sphere are 0.4771 and 0.2486 for ρ∗ = 0.34 and 0.52, respectively. The bath diffusion constants begin to approach these smooth sphere limits when the tracer becomes much smaller in size than the bath. Now consider the tracer diffusion constants. Again, for the “fixed I” calculations, the values of αtracer vary with tracer size, and as seen in Table 5.3 approach close to zero for the largest tracers. As seen in Fig. 5.4, as a function of R, these diffusion constants begin to approach the behaviour predicted by hydrodynamic theory with “slip” boundary conditions. This is not unexpected, as discussed in Section 5.3.3. In the present case the αtracer = 0 limit corresponds to a smooth hard sphere moving in a rough hard sphere bath. It has already been shown that a smooth hard sphere moving in a smooth hard sphere bath adopts “slip” boundary conditions in the large tracer size limit [23]. In this limit, the effect of the bath should be solely incorporated in the viscosity coefficient, so the large tracer size behaviour should depend only upon the bath-tracer collision dynamics. As such, it is expected that the αtracer = 0 limiting behaviour in this case should also approach the “slip” hydrodynamic limit. The curves in Fig. 5.4 show this consistency within the rough sphere model by correctly predicting the smooth sphere limit. The decrease in the diffusion coefficient may be attributed to faster decay of the velocity autocorrelation function that occurs with increase of α, that in turn destroys correlations in the fluid. For the “fixed α” class, calculations were performed at reduced densities of ρ∗ = 0.34 and 0.52 for a range of tracer sizes and for α = 2/5 and 2/3. Reduced diffusion constants for the tracer using different numbers of particles are shown in Tables 5.4 and 5.5. Also included are tracer values extrapolated to the infinite volume limit. Note that the bath diffusion constants are not included here because in the “fixed α” calculations, the values of αbath are the same regardless of the tracer size (unlike the “fixed I” calculations) so the bath diffusion constants show the same trends in all cases, being smaller than the pure bath values due to excluded volume effects and increasing towards the pure bath values as N increases. The tracer data are plotted in Fig. 5.4 where the trend against tracer size can be compared with the hydrodynamic “slip” and “stick” values. 46 Table 5.4: For the “fixed α” class with α = 2/5, reduced diffusion constants of the tracer for different tracer sizes Σ/σ and total number of particles N for two reduced densities ρ∗ . Values for N = ∞ are obtained by extrapolating to the infinite volume limit, as described in the text. ρ∗ = 0.34 ρ∗ = 0.52 3 ∗ ∗ (Σ/σ) N DN DN 1/512 128 1.15 0.635 1/512 256 1.16 0.640 512 1.16 0.642 1/512 1/512 1024 1.17 0.645 1/512 4096 1.17 N/A 1/512 ∞ 1.18 0.657 1/8 128 0.579 0.293 1/8 256 0.588 0.298 1/8 512 0.593 0.302 1/8 1024 0.606 0.305 1/8 4096 0.606 N/A 1/8 ∞ 0.620 0.317 1 128 0.306 0.144 1 256 0.314 0.149 1 512 0.322 0.153 1 1024 0.324 0.157 1 4096 0.333 N/A 1 ∞ 0.348 0.170 8 128 N/A 0.0553 8 256 N/A 0.0614 8 512 N/A 0.0661 1024 N/A 0.0698 8 8 4096 N/A 0.0743 ∞ N/A 0.0844 8 32 128 0.0613 0.0227 32 256 0.0710 0.0289 32 512 0.0789 0.0337 32 1024 0.0843 0.0374 32 4096 0.0920 0.0426 32 ∞ 0.108 0.0525 1024 1024 0.00914 0.00172 1024 4096 0.0165 0.00649 1024 13824 0.0201 0.00905 1024 ∞ 0.0290 0.0146 4096 4096 0.00529 0.00107 4096 13824 0.00929 0.00350 4096 ∞ 0.0173 0.00835 47 Table 5.5: For the “fixed α” class with α = 2/3, reduced diffusion constants of the tracer for different tracer sizes Σ/σ and total number of particles N for the reduced density ρ∗ = 0.52. Values for N = ∞ are obtained by extrapolating to the infinite volume limit, as described in the text. ∗ (Σ/σ)3 N DN 1/512 128 0.571 256 0.575 1/512 1/512 512 0.580 1/512 1024 0.583 1/512 ∞ 0.594 1/8 128 0.262 1/8 256 0.267 1/8 512 0.270 1/8 1024 0.273 1/8 ∞ 0.284 1 128 0.127 1 256 0.132 1 512 0.136 1 1024 0.139 1 ∞ 0.151 8 128 0.0482 8 256 0.0536 8 512 0.0570 8 1024 0.0608 8 4096 0.0654 ∞ 0.0734 8 32 128 0.0196 32 256 0.0249 32 512 0.0292 32 1024 0.0323 32 4096 0.0365 32 ∞ 0.0453 1024 1024 0.00143 1024 4096 0.00548 1024 13824 0.00778 1024 ∞ 0.0124 4096 4096 0.000885 4096 13824 0.00299 4096 ∞ 0.0072 48 Unlike the “fixed I” case, the “fixed α” curves in Fig. 5.4 do not approach the hydrodynamic “slip” lines at large tracer sizes. Rather, the curves for α = 2/5 appear to adopt a slope between the “slip” and “stick” values while the curves for α = 2/3 appear to approach the “stick” line. This latter case represents spheres with α at its largest allowed value. This figure shows a clear correlation between the degree of translational-rotational energy exchange, as quantified by α, and the coefficient appearing in the hydrodynamic description. This data can be used in another manner to explicitly calculate ζ as a function of R from the Stokes-Einstein equation by using the diffusion constants (extrapolated to infinite volume) from Tables 5.2, 5.4 and 5.5. Such data are plotted in Fig. 5.6 for all the calculations performed, as well as for the smooth hard sphere data of Ref. [23]. The trends are apparent. The smooth hard sphere data and the “fixed I” data both approach ζ = 4 for large R – the “slip” limit. The “fixed α” data with α = 2/3 approach ζ ≈ 6 for large R – the “stick” limit – while those with α = 2/5 approach ζ ≈ 5 for large R. This latter value is between the “slip” and “stick” limits. These limiting values of ζ for large R are plotted as a function of α in the top panel of Fig. 5.8. While only three points are present, the dependence looks fairly linear, as seen by the best fit line in the figure. This suggests that the boundary condition ultimately approached for any α is a simple linear combination of the “slip” and “stick” limits. In fluid dynamics, surfaces are often characterized by precisely this approach using an accommodation coefficient. The rough sphere model provides a direct link between accommodation coefficients and parameters related to microscopic scattering events. We are not aware of any theoretical predictions of this dependence for the rough sphere model. These results are also consistent with the results of Schmidt and Skinner[13] who examined the limiting behaviour of an isotropic tracer surrounded by a strongly attractive potential; strong enough to create a first solvation shell of strongly bound bath particles with long residence times (relative to hydrodynamic timescales). They found the system adopted “stick” boundary conditions in this limit. Effectively, the attractive potential produced a shell of particles around the tracer but because the potential was isotropic, there was no mechanism to transfer rotational energy or momentum from this shell to the tracer particle inside. In essence, the “stick” boundary condition arose solely from the motion of the shell of atoms rotating freely around the tracer core. This situation is very similar to the rough hard sphere model with α = 2/3, a case for which “stick” boundary conditions are also observed in the present case. In the low density limit, when the masses and values of α for the bath and tracer are both ∗ (α)/D ∗ (α = 0) = (1+α)/(1+2α), that is the ratio kept equal and constant, µχ = α and D∞ ∞ of the extrapolated diffusion constant for fixed α compared with α = 0 should be a function ∗ (α)/D ∗ (α = 0) = ζ(α = 0)/ζ(α), of α alone. Furthermore, for large enough tracers, D∞ ∞ again when this ratio is taken for tracers with the same size and in baths of the same density. 49 6 5 ζ 4 3 2 ρ = 0.34 1 0 6 5 ζ 4 3 2 ρ = 0.52 1 0 0 1 2 3 4 5 6 7 8 9 R/σ Figure 5.5: Using the extrapolated reduced tracer diffusion constants, and employing Eq. (3.2), values of ζ were calculated and plotted as a function of tracer size, R/σ (R = (Σ + σ)/2) for the reduced densities ρ∗ = 0.34 (upper panel) and ρ∗ = 0.52 (lower panel). The solid and dashed lines denote the hydrodynamic limit with “slip” and “stick” boundary conditions, respectively. The circles, squares, and triangles denote extrapolated results from the “fixed I”, “fixed α = 2/5”, and “fixed α = 2/3” calculations, respectively. The crosses denote values calculated from the smooth hard sphere results of Sokolovskii et al. [23] 50 6.5 6 ζ(α) 5.5 5 4.5 4 0.9 0.8 0.7 * * D (α)/D (α = 0) 3.5 1 0.6 0.5 0 0.1 0.2 0.3 α 0.4 0.5 0.6 0.7 Figure 5.6: For the largest-sized tracers, estimates of ζ for α = 2/5 and α = 2/3 from Fig. 5.4 are plotted as a function of α in the upper panel. The smooth sphere result is also included for α = 0. The solid line is the best linear fit through the three points. In the lower panel, using the smooth hard spheres results of Sokolovskii et al [23] for α = 0, ∗ (α)/D ∗ (α = 0) of extrapolated reduced diffusion constants is plotted for the ratio D∞ ∞ α = 2/5 and α = 2/3 as a function of α. The open and filled symbols denote values for the reduced densities ρ∗ = 0.34 and ρ∗ = 0.52, respectively. The circles, upwards triangles, and downwards triangles denote values Σ/σ = 1, 3.17, and 10.1, respectively. The solid line is the function (1 + α)/(1 + 2α). 51 In other words, curves should exist that universally scale these diffusion constants and that depend only upon α in both the molecular and hydrodynamic limits. This behaviour can be tested using the data in Tables 5.2, 5.4 and 5.5. The bottom panel of Sokolovskii et al. [23], for a few tracer sizes, the ratio of the diffusion constants for “fixed α” relative to those for “α = 0” (these latter values are taken to be the smooth hard sphere results of Ref. [23]) as a function of α for ρ∗ = 0.52. Also on this plot is the functional dependence (1 + α)/(1 + 2α) expected in the low density limit. A number of interesting trends are observed. First, the diffusion constant ratio does decrease as a function of α but this decrease is more than suggested by the low density function. In other words, the α dependence appears to get stronger as the density increases. Second, the values tend to cluster with density, indicating that some small density dependence still remains after the ratio is taken. Third, there is a consistent trend with respect to tracer size with values of the ratio decreasing as tracer size increases, for a fixed density. If one imagines that the diffusion constants for the larger tracers are close to the hydrodynamic limit, then the curve traced out by the filled downwards triangles is ζ(α = 0)/ζ(α). In other words, the α dependence of the ratio is different in the hydrodynamic limit than it is in the low-density (small tracer size) limit. Most of the data values in Fig. 5.6 fall between these limits, that is they are in the transition region between the molecular and hydrodynamic limits. Recall that the smooth and rough hard sphere fluids both have impulsive potentials so the fluid structure is identical in both cases. Thus, the dependencies seen in Fig. 5.6 in the transition region must not be due to configurational changes in the fluid structures but instead to dynamical features related to the different collision models (smooth versus rough) and involve translational and rotational energy exchange. While the theory for describing the transition region is likely difficult, a phenomenological approach can be useful. For smooth hard spheres, Enskog theory predicts that at higher densities the binary diffusion coefficient is given by [11] DE = D12 , g12 (R) (5.2) in which g12 (R) = [Σgb (R) + σgt (R)]/(σ + Σ) with gb (R) and gt (R) being the contact values of the radial distribution function for the bath and tracer, respectively. These values can be related to the equation of state [29], which in the case for rigid spheres is well approximated by the Carnahan-Starling formula [53]. The final result is that gt (R) = 1 1Σ η Σ η 1+ +3 (1 − η) 2 σ (1 − η) σ (1 − η) , (5.3) with the packing fraction (not to be confused with shear viscosity) given by η = πρ∗ /6 = πnσ 3 /6. The equation for gb (R) is identical with that for gt (R) except with Σ/σ set equal 52 to unity. Equation (5.2) holds for the smooth hard sphere model. Unfortunately, we have not been able to find the corresponding result in the literature for the rough sphere model. As an ansatz, Eq. (5.2) shall be applied as is to the rough sphere model and describe it as an “Enskog-like” approximation. This is not wholly unreasonable since the radial distribution function is identical for the rough and smooth sphere models, and diffusion does not depend upon the transport of rotational energy or angular momentum. Thus, the density correction in Enskog theory for rough hard spheres might be very similar to that for smooth hard spheres. Using Eq. (5.2) along with the expression for D12 of Eq. (3.18) for the rough sphere model gives the curved lines in Fig. 5.4. It can be seen that the values of DE approximate quite well the diffusion values for small tracers but of course deviate substantially for larger tracers. This agreement mirrors that seen for the smooth hard sphere model[23]. It is possible to combine the predictions of Eqs. (5.2) and (3.18) to give a phenomenological fit to the diffusion data over the entire tracer size range. In particular, consider the function D= x + Dx DE H , (DE + DH )x−1 (5.4) in which the value of ζ appearing in DH , and the value of D12 appearing in DE are taken as functions of α with the appropriate values taken for each particular calculation. In Fig. 5.7 this function is plotted for the three “fixed α” calculations using the value x = 1.4. It can be seen that the fit is quite reasonable over the whole range of tracer sizes. It is also particularly interesting that this same function, with the same value of x, also gave a quite reasonable fit to the corresponding smooth hard sphere diffusion constants[23]. This is perhaps not too surprising given that the smooth and rough sphere fluids are structurally the same. However, it might also indicate that a function of the form of Eq. (5.4) might have some more general use in estimating the size dependence of diffusion constants. 5.5 Conclusions Overall, the results show that rough hard sphere tracers approach the hydrodynamic limit with boundary conditions spanning the range between “slip” and “stick”. The actual value depends essentially linearly upon the exchange between translational and rotational energy, as characterized by α. The usual qualitative explanation given for adopting the “stick” value is that particles become associated with a tracer, due to an attractive potential. This prolongs the interaction time and allows the colliding particle time to thermalize. It then leaves the tracer with a random velocity and energy determined by this thermalization process. In other words, the collisions are diffuse and this leads to “stick” boundary conditions. The present results show that such attractive interactions are not required to attain “stick” 53 1.4 ∗ ρ = 0.34 1.2 * 2 D R /σ 2 1 0.8 0.6 0.4 0.2 0 0.6 ∗ ρ = 0.52 * 2 D R /σ 2 0.5 0.4 0.3 0.2 0.1 0 1 2 3 4 5 6 7 8 9 R/σ Figure 5.7: Extrapolated reduced diffusion constants for the fixed α calculations as a function of the tracer size. Symbols have the same meaning as in Fig 5.3. The solid curves are the corresponding predictions of Eq. (5.4) with x = 1.4 54 boundary conditions. In fact, because the rough hard sphere has no attractive potential, the fluid structure around it is precisely the same as with the smooth hard sphere (which adopts “slip” boundary conditions). The change in boundary condition for the rough hard sphere results solely from the inelasticity of the collisions with the bath gas. In the low field limit, the mobility of ions is proportional to the diffusion constant so the present calculations can also be used to gauge mobility theories in this limit. For example, the theory of capillary electrophoresis generally treats molecular ions with hydrodynamic theory using “stick” boundary conditions. The present results suggest this description should be valid for large ions but not for small ones. As well, the “stick” boundary condition was found to hold only for the largest value of α, and it is unlikely that many ions adopt this large value. Thus, the use of “stick” boundary conditions for ions in solution may be suspect. Real molecules and ions have a number of general features in common. They have interaction potentials with a repulsive core and an attractive well, and they allow for energy exchange between translational and rotational/vibrational degrees of freedom. Smooth hard spheres have only a repulsive core, and this leads to “slip” boundary conditions in the hydrodynamic limit. Rough hard spheres have a repulsive core and also the ability to exchange translational and rotational energy and momentum. In other words, the rough hard sphere model isolates one important physical effect found in real molecules, the ability to have inelastic collisions. This effect alone changes the boundary condition in the hydrodynamic limit, with values ranging from the “slip” to “stick” limit. 55 Chapter 6 Transport Properties of the Rough Hard Sphere Fluid 6.1 Objectives and Background In this chapter results are presented of extensive benchmark calculations of self-diffusion, shear and bulk viscosity, and thermal conductivity coefficients of the rough hard sphere fluid. The validity of various kinetic theory equations have been examined at various levels of approximation as a function of fluid density and degree of rotational-translational coupling. In particular, expressions from Enskog theory using different numbers of basis sets in the representation of the distribution function were tested. There are several motivations for carrying out this study. The first one is to obtain a reliable set of reference data that can be used for various studies related to kinetic theory. The rough hard sphere model is one of the most common models used in kinetic theory due to it’s theoretical simplicity. It is used in many cases where the effect of inelastic collisions has to be taken into account [30, 32, 34–40, 50, 54–60]. Given such importance, it is valuable to have a set of well-converged benchmark calculations against which comparisons with theory can be made. In the past, theoretical studies of the effect of different basis sets on the values of transport coefficients using the Chapman-Enskog approximation to the Boltzmann equation have been conducted [35]. However, these results have not been compared with actual data for the rough hard sphere fluid since such calculations were too computationally intensive at that time. In this study, the effect of using different basis sets have been re-examined and compared with simulation results. Second, it is interesting to examine the accuracy of Enskog theory for the rough hard sphere model and determine whether it is similar to that observed for the smooth hard sphere case. For larger densities, Enskog theory has been typically used in kinetic theory to predict transport coefficients. The main emphasis has been on the smooth hard sphere fluid but in fact the requisite theory has also been developed for the rough hard sphere fluid [36].There have been no comparisons of the results of this theory with computational data. In this study the comparison has been conducted and various dependencies observed. Third, the comparison of the benchmark values for the rough hard sphere fluid and 56 smooth hard sphere fluid may provide insight into the effect of rotational motion on the behaviour of the fluid. A benchmark calculation for smooth hard sphere fluid has been performed by Sigurgeisson and Heyes[61]. They compare, as a function of density, transport coefficients from simulations with those predicted by both Enskog theory and by solutions of the Boltzmann equation. By changing the parameter affecting translational rotational coupling, the results for the rough hard sphere fluid can be tuned essentially from the smooth sphere limit to the very inelastic limit. By doing so, the sole effect of translational rotational coupling upon transport coefficients can be gauged since all other aspects of the system remain unchanged. By comparing simulation data directly, one is freed from any particular assumptions used in solving specific kinetic equations, and hence can determine general results. This work provides data sets that can serve as a reliable foundation for studying the behaviour of real fluids. It has been shown previously that the rough hard sphere model provides a reasonable description for fluids of small spherical molecules, such as carbon tetrachloride and methane [62–64]. Such studies could be extended to other systems and other transport properties using the present benchmark calculations since these calculations cover the full relevant range in the parameter spaces of density and translational-rotational coupling. This chapter is organized by first providing the computational details of the simulation, followed by the data analysis, with results and conclusions following at the end. 6.2 6.2.1 Methods Computational Details As in Chapter 5, results are reported in scaled values and denoted by a superscript ∗. Specifically, ρ∗ = nσ 3 , D∗ = Dt0 /σ 2 , η ∗ = ησt0 /m, ηb∗ = ηb σt0 /m, and λ∗ = λT t30 /(mσ) with t0 = σ m/kB T . Simulations are performed using the molecular dynamics program based on an event-driven mechanism described in Chapter 4. The number of particles N varies from 32 to 5300. Values of the reduced density ρ∗ range from 0.001 to 0.9512 in increments to match the values used in the smooth hard sphere calculations of Heyes et al.[61]. Calculations were performed for α = 0.004, 0.2, 0.4, and 2/3. For each combination of parameters more than one to two hundred simulation runs have been carried out. All particles were assigned random velocities and angular velocities, and these velocities were then scaled so that kB T = 1 with net linear and angular momenta set to zero. Each run was equilibrated for 1 000 000 collisions. In production runs, 40 000 000 to 160 000 000 collisions were simulated. The correlation functions of Eqs. (4.3)–(4.6) were calculated at fixed time intervals which were generally chosen to be approximately 10τF , 57 where τF is the mean time between collisions. Correlation function values were converged to total times of 700-2000τf depending upon the number of particles in a simulation. These times are sufficient for the correlation functions to approach to their long-time linear forms. The total number of jobs ran for this work was approximately 350 000, with the longest taking about 4 weeks of CPU time for the case of N = 5300 at ρ∗ = 0.001. The data occupy nearly 1TB of disk space. 6.2.2 Convergence To assure the system is well converged, the time dependence of the correlation functions is examined in a manner similar to that described in Chapter 5. Generally, the correlation functions are well converged at short times but not as well converged at longer times because the frequency of statistical sampling decreases as time increases. Figure 6.1 shows a representative sample of raw data for each transport coefficient as a function of time. Each panel contains the average root mean square difference of the G function corresponding to each transport coefficient for four different values of α at a representative density ρ∗ = 0.4298. In order to estimate the transport coefficients, the time dependence of the correlation functions is fit to rational functions which have the correct limiting behaviour at long times. Two rational functions were used for this purpose (with constants d and ci ), f1 (s) = f2 (s) = 6ds2 + c1 s + c2 , s + c3 (6ds2 + c1 s + c2 )s , s2 + c3 s + c4 (6.1) (6.2) which differ only in their overall order. Function f2 (s) is more flexible than f1 (s). After fitting, the constants d, that is equal to the value of Φ for the corresponding transport property for a given density and number of particles, was extracted. When dealing with bulk viscosity, the term ct2 was added to both functions to account for the residual quadratic dependence remaining in the correlation functions. This dependence is usually small but it must be well removed in order for the underlying linear dependence to be obtained accurately. 6.2.3 Data Analysis The data analysis consists of several steps. First, the data is fit independently four times in the following manner. A fit with function f1 (s) is performed using a short time window. This fit is repeated using a long time window. The analogous two fits are performed as well with function f2 (s) so that there are four independent estimates of Φ in total. Second, the four values are compared with one another and the best one is chosen. If the simulation 58 0 500 1000 1500 0 500 1000 1500 2000 7000 700 α = 0.004 α = 0.2 α = 2/5 α = 2/3 600 2 <∆G > 500 6000 5000 400 4000 300 3000 200 2000 GD 100 Gλ 0 1400 0 16 14 1200 12 1000 10 800 2 3 <∆G > (10 ) 1000 8 600 6 400 4 Gη Gη 2 b 200 0 0 0 500 1000 1500 2000 0 t / τf 500 1000 1500 2000 t / τf Figure 6.1: Mean squared values < ∆G2 > as a function of time scaled for ρ∗ = 0.4298. Each panel corresponds to a different transport coefficient and contains data for four different values of α. 59 data is well converged, all four fits will produce fitted parameters with similar values and small errors. In this case, the values from the fit using f2 (s) and the larger time window are chosen, since these are expected to be the most accurate. If the simulation data is not as well converged, the extra flexibility of f2 (s) sometimes causes undesired wiggles or kinks to appear in the fitting function. The function f1 (s) is less susceptible to such behaviour owing to its lower order and will often produce good fitted values in this case. Third, finite size effects are removed by extrapolating the values of Φ to the infinite size limit. This is done using linear regression on the set of Φ for simulations with different numbers of particles. For self-diffusion, this extrapolation is performed as a function of inverse simulation edgelength, as shown by Fushiki [25]. An example of extrapolation of the the diffusion coefficient is shown on Fig. 5.3 in Chapter 5. Similar extrapolations were performed for each density in the current work. For all other properties, the extrapolation is performed as a function of 1/N . This extrapolation was chosen to account for the fact that shear and bulk viscosity and thermal conductivity depend upon the size of the entire system and are calculated for the sum of all particles, i.e. N . In practice, only the simulations with N > 800 were used for this extrapolation. Fig. 6.2 is a sample extrapolation plot for thermal conductivity, shear and bulk viscosity as functions of 1/N at ρ∗ = 0.2865 and α = 0.004. Similar plots have been created for all combinations of ρ∗ and α. Since there are hundreds of fits to perform, it is not efficient to perform them manually. A program was designed that automatically goes through steps described above, suggests the best fits and flags potential errors. Results of the automated analysis are further examined manually to confirm the fit choice. Figure 6.3 illustrates the fitting process of Gn in one of the more challenging cases. In this graph, four fit functions are compared against the raw data. The higher order f2 function is more sensitive to small variations in the data points and thus produced curves that are not suitable for estimating the fit coefficients. The other two are lower order f1 functions; they perform better, albeit not ideal. In such cases, fits from the higher order are discarded, and the fit that is closest to the raw data is chosen. This analysis is performed for each set of parameters and the best values are used for the next step, which is extrapolation to the infinite size limit. The contact values of the radial distribution function, g, in Eqs. (3.23)-(3.27), were calculated for each density using the Carnahan-Starling[53] equation of state, that is g= in which 1 − /2 , (1 − )3 (6.3) = πρ∗ /6 is the packing fraction. 60 2.1 ρ = 0.2865, α = 0.004 λN 2.08 2.06 2.04 ηb N 13 12 11 0.295 ηN 0.29 0.285 0.28 0.275 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 1/N Figure 6.2: Plots of reduced constants of thermal conductivity, shear and bulk viscosity as a function of inverse particle number for ρ∗ = 0.2865 and α = 0.004. The lines are linear function fits to the data, whose y intersept gives the extrapolated value of the corresponding transport coefficient. 61 1500 ∗ 1350 Gη at ρ = 0.4298 and α = 2/5 f1(t), 10 000 samples f1(t), 90 000 samples f2(t), 90 000 samples f2(t), 10 000 samples 1200 1050 900 2 <∆ G > 750 600 450 300 150 0 -150 -300 0 200 400 600 800 1000 1200 1400 1600 1800 2000 t/τf Figure 6.3: The solid line represents the mean squared values of the off–diagonal elements of the pressure tensor plotted as a function of time at number density ρ∗ = 0.4298 and α = 2/5. Dashed and dotted lines represent the four different fits used to extract the value of the transport coefficient. 62 6.3 Results and Discussion Predicted values of transport coefficients, extrapolated to the infinite size limit, are presented in Tables 6.1–6.4. Also included in Tables 6.1–6.4 are predictions from various low density formulae for comparison with the simulation results. The data for two low densities are included for comparison, ρ∗ = 0.01 and 0.001. The change in the transport coefficients on moving from ρ∗ = 0.01 to ρ∗ = 0.001 is only a few percent, so the values for ρ∗ = 0.001 are very close to the low density limit. The actual simulation results as a function of N , from which the extrapolated results were derived, are given in Tables 6.5–6.8. Let’s start by comparing the low density predictions for the transport coefficients with simulation results. For bulk viscosity, shown in Table 6.1, the difference between values predicted by the Pidduck formula of Eq. (3.17) and simulation is about 5 − 6%. However, the Pidduck formula predictions are too large for small α but too small for large α. Bulk viscosity does not exist at low density for the smooth sphere fluid but this is not the case for the rough sphere one. As α grows small, a careful limiting procedure for Eqs. (2.13)– (2.16) shows that the smooth sphere collision conditions are not strictly obtained. Rather, exchange can still occur between the orbital angular momentum of a colliding pair and the rotational angular momentum. This effect is somewhat unphysical since molecular systems do not behave in this manner. In practice this means values of ηb for small α are difficult to predict theoretically and also difficult to estimate numerically. Examining the low density shear viscosity values in Table 6.2 shows predictions from Pidduck’s formula in Eq. (3.16) differ from the simulation values by as much as 8% for α = 2/3, decreasing slightly as alpha decreases. The values predicted by Eq. (3.21) are much better and essentially match exactly the simulation data, within numerical uncertainty. As seen in Table 6.3 for thermal conductivity, Pidduck’s formula of Eq. (3.15) predicts low density values which are too small by about 8% while Eq. (3.20) does much better having errors of only a few percent. For both thermal conductivity and shear viscosity, the inclusion of additional terms in the basis set expansions in kinetic theory leads to much improved theoretical results. For diffusion, as shown in Table 6.4, the difference between the prediction of Pidduck’s formula of Eq. (3.14) and the simulation data varies from about 4% at large α to 1% at small α. It is expected that the error should become smaller as α decreases because Pidduck’s formula works well for the smooth sphere fluid, and the diffusion constant for the rough sphere fluid approaches that of the smooth sphere one when α → 0. In the smooth sphere case only translations contribute to the diffusion coefficient. In the rough sphere case translational-rotational coupling plays a role and the extra terms in the kinetic theory that contribute to Eq. (3.19) help account for this effect. Interestingly, while the predictions of Eq. (3.14) are generally too large, those of Eq. (3.19) are generally too small. Thus, it 63 appears as if there is an overcorrection due to the added terms include in Eq. (3.19). The simulation data are very close to the averages of the predictions of Eqs. (3.14) and (3.19). Consider now the Enskog theory predictions. In order to do so plots of the density dependence of transport coefficients will be compared with theoretical predictions. In almost all cases, the data are normalized by plotting values relative to the low density values (which are taken to be the values at ρ∗ = 0.01 when simulation values are constructed) multiplied by g. With this normalization, the Enskog theory predictions are simply the functions contained within the curly braces on the right hand sides of Eqs. (3.23)–(3.27). Generally speaking, the coefficients for thermal conductivity, shear viscosity, and bulk viscosity should increase as density increases. There are two mechanisms for transport in the fluid, particle motion and collisional transfer. When the density is low, particles travel relatively long distances between collisions, and therefore particle motion is the primary means for transport. In this sense, both rough and smooth hard sphere fluids behave similarly. When the density increases, collisions occur more often, and collisional transfer begins to dominate. In this way, as density is increased the fluid moves from a regime where free flow defines property transfer to a regime where collisional transfer becomes important. This latter mechanism becomes very efficient at high densities leading to transport coefficients that increase in value with density. The strength of Enskog theory is its ability to account in part for the collisional transfer mechanism. For bulk viscosity, simulation data is compared with predictions from Enskog theory in Fig. 6.4. The Enskog predictions of Eq. (3.25) are generally good for densities up to about ρ∗ = 0.6 but begin to deviate substantially for larger densities. Specifically, the theoretical values tend to be too low and the values grow too weakly with density. This is not surprising since the level of approximation used in obtaining Eq. (3.25) allows only for a quadratic dependence upon density, and the simulation data growth at large ρ∗ is much greater than quadratic. It is possible that higher order Enskog theory may account better for this large density behaviour, although this would be quite challenging to obtain analytically. In Fig. 6.5 shear viscosity values from simulation data are compared with Enskog predictions. In this case, Enskog theory was used to give two different predictions. The first uses only the basis functions employed with Pidduck’s formula and leads to Eq. (3.23). The second, leading to Eq. (3.24) uses an additional basis function. In other words, from the point of view of basis functions included in kinetic theory, the Enskog predictions of Eqs. (3.23) and (3.24) are at the same level as the corresponding low density results of Eqs. (3.16) and (3.21), respectively. In general, Eqs. (3.23) and (3.24) predict essentially the same values for all densities. The predictions of these two formulae differ by only a few percent for α = 2/3 and by even less for lower values of α. This appears to be at odds with the fact that in the low density limit, these two equations give results that differ noticeably. It must be remembered though that the quantities plotted in Figs. 6.2-6.6 are 64 800 1400 1200 α=2/3 600 α=2/5 1000 400 300 800 600 η ∗ 500 ∗ B ∗ g / η B(ρ =0) 700 400 200 100 200 0 0 160 700 α=0.2 140 α=0.004 120 500 100 400 80 300 60 η ∗ B ∗ ∗ g / η B(ρ =0) 600 200 40 100 20 0 0 0 0.2 0.4 0.6 ∗ ρ 0.8 0 0.2 0.4 0.6 0.8 1 ∗ ρ Figure 6.4: Values of the bulk viscosity coefficient ηb∗ multiplied by the radial distribution function at contact g and divided by the value of ηb∗ at ρ∗ = 0 plotted as a function of number density ρ∗ . Each panel corresponds to a different value of α. For the simulation data, plotted as open circles, the values for ρ∗ = 0.01 are taken to be the values for ρ∗ = 0. The solid lines are Enskog theory predictions of Eq. 3.25. scaled relative to the values at low density. Any differences in the low density values are automatically scaled to unity. Thus, the essentially identical curves seen in Fig. 6.5 do not imply that the absolute coefficients are predicted equally well by Eqs. (3.23) and (3.24). In truth, Eq. (3.24) is performing much better in this regard than Eq. (3.23). Rather, the identical curves seen in Fig. 6.4 imply that the main effect from including additional basis 65 functions in the kinetic theory leading to these two equations is to correct the low density values. The additional basis functions do very little to change the density dependence of the Enskog result. As seen with the bulk viscosity, the Enskog predictions fall below the simulation data at large densities and grow less rapidly. However for ρ∗ < 0.4 the Enskog values differ by only a few percent from the simulation data. 400 600 α=2/3 α=2/5 300 400 ∗ 300 ∗ ∗ η g / η (ρ =0) 500 200 200 100 100 0 400 0 300 α=0.2 200 200 150 ∗ ∗ ∗ η g / η (ρ =0) 300 250 α=0.004 100 100 50 0 0 0 0.2 0.4 0.6 ∗ ρ 0.8 01 0.2 0.4 0.6 0.8 1 ∗ ρ Figure 6.5: Values of the shear viscosity coefficient η ∗ multiplied by the radial distribution function at contact g and divided by the value of η ∗ at ρ∗ = 0 plotted as a function of number density ρ∗ . Each panel corresponds to a different value of α. For the simulation data, plotted as open circles, the values for ρ∗ = 0.01 are taken to be the values for ρ∗ = 0. The solid and dashed lines are Enskog theory predictions of Eqs. (3.23) and (3.24), respectively. 66 As with shear viscosity, there are also two different Enskog predictions for thermal conductivity, depending upon the number of basis functions retained in the corresponding kinetic theory. Equations (3.26) and (3.27) are the Enskog predictions keeping the same number of basis functions as used to obtain the low density results in Eqs. (3.15) and (3.20), respectively. These are compared with simulation data in Fig. 6.6. For all densities, the predictions of the two equations are almost identical, varying by a few per cent at most from each other. This behaviour is analogous to that seen for shear viscosity and implies that the higher order Enskog expression effectively corrects only the low density values but leaves the density dependence unaltered. Of all the transport properties, the thermal conductivity is predicted the best by Enskog theory, even at higher densities. Overall, the agreement seen in Fig. 6.6 is quite reasonable. In fact, for ρ∗ < 0.5 the Enskog and simulation data differ by only a few percent. The diffusion coefficient data is plotted in Fig. 6.7. We are unaware of any Enskog prediction for this coefficient for rough spheres. For smooth spheres, Enskog theory gives DE = D/g [31], which if applied with the normalization in Fig. 6.7 would predict, for the data plotted as squares, a constant value of one. Clearly, the simulation data are not constant and deviate from one by up to 50%. This behaviour is distinct from the other transport coefficients but is consistent with that observed for the smooth hard sphere fluid [65, 66] (which has the same structure as the rough hard sphere fluid). We will argue below this curved shape of the data is an artefact of the way the data is plotted, and results from the factor of g. When the data is plotted as the diffusion ratio without this factor, ρ∗ D∗ /ρ∗ D∗ (ρ∗ = 0), shown as circles in Fig. 6.7, the peak disappears. For each value of α, the diffusion ratios plotted as circles in Fig. 6.5 show the same trend, being composed of two linear regions. The first region extends from a reduced density of zero to about 0.4, the second from about 0.5 to 0.8. The data is quite linear in each region but the magnitude of the slope in the first region is smaller than in the second region, and the transition between the two occurs at approximately ρ∗ = 0.5. The velocity autocorrelation function for smooth hard spheres [52] shows that for reduced densities greater than 0.8, caging effects become important and the diffusion coefficient begins to decrease even more rapidly. The decrease in the diffusion constant at lower densities is due to increased density around a particle. Radial distribution functions were calculated for the systems shown in Fig. 6.7 and one case in shown in Fig. 6.8. For small ρ, the values show a single monotonic decay with contact values increasing with ρ∗ . For ρ ≈ 0.5, a peak in the radial distribution function begins to grow at approximately r = 2σ. For ρ∗ > 0.8 an additional peak begins to grow at approximately r = 3σ. This behaviour correlates well with the dependence plotted as circles in Fig. 6.7. In other words, the linear decrease in the diffusion ratio from ρ∗ = 0 to 0.4 results from a buildup of density near the particle surface that reduces the diffusion coefficient beyond 67 120 120 100 α=2/5 80 60 60 40 40 20 20 0 120 0 140 ∗ 80 ∗ α=2/3 ∗ λ g / λ (ρ =0) 100 α=0.2 ∗ λ g / λ (ρ =0) 100 α=0.004 120 100 80 80 60 60 40 40 20 20 0 0 0 0.2 0.4 0.6 ∗ ρ 0.8 0 0.2 0.4 0.6 0.8 1 ∗ ρ Figure 6.6: Values of the thermal conductivity coefficient λ∗ multiplied by the radial distribution function at contact g and divided by the value of λ∗ at ρ∗ = 0 plotted as a function of number density ρ∗ . Each panel corresponds to a different value of α. For the simulation data, plotted as open circles, the values for ρ∗ = 0.01 are taken to be the values for ρ∗ = 0. The solid and dashed lines are Enskog theory predictions of Eqs. (3.23) and (3.24), respectively. the value expected by scaling against the bulk density. This decrease though is less than the rate of increase of g so that when the factor of g is included in the diffusion ratio data, plotted as squares in Fig. 6.7, the resulting values increase linearly. At ρ∗ ≈ 0.5, the second solvation shell begins to form around a particle, and this causes the diffusion coefficient to 68 1.4 1.4 α=2/3 1.2 α=2/5 1 1 ∗ ρ D g / ρ D (ρ =0) 1.2 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 1.6 0 1.6 ∗ ∗ 0.8 α=0.2 α=0.004 1.4 1.2 1.2 1 1 ∗ ∗ ρ D g / ρ D (ρ =0) 1.4 0.8 0.6 0.6 0.4 0.4 0.2 0.2 ∗ 0.8 0 0 0 0.2 0.4 0.6 0.8 0 ∗ ρ 0.2 0.4 0.6 0.8 1 ∗ ρ Figure 6.7: Squares represent the values of ρ∗ D∗ g divided by ρ∗ D∗ at ρ∗ = 0 plotted as a function of number density ρ∗ . Circles represent the values of ρ∗ D∗ divided by ρ∗ D∗ at ρ∗ = 0 plotted as a function of number density ρ∗ . Each panel corresponds to a different value of α, as indicated. The values for ρ∗ = 0.01 are taken to be the values for ρ∗ = 0. decrease even faster, giving a linear curve of greater negative slope. This slope exceeds the rate of growth of g so the diffusion ratio data plotted as squares begins to decrease with increasing density, producing the peaked shape. In some cases it has been argued that the peaked shape in Fig. 6.7 is the result of vortex velocity fields forming around a particle, leading to a long-time tail in the correlation function decaying as t−3/2 [67]. The contribution from this tail results in a larger diffusion 69 3 * ρ = 0.191 * ρ = 0.382 2.5 * ρ = 0.5253 * ρ = 0.7163 * ρ = 0.8404 g (r) 2 1.5 1 0.5 0 1 1.5 2 2.5 3 3.5 r/σ Figure 6.8: Values of the radial distribution function plotted as a function of r/σ for α = 0.2 at various reduced densities. These values were calculated from the simulation with N = 1360 particles. constant than otherwise expected. While it is likely such tails are contributing some fraction to the diffusion constants at higher reduced densities, the general trends do not appear to be affected by them. For example, diffusion ratios plotted as circles in Fig. 6.7 begin to decrease linearly as soon as ρ∗ exceeds zero. At these low densities, the contributions from the long-time tails are negligible yet the data plotted as squares and circles in Fig. 6.7 begin to diverge immediately. The simplest explanation is that Enskog theory has a fundamental difficulty predicting diffusion coefficients accurately, and this likely results from the molecu70 lar chaos approximation inherent within it. As density increases, diffusion will be sensitive to correlations in the bath, and such correlations are not accounted for in Enskog theory. In other words, the value of g alone cannot account for the structure of the fluid away from a particle, an area that significantly affects diffusion. Consider now the effect of translational-rotational coupling on the transport coefficients by plotting data as a function of α for fixed densities. Specifically, the ratio of the transport coefficient to its value at α = 0 is plotted where in the present case the values for α = 0.004 are taken to be the values for α = 0 when plotting the simulation data. Analogously, the theoretical predictions for the transport coefficients both at low density and at the level of Enskog theory are plotted using the same ratio. In essence, the α dependence of these expressions is being tested. The only difficulty occurs with the bulk viscosity for which the α = 0 value diverges. For this specific case, the theoretical values are normalized to their values at α = 0.004 in precisely the same manner as for the simulation data. The results are plotted in Figs. 6.9 – 6.12. The bulk viscosity, as seen in Fig. 6.9, shows a strong α dependence especially for small α where it diverges. For α > 0.2 the dependence is much less varying. The lines in Fig. 6.9 show the theoretical predictions of Eq. (3.25) at the Enskog level for the α dependence. Generally the theoretical curves give a good approximation to the simulation data at low densities. For intermediate and larger densities the deviations are more significant. In all cases the theoretical values are larger than the simulation data. For the shear viscosity, as seen in Fig. 6.10, an even stronger increase is observed with α. At the highest density, the coefficient increases by almost 2.5 times in comparison with the smallest value of α. This shows that translational-rotational coupling greatly affects the shear viscosity coefficient. Theoretical predictions from both Eqs. (3.23) and (3.24) are also shown. The Enskog expression of Eq. (3.24) includes more basis functions in its kinetic theory derivation, and predicts values slight larger than those of Eq. (3.23) which includes less basis functions. At low densities, Eq. (3.24) gives quite good results but this agreement begins to deteriorate at intermediate densities. The theoretical curves approach a limiting value at high density which is significantly smaller than the simulation data. This undervaluation is analogous to that seen in Fig. 6.5 at higher densities. Overall Enskog theory is not performing well for high densities either as a function of density or of α. Translational-rotational coupling does not affect the thermal conductivity coefficient as much as the shear viscosity coefficient, as shown in Fig. 6.11. When looking at the full range of α values, the coefficient changes by about 5% at the lowest density and by about 15% for highest density. Interestingly the trend changes noticeably with density. At low densities the coefficient increases with α but this trend slowly and continuously changes until at high densities it decreases with α. At higher densities, collisional transfer of energy becomes important in determining the thermal conductivity coefficient, and the data indicates that 71 0.12 1 0.11 ρ∗=0.01 ρ∗=0.2865 ρ∗=0.4298 ρ∗=0.6685 ρ∗=0.8404 0.8 0.6 0.1 0.4 0.2 0.09 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.07 ∗ η b / η b(α=0.004) 0.08 ∗ 0.06 0.05 0.04 0.03 0.02 0.01 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 α Figure 6.9: Values of the bulk viscosity coefficient ηb∗ divided by the value of ηb∗ for α = 0.004 plotted as a function of α. Symbols represent the values corresponding to different reduced densities, as shown in figure legend. Solid lines represent the Enskog predictions of Eq. (3.25). this process becomes less efficient as the degree of translational-rotational energy exchange increases. The Enskog predictions of Eqs. (3.26) and (3.27) both fail to capture this change. This is a bit surprising since this same theory gave quite reasonable agreement for thermal conductivity coefficients as a function of density. The scale in Fig. 6.11 is much smaller though compared with Fig. 6.6 so in a sense, the prediction of the α dependence is more difficult. All the theoretical curves predict coefficients that increase as a function of α for 72 2.2 2.1 * ρ = 0.01 * ρ = 0.2865 * ρ = 0.4298 * ρ = 0.6685 * ρ = 0.8404 2 1.9 1.8 1.6 ∗ ∗ η / η (α=0) 1.7 1.5 1.4 1.3 1.2 1.1 1 0.9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 α Figure 6.10: Values of the shear viscosity coefficient η ∗ divided by the value of η ∗ for α = 0 plotted as a function of α. For the simulation data, plotted as symbols, values for α = 0.004 are taken to be the values at α = 0. Symbols represent the values corresponding to different reduced densities, as shown in the figure legend. Solid and dashed lines represent the Enskog predictions of Eqs. (3.23) and (3.24), respectively. all densities. It is likely that many basis functions would have to be included in the kinetic theory formulation of the Enskog result in order to predict the correct α dependency. The diffusion coefficient is shown in Fig. 6.12 where it shows a fairly strong dependence upon α, decreasing by a factor of 2 for the largest α and density. The trends are smooth and show that translational-rotational coupling affects the diffusion constant more at high 73 1.1 1 ∗ ∗ λ / λ (α=0) 1.05 0.95 * ρ = 0.01 * ρ = 0.2865 * ρ = 0.4298 * ρ = 0.6685 * ρ = 0.8404 0.9 0 0.1 0.2 0.3 α 0.4 0.5 0.6 0.7 Figure 6.11: Values of the thermal conductivity coefficient λ∗ divided by the value of λ∗ for α = 0 plotted as a function of α. For the simulation data, plotted as symbols, the values for α = 0.004 are taken to be the values at α = 0. Symbols represent the values corresponding to different reduced densities, as shown in the figure legend. Solid and dashed lines represent the Enskog predictions of Eqs. (3.26) and (3.27), respectively. densities than at low densities. This is not surprising because the increase in translationalrotational coupling at higher α is expected to change the velocities more, causing a faster decay of the velocity auto-correlation function and hence a smaller diffusion constant. 74 1.05 1 * ρ = 0.01 * ρ = 0.2865 * ρ = 0.4298 * ρ = 0.6685 * ρ = 0.8404 0.95 0.9 0.85 0.75 * * D / D (α=0) 0.8 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0 0.1 0.2 0.3 0.4 0.5 0.6 α Figure 6.12: Values of the self-diffusion coefficient D∗ divided by the value of D∗ for α = 0 plotted as a function of α. For the simulation data, plotted as symbols, the values for α = 0 are taken to be the values at α = 0. Symbols represent the values corresponding to different reduced densities, as shown in the figure legend. 6.4 Summary The validity of various kinetic theory relations in the rough hard sphere fluid have been investigated for a broad range of densities and degrees of translational-rotational coupling. At low density the simulation results are closer to theoretical predictions that include higher levels of approximation. In the cases of shear viscosity, thermal conductivity and self-diffusion, 75 the agreement improves by a few percent for each transport coefficient in comparison with estimates given by Pidduck’s formulae. The quality of different Enskog predictions was also examined. As a function of density, Enskog theory predicts the thermal conductivity coefficient the best. The agreement between simulation data and theoretical predictions is good over the whole range of densities studied. This is also consistent with the analogous results for smooth hard spheres[61] where the same level of agreement is observed. This is likely because the transport of energy is fairly insensitive to correlations in the fluid. For shear and bulk viscosity, the agreement between theory and simulation data is worse. Essentially, the Enskog predictions are reasonable only for reduced densities below about 0.5 where they differ from simulation data by only few percent. This is in contrast to behaviour seen with smooth hard spheres[61] in which agreement between Enskog theory and simulation data is good to much higher densities. The collisional transfer of momentum, important at high densities, is not being well described by the Enskog theory. The low density predictions though are very good. For some of the transport coefficients, it was possible to compare the predictions of Enskog theory using different numbers of basis functions included in the kinetic theory description of the distribution function. The Pidduck formulae include the same basis functions commonly used for the smooth sphere fluid so the accuracy of these formulae gives a measure of the degree to which only translational effects are important in determining transport coefficients. The higher order Enskog formulae used in this work include additional terms arising from translational-rotational polarization effects. In principle, these higher order formulae should be more accurate than the Pidduck ones (and this was generally found to be the case) because they incorporate effects arising from the exchange of rotational and translational energy and momentum upon collision. Interestingly, it was found that the density dependence (when normalized against the value of a coefficient at low density) of the coefficients was virtually unaffected by these additional basis functions, with the Pidduck and higher order expressions predicting essentially the same result. The main effect seems to be the correcting the low density values of the coefficients. The more stringent test of Enskog theory involves prediction of the coefficients as a function of α rather than density. This is the case for two reasons. First, the coefficients change much less as a function of α than of density, so accurate predictions are needed to capture them. Second, varying α alters the degree of translational-rotational coupling. Energy exchange is zero when α = 0 and a maximum when α = 2/3. Therefore, the entire change in the coefficients relies upon an accurate representation, in the kinetic theory, of the exchange between these two degrees of freedom. This is more challenging than determining the effects arising from translational motions alone (as is needed in describing the smooth sphere fluid). Generally, as a function of α, the Enskog predictions are poor. In some cases, as with the bulk and shear viscosity, the qualitative trends are good but the quantitative 76 values can be poor. In other cases, as with thermal conductivity, even the qualitative trends can be incorrect. In principle, Enskog theory might perform much better were more basis functions included in the expansion of the distribution function. The theoretical framework to do this has already been presented[36] but the mathematical effort required in doing so is formidable. This would also necessitate a purely numerical approach as evidenced by the complexity of the higher order Enskog formulae used in this work. It is likely that many basis functions would be required in order to converge the results to within the tolerance necessary to predict the dependencies on α. Overall, a benchmark calculation has been performed for the rough hard sphere fluid. A series of transport coefficients covering a wide range of density and reduced moment of inertia have been reported. These values can be used for comparison with different theoretical predictions, and to gauge the effect of translational-rotational coupling on transport coefficients. 77 Table 6.1: Bulk viscosity coefficients ηb∗ for a range of reduced densities ρ∗ and corresponding packing fraction at different values of the dimensionless radius of gyration α. The reported values were calculated using the Eq. (4.10), as described in the text. The values have been extrapolated to the infinite size limit. The low density values denoted by (a) are predictions of Eq. (3.17). ρ∗ α = 2/3 α = 2/5 α = 0.2 α = 0.004 (a) 0.0735 0.0864 0.127 4.44 0.0005 0.001 0.0779 0.0928 0.121 4.24 0.0922 0.133 5.39 0.005 0.0100 0.0779 0.100 0.1910 0.176 0.197 0.269 9.19 0.150 0.2865 0.292 0.308 0.410 11.5 0.175 0.3343 0.355 0.390 0.498 13.4 0.200 0.3820 0.464 0.505 0.622 16.5 0.225 0.4298 0.583 0.637 0.780 19.1 0.745 0.790 0.964 21.3 0.250 0.4775 0.275 0.5253 0.924 0.980 1.18 28.7 0.300 0.5730 1.17 1.21 1.45 30.1 0.325 0.6208 1.51 1.56 1.86 42.4 1.91 1.97 2.28 42.9 0.350 0.6685 0.375 0.7163 2.63 2.50 2.83 47.0 0.400 0.7640 3.66 3.52 3.69 75.6 0.425 0.8118 5.22 4.85 5.05 91.6 0.440 0.8404 6.90 6.12 6.04 95.5 8.37 7.34 7.39 99.3 0.450 0.8595 0.480 0.9168 17.6 15.2 13.6 136 0.490 0.9359 22.8 20.5 16.9 152 78 Table 6.2: Shear viscosity coefficients η ∗ for a range of reduced densities ρ∗ and corresponding packing fraction at different values of the dimensionless radius of gyration α. The values have been extrapolated to the infinite size limit. The low density values denoted by (a) and (b) are predictions of Eqs. (3.16) and (3.21), respectively. ρ∗ α = 2/3 α = 2/5 α = 0.2 α = 0.004 (a) 0.200 0.185 0.177 0.176 (b) 0.217 0.195 0.181 0.176 0.0005 0.0010 0.216 0.194 0.184 0.177 0.219 0.199 0.185 0.180 0.005 0.0100 0.100 0.1910 0.327 0.286 0.255 0.229 0.150 0.2865 0.427 0.375 0.327 0.287 0.175 0.3343 0.511 0.438 0.379 0.331 0.200 0.3820 0.602 0.507 0.443 0.376 0.225 0.4298 0.723 0.606 0.520 0.437 0.860 0.727 0.623 0.515 0.250 0.4775 0.275 0.5253 1.05 0.898 0.752 0.615 0.300 0.5730 1.30 1.08 0.920 0.734 0.325 0.6208 1.59 1.36 1.14 0.897 2.03 1.70 1.41 1.10 0.350 0.6685 0.375 0.7163 2.59 2.19 1.84 1.37 0.400 0.7640 3.41 2.87 2.39 1.76 0.425 0.8118 4.74 3.96 3.22 2.32 0.440 0.8404 5.87 4.93 4.03 2.83 6.90 5.78 4.61 3.23 0.450 0.8595 0.480 0.9168 12.1 10.1 7.99 5.26 0.490 0.9359 15.9 12.8 10.2 6.61 79 Table 6.3: Thermal conductivity coefficients λ∗ for a range of reduced densities ρ∗ and corresponding packing fraction at different values of the dimensionless radius of gyration α. The values have been extrapolated to the infinite size limit. The low density values denoted by (a) and (b) are predictions of Eqs. (3.15) and (3.20), respectively. ρ∗ α = 2/3 α = 2/5 α = 0.2 α = 0.004 (a) 1.03 1.01 0.993 0.978 (b) 1.12 1.12 1.09 1.02 0.0005 0.0010 1.11 1.09 1.08 1.06 1.11 1.12 1.09 1.06 0.005 0.0100 0.100 0.1910 1.60 1.59 1.59 1.60 0.150 0.2865 2.07 2.08 2.08 2.09 0.175 0.3343 2.39 2.39 2.38 2.43 0.200 0.3820 2.75 2.81 2.80 2.86 0.225 0.4298 3.22 3.20 3.27 3.39 3.79 3.85 3.86 4.09 0.250 0.4775 0.275 0.5253 4.48 4.53 4.62 4.82 0.300 0.5730 5.30 5.39 5.59 5.76 0.325 0.6208 6.37 6.40 6.53 6.95 7.62 7.68 7.92 8.36 0.350 0.6685 0.375 0.7163 9.12 9.20 9.47 10.2 0.400 0.7640 11.0 11.1 11.5 12.3 0.425 0.8118 13.3 13.5 13.9 15.0 0.440 0.8404 15.0 15.2 15.6 16.7 16.1 16.3 17.0 18.3 0.450 0.8595 0.480 0.9168 20.5 20.6 21.3 23.5 0.490 0.9359 22.2 22.3 23.2 25.2 80 Table 6.4: Products of the density and diffusion coefficient ρ∗ D∗ for a range of reduced densities ρ∗ and corresponding packing fraction at different values of the dimensionless radius of gyration α. The values have been extrapolated to the infinite size limit. The low density values denoted by (a) and (b) are predictions of Eqs. (3.14) and (3.19), respectively. ρ∗ α = 2/3 α = 2/5 α = 0.2 α = 0.004 (a) 0.1511 0.1646 0.1814 0.2107 (b) 0.1438 0.1568 0.1747 0.2105 0.0005 0.0010 0.1470 0.1604 0.1786 0.2146 0.1456 0.1586 0.1766 0.2122 0.005 0.0100 0.100 0.1910 0.1238 0.1355 0.1509 0.1807 0.150 0.2865 0.1131 0.1250 0.1399 0.1680 0.175 0.3343 0.1071 0.1188 0.1336 0.1613 0.200 0.3820 0.1003 0.1118 0.1267 0.1540 0.225 0.4298 0.09298 0.1042 0.1187 0.1458 0.1366 0.250 0.4775 0.08515 0.09598 0.1102 0.275 0.5253 0.07708 0.08718 0.1008 0.1266 0.300 0.5730 0.06845 0.07801 0.0908 0.1159 0.325 0.6208 0.05993 0.06855 0.0804 0.1044 0.09236 0.350 0.6685 0.05137 0.05898 0.06988 0.375 0.7163 0.04283 0.04949 0.05893 0.07994 0.400 0.7640 0.03461 0.04028 0.04852 0.06744 0.425 0.8118 0.02682 0.03142 0.03829 0.05475 0.440 0.8404 0.02245 0.02642 0.03243 0.04729 0.04236 0.450 0.8595 0.01970 0.02319 0.02864 0.480 0.9168 0.01216 0.01446 0.01813 0.02822 0.490 0.9359 0.009988 0.01189 0.01500 0.02386 81 Table 6.5: Values of densities ρ∗ , 1 V kB T Φ(Gηb ) from Eq. 4.6 in reduced units for a range of reduced dimensionless radii of gyration α and number of particles N spanning from 32 to 5300. In order to obtain the reduced bulk viscosity coefficients ηb∗ , 34 η ∗ should be subtracted from each number, where η ∗ is the corresponding reduced shear viscosity coefficient from Table 6.2. ρ∗ 0.0010 α 0.2865 0.3343 0.3820 250 0.3754 490 850 1360 4.6504 4.6420 4.5175 2900 5300 0.3774 0.3796 0.3728 0.3761 0.3764 0.3698 0.3729 2/5 0.3518 0.3517 0.3524 0.3514 0.3519 2/3 0.3649 0.3648 0.3650 0.3659 0.004 4.803 4.732 4.7706 4.8951 4.6380 5.484 5.743 0.3822 0.3803 0.3672 0.3815 0.3800 0.3805 0.3805 2/5 0.3566 0.3568 0.3569 0.3530 0.3593 0.3581 2/3 0.3692 0.3715 0.3722 0.3714 0.3736 0.3687 0.3716 0.2 0.1910 100 0.004 0.2 0.0100 32 0.3789 0.004 8.260 8.119 8.033 9.419 9.330 9.448 9.562 9.457 0.2 0.5932 0.5975 0.5885 0.6038 0.5980 0.6061 0.6096 0.6070 2/5 0.3566 0.3568 0.3569 0.3529 0.3593 0.3581 0.3542 2/3 0.5954 0.6019 0.6060 0.6058 0.6083 0.6099 0.6055 0.6119 0.004 11.17 11.04 10.93 12.92 10.88 12.93 12.58 12.78 0.2 0.8259 0.8391 0.8355 0.8384 0.8429 0.8392 0.8444 0.8462 2/5 0.5617 0.5690 0.5702 0.5670 0.5729 0.5768 0.5762 0.5746 2/3 0.8480 0.8570 0.8581 0.8603 0.8603 0.8623 0.8615 0.8617 0.004 13.17 12.92 12.59 15.23 12.57 15.06 15.00 15.11 0.2 0.9939 0.9953 0.9993 1.003 1.000 0.9982 1.007 1.005 2/5 0.7932 0.8060 0.8035 0.8024 0.8042 0.8092 0.8035 0.8083 2/3 1.025 1.030 1.029 1.033 1.036 1.037 1.039 1.036 0.004 15.44 15.30 15.31 17.66 14.97 17.48 18.19 17.71 0.2 1.204 1.203 1.195 1.224 1.205 1.208 1.217 1.213 2/5 0.9565 0.9609 0.9710 0.9528 0.9712 0.9729 0.9709 0.9791 82 Table 6.5 – continued from previous page ρ∗ 0.4298 0.4775 0.5253 0.5730 0.6208 0.6685 0.7163 0.7640 α 32 100 250 490 850 1360 2900 5300 2/3 1.250 1.256 1.246 1.250 1.255 1.260 1.260 1.268 0.004 18.27 17.56 17.89 20.96 17.68 20.95 21.04 20.66 0.2 1.464 1.461 1.469 1.454 1.468 1.446 1.474 1.467 2/5 1.162 1.171 1.173 1.183 1.172 1.178 1.174 1.187 2/3 1.530 1.527 1.515 1.561 1.534 1.542 1.541 1.537 0.004 21.63 21.07 21.31 20.15 20.80 25.00 18.47 25.25 0.2 1.790 1.777 1.765 1.791 1.778 1.799 1.802 1.793 2/5 1.425 1.433 1.429 1.424 1.440 1.437 1.431 1.444 2/3 1.887 1.873 1.878 1.871 1.871 1.887 1.885 1.901 0.004 25.62 25.38 23.87 22.46 24.35 29.41 29.14 29.87 0.2 2.206 2.173 2.175 2.175 2.176 2.189 2.189 2.171 2/5 1.758 1.745 1.738 1.751 1.764 1.759 1.745 1.758 2/3 2.361 2.321 2.321 2.330 2.322 2.322 2.330 2.316 0.004 31.20 29.47 29.37 34.87 29.32 34.88 25.64 35.82 0.2 2.719 2.693 2.626 2.681 2.701 2.671 0.9193 2.672 2/5 2.171 2.114 2.155 2.127 2.149 2.149 2.162 2.181 2/3 2.955 2.888 2.921 2.895 2.892 2.864 2.907 2.897 0.004 37.11 35.46 35.83 33.11 36.09 42.24 43.29 42.41 0.2 3.387 3.350 3.306 3.284 3.315 3.350 3.392 3.335 2/5 2.715 2.665 2.692 2.605 2.660 2.660 2.653 2.669 2/3 3.750 3.630 3.632 3.605 3.642 3.645 3.637 3.649 0.004 44.60 42.73 42.56 49.92 42.57 37.45 50.33 53.09 0.2 4.272 4.160 4.173 4.057 4.149 4.158 4.156 4.156 2/5 3.431 3.337 3.344 3.330 3.362 3.353 3.347 3.248 2/3 4.866 4.620 4.643 4.633 4.617 4.605 4.639 4.640 0.004 54.01 52.49 51.73 60.88 50.75 60.90 61.91 45.11 0.2 5.550 5.221 5.298 5.255 5.281 5.230 5.259 5.298 2/5 4.414 4.210 4.237 4.243 4.231 4.230 4.226 4.253 2/3 6.520 6.016 6.038 6.083 6.083 6.084 6.070 6.081 0.004 65.72 62.16 62.08 73.81 61.80 106.2 74.79 72.03 83 Table 6.5 – concluded from previous page ρ∗ 0.8118 0.8404 0.8595 0.9168 0.9359 α 32 100 250 490 850 1360 2900 5300 0.2 7.480 6.832 6.848 6.842 6.784 6.892 6.880 6.747 2/5 5.844 5.462 5.463 5.515 5.472 5.442 5.458 5.459 2/3 9.255 8.055 8.147 8.163 8.170 8.259 8.178 8.143 0.004 79.84 77.02 73.99 89.50 75.52 133.4 91.00 91.52 0.2 10.79 9.144 9.381 9.281 9.342 9.370 9.289 9.333 2/5 8.138 7.177 7.366 7.483 7.276 7.319 7.329 7.322 2/3 14.16 11.44 11.54 11.59 11.61 11.57 11.51 11.63 0.004 91.42 85.23 84.73 99.53 85.58 103.5 102.4 103.5 0.2 13.92 11.22 11.30 11.45 11.43 11.47 11.39 11.47 2/5 12.24 9.979 10.05 10.20 10.19 10.17 10.14 10.12 2/3 18.90 14.39 14.59 14.71 14.69 14.55 14.75 14.75 0.004 98.24 92.57 91.72 110.5 92.00 109.0 109.7 109.5 0.2 16.79 13.06 13.34 13.36 13.46 13.51 13.46 13.47 2/5 16.21 12.54 12.74 12.80 12.77 12.70 12.74 12.70 2/3 23.44 16.99 17.39 17.57 17.48 17.40 17.64 17.54 0.004 613.3 132.6 133.2 142.9 131.1 144.8 143.7 145.2 0.2 809.4 23.32 24.25 24.31 23.90 24.39 23.96 24.30 2/5 19.95 14.78 15.03 15.25 15.22 15.03 15.21 15.08 2/3 1041. 33.83 32.95 34.48 33.80 33.81 33.98 33.96 148.8 136.0 147.1 158.3 152.5 159.7 0.004 0.2 31.59 30.60 30.87 30.76 30.68 30.49 30.64 2/5 28.61 28.26 28.26 28.42 27.75 28.85 28.80 2/3 45.84 44.89 45.67 44.12 44.75 44.75 44.65 84 Table 6.6: Reduced shear viscosity coefficients η ∗ for a range of reduced densities ρ∗ , dimensionless radii of gyration α, and number of particles N spanning from 32 to 5300. ρ∗ α 32 100 250 490 850 1360 2900 5300 0.0010 0.004 0.2 0.1819 0.1827 0.1839 2/5 0.0100 0.2865 0.3343 0.3820 0.4298 0.1784 0.1718 0.1829 0.1829 0.1861 0.1971 0.1967 0.1979 0.2173 0.2126 0.2167 0.1854 0.1813 2/3 0.2149 0.004 0.1772 0.1744 0.1796 0.1797 0.1795 0.1806 0.1799 0.1841 0.1846 0.1859 0.1884 0.1852 0.1854 0.1880 2/5 0.1974 0.2012 0.1967 0.1964 0.1981 0.1996 0.1981 2/3 0.2176 0.2182 0.2174 0.2169 0.2197 0.2188 0.2176 0.2 0.1910 0.1910 0.1807 0.004 0.2142 0.2220 0.2243 0.2268 0.2256 0.2279 0.2296 0.2232 0.2 0.2395 0.2491 0.2474 0.2484 0.2506 0.2557 0.2530 0.2539 2/5 0.2712 0.2783 0.2829 0.2813 0.2828 0.2843 0.2833 0.2861 2/3 0.3133 0.3208 0.3220 0.3236 0.3233 0.3273 0.3268 0.3263 0.004 0.2703 0.2794 0.2803 0.2848 0.2775 0.2836 0.2844 0.2894 0.2 0.3137 0.3205 0.3196 0.3226 0.3296 0.3236 0.3315 0.3255 2/5 0.3615 0.3674 0.3695 0.3729 0.3749 0.3744 0.3750 0.3752 2/3 0.4227 0.4289 0.4294 0.4318 0.4335 0.3382 0.4333 0.4330 0.004 0.3111 0.3169 0.3196 0.3173 0.3232 0.3233 0.3299 0.3312 0.2 0.3665 0.3718 0.3777 0.3760 0.3788 0.3775 0.3784 0.3765 2/5 0.4303 0.4286 0.4355 0.4343 0.4360 0.4378 0.4385 0.4368 2/3 0.4981 0.5038 0.5015 0.5073 0.5099 0.5092 0.5049 0.5112 0.004 0.3615 0.3711 0.3666 0.3719 0.3647 0.3751 0.3753 0.3762 0.2 0.4288 0.4380 0.4317 0.4304 0.4456 0.4463 0.4432 0.4431 2/5 0.5033 0.5068 0.5084 0.5123 0.5122 0.5107 0.5073 0.5073 2/3 0.5947 0.5978 0.5971 0.6111 0.6034 0.5944 0.6041 0.6024 0.004 0.4248 0.4300 0.4322 0.4347 0.4329 0.4371 0.4366 0.4373 Continued on next page 85 Table 6.6 – continued from previous page ρ∗ 0.4775 0.5253 0.5730 0.6208 0.6685 0.7163 0.7640 α 32 100 250 490 850 1360 2900 5300 0.2 0.5132 0.5182 0.5206 0.5230 0.5308 0.5197 0.5151 0.5237 2/5 0.6005 0.6057 0.6079 0.6086 0.6112 0.5963 0.6018 0.6102 2/3 0.7107 0.7137 0.7142 0.7166 0.7167 0.7119 0.7175 0.7226 0.004 0.5049 0.5096 0.5153 0.5128 0.5117 0.4998 0.5181 0.5146 0.2 0.6164 0.6224 0.6230 0.6186 0.6199 0.6201 0.6253 0.6220 2/5 0.7249 0.7261 0.7277 0.7310 0.7207 0.7313 0.7273 0.7299 2/3 0.8535 0.8582 0.8588 0.8623 0.8623 0.8614 0.8614 0.8598 0.004 0.6067 0.6076 0.5956 0.5982 0.6135 0.6146 0.6109 0.6160 0.2 0.7462 0.7461 0.7560 0.7540 0.7515 0.7519 0.7524 0.7531 2/5 0.8803 0.8825 0.9102 0.8819 0.8890 0.8921 0.8764 0.9018 2/3 1.037 1.044 1.051 1.043 1.004 1.038 1.048 1.050 0.004 0.7342 0.7341 0.7350 0.7492 0.7400 0.7340 0.7347 0.7355 0.2 0.9169 0.9101 0.9111 0.8939 0.9160 0.9175 0.9193 0.9182 2/5 1.080 1.083 1.082 1.084 1.087 1.087 1.081 1.084 2/3 1.275 1.246 1.283 1.250 1.277 1.283 1.284 1.295 0.004 0.8949 0.8837 0.8895 0.8935 0.8962 0.8981 0.8928 0.8947 0.2 1.126 1.128 1.127 1.157 1.122 1.128 1.146 1.137 2/5 1.349 1.332 1.341 1.327 1.340 1.347 1.354 1.348 2/3 1.590 1.587 1.636 1.598 1.574 1.601 1.619 1.600 0.004 1.104 1.097 1.092 1.106 1.103 1.107 1.103 1.097 0.2 1.409 1.411 1.407 1.395 1.410 1.426 1.410 1.414 2/5 1.688 1.693 1.765 1.689 1.715 1.701 1.691 1.694 2/3 1.999 2.012 2.026 2.043 2.008 2.014 2.028 2.028 0.004 1.378 1.359 1.383 1.383 1.378 1.373 1.371 1.363 0.2 1.803 1.786 1.816 1.800 1.814 1.813 1.818 1.813 2/5 2.167 2.112 2.189 2.165 2.172 2.178 2.201 2.183 2/3 2.577 2.575 2.628 2.576 2.561 2.605 2.601 2.585 0.004 1.757 1.754 1.769 1.768 1.743 1.778 1.755 1.767 0.2 2.337 2.344 2.368 2.359 2.390 2.375 2.375 2.387 Continued on next page 86 Table 6.6 – concluded from previous page ρ∗ 0.8118 0.8404 0.8595 0.9168 0.9359 0.9512 α 32 100 250 490 850 1360 2900 5300 2/5 2.845 2.920 2.879 2.879 2.883 2.870 2.879 2.878 2/3 3.387 3.434 3.434 3.445 3.419 3.447 3.415 3.427 0.004 2.297 2.309 2.327 2.337 2.325 2.340 2.330 2.313 0.2 3.117 3.174 3.251 3.262 3.226 3.221 3.229 3.237 2/5 3.781 3.874 3.973 3.952 3.961 3.964 3.940 3.986 2/3 4.518 4.668 4.718 4.746 4.730 4.718 4.756 4.752 0.004 2.729 2.759 2.817 2.825 2.811 2.824 2.810 2.847 0.2 3.732 3.856 3.996 3.965 3.982 4.013 4.012 4.012 2/5 4.602 4.811 4.927 4.972 4.919 4.910 4.938 4.963 2/3 5.511 5.779 5.915 5.869 5.902 5.887 5.891 5.832 0.004 3.101 3.139 3.238 3.283 3.245 3.228 3.243 3.222 0.2 4.298 4.547 4.500 4.678 4.638 4.623 4.626 4.630 2/5 5.304 5.606 5.804 5.784 5.772 5.775 5.700 5.789 2/3 6.356 6.766 6.921 6.887 6.891 6.900 6.885 6.907 0.004 4.970 5.200 5.290 5.332 5.279 5.285 5.275 5.274 0.2 7.318 8.056 8.075 8.076 8.070 7.937 7.937 8.054 2/5 9.024 10.24 10.22 10.10 10.25 10.11 10.20 10.04 2/3 10.87 12.37 12.06 12.25 12.10 12.18 12.14 12.20 0.004 5.527 6.742 6.468 6.444 6.490 6.559 6.586 0.2 8.000 10.59 10.18 10.13 10.19 10.28 10.30 10.20 2/5 9.583 13.43 12.78 12.85 12.66 12.77 12.84 12.77 2/3 11.96 16.09 15.76 15.71 15.49 15.96 15.75 15.97 0.004 4.851 8.841 7.805 7.708 7.657 7.660 7.743 0.2 6.709 14.05 12.27 12.11 12.02 11.10 12.29 12.35 2/5 8.179 18.05 15.67 15.34 2/3 9.823 21.82 19.29 18.57 18.55 22.42 19.44 19.69 87 Table 6.7: Reduced thermal conductivity coefficients λ∗ for a range of reduced densities ρ∗ , dimensionless radii of gyration α, and number of particles N spanning from 32 to 5300. ρ∗ α 32 100 250 490 850 1360 2900 5300 0.0010 0.0100 0.004 0.2865 0.3343 0.3820 0.4298 1.049 1.049 1.082 1.150 0.2 1.028 1.065 1.076 1.081 1.080 2/5 1.088 1.097 1.098 1.093 1.083 1.100 1.109 1.103 1.073 1.067 2/3 1.108 0.004 1.042 1.053 1.057 1.059 1.060 1.062 1.059 1.077 1.088 1.098 1.089 1.088 1.089 1.108 2/5 1.095 1.112 1.117 1.128 1.114 1.115 1.129 2/3 1.102 1.119 1.121 1.103 1.092 1.123 0.2 0.1910 1.040 1.035 0.004 1.447 1.525 1.513 1.536 1.562 1.577 1.563 1.600 0.2 1.462 1.535 1.565 1.564 1.587 1.581 1.589 1.595 2/5 1.465 1.536 1.575 1.586 1.579 1.590 1.601 1.594 2/3 1.454 1.540 1.563 1.512 1.578 1.593 1.594 1.581 0.004 1.930 2.026 2.103 2.068 2.059 2.064 2.088 2.087 0.2 1.915 2.006 2.040 2.043 2.062 2.064 2.070 2.086 2/5 1.898 2.008 2.039 2.062 2.065 2.070 2.073 2.073 2/3 1.892 1.994 2.057 2.046 2.049 2.055 2.070 2.059 0.004 2.259 2.355 2.404 2.426 2.414 2.425 2.453 2.441 0.2 2.225 2.370 2.391 2.380 2.380 2.382 2.382 2.372 2/5 2.203 2.268 2.346 2.368 2.386 2.387 2.384 2.383 2/3 2.181 2.302 2.316 2.345 2.368 2.367 2.369 2.382 0.004 2.668 2.752 2.809 2.862 2.845 2.858 2.866 2.860 0.2 2.605 2.731 2.769 2.772 2.791 2.786 2.823 2.813 2/5 2.572 2.682 2.727 2.733 2.794 2.787 2.794 2.824 2/3 2.544 2.671 2.698 2.716 2.756 2.750 2.734 2.749 0.004 3.172 3.291 3.338 3.354 3.372 3.375 3.387 3.388 Continued on next page 88 Table 6.7 – continued from previous page ρ∗ 0.4775 0.5253 0.5730 0.6208 0.6685 0.7163 0.7640 α 32 100 250 490 850 1360 2900 5300 0.2 3.077 3.170 3.246 3.251 3.271 3.275 3.288 3.263 2/5 3.020 3.156 3.205 3.245 3.216 3.222 3.232 3.183 2/3 2.998 3.133 3.181 3.194 3.207 3.211 3.212 3.233 0.004 3.796 3.917 4.089 3.892 3.987 3.978 4.068 4.072 0.2 3.664 3.787 3.842 3.867 3.866 3.903 3.887 3.851 2/5 3.578 3.712 3.789 3.794 3.801 3.878 3.781 3.833 2/3 3.525 3.690 3.663 3.907 3.806 3.781 3.798 3.796 0.004 4.570 4.696 4.763 4.764 4.803 4.793 4.810 4.856 0.2 4.363 4.497 4.539 4.575 4.594 4.588 4.625 4.612 2/5 4.227 4.389 4.454 4.534 4.492 4.509 4.526 4.520 2/3 4.195 4.376 4.402 4.488 4.453 4.475 4.478 4.481 0.004 5.500 5.649 5.719 5.713 5.752 5.793 5.741 5.756 0.2 5.218 5.356 5.454 5.475 5.535 5.446 2/5 5.046 5.268 5.031 5.308 5.325 5.389 5.386 5.358 2/3 5.003 5.175 5.246 5.306 5.287 5.313 5.322 5.225 0.004 6.633 6.801 6.884 6.964 6.851 6.872 6.950 6.945 0.2 6.263 6.406 6.506 6.642 6.566 6.562 6.545 6.581 2/5 6.073 6.254 6.125 6.383 6.380 6.372 6.404 6.388 2/3 5.978 6.192 6.253 6.381 6.346 6.352 6.362 6.371 0.004 8.027 8.158 8.225 8.306 8.363 8.390 8.253 8.341 0.2 7.539 7.697 7.782 7.883 7.907 7.915 7.894 7.910 2/5 7.240 7.460 7.613 7.494 7.577 7.677 7.651 7.675 2/3 7.163 7.357 7.496 7.402 7.613 7.584 7.624 7.644 0.004 9.764 9.958 10.06 9.920 10.11 10.10 10.09 10.16 0.2 9.107 9.310 9.437 9.412 9.511 9.512 9.497 9.316 2/5 8.781 9.378 9.193 9.111 9.196 9.103 9.260 9.223 2/3 8.621 8.780 8.935 9.076 9.100 9.036 9.124 9.123 0.004 11.84 11.99 12.18 12.28 12.25 12.27 12.26 12.42 0.2 10.98 11.25 11.29 11.38 11.47 11.48 11.51 11.52 5.515 Continued on next page 89 Table 6.7 – concluded from previous page ρ∗ 0.8118 0.8404 0.8595 0.9168 0.9359 0.9512 α 32 100 250 490 850 1360 2900 5300 2/5 10.58 10.84 10.99 11.12 10.99 11.06 11.09 11.08 2/3 10.40 10.75 10.85 10.92 10.94 11.02 11.05 11.04 0.004 14.35 14.63 14.42 14.84 14.96 14.95 14.95 14.96 0.2 13.30 13.60 13.75 13.85 13.88 13.81 13.92 13.87 2/5 12.69 13.12 13.33 13.30 13.49 13.39 13.45 13.54 2/3 12.54 12.93 13.10 13.24 13.29 13.31 13.05 13.43 0.004 16.07 16.31 16.72 16.52 16.75 16.64 16.81 16.64 0.2 14.77 15.29 15.46 16.49 15.66 15.60 15.61 15.63 2/5 14.22 14.65 14.97 14.89 14.99 15.08 15.11 15.11 2/3 14.04 14.42 14.78 14.97 14.82 14.88 15.02 14.88 0.004 17.29 17.89 18.16 18.28 18.15 18.25 18.36 18.25 0.2 15.86 16.47 16.70 16.95 17.01 17.06 16.94 17.01 2/5 15.26 15.82 16.13 16.28 16.27 16.28 16.26 16.33 2/3 15.01 15.70 15.40 16.13 15.94 16.11 16.30 16.04 0.004 21.24 22.54 23.10 24.35 23.25 23.32 23.37 23.45 0.2 19.38 20.43 21.33 21.35 21.34 21.49 21.33 21.19 2/5 18.65 19.90 20.38 19.87 20.72 21.02 20.58 20.55 2/3 18.41 19.76 20.28 20.45 20.55 20.46 20.46 20.66 0.004 21.42 24.43 24.89 25.23 25.42 25.24 25.16 0.2 19.80 22.43 23.00 22.76 22.83 23.38 23.03 23.12 2/5 18.98 21.82 22.31 22.33 22.06 22.47 22.31 22.12 2/3 18.83 21.42 21.94 22.06 22.19 22.12 22.17 22.11 0.004 19.89 26.03 26.80 26.85 26.83 27.02 26.79 0.2 19.32 24.07 24.54 24.79 24.71 23.59 24.70 24.90 2/5 18.17 22.99 23.51 23.60 2/3 17.99 22.90 23.44 23.31 23.46 22.42 23.54 23.86 90 Table 6.8: Reduced self diffusion coefficients D∗ for a range of reduced densities ρ∗ , dimensionless radii of gyration α and number of particles N spanning from 32 to 5300. ρ∗ α 32 100 250 490 850 1360 2900 5300 0.0010 0.004 0.2 0.0100 0.1910 0.2865 0.3343 0.3820 0.4298 213.93 214.12 176.69 177.61 177.86 178.05 178.06 2/5 158.58 159.49 159.78 159.90 159.96 2/3 145.56 146.61 146.69 146.74 146.83 0.004 20.979 21.082 21.137 21.161 21.169 21.180 21.185 17.419 17.525 17.571 17.586 17.600 17.610 17.620 2/5 15.637 15.746 15.782 15.802 15.813 15.824 15.828 2/3 14.336 14.437 14.473 14.495 14.503 14.514 14.517 0.2 173.52 213.85 17.053 178.17 178.15 0.004 0.8435 0.8855 0.9048 0.9137 0.9200 0.9235 0.9281 0.9314 0.2 0.6918 0.7314 0.7495 0.7587 0.7641 0.7684 0.7729 0.7758 2/5 0.6152 0.6532 0.6706 0.6799 0.6847 0.6885 0.6931 0.6960 2/3 0.5589 0.5949 0.6115 0.6196 0.6247 0.6282 0.6325 0.6350 0.004 0.4983 0.5312 0.5471 0.5556 0.5608 0.5645 0.5694 0.5724 0.2 0.4046 0.4353 0.4506 0.4587 0.4637 0.4672 0.4718 0.4748 2/5 0.3580 0.3865 0.4009 0.4080 0.4131 0.4162 0.4205 0.4233 2/3 0.3232 0.3496 0.3627 0.3695 0.3739 0.3767 0.3808 0.3833 0.004 0.4005 0.4298 0.4447 0.4529 0.4576 0.4613 0.4661 0.4691 0.2 0.3231 0.3503 0.3644 0.3717 0.3765 0.3797 0.3843 0.3869 2/5 0.2845 0.3097 0.3226 0.3293 0.3338 0.3368 0.3410 0.3435 2/3 0.2560 0.2789 0.2908 0.2968 0.3007 0.3035 0.3073 0.3095 0.004 0.3269 0.3535 0.3673 0.3746 0.3796 0.3829 0.3874 0.3900 0.2 0.2620 0.2861 0.2988 0.3054 0.3097 0.3131 0.3171 0.3197 2/5 0.2296 0.2516 0.2631 0.2694 0.2733 0.2761 0.2798 0.2821 2/3 0.2058 0.2258 0.2360 0.2414 0.2450 0.2475 0.2509 0.2530 0.004 0.2697 0.2932 0.3059 0.3127 0.3171 0.3205 0.3244 0.3271 Continued on next page 91 Table 6.8 – continued from previous page ρ∗ 0.4775 0.5253 0.5730 0.6208 0.6685 0.7163 0.7640 α 32 100 250 490 850 1360 2900 5300 0.2 0.2140 0.2353 0.2463 0.2525 0.2565 0.2594 0.2632 0.2655 2/5 0.1866 0.2057 0.2159 0.2214 0.2250 0.2275 0.2308 0.2330 2/3 0.1665 0.1838 0.1928 0.1977 0.2008 0.2030 0.2060 0.2079 0.004 0.2232 0.2445 0.2558 0.2620 0.2660 0.2688 0.2729 0.2752 0.2 0.1753 0.1940 0.2039 0.2092 0.2131 0.2154 0.2189 0.2210 2/5 0.1521 0.1687 0.1775 0.1824 0.1854 0.1878 0.1908 0.1926 2/3 0.1352 0.1500 0.1578 0.1621 0.1648 0.1667 0.1694 0.1710 0.004 0.1850 0.2036 0.2138 0.2193 0.2229 0.2256 0.2291 0.2312 0.2 0.1434 0.1597 0.1683 0.1730 0.1762 0.1785 0.1815 0.1833 2/5 0.1235 0.1380 0.1456 0.1497 0.1525 0.1544 0.1571 0.1587 2/3 0.1094 0.1221 0.1288 0.1324 0.1348 0.1365 0.1387 0.1402 0.004 0.1527 0.1692 0.1780 0.1829 0.1862 0.1885 0.1916 0.1935 0.2 0.1166 0.1307 0.1381 0.1423 0.1449 0.1470 0.1492 0.1512 2/5 0.09975 0.1122 0.1187 0.1222 0.1245 0.1262 0.1285 0.1299 2/3 0.08790 0.09887 0.1045 0.1076 0.1096 0.1110 0.1130 0.1142 0.004 0.1252 0.1395 0.1472 0.1513 0.1542 0.1562 0.1589 0.1606 0.2 0.09387 0.1060 0.1123 0.1157 0.1181 0.1197 0.1219 0.1233 2/5 0.07979 0.09028 0.09570 0.09874 0.1007 0.1021 0.1039 0.1052 2/3 0.06992 0.07913 0.08386 0.08647 0.08816 0.08937 0.09094 0.09199 0.004 0.1013 0.1138 0.1202 0.1239 0.1262 0.1280 0.1303 0.1318 0.2 0.07455 0.08483 0.09008 0.09291 0.09487 0.09625 0.09816 0.09926 2/5 0.06280 0.07167 0.07615 0.07862 0.08020 0.08141 0.08295 0.08395 2/3 0.05476 0.06252 0.06641 0.06847 0.06987 0.07092 0.07216 0.07307 0.004 0.08053 0.09131 0.09668 0.09971 0.1017 0.1031 0.1051 0.1062 0.2 0.05796 0.06663 0.07087 0.07332 0.07479 0.07590 0.07741 0.07836 2/5 0.04840 0.05579 0.05941 0.06136 0.06268 0.06359 0.06476 0.06567 2/3 0.04195 0.04837 0.05152 0.05316 0.05430 0.05508 0.05612 0.05676 0.004 0.06237 0.07166 0.07609 0.07855 0.08015 0.08130 0.08292 0.08389 0.2 0.04376 0.05102 0.05443 0.05625 0.05751 0.05836 0.05953 0.06029 Continued on next page 92 Table 6.8 – concluded from previous page ρ∗ 0.8118 0.8404 0.8595 0.9168 0.9359 0.9512 α 32 100 250 490 850 1360 2900 5300 2/5 0.03614 0.04226 0.04514 0.04665 0.04769 0.04841 0.04936 0.04995 2/3 0.03110 0.03644 0.03894 0.04022 0.04106 0.04169 0.04250 0.04299 0.004 0.04629 0.05450 0.05804 0.05995 0.06119 0.06206 0.06332 0.06408 0.2 0.03150 0.03770 0.04031 0.04171 0.04263 0.04330 0.04414 0.04470 2/5 0.02572 0.03088 0.03306 0.03421 0.03496 0.03550 0.03623 0.03667 2/3 0.02199 0.02644 0.02831 0.02932 0.02993 0.03039 0.03099 0.03137 0.004 0.03764 0.04532 0.04830 0.04997 0.05101 0.05177 0.05276 0.05346 0.2 0.02513 0.03073 0.03289 0.03409 0.03484 0.03539 0.03611 0.03655 2/5 0.02034 0.02498 0.02679 0.02775 0.02836 0.02881 0.02940 0.02977 2/3 0.01732 0.02130 0.02284 0.02367 0.02419 0.02454 0.02503 0.02536 0.004 0.03232 0.03964 0.04232 0.04374 0.04469 0.04534 0.04622 0.04678 0.2 0.02125 0.02647 0.02840 0.02940 0.03007 0.03053 0.03115 0.03155 2/5 0.01710 0.02140 0.02300 0.02383 0.02437 0.02476 0.02525 0.02556 2/3 0.01452 0.01820 0.01956 0.02026 0.02072 0.02104 0.02144 0.02173 0.004 .01783 .02460 .02637 .02737 .02793 .02835 .02890 .02924 0.2 .01113 .01558 .01678 .01745 .01788 .01813 .01849 .01873 2/5 .008777 .01238 .01335 .01390 .01423 .01444 .01475 .01491 2/3 .007396 .01041 .01124 .01171 .01197 .01215 .01241 .01256 0.004 .01120 .02023 .02175 .02311 .02344 .02391 .02420 0.2 .006727 .01257 .01357 .01415 .01447 .01470 .01501 .01517 2/5 .005317 .009911 .01072 .01118 .01144 .01163 .01187 .01201 2/3 .004428 .008312 .008991 .009380 .009600 .009760 .009960 .01008 0.004 .004498 .01700 .01830 .01943 .01913 .02019 .02043 0.2 .002738 .01036 .01119 .01149 .01159 .01002 .01244 .01259 2/5 .002034 .008130 .008794 .008937 2/3 .001726 .006793 .007345 .007471 .007756 .006062 .008198 .008298 93 Chapter 7 Applications 7.1 Motivation This chapter will illustrate how the theoretical and computational findings presented in the previous chapters can be applied to analysing real experimental data. The main question being asked is whether it is possible to establish simple criteria based by which the fluid regime can be identified and described with appropriate theory. For example, can the Knudsen number Kn, which is usually used to describe rarefied fluids, serve as an identifier for the nature of diffusion in any fluid? Are there other parameters that play roles in predicting diffusion dynamics? What are the ranges in which these parameters give best prediction? It is particularly interesting to look at the diffusion coefficient for the following reason. As described in Chapter 1, the diffusion coefficient is related to ion mobility in the weak electric field limit and is often sought in experiments with aerosols, colloids and with electrophoresis techniques. The ion mobility can be measured quite easily, but it is difficult to interpret the size of the particles correctly. The drag force acting on a large spherical particle of radius R moving through a fluid of viscosity η with a constant relative speed V is given by the Stokes equation [29, 68] F0 = −6πηRV . (7.1) In this equation, the coefficient 6, correponding to “stick” boundary conditions, is typically used. The drag coefficient kD = − F0 /V is related to the diffusion coefficient via expression D= kT . kD (7.2) In the hydrodynamic limit, the diffusion coefficient is given by Eq. (1.2). The Cunningham correction factor was introduced to account for deviation from the hydrodynamic diffusion coefficient, caused by drag force acting on smaller particles when the surrounding fluid is not continuous. The Cunningham correction factor is given by C(Kn) = 1 + A Kn, (7.3) 94 in which A is a parameter [69]. With the Cunningham correction factor, the expression for the drag force becomes F =− 6πηRV . C (7.4) The form of Eq. (7.3) was chosen according to the following considerations. In the hydrodynamic limit Kn 1 making the term A Kn negligible and C(Kn) = 1. Equation (7.4) then takes the form of Eq. (7.1), that has a linear dependence upon R, as expected for the hydrodynamic limit. For large Kn, in the molecular limit, the term A Kn dominates and leads to R2 dependence in Eq. (7.1) that is expected in the molecular regime. Thus Eq. (7.3) approaches theoretically predicted expressions for both limiting cases, regardless of the way the system approaches the limits. Modifying the coefficient A does not change the two limits even though it changes the way the system moves between them. Knudsen and Weber [70] have noticed that A is constant only for Kn 1 and suggested another correction, that would include all possible Kn values. With their suggestions expression (7.1) takes form F =− 6πµRV , 1 + Kn [A + B exp (−E/Kn)] (7.5) where A, B and E are constants. Equation (7.5) is known as the Stokes-Cunningham equation. Similarly, the form of this expression was chosen such that F has R dependence in the hydrodynamic limit and R2 in the molecular limit. Combined with Eq. (1.2) Stokes-Cunningham equation gives the following expression for the diffusion coefficient D= kB T [1 + Kn (A + B exp (−E/Kn))] . 6πηR (7.6) The constants A, B and E can be determined by fitting experimental data. One of the first sets of diffusion data was obtained by Millikan in experiments done with oil droplets. The Stokes-Cunningham equation fits the Millikan data very well and parameters obtained from the Millikan data are often used to analyze various other experimental data sets [70–78]. Based on the background information above the following specific questions arise, that are testable: (i) Can the diffusion expression given by Eq. (5.4) suggested by R. Sokolowskii et al. [23] be a good predictor for the diffusion coefficient for a broad range of Kn? (ii) Can Eq. (7.6) be used for obtaining the diffusion coefficient and what is its’ range of applicability? (iii) What are the boundary conditions to be used in the hydrodynamic formula for various types of fluids? In order to answer these questions several experimental data sets obtained from other published sources have been analysed and a series of plots produced as a function of Kn for comparison with the computational results for the rough and smooth hard sphere fluids. 95 7.2 Data Analysis A number of studies where diffusion or particle size are obtained at various Kn values have been done, but a data set that covers a wide range of nanoscale particle sizes, where various values of Kn are achieved by changing the size of the particle relative to the size of the bath, while keeping the density of the bath constant, is difficult to find. The data sets chosen for this analysis cover a range of tracer sizes at low density. The experimental particle data used in this analysis are presented in Table 7.1. Several studies have been done in which particle sizes determined by a differential mobility analyser(DMA) and measured by transmission electronic microscopy (TEM) have been compared. The tracer particles used in these experiments are Ag [79, 80] and Cu2 O [81] in a nitrogen bath. In all cases the tracer particles are nearly spherical with the smallest radius of approximately 1.3 nm. The radii data used in this analysis come from the TEM measurement due to fact that the radius obtained via TEM is more precisely defined and therefore better suited for our purpose of comparing with the hard sphere model. Based on experimental conditions, the reduced density ρ∗ in all three experiments is in the range of 10−5 − 10−4 . In the first column of Table 7.1 the TEM radii are converted into relative units used throughout this thesis by dividing the measured value by σ, the diameter of the nitrogen bath particles, that has the value of 0.365 nm [16]. The second column provides value for the corresponding Knudsen number, calculated according to Eq. (3.1), with the mean free √ −1 path λ = 2πσ 2 ρ∗ , where σ us the diameter of the nitrogen bath and ρ∗ is the reduced density of the bath. The third column represents the distance between the centers of the tracer and bath particles upon closest approach, i.e. R∗ = (Σ + σ)/2σ, where Σ = RT EM . The last column shows the scaled value of the diffusion coefficient D∗ . Consider experimental values of the diffusion coefficient as a function of the distance between the centres of the particles, as shown in Fig. 7.1. In the upper panel the values D∗ R∗ are shown. If the fluid were in the hydrodynamic regime where D∗ ∼ 1/R∗ , the plot would stay constant as R∗ increases. Instead, decaying curves are observed, that is consistent with D∗ ∼ 1/R∗ 2 . In the lower panel, the quantity D∗ R∗ 2 is plotted. If the fluid were described by the Boltzmann equation the curve would be expected to fluctuate around a constant value, while it should increase linearly if the hydrodynamic regime is in place. Unfortunately, it is difficult to obtain a clear trend in the representation shown on the lower panel due to data scatter, but it can be safely presumed that there is no clear increase. Figure 7.2 was created to compare fit done with Eq. (5.4), suggested by Sokolovskii et. al. [23] and Eq. (7.6) to the Millikan data. The Millikan data is represented by open circles. On the vertical axis the quantity A+B exp(−E/Kn) = (D/DH −1)/Kn is shown. The reasons for representing Millikan’s data in this form are the following. The ratio D/DH shows how much the actual diffusion coefficient deviates from the diffusion coefficient predicted for the hydrodynamic regime. Subtracting 1 from this ratio removes the integer part and the 96 Table 7.1: Experimental data for particle sizes measured by TEM with corresponding diffusion coefficients. Numbers in square parentheses are references to published sources. Ag [79] RT EM /σ Kn R∗ D∗ 30.5 6.15 31.0 0.073 25.9 7.25 26.4 0.111 19.7 9.52 20.2 0.225 18.8 10.0 19.3 0.204 18.2 10.3 18.7 0.142 17.7 10.6 18.2 0.204 16.0 11.7 16.5 0.304 13.2 14.7 0.304 14.2 13.6 13.8 14.1 0.455 13.6 13.8 14.1 0.283 13.5 13.9 14.0 0.379 11.7 16.1 12.2 0.455 10.9 17.3 11.4 0.742 10.8 17.4 11.3 1.121 9.37 20.0 9.87 0.742 20.8 9.54 0.568 9.04 8.43 22.5 8.94 0.742 7.42 25.3 7.92 1.121 7.17 26.2 7.67 1.121 5.29 35.5 5.79 2.197 4.55 41.3 5.05 2.252 Ag [80] 29.4 6.38 29.9 0.0907 8.52 22.5 0.135 22.0 21.3 8.81 21.8 0.202 16.8 11.2 17.3 0.300 10.3 18.2 10.8 0.455 10.0 18.7 10.5 0.677 8.13 23.1 8.63 1.001 Cu2 O [81] 12.2 15.3 12.7 0.416 18.2 10.8 0.623 10.3 8.65 21.7 9.15 0.934 8.00 23.5 8.50 1.141 7.34 25.6 7.84 1.399 6.57 28.6 7.07 1.706 6.19 30.3 6.69 2.093 5.31 35.3 5.81 2.560 4.57 41.0 5.07 3.142 4.03 46.6 4.53 3.821 53.5 4.00 4.690 3.51 97 20 16 DR * * 12 8 4 0 140 DR * *2 120 100 80 60 40 0 4 8 12 16 * R 20 24 28 32 Figure 7.1: Experimental values for the reduced diffusion coefficient multiplied by R∗ in 2 the upper panel and by R∗ in the lower panel as a function of the distance between the centers of the particles, drawn for two different representation of the particle’s size. Circles, squares and crosses represent Ag [79], Ag [80] and Cu2 O [81], respectively. remaining expression is then expected to behave linearly with Kn. Dividing by Kn removes this linear dependence which may overpower the trend and allows one to see a more detailed picture of the variation in the data. The solid line is produced by using the Stokes-Cunningham formula with coefficients estimated by Allen and Raabe [77]. The formula fits the Millikan data well, as expected. The dashed and dashed-dotted lines represents the fit given by Eq. (5.4), that combines both molecular and hydrodynamic diffusion coefficients. The dashed dotted line represents 98 the fit by Eq. (5.4) with χ = 1.4 obtained by varying the ratio DE /DH , as suggested by Sokolovskii et. al. [23] for smooth hard sphere fluid [23]. For the dashed line, χ was treated as a fit variable, and was determined to have the value of approximately 1.1. As can be seen, the formula fits the Millikan data quite well in the case of χ = 1.1, though deviates for χ = 1.4. Thus, Eq. (5.4) can fit a wider range of Kn. 2 1.9 A+B exp(-E/Kn) 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1 0 10 20 30 40 50 60 70 80 90 Kn Figure 7.2: Millikan diffusion data and the fitted curves plotted as a function of the Knudsen number Kn. The circles represent the values derived from the raw Millikan data, the solid line represents Eq. (7.6) with parameters obtained by Allen and Raabe [77], the dashed line represents Eq. (5.4) with x = 1.1 and the dashed-dotted line represents Eq. (5.4) with x = 1.4. Figure 7.3 shows computational results for the rough hard sphere fluid at different values of α, as well as the smooth hard sphere fluid obtained by Sokolovskii [23] and Heyes [61] as a function of Kn in order to estimate the value for the boundary conditions. The hydrodynamic equation is intentionally applied outside the range of theoretical applicability, in order to observe, whether it will predict diffusion correctly, and if so, to estimate the the values for the boundary conditions. By plotting the quantity 1/(πD∗ η ∗ ), one can determine which values for the boundary conditions are adopted at different Kn, if the hydrodynamic description is applicable. To understand what the expected behaviour might be, let’s consider two limiting cases. In the hydrodynamic regime, when Kn 1, Eq. (7.6) transforms into the Stokes-Einstein formula given by Eq. (1.2). In reduced units and after simple transformation this expression 99 can be written as 1 D∗ η∗ π = R∗ ζ. (7.7) In the case of self-diffusion, when the size of the bath particle is taken into account, R∗ = 1. When the size of the bath particle is not considered, R∗ = 1/2. Plotting the left side of Eq. (7.7) allows one to estimate the value of ζ in the hydrodynamic regime. There is no Kn dependence in this expression, therefore, if the hydrodynamic regime is reached, the value for ζ will remain unchanged and the plot will result in a straight line. On the other hand, when Kn is large, the fluid will be in the Boltzmann regime, where diffusion for the rough hard sphere fluid is described by Eq. (3.12). Rewriting this expression in reduced units and transforming it in a similar manner to Eq. (7.7) leads to 1 D∗ η∗ π = 64 √ 15 2πR∗ 1 (1 + 2α)(6 + 13α) , Kn (1 + α)3 (7.8) where R∗ = 1 when the size of the bath is taken into account, and R∗ = 1/2, when it’s not. Note that for the smooth hard sphere fluid Eq. (7.8) will be the same except for the part in square brackets, which for smooth spheres is equal to 1 due to α = 0. From this expression it follows that as Kn increases, 1/(D∗ η ∗ π) should decrease. As seen in Fig. 7.3 the values of 1/(D∗ η ∗ π) decrease as Kn increases, but reache somewhat stable values of 2.1 ± 0.1 near Kn ≈ 0.6. If Eq. (7.7) is applied for Kn ≈ 0.6 then the values of ζ is either 2.1 if R∗ = 1 or 4.2 if R∗ = 1/2. There is some variation in this value with α, although no simple trends are observed. This suggests that if one is to use the hydrodynamic equation near this range of Kn, the “slip” value for the boundary conditions may be more appropriate in cases when the size of the bath particles is ignored. The purpose of Figs. 7.4 and 7.5 is to illustrate that defining the representative size of the particle in different ways may affect the interpretation of results. In Fig. 7.4 the deviation from the hydrodynamic behaviour is shown. The difference between the top and bottom plots is in the definition of the radius of the particles in the expressions for the diffusion coefficient. In the top panel the radius is defined as the diameter of the tracer particle divided by 2, that is R∗ = Σ/2, i.e. the size of the bath particles is not taken into account. In the bottom panel, the distance between the particles upon closest approach is considered, and thus the radius of the bath is taken into account. With the density unchanged, for larger Kn and smaller particles the fluid is expected to be described by the Boltzmann theory, i.e. the diffusion coefficient is proportional to the inverse of the R2 , and thus the ratio D/Dh , where D = DB should be linear. This is observed for the lowest density in the bottom panel, where the size of the bath is included, but not for the top panel, where only the radius of the tracer is considered. The curves in the inserts on both panels represent higher densities. If the hydrodynamic regime were 100 2.5 2.4 2.3 ∗ 2.1 * 1/(πD η ) 2.2 2 1.9 1.8 α = 2/3 α = 2/5 α = 0.2 α = 0.004 smooth 1.7 1.6 1.5 1.4 0 0.5 1 1.5 2 2.5 3 Kn Figure 7.3: Values 1/πD∗ η ∗ for the rough hard sphere data given in Chapter 5 and smooth hard sphere data given by Sokolovskii at. al. [23] and Heyes [61] are plotted as a function of Knudsen number Kn. in place, the curves would have been constant and originate at 1. In both cases they do originate at 1 but deviate quickly from the hydrodynamic behaviour. Figure 7.5, again, shows the deviation from the hydrodynamic behavior while employing a more sensitive functional form, as in Fig. 7.2. The dotted line in the panels represents the analytic fit using Eq. (7.6) with Allen and Raabe coefficients [77]. When the size of the bath is taken into account, the analytic fit line matches the low density smooth hard sphere data almost exactly, although the analytic parameters were obtained for a different density range. On both Figs. 7.4 and 7.5 a density dependence is observed. The data obtained at higher densities do not fall directly onto the lowest density curve. It means that Eq. (7.6) performs well only at low density, i.e. the parameters A, B and E are density-dependent and Eq. (7.6) will break down at higher densities if used with parameters obtained at lower densities. Therefore, additional sets of these parameters should be obtained for higher densities. 101 12 * ρ = 0.02 10 * ρ = 0.18 D / Dh 8 * ρ = 0.34 6 2 1.5 1 0.5 0 4 2 0 * ρ = 0.34, α = 2/5 * ρ = 0.52, α = 2/5 * ρ = 0.52, α = 2/3 0 1 2 3 4 5 6 4 50 3 2 40 D / Dh * ρ = 0.52 1 30 0 0 0.4 0.8 1.2 20 10 0 0 5 10 15 20 25 Kn Figure 7.4: Deviation from the hydrodynamic diffusion coefficient as a function of Kn at various densities, represented by D/Dh , for the rough and smooth hard sphere particles at various densities and degree of rotational-translational coupling. In the top panel the radius of the bath particle is not taken into account, while in the bottom panel the radius of the bath particle is accounted for. 7.3 Summary Several sets of experimental data have been analysed and compared with the computational sets of rough and smooth hard sphere data in order to determine the applicability of the hydrodynamic or molecular expressions at different values Kn. It has been found, that in the case of the nanoparticles Ag and Cu2 O in a nitrogen bath, hydrodynamic behaviour does not manifest itself, and therefore molecular (Boltzmann) expressions should be employed 102 3 * ρ = 0.02 analytic fit (D / Dh-1)/Kn 2.5 * 2 ρ = 0.34 1.5 ρ = 0.52 * * ρ = 0.34, α = 2/5 1 * ρ = 0.52, α = 2/5 * ρ = 0.52, α = 2/3 0.5 (D / Dh-1)/Kn 06 3 5 2 4 1 3 0 0 0.2 0.4 0.6 0.8 4 6 1 1.2 2 1 0 0 2 8 10 12 14 16 18 20 Kn Figure 7.5: Deviation from the hydrodynamic diffusion coefficient as a function of Kn at various densities, represented by (D/Dh − 1) /Kn , for the rough and smooth hard sphere particles at variuos densities and degree of rotational-translational coupling. In the top panel the radius of the bath particle is not taken into account, while in the bottom panel the radius of the bath particle is accounted for. for interpreting the radius of the particles. Comparing the applicability of Eq. (5.4) and (7.6) showed that both expressions perform well in the hydrodynamic regime. It was already known that Eq. (5.4) can be applied in the molecular regime, and therefore the larger range of applicability of this formula is established. Various ways of representing the diffusion coefficient have been explored. It is found, that the ratio D/Dh is less sensitive towards deviation from the hydrodynamic diffusion and 103 therefore has to be used more carefully. Instead, the expression (D/DH − 1)/Kn can be employed for more detailed analysis. It has been found that for hydrodynamic behaviour the “slip” value for the surface boundary conditions near Kn = 0.6 is more appropriate rather than the “stick” value that is commonly used. The difference between the definition of the representative size of the particle was explored. It was found that including the size of the bath particles in the expression for the diffusion coefficient produces a trend that corresponds to expected behaviour for the molecular and hydrodynamic cases, while excluding the size of the bath does not show the correct behaviour for these limits. A density dependence has been observed at higher densities and therefore Kn can not be the only identifier for the behaviour of the fluid at higher densities. Overall, more experimental data on diffusion as a function of particles size at different densities is needed to study a wider range of Kn. 104 Chapter 8 Summary 8.1 Computational Results In this thesis the results of a systematic study of the transport properties of the rough hard sphere fluid are presented. Molecular dynamics was employed to perform simulations and obtain transport coefficients. Diffusion coefficients were studied as a function of particle size at different densities and with different degrees of translational-rotational coupling. Self-diffusion, shear and bulk viscosity and thermal conductivity were studied over a wide range of densities and degrees of rotational-translational coupling. Even though the rough and smooth hard sphere models don’t account for all physical effects in real molecules, their simplicity makes them powerful tools for the qualitative study of the properties of the real fluids. The smooth hard sphere model does not presume rotation of the individual molecules and therefore no exchange of rotational energy and momenta, while the rough hard sphere fluid differs from the smooth hard sphere fluid by the presence of rotation. The smooth hard sphere particles are characterized by “slip” boundary conditions at the surface in the hydrodynamic limit. Comparing the two models allows one to isolate and test the physical effect found in real fluids, that is the ability to collide inelastically. Introducing this effect alone changes boundary conditions from “slip” to “stick” with a range of values in between. In Chapter 5, results of tracer diffusion in the rough hard sphere fluid are reported. It is found that the rough hard sphere fluid can adopt more than two limiting values for hydrodynamic boundary conditions, spanning between “slip” and “stick” values, depending on the degree of translational - rotational coupling between the particles during collisions. The value of the boundary condition is found to have a linear dependence upon α, a parameter that characterizes translational and rotational energy exchange. It is commonly believed that particles adopt “stick” boundary condition upon collision due to the presence of an attractive potential. When an attractive potential is present, the colliding particle is “trapped” in the attractive well of the other particle. After it thermalizes it leaves with random velocity, consistent with the Boltzmann distribution, making such collisions so-called “diffuse”. The result of this work shows that in fact an attractive potential is not necessary to reach the “stick” boundary condition. The “stick” boundary condition may also be reached by simply “tuning” the degree of rotational - translational coupling. The absence of an attractive potential in the current work keeps the structure of the fluid unchanged and 105 allows one to test the mere effect of rotational-translational coupling upon boundary conditions on the surface of the tracer. It has been shown that the “stick” boundary conditions can result from inelasticity in collisions alone. As was mentioned in Chapter 1, the ion mobility coefficient is proportional to the diffusion coefficient in the weak electric field limit. The findings of this thesis may provide a more precise interpretation of data in cases such as capillary electrophoresis, where molecular ions are usually treated with hydrodynamic theory employing “stick” boundary conditions. According to this work, “stick” boundary conditions are only applicable in the case of large ions and for the case of largest α. Therefore, one should be careful applying “stick” boundary conditions for smaller ions with a lesser degree of translational-rotational coupling. As well, in Chapter 6, the results of a thorough molecular dynamics benchmark study of the self-diffusion, shear and bulk viscosities and thermal conductivity coefficients for a broad range of densities and different degrees of translational-rotational coupling are presented. The values can be used to improve the treatment of various experimental and computational data and for comparison with theoretical predictions. The validity of several kinetic theory expressions is investigated at different levels of approximation. It was found that at low density the simulation results are closer to the theoretical expressions of the higher level approximation. In the case of self diffusion, shear viscosity and thermal conductivity the accuracy of the result is improved by a few percent in comparison with the estimates given by the formulas with the lower level approximation (Pidduck’s results). Various Enskog predictions have been studied as a function of density and α and compared with the simulation results. The best agreement with the Enskog predictions was observed for the thermal conductivity coefficient as a function of density. Over the whole range of densities the computational data were consistent with theoretical predictions. This can be explained by the fact that the transport of energy depends only weakly upon correlations in the fluid. In the case of shear and bulk viscosity the Enskog predictions work quite well for reduced densities below 0.5, where the difference is a few percent, but at higher densities the discrepancy grows. This is caused by the fact that Enskog theory doesn’t describe collisional momentum transfer well enough and therefore the difference becomes more pronounced at higher densities where collisional transfer becomes more important. For some of the transport coefficients, differing numbers of basis functions included in the kinetic theory descriptions have been tested. Pidduck’s formulae includes rotational effects in its’ tensor basis set but it has the same number of basis functions as usually used to describe the smooth hard sphere fluid, and therefore the degree of accuracy it provides is suited for systems where translational effects prevail. The higher order Enskog expressions include more terms that account for rotational-translational exchange and generally give more accurate results in comparison with the simulations, as found to be the case in this work. It was also found that the density dependence is not affected by additional basis 106 functions but including them allows one to attain better accuracy at low densities. Enskog theory has been tested at different values of α. Varying the value of α allows one to change the degree of rotational-translational energy exchange and thus accurately measure the change in transport coefficients due entirely to this exchange. In the smooth sphere case, translational motion defines the values for the transport coefficient, while in the rough hard sphere case, the translational-rotational exchange is non-zero. It was found that Enskog theory provides good qualitative estimates as a function of α in some cases, but quantitative estimates are generally poor. For shear and bulk viscosity the qualitative predictions are good, while quantitative ones are poor. For thermal conductivity both qualitative and quantitative predictions produce poor results. It is reasonable to expect that including more basis functions in the kinetic theory expressions might improve the agreement, but testing it computationally is very difficult, although the theoretical formulation has already been presented [36]. 8.2 Experimental Data Analysis To establish the framework for analysing experimental data, the findings of this work have been applied to four sets of experimentally measured tracer sizes in Chapter 7. Three sets were obtained via TEM measurements, out of which two sets consisted of nearly spherical silver particles in a nitrogen bath, and one set of copper oxide particles in a nitrogen bath. In all three sets the sizes of the tracer particles were only a few times bigger than the sizes of the bath particles. The fourth set was the Millikan oil droplet experiment data, that is widely used as a reference for hydrodynamic diffusion. As well, smooth and rough hard sphere diffusion data have been examined as a function of Kn to test whether knowing Kn alone is enough to predict diffusion of a particle in a fluid, and to find other factors that may affect this prediction. The following observations have been made. When diffusion was plotted as a function of tracer size, the hydrodynamic 1/R∗ dependence was not observed for either one of the first three sets mentioned above. Even though experimental data was quite scattered, the curves clearly did not show hydrodynamic trends. Instead, they were more consistent with molecular 1/R∗ 2 dependence, suggesting that the Boltzmann equation may be more appropriate for describing the system with a particle size distribution similar to these three sets. The Millikan data set for hydrodynamic diffusion was used to test and compare the applicability of two expression for diffusion, one obtained from the Stokes-Cunningham formula and given by Eq. (7.6), and the other one suggested by Sokolovskii et. al. for smooth hard spheres and given by Eq. (5.4) [23]. It was shown that both expressions perform well in the hydrodynamic regime. Equation (7.4) follows the Millikan data very closely, as expected. Equation (5.4), when used with χ = 1.4, which was the case in the smooth hard sphere analysis by Sokolovskii et. al., works best for Kn < 10; for all other 107 Kn it follows the general trend in the Millikan data with the deviation of about 10% by the time Kn = 90. When χ = 1.1 the fit mimics the one performed with Eq. (7.6). Earlier it has been shown that Eq. (5.4) performs well in the molecular regime [23], where Eq. (7.4) is not applicable. Thus Eq. (5.4) covers wider range of Kn than Eq. (7.4). Smooth and rough hard sphere self diffusion values [23, 41, 42] have been analysed by applying the hydrodynamic formula outside of its’ range of applicability in order to illustrate deviations of the boundary conditions from the hydrodynamically predicted at different values of Kn. It was found, that in cases when the size of the bath particles can be ignored, the fluid shows hydrodynamic-like behaviour for a short range of Kn values, with a “slip” value for the boundary conditions near Kn = 0.6. Smooth and rough hard sphere diffusion data have been analysed as a function of Kn using two different definitions for the radius of the tracer: one that accounts for the size of the bath particles and one that neglects it. It was found that the definition of the radius has a great effect on the interpretation of the data. Both in the molecular and hydrodynamic regimes, including the radius of the bath particles in the expression for the diffusion coefficient led to theoretically predicted results, while excluding the radius gave qualitatively different result, that did not follow theoretically predicted trends. This suggests that the radius of the bath particles should be taken into account when interpreting experimental data, especially for small particles. Lastly, the density dependence was studied for smooth and rough hard sphere diffusion. When a more sensitive way of representing the deviation from hydrodynamic diffusion was employed, the density dependence was observed to become stronger with increase in density. It means that Kn alone is not enough to correctly predict the behaviour of the fluid, especially at higher densities. 8.3 Future Work The results presented in this thesis suggest several directions for future work. First, a benchmark study of the transport coefficients at a wider range of Kn for the rough hard sphere fluid is needed. The results of such simulations can be used for interpreting experimental data obtained in the drift tube experiments with aerosols, colloids and biological molecules, as well as capillary electrophoresis. The transition region between molecular and hydrodynamic regime should be investigated thoroughly to have a clear picture of the fluid behaviour at a given set of parameters. Mass dependence of the transport coefficients can also be investigated. The low density region requires special attention. Studying this region computationally is difficult with the current method. Because collisions are so rare, performing simulations of a desirable quality, defined by a high level of convergence, takes a significant amount of time, that is on the order of many months, if not years. More effective algorithms for 108 performing such simulations may be beneficial. Another direction to explore is the effect of attractive potential upon the boundary conditions on the surface of the tracer in the hydrodynamic regime. Particles of real fluids interact via complex potentials that have both repulsive and attractive parts. A simple square well attractive potential may be introduced into the rough hard sphere fluid, that already assumes a repulsive core defined by the radius, in order to study how it affects the boundary conditions on the surface of the tracer. The attractive potential will change the distribution of the particles around the tracer, which eventually may lead to changes in the boundary conditions. Yet this question should be examined to obtain a definitive answer about the role of attractive potentials. Last, exploring various ways to analyse experimental data in a similar way as presented in Chapter 7 would be most beneficial. More extensive experimental data is required in order to accomplish that. 109 Bibliography [1] B. Carrasco and G. Garcia de la Torre, Biophys. J. 75, 3044 (1999) [2] M. X. Fernandez and G. Garcia de la Torre, Biophys. J., 83, 3039 (2002) [3] G. R. McMeeking, S. M. Kreidenweis, C. M. Carrico, T. Lee, J. L. Collett, and W. C. Malm, J. Geophys. Res., 110, D09206 (2005) [4] R. Baranowski and M. Thachuk, J. Chem. Phys 110, 11383 (1999); 111, 10061 (1999); Phys Rev. A 63, 032503 (2001) [5] R. Baranowski, B. Wagner and M. Thahcuk, J. Chem. Phys, 114, 6662 (2001) [6] C. S. Hoaglund-Hyzer and D.E. Clemmer, Anal. Chem. 73, 177 (2001) [7] C. A. Srebalus Barnes and D. E. Clemmer, Anal. Chem., 73, 424 (2001) [8] G. R. McMeeking, S. M. Kreidenweis, C. M. Carrico, T. Lee, J. L. Collett, and W. C. Malm., J. Geophys. Res. 110 D09206 (2005) [9] J. J. Thomas, B. Bothner, J. Traina, W. H. Benner, and G. Siuzdak, Spectroscopy 18, 31 (2004) [10] E. A. Mason and E. W. McDaniel, Transport Properties of Ions in Gases, (WileyInterscience) (1988) [11] J. O. Hirschfelder, C. F. Curtiss and R. B. Bird, Molecular Theory of Gases and Liquids, (John Wiley, New York) (1964) [12] J.R. Schmidt and J.L., Skinner J. Chem. Phys., 119, 15 (2003) [13] J. R. Schmidt and J. L. Skinner. J. Phys. Chem. B 108, 6767 (2004) [14] A. J. Masters and P. A. Madden J. Chem. Phys., 74, 2450 (1981) [15] S. Bhattacharyya and B. Bagchi, J. Chem. Phys., 106 1757 (1997) [16] Z. Li, Phys. Rev. E, 80, 061204 (2009) [17] A. Uvarov and S. Fritzsche, Chem. Phys. Lett. 401 296 (2005) [18] A. Uvarov and S. Fritzsche, Europhys. Lett. 79 68001 (2007) 110 [19] S. M. Bhattacharyya, Chem. Phys. Lett. 386, 83 (2004) [20] F. Ould-Kaddour and D. Levesque, J. Chem. Phys. 127, 154514 (2007) [21] W. Fenz, I. M. Mryglod, O. Prytula, and R. Folk, Phys. Rev. E 80, 021202 (2009) [22] S. H. Lee and R. Kapral, J. Chem. Phys. 121, 11163 (2004) [23] R. O. Sokolovskii, M. Thachuk, and G. N. Patey, J. Chem. Phys. , 125, 204502 (2006) [24] M. Cappelezzo, C. A. Capellari, S. H. Pezzin, and L. A. F. Coelho, J. Chem. Phys. 126, 224516 (2007) [25] M. Fushiki, Phys. Rev. E 021203 (2003) [26] J. R. Mehaffey and R. I. Cukier. Phys. Rev. Lett. 38, 1039 (1977) [27] H. van Beijeren and J. R. Dorfman, Phys. Rev. A 19, 416-418 (1979) [28] A. Mulero, Theory and Simulations of Hard-Sphere Fluids and Related Systems, (Springer. Berlin Heidelberg) (2008) [29] J. P. Hansen and I. R. McDonald, Theory of Simple Liquids, 2nd ed. (Academic, London) (1986) [30] J. T. Hynes, Ann. Rev. Phys. Chem. 28, 301 (1977) [31] S. Chapman and T. G. Cowling, The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases, (Cambridge University Press) (1952) [32] G. H. Bryan, Br. Assoc. Rep., 83 (1894) [33] F. B. Pidduck, Proc. R. Soc. Lond. A A101, 101 (1922) [34] M. Theodosopulu and J. Dahler, J. Chem. Phys. 60, 4048 (1974) [35] D. W. Condiff, W.-K. Lu, and J. S. Dahler, J. Chem. Phys. 42, 3445 (1965) [36] B. J. McCoy, S. I. Sandler, and J. S. Dahler, J. Chem. Phys. 45, 3485 (1966) [37] J. Montgomery and B. Berne, J. Chem. Phys. 67, 4580 (1977) [38] J. O’Dell and B. Berne, J. Chem. Phys., 63, 2376 (1975) [39] B. Berne, J. Chem. Phys. 66, 2821 (1977) [40] C. Pangali and B. Berne, J. Chem. Phys. 67, 4571 (1977) [41] O. Kravchenko and M. Thachuk, J. Chem. Phys. 134, 114310 (2011) 111 [42] O. Kravchenko and M. Thachuk, J. Chem. Phys. 136, 044520 (2012) [43] G. Subramanian, H.T. Davis, Phys. Rev. A 11, 1430 (1975) [44] J. Schofield and I. Oppenheim Physica A, 187, 210 (1992) [45] R. I. Cukier, R. Kapral, J. R. Lebenhaft and J. R. Mehaffey, J. Chem. Phys., 73, 5244 (1980) [46] A. J. Masters and T. E. Keyes, J. Stat. Phys., 36, 401 (1984) [47] T. Keyes and I. Oppenheim. Phys. Rev. A. 8, 937 (1973) [48] J. R. Mehaffey and R. I. Cukier. Phys. Rev. A. 17, 1181 (1978) [49] L. Waldmann, S. Naturforsch. 18a, 1033 (1963) [50] W. M. Klein, D. K. Hoffman, and J. S. Dahler, J. Chem. Phys. 49, 2321 (1968) [51] M. P. Allen and D. Tildesley, Computer simulation of liquids (Clarendon Press, Oxford) (1987) [52] B. J. Alder, D. M. Gass, and T. E. Wainwright, J. Chem. Phys. 53, 3813 (1970) [53] N. F. Carnahan and K. E. Starling. J. Chem. Phys. 51, 635 (1969) [54] K. L. Rider and M. Fixman, J. Chem. Phys. 57, 2548 (1972) [55] J. R. Mehaffey and R. C. Desai, J. Chem. Phys. 66, 5489 (1977) [56] J. T. Hynes, R. Kapral, and M. Weinberg, J. Chem. Phys. 67, 3256 (1977). [57] R. Kapral and R. C. Desai, J. Chem. Phys. 67, 5645 (1978) [58] J. R. Lebenhaft and R. Kapral, J. Chem. Phys. 74, 6888 (1981) [59] B. J. McCoy and J. S. Dahler, J. Chem. Phys. 50, 2411 (1969) [60] C. Cercignani and M. Lampis, J. Stat. Phys. 53, 655 (1988) [61] H. Sigurgeirsson and D. M. Heyes, Mol. Phys. 101, 469 (2003) [62] D. Chandler, J. Chem. Phys., 62, 1358 (1975) [63] K. Harris, Mol. Phys., 77, 1153 (1992) [64] A. Eastel and L. Woolf, JCS Faraday Trans, I, 80, 549 (1984) [65] R. J. Speedy, Mol. Phys. 62, 509 (1987) [66] M. J. Assael, J. H. Dymond and P. M. Patterson, Int. J. Thermophys. 13, 729 (1992) 112 [67] J. J. Erpenbeck, W. W. Wood, Phys. Rev. A, 43, 4254 (1991) [68] G. G. Stokes, Proc. Cambridge Philos. Soc. 9, 8 (1851) [69] E. Cunningham, Proc. R. Soc. London, Ser. A 83, 357 (1910) [70] M. Knudsen and S. Weber, Ann. Phys. (Leipzig) 36, 981 (1911) [71] I. Langmuir, Filtration of Aerosoles and Development of Filter Materials, OSRD-865 (Office of Scientific Research and Development, Washington, D.C.) (1942) [72] C. N. Davies, Proc. Phys. Soc. London 57, 259 (1945) [73] W. DeMarcus and J. W. Thomas, Theory of a Diffusion Battery, ORNL-1413 (Oak Ridge National Lab, Oak Ridge, TN) (1952) [74] A. E. Reif, Aviation Medicine Selected Reviews (Pergamon Press, Oxford) (1958) [75] N. A. Fuchs, The Mechanics of Aerosols (Pergamon Press, Oxford) (1964) [76] B. E. Dahneke, J. Aerosol Sci. 4, 139 (1972) [77] M. D. Allen and O. G. Raabe, J. Aerosol Sci. 13, 537 (1982) [78] R. L. Buckley and S. K. Loyalka, J. Aerosol Sci. 20, 347 (1988) [79] H. G. Scheibel and J. Porstendorfer, J. Aerosol Sci. 14 113 (1983) [80] Y. Kuga, K. Okauchi, D. Takeda, Y. Ohira, and K. Ando, J. Nanopart. Res. 3, 175 (2001) [81] V. Ya. Rudyak and S. L. Krasnolutsky, Dokl. Phys. 47, 758 (2002) 113
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Transport properties of the rough hard sphere fluid
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Transport properties of the rough hard sphere fluid Kravchenko, Olga 2013
pdf
Page Metadata
Item Metadata
Title | Transport properties of the rough hard sphere fluid |
Creator |
Kravchenko, Olga |
Publisher | University of British Columbia |
Date Issued | 2013 |
Description | This thesis examines the dynamics and physical properties of the rough hard sphere fluid (RHSF). The RHSF model consists of spherical particles with well-defined radii that exchange linear and angular momenta upon collision. The simplicity of this model allows for a precise theoretical description that provides a basis for studying fluid properties on the most fundamental level. Extensive molecular dynamics calculations were made of transport properties as functions of density, tracer particle size, and degree of rotational-translational coupling. Results were compared with the smooth hard sphere case and it was found that transport coefficients change significantly due to rotational-translational coupling, becoming stronger with an increase in coupling. Tracer diffusion coefficients were examined for a range of tracer sizes and at various densities. As tracer particles become larger, their diffusion coefficient moves from an Enskog (molecular) to a Stokes-Einstein (hydrodynamic) functional form; the latter depends upon the boundary condition at the surface of the tracer. These boundary conditions for the RHSF are directly proportional to the degree of rotational-translational energy exchange, and can be tuned from "slip" to "stick" values. The validity of several kinetic theory equations have been examined as functions of density and translational-rotational coupling. At very low densities, Boltzmann theory was accurate even at low order except for describing the dependence upon rotational-translational coupling, where low order expansions are less accurate. Enskog theory performed well at low and moderate densities but deviated at larger densities, as expected. For thermal conductivity as a function of translational-rotational coupling even the qualitative behavior was incorrect. The Enskog predictions for diffusion were also found to be quite poor at low order. Finally, motivated by the results of the thesis, experimental diffusion coefficient data were analyzed, especially for nanoparticles. It was shown that defining the correct radius is crucial for describing such systems. In addition, a new formula for predicting tracer diffusion was tested. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2013-04-12 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NoDerivs 3.0 Unported |
DOI | 10.14288/1.0071935 |
URI | http://hdl.handle.net/2429/44173 |
Degree |
Doctor of Philosophy - PhD |
Program |
Chemistry |
Affiliation |
Science, Faculty of Chemistry, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2013-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nd/3.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2013_spring_kravchenko_olga.pdf [ 769.71kB ]
- Metadata
- JSON: 24-1.0071935.json
- JSON-LD: 24-1.0071935-ld.json
- RDF/XML (Pretty): 24-1.0071935-rdf.xml
- RDF/JSON: 24-1.0071935-rdf.json
- Turtle: 24-1.0071935-turtle.txt
- N-Triples: 24-1.0071935-rdf-ntriples.txt
- Original Record: 24-1.0071935-source.json
- Full Text
- 24-1.0071935-fulltext.txt
- Citation
- 24-1.0071935.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0071935/manifest