A smoothed dissipative particle dynamics methodology for wall-bounded domains
A SMOOTHED
DISSIPATIVE
PARTICLE
DYNAMICS
METHODOLOGY
FOR
WALL-BOUNDED
DOMAINS by Jun Yang A Dissertation Submitted to the Faculty of WORCESTER POLYTECHNIC INSTITUTE in partial fulfillment of the requirement for the degree of Doctor of Philosophy in Mechanical Engineering ____________________________________ April 25, 2013 APPROVED: ______________________________________________ Dr. Nikolaos A. Gatsonis, Advisor Professor, Mechanical Engineering Department, WPI ______________________________________________ Dr. David J. Olinger, Committee Member Associate Professor, Mechanical Engineering Department, WPI ______________________________________________ Dr. Simon W. Evans, Committee Member Assistant Professor, Mechanical Engineering Department, WPI ______________________________________________ Dr. Marcus Sarkis, Committee Member Professor, Department of Mathematics, WPI ______________________________________________ Dr. George E. Karniadakis, Committee Member Professor, Division of Applied Mathematics, Brown University ______________________________________________ Dr. Mark W. Richman, Graduate Committee Representative Associate Professor, Mechanical Engineering Department, WPI 2 ABSTRACT
This work presents the mathematical and computational aspects of a smooth dissipative particle dynamics with dynamic virtual particle allocation method (SDPD-DV) for modeling and simulation of mesoscopic fluids in wall-bounded domains. The SDPD-DV method is realized with fluid particles, boundary particles and dynamically allocated virtual particles near solid boundaries. The physical domain in SDPD-DV contains external and internal solid boundaries, periodic inlets and outlets, and the fluid region. The solid boundaries of the domain are represented with boundary particles which have an assigned position, wall velocity, and temperature upon initialization. The fluid domain is discretized with fluid particles placed in a global index. The algorithm for nearest neighbor particle search is based on a combination of the linked-cell and Verlet-list approaches and utilizes large rectangular cells for computational efficiency. The density model of a fluid particle in the proximity of a solid boundary includes the contribution from the virtual particles in its truncated support domain. The thermodynamic properties of a virtual particle are identical to those of the corresponding fluid particle. A periodic boundary particle allocation method is used at periodic inlets and outlets. Models for the conservative and dissipative forces on a fluid particle in the proximity of a solid boundary are presented and include the contributions of the virtual particles in its truncated support domain. The integration of the fluid particle position and momentum equations is accomplished with an implementation of the velocity-Verlet algorithm. The integration is supplemented by a bounce-forward algorithm in cases where the virtual particle force model is not able to prevent particle penetration. The integration of the entropy equation is based on the Runge-Kutta scheme. In isothermal simulations, the pressure of a fluid particle is obtained by an artificial compressibility formulation for liquids and the ideal gas law for compressible fluids. Sampling methods used for 4 particle properties and transport coefficients in SDPD-DV are presented. The self-diffusion coefficient is obtained by an implementation of the generalized Einstein and the Green-Kubo relations. Field properties are obtained by sampling SDPD-DV outputs on a post-processing grid that allows harnessing the particle information on desired spatio-temporal scales. The isothermal (without the entropy equation) SDPD-DV method is verified and validated with simulations in bounded and periodic domains that cover the hydrodynamic and mesoscopic regimes. Verification is achieved with SDPD-DV simulations of transient, Poiseuille, body-force driven flow of liquid water between plates separated by 10 -3 m. The velocity profiles from the SDPD-DV simulations are in very good agreement with analytical estimates and the field density fluctuation near solid boundaries is shown to be below 5%. Additional verification involves SDPD-DV simulations of transient, planar, Couette liquid water flow. The top plate is moving at and separated by 10 -3 m from the bottom stationary plate. The numerical results are in very good agreement with the analytical solutions. Additional SDPD-DV verification is accomplished with the simulation of a body-force driven, low-Reynolds number flow of water over a cylinder of radius . The SDPD-DV field velocity and pressure are compared with those obtained by FLUENT. An extensive set of SDPD-DV simulations of liquid water and gaseous nitrogen in mesoscopic periodic domains is presented. For the SDPD-DV simulations of liquid water the mass of the fluid particles is varied between 1.24 and 3.3 × real molecular masses and their corresponding size is between 1.08 and 323 physical length scales. For SDPD-DV simulations of gaseous nitrogen the mass of the fluid particles is varied between 6.37×10 and 6.37×10 real molecular masses and their corresponding size is between 2.2×10 and 2.2 × physical length scales. The equilibrium states are obtained and show that the particle speeds scale inversely with 5 particle mass (or size) and that the translational temperature is scale-free. The self-diffusion coefficient for liquid water is obtained through the mean-square displacement and the velocity auto-correlation methods for the range of fluid particle masses (or sizes) considered. Various analytical expressions for the self-diffusivity of the SDPD fluid are developed in analogy to the real fluid. The numerical results are in very good agreement with the SDPD-fluid analytical expressions. The numerical self-diffusivity is shown to be scale dependent. For fluid particles approaching asymptotically the mass of the real particle the self-diffusivity is shown to approach the experimental value. The Schmidt numbers obtained from the SDPD-DV simulations are within the range expected for liquid water. The SDPD-DV method (with entropy) is verified and validated with simulations with an extensive set of simulations of gaseous nitrogen in mesoscopic, periodic domains in equilibrium. The simulations of N (g) are performed in rectangular domains with in the range 0.25×10 -6 ~10×10 -6 m, with fluid mass in the range 1.85×10 -20 ~1.184×10 -15 kg. The mass of the fluid particles is varied between 118 and 7.5×10 real molecular masses. The self-diffusion coefficient for N (g) at equilibrium states is obtained through the mean-square displacement for the range of fluid particle masses (or sizes) considered. The numerical self-diffusion is shown to be scale dependent. The simulations show that self-diffusion decreases with increasing mass ratio. For a given mass ratio, increasing the smoothing length, increases the self-diffusion coefficient. The shear viscosity obtained from SDPD-DV is shown to be scale free and in good agreement with the real value. We examine also the effects of timestep in SDPD-DV simulations by examining thermodynamic parameters at equilibrium. These results show that the time step can lead to a significant error depending on the fluid particle mass and smoothing 6 length. Fluctuations in thermodynamic variables obtained from SDPD-DV are compared with analytical estimates. Additional verification involves SDPD-DV simulations of steady planar thermal Couette flow of N (g). The top plate at temperature T =330K is moving at V xw =30m/s and is separated by 10 -4 m from the bottom stationary plate at T =300K. The SDPD-DV velocity and temperature fields are in excellent agreement with those obtained by FLUENT. 7 ACKNOWLEDGEMENTS
I would like to express the deepest appreciation to my advisor, Prof. Gatsonis, who continually and convincingly conveyed a spirit of adventure in regard to research and scholarship, and an excitement in regard to teaching. You have shaped the professional that I have become and further defined my sense of scientific rigor and engineering creativity. Thank you for the years of support and inspiration. Great thanks are given to Dr. Raffaele Potami whose initial work on SDPD-DV has been the basis for the current implementation of SDPD-DV. His assistance was essential and invaluable to the advancements of SDPD-DV pursued by me. I am very appreciative for the interest that Professor Karniadakis expressed in this work. His inputs strengthened this work and deepened my understanding. I would like to thank also all the members of my committee for their time, patience, and expertise. Several of the WPI faculty provided invaluable expertise on many aspects of this research. However, I would like to specifically acknowledge the efforts of Sia Najafi not only for his expertise, continued support and endless patience with the Linux uninitiated, but also regarding life overall. I would like to thank all of the teachers, coaches and professors who have collectively molded the person that I am today. I would like to specifically thank Prof. Dapeng Hu of the Dalian Univ. of Technology in China, it is your discipline and scrutiny of detail that is instilled in all the work that I do, and at times of need, it is your teachings that resound in my mind and enforce the purity of science and the rigor that must define all that it touches. 8 I would like to extend the sincerest thanks and appreciation to all of the friends that I have made here at WPI throughout the years. You are too numerous to name yet too dear to forget. Gratitude is extended to the Mechanical Engineering Department, and especially to Gloria Boudreau, Statia Canning, Barbara Edilberti, Barbara Furhman and Patricia Rheaume. Most of important, I would like to thank my parents. Their continued love, support, and admiration are empowering. From my parents, I learned the value of education both from their instruction and influence. My success is a reflection of their hard work and sacrifice. I thank them for all that I am and all that they have done for me. Special thanks go to my friends I met in Christian Gospel Church of Worcester for their continuous help and support throughout the years that I had spent at WPI. This work was partially supported by AFOSR’s Computational Mathematics Program under Grant FA9550-10-1-0206. I would also like to acknowledge the support obtained through Teaching Assistantships obtained from the Mechanical Engineering Department at WPI.
This work is dedicated to my mother and father. TABLE OF CONTENTS
ABSTRACT .................................................................................................................................... 3
ACKNOWLEDGEMENTS ............................................................................................................ 7
TABLE OF CONTENTS ................................................................................................................ 9
LIST OF FIGURES ...................................................................................................................... 11
LIST OF TABLES ........................................................................................................................ 16
NOMENCLATURE ..................................................................................................................... 17
1. INTRODUCTION ................................................................................................................... 18
2. SDPD-DV METHODOLOGY AND IMPLEMENTATION .................................................. 45
10 2.2.7 Integration of Fluid Particle Position, Momentum and Entropy Equations ................ 67
3. VERIFICATION, VALIDATION, AND ERROR OF THE ISOTHERMAL SDPD-DV ...... 75
4. VERIFICATION, VALIDATION AND ERROR ANALYSIS OF SDPD-DV ..................... 97
5. SUMMARY, CONCLUSIONS AND RECOMMENDATIONS ......................................... 110
REFERENCES ........................................................................................................................... 117
APPENDIX A. SDPD Analytic Self-diffusion Coefficient ....................................................... 128
APPENDIX B. Boundary Method in 3D with Lucy’s Smoothing Function ............................. 131 LIST OF FIGURES
Figure 1. The DPD cutoff radius used in calculating interparticle forces. ................................... 27
Figure 2. Discrete SPH particles representing the fluid and the support domain of the interpolating function around an SPH particle i . .................................................................. 31 Figure 3. Contribution of neighboring particles in particle i support domain to the number density evaluation. ................................................................................................................ 46 Figure 4. General flow chart of SDPD-DV and post-processing.................................................. 52
Figure 5. Physical wall-bounded domain with an exterior wall and an interior solid body showing the fluid particles (FPs), boundary particles (BPs), and virtual particles (VPs) used in the SDPD-DV. Inlets and outlets are considered periodic. The large rectangular cells are used for nearest neighbor particle search (NNPS). ....................................................................... 53
Figure 6. Normal vector for a BP on a solid boundary. ................................................................ 54
Figure 7. Rectangular cells used for NNPS in the interior, near a solid and near a periodic boundary. .............................................................................................................................. 56
Figure 8. Support domain for a fluid particle (FP) for three cases considered in SDPD-DV. ..... 57
Figure 9. Dynamic virtual particle allocation (DVPA) method used in SDPD-DV for density and force evaluation. A virtual particle (VP) is generated for a fluid particle (FP) near a solid boundary represented by boundary particles (BP). ............................................................... 59
Figure 10. Reflective bounce-forward method in SDPD-DV applied in a case where a fluid particle (FP) penetrates a solid boundary represented by boundary particles (BP). ............. 69
Figure 11. Sampling used inn SDPD-DV for evaluation of the self-diffusion coefficient. The delay time is τ . The total number of samples is M , and K is the number of samples in delay time . The time origins are indexed from m= M-K+
12 Figure 12. Post-processing grid used for evaluation of field (Eulerian) instantaneous and time-averaged properties from SDPD-DV particle samples. ........................................................ 74
Figure 13. SDPD-DV simulations of transient, body-force driven, planar Poiseuille flow. Physical domain showing the BPs and FPs. ......................................................................... 76
Figure 14. Sample-averaged fluid density ρ(Z d ) and velocity Vx ( Z d ,t) from SDPD-DV simulations of transient, body-force driven, planar Poiseuille flow. The domain has L X =10 -3 m, L Y =3×10 -4 m, L Z =10 -3 m, the fluid is H O(l) with ρ a =1,000 kg/m , T i =300 K, η =10 -3 kg·m -1 s -1 , and f x =10 -4 m/s . The analytical profiles are plotted for verification (Morris et al. 1997). .................................................................................................................................... 79 Figure 15. Forces due to virtual particles from SDPD-DV simulations of transient, body-force driven, planar Poiseuille flow. The domain has L X =10 -3 m, L Y =3×10 -4 m, L Z =10 -3 m, the fluid is H O(l) with ρ a =1,000 kg/m , η =10 -3 kg·m -1 s -1 , and f x =10 -4 m/s . ........................... 80 Figure 16. SDPD-DV simulations of transient planar Couette flow. Physical domain showing the BPs and FPs. ......................................................................................................................... 82
Figure 17. Sample-averaged fluid density ρ(Z d )and velocity V x ( Z d, t ) from SDPD-DV simulations of transient, planar, Couette flow. The domain has L X =10 -3 m, L Y =3×10 -4 m, L Z =10 -3 m, the fluid is H O(l) with ρ a =1000 kg/m , T i =300 K, η =10 -3 kg·m -1 s -1 , and V xw =1.25x10 -5 m/s. The analytical profiles are plotted for verification (Morris, Fox and Zhu, 1997). ............... 83 Figure 18. Computational domains used in SDPD-DV and FLUENT simulations of steady, low-Reynolds number water flow with ρ a =1,000 kg/m , T i = 300 K, η =10 -3 kg·m -1 s -1 . The cylinder has R =0.02m and the body force per unit mass is f x = g x =1.5x10 -7 m/s . (a) SDPD-DV domain with L X =0.1 m, L Y =0.015 m, L Z =0.1 m showing the BPs on the cylinder and FPs in the domain. (b) FLUENT domain L X =0.3 m, L Y =0.015 m, L Z =0.3 m includes an array of 13 cylinders to simulate the periodic flow. The insert shows the extend of the SDPD-DV physical domain. (c) Planes P1,P2,P3 and contours C1,C2,C3,C4 used for comparison of SDPD-DV and FLUENT results. .......................................................................................... 85 Figure 19. Sample-averaged V x (r) and P (r) on y = 0 m, y = 7.5×10 -3 m, y = 1.5×10 -2 m planes from SDPD-DV simulations. V x ( x,z ) and P (r) from FLUENT simulation. Steady, low-Reynolds number water flow with ρ a =1,000 kg/m , T i = 300 K, η =10 -3 kg·m -1 s -1 . The cylinder has R =0.02m and the body force per unit mass is f x = g x =1.5x10 -7 m/s . ................ 86 Figure 20. Sample-averaged V x (r) and P (r) on y = 7.5×10 -3 m plane and along C1,C2,C3,C4 from SDPD-DV simulations. V x ( x,z ) and P (r) from FLUENT simulation. Steady, low-Reynolds number water flow with ρ a =1,000 kg/m , T i = 300 K, η =10 -3 kg·m -1 s -1 . The cylinder has R =0.02m and the body force per unit mass is f x = g x =1.5x10 -7 m/s . ................ 88 Figure 21. Phase plot v xi -v yi from SDPD-DV simulations of H O(l) with ρ a =1,000 kg/m , T i = 300 K, η =10 -3 kg·m -1 s -1 and N (g) with ρ a =1,184 kg/m , T i = 300 K, η =10 -5 kg·m -1 s -1 . Results show the scale effects of fluid particle size on velocity. .......................................... 93 Figure 22. Average translational temperature from SDPD-DV simulations of H O(l) with ρ a =1,000 kg/m , T i = 300 K, η =10 -3 kg·m -1 s -1 and N (g) with ρ a =1,184 kg/m , T i = 300 K, η =10 -5 kg·m -1 s -1 . (a) Average translational temperature as a function of time (Case 1, H O(l)). (b) Average translational temperature and standard deviation as a function of m i and D i for H O(l). (c) Average translational temperature and standard deviation as a function of m i and D i for N (g). ............................................................................................. 94 Figure 23. Self-diffusion coefficient and Schmidt number from SDPD-DV simulations of H O(l) , ρ a =1,000 kg/m , η =10 -3 kg·m -1 s -1 , T =300K for various sizes of the SDPD-DV fluid 14 particles. Experimental and analytical estimates, as well as analytical estimates using SDPD-DV parameters are shown for validation and verification. ....................................... 95 Figure 24: Self-diffusion coefficient from SDPD-DV simulations of N (g) for values of h and m i / m N2(g) used in Cases 1-20. Analytical estimates are from Eq.(4.4) and data by Holz and Sacco (2000). ...................................................................................................................... 101
Figure 25 : Shear viscosity from SDPD-DV simulations of N (g) for value of m i / m N2(g) used in Cases 1-20. Initial input viscosity η ( t =0) is plotted for verification. ................................. 102 Figure 26: Effects of time step on equilibrium temperature , translational temperature and pressure from SDPD-DV simulations of N (g). Results are for Case 9, 10, and 11. ...... 103 Figure 27: Effects of time step on equilibrium temperature , translational temperature and pressure from SDPD-DV simulations of N (g). Results are for Case 17, 18, and 19. ... 103 Figure 28: Computational domains used in SDPD-DV and FLUENT simulations of steady N (g) Couette flow with ρ a =1.184 kg/m , η =1.79×10 -5 kg·m -1 s -1 , κ =0.026 W·m -1 K -1 . The upper wall T = 300 K, lower wall T = 300 K and V xw =30m/s. ..................................................... 107 Figure 29: Sample-averaged V x (r) and T (r) on y = 0 m, y = 1.5×10 -5 m, y = 3×10 -5 m planes from SDPD-DV and FLUENT simulations of steady state Couette flow. The domain has L X =10 -4 m, L Y =3×10 -5 m, L Z =10 -4 m, upper wall T = 300 K, lower wall T = 300 K, and V xw =30 m/s. The fluid is N (g) with ρ a =1.184 kg/m , η =1.79×10 -5 kg·m -1 s -1 and , κ =0.026 W·m -1 K -1 .............................................................................................................................................. 108 Figure 30: Sample-averaged fluid temperature T (r d ) and velocity V x (r d ) of steady-state Couette flow from SDPD-DV and FLUENT simulations. The domain has L X =10 -4 m, L Y =3×10 -5 m, L Z =10 -4 m, upper wall T = 300 K, lower wall T = 300 K and V xw =30 m/s. The fluid is N (g) with ρ a =1.184 kg/m , η =1.79×10 -5 kg·m -1 s -1 and , κ =0.026 W·m -1 K -1 . .......................... 109
15 Figure 31: Function i h for Lucy kernel in 2-D and 3-D. ...................................................... 134 Figure 32: Function i h for Lucy kernel in 2-D and 3-D. ..................................................... 141 Figure 33: Function i h for Lucy kernel in 2-D and 3-D. .................................................... 141 Figure 34: Function i h for Lucy kernel in 2-D (Vazquez-Quesada formulation) and in 3-D.............................................................................................................................................. 142 Figure 35: Function i h for Lucy kernel in 2-D (corrected) and 3-D. .................................. 145 LIST OF TABLES
Table 1. Properties of gaseous nitrogen and liquid water at STP. ................................................ 21
Table 2. Nearest Neighboring Particle Search Algorithm used in SDPD-DV ............................. 56
Table 3. Dynamic Virtual Particle Allocation and Density Evaluation Algorithm used in SDPD-DV. ........................................................................................................................................ 60
Table 4. Time integration in SDPD-DV. Velocity-Verlet is used for the particle position and momentum equation, and Runge-Kutta is used for the entropy equation. ............................ 68
Table 5. Bounce-forward Method used in SDPD-DV in case of particle penetration. ................. 69
Table 6: Input parameters used in SDPD-DV simulations of transient planar Poiseuille flow, transient planar Couette flow, and flow over a cylinder. ...................................................... 77
Table 7. Input and derived parameters in SDPD-DV simulations of mesoscale flows of liquid water at equilibrium states in rectangular domains with periodic boundaries. ..................... 90
Table 8. Input and derived parameters in SDPD-DV simulations of mesoscale flows of N (g) at equilibrium states in rectangular domains with periodic boundaries. ................................... 91 Table 9: Input parameters in full-set SDPD-DV simulations of mesoscale flows of N (g) at equilibrium states in rectangular domains with periodic boundaries. ................................... 98 Table 10: Derived parameters in SDPD-DV simulations of mesoscale flows of N (g) at equilibrium states in rectangular domains with periodic boundaries. ................................... 99 Table 11: Fluctuations in temperature, density, pressure and velocity, from SDPD-DV and analytical expressions. ........................................................................................................ 104
Table 12: Input parameters used in SDPD-DV non-equilibrium simulations of Couette flow. . 106 NOMENCLATURE
Boldface denotes a vector. The magnitude of a vector is denoted using the same symbol as the vector, but without boldface. Duplicate use of a symbol, or usage not defined below, will be clarified within the text. c Sound speed V c Specific heat capacity at constant V C Heat capacity at constant volume i D Individual particle diameter Density E Energy f Body Force C F Total conservative force D F Total dissipative force R F Total random force h Smoothing length Thermal conductivity B k Bolztmann constant Kn Knudsen Number m Mass Mean free path N Number of particles Viscosity P Pressure r position S Entropy t Time step T Temperature v velocity V Volume Bulk viscosity 18
1. INTRODUCTION
The numerical modeling and simulation of gases, liquids and multi-phase systems above the atomic spatiotemporal scales and below macroscopic scales has become increasingly important due to numerous applications in physical, biological and engineering systems. This regime can be characterized as mesoscopic and requires new mathematical and computational methods because the traditional atomistic (microscopic) and hydrodynamic (macroscopic) descriptions are not valid for either computational or theoretical reasons. One of the key characteristics of the mesoscopic flow regime is the presence of thermal fluctuations. The Smooth Dissipative Particle Dynamics (SDPD) developed by Espanol and Revenga (2003) has been proposed as a method appropriate for mescoscopic flows with fluctuations. The SDPD invokes the Smoothed Particle Hydrodynamics (SPH) which is a well-developed method for the Navier-Stokes equations (Liu and Liu, 2003). Using the GENERIC framework developed by Ottinger (2005) to describe hydrodynamic fluctuation, Espanol and Revenga arrived to SDPD discrete equations, which include thermal fluctuations. A three-dimensional implementation of the SDPD method in unbounded domains with dynamic virtual particle allocation (SDPD-DV) has been pursued in the Computational Gas&Plasma Dynamics Lab (CGPL) at WPI (Yang et al., 2011; Yang et al., 2012; Yang et al., 2013; Gatsonis et al., 2013). This work on SDPD-DV is part of the effort at the CGPL to develop simulation methods that apply to the investigation of micro- and nano-scale flows and involves the Direct Simulation Monte Carlo (DSMC) method for gases (Gatsonis et al., 2010; 2013) and particle-in-cell for plasmas (Gatsonis et al., 2009). There are three goals in this work: Revise and further implement existing algorithms in SDPD-DV, as well as develop and implement new algorithms. 19 Validate and verify the algorithms of the SDPD-DV code by comparisons with experimental, analytical and numerical solutions. Investigate thermodynamic properties and transport coefficients with applications of the SDPD-DV code to mesoscopic systems in equilibrium. We review below characteristic length, time scales and non-dimensional parameters used to characterize liquids and gases in order to establish notions of importance to mesoscale, those of continuum, local thermodynamic equilibrium, and fluctuations . The molecular interaction featured in all states of matter (solids, liquids, and gases) can be represented by a force due to the Lennard-Jones potential (Allen and Tildesley, 1987) given by
12 6 ij ij ij r rV r c d , (1.1) where, r is the distance separating the molecules i and j , ij c and ij d are parameters particular to the pair of interacting molecules, is a characteristic energy scale, and is a characteristic length scale. The force between two molecules can be achieved by derivative of Eq.(1.1)
13 7
48 2 ij ijij ij
V r dr rF r cr . (1.2) The characteristic atomic time scale associated with this Lennard-Jones molecular interaction is m , (1.3) where m is the mass of the individual molecule. The mean molecular spacing assuming that molecules are in contact with one another is defined as , (1.4) 20 where the macroscopic density is . For H O(l) typical scales are shown in Table 1. For gases additional physical scales can be defined (Vincenti and Kruger, 1975; Gombosi, 1994). The number of molecules in one mole of gas is the Avogadro number
23 1 A N . The volume occupied by 1 mole of gas at a given temperature and pressure is constant m V at STP. Gases obey the perfect gas law B P nk T (1.5) where, n is the number density,
23 2 2 B k is the Bolztmann constant, T (K) is the temperature. An important length scale in a gas is the mean molecular spacing defined as n . (1.6) The molecular diameter of a gas molecule i D is another physical scale which cannot be precisely defined but can be derived from viscosity measurements. Table 1 provides estimates based on the hard-sphere model. For gases undergoing binary collisions (dilute gases), it is required that i D . (1.7) Transport properties (such as viscosity and thermal conductivity) in a gas relate to collisions. The characteristic collision time in a dilute gas derived from the hard-sphere collision frequency is c i D nv (1.8) where v is an average relative speed. The average distance between collisions, called the mean free path, is given by i D n (1.9) 21 The mean-square thermal molecular speed of a molecule for a gas in equilibrium is given B k Tc m (1.10) The properties of a typical gas and liquid at standard conditions are shown in Table 1. Table 1. Properties of gaseous nitrogen and liquid water at STP.
Property N (g) H O (l) Molecular diameter, D i -10 m 3×10 -10 m Number density, n m -3 m -3 Intermolecular spacing, δ -9 m 4×10 -10 m Molecular speed, c
500 m/s 1,000 m/s Mean free path, λ 6.54×10 -8 m 2.5×10 -10 m Collision time, c -10 s 1.25×10 -13 s While all fluid systems consist of molecules, macroscopic variables are obtained from averages of properties that depend on molecular velocities. These macroscopic averages can be functions of time and space. The continuum hypothesis allows us to describe the state of the fluid using a number of thermodynamic variables (or fields) that depend on position r and time t , for example, , t v r , , t r , , T t r , , P t r . One must be able to define sampling volumes, infinitesimal subsystems, with enough molecules so that statistical averages can be performed and provide a homogenous thermodynamic field. The local thermodynamic equilibrium hypothesis implies that the thermodynamic variables for these infinitesimal subsystems vary with time and space but satisfy the same relations as the equilibrium thermodynamics properties. Once these fields (or point variables) are defined conservation 22 equations can be obtained which contain dissipative fluxes (the stress tensor, the heat flux, and diffusion flux). Phenomenological expressions are used to relate the dissipative fluxes , J t r with corresponding conjugate thermodynamic forces , X t r in the form given [Ortiz De Zarate and Sengers, 2006] , , , J t M t X t r r r , (1.11) The functional derivatives of the dissipative fluxes , M t r are phenomenological coefficients, commonly referred to as the Onsager coefficients, and must be symmetric which implies , , M t M t r r . The dependence of the Onsager coefficients, in principle, on space and time through the local state variables, in most practical applications can be neglected. For example, the heat flux , t Q r is given by Fourier’s law with thermal conductivity ( , ) ( , ) t T t Q r r . (1.12) With the phenomenological models defined, conservation laws, such as the Navier-Stokes equations are derived. It is well known that fluctuations exist in fluids in thermodynamic equilibrium as introduced by Landau and Lifshitz (Statistical Physics, Part 1, Ch XII, 1980) as well as those in non-equilibrium (Ortiz De Zarate and Sengers, 2006). These fluctuations are spontaneous variations of thermodynamic variables about their mean values and are characterized by correlation functions of the fluctuating properties. Fluctuations in systems in thermal equilibrium can be derived from statistical physics (Landau and Lifshitz, Statistcial Physics, Part 1, Ch XII, 1980) and kinetic theory in case of gases (Groot and Mazur, Ch IX, 1962). The fluctuation for a thermodynamic quantity f about its mean f measured in a volume V is given as 23 f f f , (1.13) where The root-mean-square fluctuation is given by the variance ( ) f f f . (1.14) and the relative fluctuation is inversely proportional to N ( ) 1 ff N . (1.15) where the denote the mean values of lengthy expressions. For example, in an infinite fluid in equilibrium the variance in number of particles in a volume V with temperature T is given by (Landau and Lifshitz, 1980, Ch XII; Ortiz de Zarate and Sengers, 2006, Ch.3) B B TT k TN V k TNN V P V , (1.16) where the isothermal compressibility is T T
VV P . (1.17) In a given volume element V , the variance in macroscopic density is independent of time t and equals to (Ortiz de Zarate and Sengers, 2006, Ch.3) E TB
Sm k TV V , (1.18) where
E T B mS k T . Similarly, the variance in pressure is ( ) BB S T
P k TP k T V V . (1.19) The variance in volume is given by B B TT
VV k T k T VP . (1.20) 24 For dilute gas T P and the above expressions become N N , (1.21) n , (1.22) ( ) B PP k T V , (1.23) VV n . (1.24) Meanwhile in general, the variance in temperature is B V k TT C . (1.25) The variance in entropy is P S C . (1.26) where V C and P C are the heat capacity of the body as a whole at constant volume and constant pressure respectively. In general fluctuations increase with decreasing volume V , where the property f is measured. It should be noted that the fluctuations of velocity are statistically independent of those of the other thermodynamic quantities. The variance of each Cartesian component of the velocity is equal to B k Tm V (1.27) Fluctuating hydrodynamics for fluids in thermodynamic equilibrium can be described by the usual hydrodynamic equations (e.g. Navier-Stokes) with introduction of random noise terms. One approach introduced by Landau and Lifshitz (Statistical Physics, Part 2, Ch IX, 1980) is 25 referred to as stochastic forcing, and treats the dissipative fluxes as the sum of an average plus a fluctuation term, which is defined as the fluctuating dissipative flux. ( , ) ( , ) ( , ) ( , ) a J t M t X t J t r r r r (1.28) The fluctuation term , J t r should satisfy that ( , ) 0, J t r (1.29) ( , ) ( , ) , J t J t C t t r r r r , (1.30) where , 2 . B C k TM r r r r (1.31) The expressions in Eq. (1.31) is the so called fluctuation-dissipation theorem (FDT). Various extensions of fluctuating hydrodynamics to systems of thermal non-equilibrium under the local equilibrium assumption have been proposed. (de Groot and Mazur, 1962; Ottinger, 2005). A non-dimensional parameter used to describe the applicability of the continuum approach in a gaseous flow is the local Knudsen Number given by
Kn L . (1.32) In this definition L is the local length scale of the gradient of a macroscopic quantity , f t r fL f , (1.33) (Bird, 2007; Gombosi, 1994). When the Knudsen number is used to define the validity of the Navier-Stokes equations , it is often required that Kn . In order to simulate mesoscopic flows, new mathematical and computational methods are required because the traditional atomistic (microscopic) and hydrodynamic (macroscopic) 26 descriptions are not valid for either computational or theoretical reasons. The fundamental method appropriate for molecular (atomistic) scales is Molecular Dynamics (MD) (Allen and Tildesley, 1987; Haile, 1997). For example, a simulation of 1 micron volume at STP of N (g) and H O(l) contains 3×10 and 2×10 molecules respectively. A simulation based on MD is certainly not feasible. Various particle and hybrid (particle-fluid) simulation methods covering from atomistic to hydrodynamic scales are reviewed in Koumoutsakos (2005). From the continuum approach, the Navier-Stokes cannot represent the thermal fluctuation in mesoscale. Several models and numerical methods have been proposed for mesoscopic fluid dynamics, derived from either bottom-up (molecular) and from top-down (hydrodynamic) scale. We review below the Dissipative Particle Dynamics (DPD) and Smoothed Particle Hydrodynamics (SPH). A fundamental method for mesoscopic domains is the Dissipative Particle Dynamics (DPD) method introduced by Hoogerbrugge and Koelman (1992). The method can be interpreted as a coarsening approach, where the DPD particles represent clusters of molecules interacting by means of repulsive, dissipative and random forces. This clustering permits use of larger integration time steps and allows simulation of spatial scales much larger than those covered by MD (Keaveny et al., 2005). The statistical mechanics context behind DPD is presented by Espanol and Warren (1995). In DPD the fluid region is discretized by a number of fluid particles each having mass, and gorvened by Newton’s equation of motion (Hoogerbrugge and Koelman, 1992) ii ddt rv , (1.34) 27 ii i dm dt vF . (1.35) where i v is its velocity, i r is its position. i F is the total force exerted on particle i by particle j , given as the sum of interparticle forces, consisting of a conservative component Cij F , a dissipative component Dij F , and a random component Rij F (Hoogerbrugge and Koelman, 1992) C C ijij ij rF F e (1.36) // ,
D D ij ij ijij ijij ij r rr rrF v (1.37) / .
R R ij ij ijij R ij r r rF (1.38) All particles j in a sphere of radius c r , which is called cutoff radius, are interacting with particle i as shown in Figure 1. The conservative force is usually a soft repulsion given by /max 1 , 0 C ij ij cij r r rF a . The strength
D ij r and R ij r in the dissipative force and random force, are coupled by /max , 01 D Rij ij ij c r r r r . The coefficients and R are coupled by R B k T . The ij in Eq. (1.38) is a symmetric Gaussian random variable with zero mean and unit variance. Figure 1. The DPD cutoff radius used in calculating interparticle forces. c r
28 The total force exerted on particle i is given as (Espanol and Warren, 1995) C D Ri ij ij ijj i j i j i
F F F F . (1.39) The time evolution of the position and momentum equation of a DPD particle using a timestep dt is given by , i d dt r v (1.40) C D Ri i i i d dt dt dtm v F F F (1.41) There have been many improvements to DPD since its introduction and the method has been applied to a variety of mesoscale systems, including binary immiscible fluids (Novik and Coveney, 1997), colloidal behavior (Dzwinel et al., 2006), DNA in microchannels (Symeonidis et al., 2006), two phase flows (Tiwari and Abraham, 2006a), flow over rotating cylinder (Haber et al., 2006), nanojet breakup (Tiwari and Abraham, 2006b), polymers (Symenodis et al., 2006; Symeonidis and Karniadakis, 2006; Pan et al., 2008), and water in microchannels (Kumar et al., 2009), evaluation of transport properties (Ripoll et al., 2001), and fluids out of equilibrium (Ripoll and Ernst, 2004). A considerable amount of effort has been devoted to boundary conditions in DPD. These studies include the non-slip boundary condition (Hoogerbrugge and Koelman, 1992; Espanol and Warren, 1995; Willemsen et al., 2000; Xu and Meakin, 2009; Wang et al., 2006; Duong-Hong et al., 2004; Pivkin and Karniadakis, 2005), slip boundary conditions (Smiatek et al., 2008), periodic boundary conditions (Chatterjee, 2007), and wall reflection laws (Revenga et al., 1999). Attempts to arrive at an energy-conserving DPD have appeared by adding a random heat term (Avalos and Mackie, 1997; Mackie et al., 1999) or mechanical energy (Chaudhri and 29 Lukes, 2009). A generalized form of DPD incorporating an internal energy and a temperature variable for each particle was presented by Espanol, (1997) and Ripoll et al. (1998). The energy-conserving DPD model was used to investigate the heat conduction in nanofluids by He and Qiao (2008).
The SDPD method has its origins on Smoothed Particle Hydrodynamics (SPH), which was originally developed for modeling of astrophysical phenomena (Lucy, 1977; Benz et al., 1989). SPH was extended to simulate problems of continuum solid and fluid mechanics. Several forms of SPH equations can be derived based on the form of equations and the particle approximation involved. We review the SPH derivation for the compressible, viscous Navier-Stokes equations written in terms of field variables density , t r , velocity , t v r and internal energy . , e t r ., in the case of free external force (Liu and Liu, 2003) , ,, , D t ttDt r v rr x (1.42) , 1 ,,
D t PDt t v r r x x (1.43) , ,, 2 ,
De t P tDt t t r v rr x r (1.44) The superscripts and denote the coordinate directions. The viscous stress is proportional to shear stress with dynamic viscosity by
2( )3 v v vx x (1.45) where the shear stress v v vx x (1.46) The internal energy per unit mass ( , ) e t r for an ideal gas is given in terms of specific heat capacity V c by , V e t c T r (1.47) and the pressure can be expressed as ( , ) 1 ( , ) ( , ) P t t e t r r r (1.48) The energy Eq. (1.44) assumes that heat flux and the body forces are neglected. The derivation of SPH equations encompasses two steps: (a) the integral approximation of fluid fields such as , t r , , t v r , , e t r and their derivatives , and (b) the particle approximation of these fields . In SPH the fluid is represented by a finite number of particles each with mass i m and volume i V , as shown in Figure 2(a). The particle approximation of any field i f x is given by , N ji j i jj j mf x f x W x x h (1.49) where , i j W x x h is the smoothing kernel function or smoothing function, h is the smoothing length defining the influence volume of the smoothing function , i j W x x h as shown in Figure 2(b). The smoothing function , W x x h in SPH must be normalized over its support domain as follows , 1. W x x h dx (1.50) and be symmetric 31 , 0. x x W x x h dx (1.51) (a) (b) Figure 2. Discrete SPH particles representing the fluid and the support domain of the interpolating function around an SPH particle i . The discrete counterparts of the constant and linear consistency conditions as expressed in Eq. (1.50) and Eq. (1.51) are , 1, N i j jj
W x x h x (1.52) , 0. N i j i j jj x x W x x h x (1.53) As discussed by Liu and Liu (2003), the discretized consistency conditions are not always satisfied due to the unbalanced particle distribution in cases where the support domain intersects with the boundary or support domain is irregularly distributed. For the Navier-Stokes equations, Eq. (1.42)-(1.44) can be expressed in form of SPH equations (Liu and Liu, 2003) given as , i j W x x h j i h , Ni j ijj m W (1.54) N Nj ij j j iji i i ij jj ji j i i j i p W Wd pm mdt x x v (1.55) N i j iji ij ij i ij i j i i p p WdE m vdt x (1.56) SPH is a well-developed method and a plethora of review papers and books have appeared including, Monaghan (1994), Randles and Libersky (1996), Liu and Liu (2003), Monaghan (2005), Li and Liu (2007), Monaghan (2009). There has benn a wide range of SPH applications and a comprehensive review is outside of the scope of this work. Studies include, incompressible fluids (Morris et al., 1997; Ellero et al., 2007), free surface flow (Monaghan, 1994; Fang et al., 2009), viscoelastic flow (Vazquez-Quesada and Ellero, 2012; Ellero et al., 2006), error estimation (Amicarelli et al., 2011; Fatehi and Manzari, 2011), immersed boundaries (Hieber and Koumoutsakos, 2008), multi-phase multiscale flow (Hu and Adams, 2006), thermal fluctuations in viscoelastic fluid (Vazquez-Quesada et al., 2009 c), and phase separating fluid mixture (Thieulot et al., 2005). The SDPD was developed by Espanol and Revenga (2003) from a top-bottom approach from SPH as the thermodynamically consistent alternative to DPD. The SDPD invokes the Smoothed Particle Hydrodynamics (SPH), applied to the discretization of the compressible, viscous Navier-Stokes equations written in terms of field variables density , t r , velocity , t v r and entropy , S t r which is given by Batchelor (1967) , , D t tDt r r v , (1.57) 33 , 1, , , 3 D t P FDt t t t v r v vr r r , (1.58) , , , DS tT TDt t t r r r . (1.59) The material derivative, following a fluid particle, is
D Dt t v . In the above system P is pressure, F is the external force exerted per unit mass, is the shear viscosity, is the bulk viscosity, and is the thermal conductivity. The viscous heating field ( , ) t r is defined by v v v (1.60) where the traceless symmetric part of the velocity gradient tensor is T v v v v (1.61) Through the use of the GENERIC framework developed by Ottinger (2005) to describe hydrodynamic fluctuations, Espanol and Revenga arrived to the SDPD discrete equations that include thermal fluctuations. In SDPD the independent variables are position t r , velocity t v and entropy S t of each fluid particle, a deviation from most mesoscopic and hydrodynamic models that involve the energy of the particle. i i d dt r v , (1.62) d = d d i i C D R m t t v F F F , (1.63) i i Vi C R T dS E dt E dt E . (1.64) where C F , D F and R F are the conservative term, dissipative term and random terms; Vi E , C E and R E are viscous, conductive and random terms. Details of the formulation are presented in Ch. 2. SDPD requires a state equation to be expressed in the form , , E N V S . Then, the pressure , P t r and temperature , T t r are obtained through the Maxwell relations 34 , ,, E N V SP t V r , (1.65) , ,, , E N V ST t S t r r (1.66) The work done by Espanol and Revenga (2003) follows the idea introduced by Espanol and Serrano (1999) about treating an SPH or a DPD particle as moving thermodynamic subsystems. In both of their works, the Sackur-Tetrode equation for internal energy , ,
E S N V is employed to close the system. However, the testing case by Espanol and Revenga (2003) only involved the deterministic part of system equation. Espanol and Serrano (1999) perform test with a full model to show the fluctuation in entropy within equilibrium system are in the order of Boltzmann constant and around a constant value. All of the testing cases are under reduced unit. The SDPD investigations so far addressed fundamental issues of the method. Serrano (2006) compared the SDPD and the Voronoi fluid particle model in a shear stationary flow. He found the efficiency of the two methods comparable. The accuracy of the Voronoi approach was found superior for regular ordinate systems while SDPD produced more accurate results for arbitrary, disordered configurations. An incompressible, isothermal SDPD model was used by Litvinov et al. (2008) to model polymer molecules in suspension. The entropy equation was decoupled from the governing equations. A quantic spline kernel function was applied to their model. Periodic boundary conditions are considered to model the bulk fluid. Virtual particles mirrored at the surface of solid boundary are implemented when the support domain of a fluid particle overlaps with the solid wall surface. They concluded that the virtual particles would introduce errors on curved surfaces as Morris et al. (1997) pointed out. In their model, the SDPD interact hydrodynamically as well as by an additionally finitely extendable nonlinear elastic springs. In the test cases, the 35 radius of gyration and the end-to-end radius are compared with analytical solutions at different chain lengths. A static structure factor is employed to extract the static factor exponent. A further test involved the Rouse mode. A form of diffusion coefficient evaluated from the mean square displacement of the center of mass is introduced and compared to an analytical formula for the 2D case. All parameters in these cases are in reduced units. Litvinov et al. (2008) concluded that the confinement of solid boundary affects the polymer configuration statistics and induces anisotropic effects. Their model is a guideline for further handling realistic microfluidic applications. Vazquez-Quesada et al. (2009) investigated the consistent scaling of thermal fluctuations in SDPD by taking the point of view that the SDPD particles are real portions of fluid material instead of only just simple Lagrangian moving nodes. They used the isothermal SDPD model and normalized units in the simulations. They pointed out that the deterministic part of SDPD governing equation is scale free while the velocity variance from the stochastic part shows dependence on the physical length of SDPD particle. The tests are performed on colloidal particles (radius R in reduced unit) suspended in a Newtonian fluid (viscosity , mass density , temperature T ) with three different particle resolutions (number of particle N ) in a box size L ). Their results show that the velocity variance vary with the resolution of solvent, however it does not affect the velocity variance of the colloid particle. As a further proof of their result, they examined a polymer molecule in a suspension with different resolution of solvent ( N ), and added finitely extensible nonlinear elastic spring forces between the polymer molecules. The model utilized by Vazquez-Quesada et al. (2009) for viscoelastic fluids is identified as SPH with thermal fluctuations, which is actually the SDPD form added by a dimensionless 36 conformation tensor to characterize the elongation of the polymer molecules within the fluid particles. They provided a full scheme of SDPD with closed equations such as total energy, entropy, conformation tensor. The entropy function is given by the logarithm of the number of microstates and coupled with total energy and conformation tensor. Their simulations are performed in a 2D periodic square box of length L and N fluid particles. They considered uniform shear flow with density , kinematic viscosity , and Kolmogorov flow at Reynolds number Re 1 . The dynamics of the conformation tensor are studied in terms of the dynamics of its eigenvalues and eigenvectors and have a non-eligible contribution to thermal fluctuations. They state that the introduction of thermal fluctuations in the viscoelastic fluid particle model is analogous to the stochastic contribution introduced by Landau and Lifshitz (Statistical Physics, part 2, Ch. IX, 1980) in fluctuating hydrodynamics. Litvinov et al. (2009) provided an analytical expression of the self-diffusion coefficient for isothermal incompressible SDPD model. The coefficient in the self-diffusion coefficient formula depends on the quantic spline smoothing function. Their expression was tested by several simulations performed in a 3D periodic box with N SDPD particles. The simulations assumed , B k T , m L N , L , / 15 h L with resulting self-diffusion D . The simulations varied the dynamic viscosity . The actual self-diffusion coefficient and Schmidt number defined as Sc D was calculated from the mean-square displacement and was compared with the analytical value under different dynamic viscosity. Schimdt number that depends on fluid particle size does not always agree with the theoretical predictions, and extensive preliminary computations are needed to characterize the diffusion properties of the solvent. They point out that the diffusion properties depend on the choice of the smoothing length. They also state that their simulations do not show the solid-like structures 37 found in DPD simulations at high coarse-graining levels as discussed by Pivkin and Karniadakis (2006). Bian et al. (2012) utilized the SDPD method to model solid particles in suspension. They start from a reduced set of Navier-Stokes equation for incompressible flow to obtain an SDPD model. Closure is obtained with a state equation in the form of artificial compressibility given by Batchelor (1967). The non-slip velocity boundary condition is embedded on all solid-liquid interface by applying frozen particle as described by Morris et al. (1997). The velocity dynamically assigned to a frozen particle is calculated based on the distance to the tangent plane of the closest point to fluid particle. Therefore, it requires that the cutoff radius of support domain should be smaller than the smallest surface (both of convex and concave) curvature radius. The total force exerted on solid particle would be the summation of all forces exerted on each frozen particles in solid particle. They use the Velocity Verlet algorithm as the time integrator. Their numerical simulations are carried out with reduced units, including cases under non-Brownian and Brownian conditions. Several tests under non-Brownian condition involved flow through a fixed circular or spherical object in a periodic array, a particle moving in a Newtonian fluid under unsteady situations, a particle rotating under shear flow, and hydrodynamic interactions between two approaching spheres. The Brownian disk (in 2D) and Brownian sphere (3D) are investigated with thermal fluctuations producing its ultimate Brownian diffusive dynamics. The mean square displacement in a 2D domain is analyzed, and the diffusional behaviors are studied for both cases. A corrected form of diffusion coefficient related to drag force is provided and is related to the Einstein-Stokes equation for the actual diffusion coefficient from mean square displacement. A more complicated simulation involves a 38 colloidal particle in the vicinity of an external boundary. The split form of the diffusion coefficient is introduced which depends on the distance to the boundary. Among the outstanding theoretical and computational issues in SDPD is its implementation in domains with solid boundaries of arbitrary geometry. We address this issue in this work by developing an SDPD method for wall-bounded domains. We review boundary condition approaches in DPD and SPH in order to provide the necessary background for our method. In mesh-free methods the wall and its effects on the fluid are modeled using various types and layers of ghost particles or frozen particles. Ghost particles can be loaded initially with static properties (Morris et al., 1997) or can be dynamically generated with properties updated during the simulation (Randles and Libersky, 1996). Static ghost or virtual particles are preloaded as uniformly distributed layers (Duong-Hon et al., 2004) or as interacting particles in the flow and solid boundary regions loaded with the same density as the fluid particles. (Hoogerbrugge and Koelman, 1992; Boek et al., 1997; Revenga et al., 1999; Liu and Liu, 2003). Dynamically allocated ghost or virtual particles are generated by reflecting neighboring fluid particles which may lead to an imperfect representation of a curved boundary at low resolution (Morris et al., 1997). The difficulty arises in assigning the physical properties of such ghost particles and in defining the interaction forces between them and the fluid particles in the proximity of the wall. Such forces must ensure that a fluid particle does not penetrate the wall and that the no-slip condition is enforced. In general, the fluid particle force is composed of a repulsive and a dissipative term. The soft repulsive force acting on a fluid particle near a wall from the neighboring fluid particles may not be sufficiently strong to prevent wall penetration. To 39 overcome this problem a stronger repulsive force between the fluid and the wall ghost particles, in the Lennard-Jones form, has been used in SPH (Monaghan, 1994). Different solutions have been proposed in DPD for imposing wall conditions, such as increasing the repulsive force coefficient for wall/fluid interaction (Pivkin and Karniadakis, 2006), increasing the wall particle density (Fedosov et al., 2008) or imposing bounce-back and bounce-forward boundary conditions (Revenga et al., 1999; Pivkin and Karniadakis, 2005). Some implementations of SDPD used dynamical ghost or virtual particle, including Hu and Adams (2006), Litvinov et al. (2008). A Lees–Edwards boundary condition is applied for solid boundary by Litvinov et al. (2010). Quesada (2009c) used the distance between a fluid particle and the solid wall in order to evaluate the truncated area of the support domain. This area is then used to evaluate the density and force on the fluid particle. This method may lead to an overestimation of the fluid particle density compared to the interior region. The effectiveness of wall reflections was discussed by Revenga et al. (1999), including specular reflection, bounce-forward reflection, Maxwellian reflection, and bounce-back reflection. Pivkin and Karniadakis (2005) placed an extra thin layer of DPD particles inside the domain and adjacent to the solid boundary with an adjusted wall-fluid conservative force parameter, which is estimated according to the fluid density, to hold the no-slip boundary conditions. A similar model was developed by Pivkin and Karniadakis (2006) to control the density fluctuations near a solid boundary by applying an adaptive force directed perpendicular to the wall. The force is also adjusted according to the density and bins adjacent to the solid boundary. These adapted models are verified and compared by Fedosov et al. (2008). Another algorithm that modifies the boundary force is implemented by Altenhoff et al. (2007), in which the boundary force is estimated by the probability density function of the force contributions in 40 the bin adjacent to solid wall. A phase-field interface representation to DPD for imposing the no-slip boundary condition was proposed by Xu and Meakin (2009), which considers the solid boundary as a phase indicated by a variable.
In this work, we use the summation form in Eq. (1.54) for evaluation of density which is a direct way of the approximation of SPH to the density itself. This form involves the contribution from neighboring particles in support domain by smoothing function. Flebbe et al. (1994) suggested that the contribution of the self-density should be discounted since overestimation arises when self-density in involved. Whitworth et al. (1995) suggested that although an overestimate in density arises from the thermal fluctuation initially with randomly distributed particles, it is eliminated after the system equilibrates.
The SDPD investigations so far addressed fundamental issues of the method but considered only the isothermal formulation, which involves continuity and momentum equations alone (Litvinov et al., 2008; Litvinov et al. 2010, Vazquez-Quesada, 2009 b; Bian et al., 2012). One of the difficulties associated with the full SDPD is that it requires a formulation of the state equation indicated in Eq. (1.64)-(1.66). Espanol and Revenga (2003) performed SDPD simulations for an ideal gas and used the Sackur-Tetrode equation for , ,
E S N V with only the deterministic part in Eq.(1.63) and (1.64). The full thermodynamically consistent SDPD model has been validated by Vazquez-Quesada (2009 a) based on a Fourier problem that involves a fluid between two wall at different temperature. For liquids, Vazquez-Quesda (2009) derived a 41 state equation through a second-order Taylor expansion around the equilibrium state and applied it to viscoelastic fluid and colloidal suspension flow.
This work is part of a research effort at the CGPL at WPI to develop a SDPD-DV methodology and apply to the investigation of mesoscopic wall-bounded flows. There are three goals: First, to revise and further implement existing algorithms in SDPD-DV, as well as develop and implement new algorithms. Second, to validate and verify the algorithms of the SDPD-DV code by comparisons with experimental, analytical and numerical solutions. Third, to investigate thermodynamic properties and transport coefficients with applications of the SDPD-DV code to mesoscopic systems in equilibrium. The objectives and approaches are divided in three main categories. 1.
Revise and further implement existing algorithms in SDPD-DV, as well as develop and implement new algorithms in order to achieve a fully functional SDPD-DV methodology: a.
Revise the periodic boundary conditions algorithm and the implementation of the periodic boundary cells searching list. b.
Modify evaluation of smoothing function in order to include the contribution of the self-density for each fluid particle. c.
Modify and further implement the dynamic virtual particle allocation algorithm for the modeling of solid boundary in order to minimize the truncation error of density. d.
Develop and implement the boundary normal vector algorithm used in the reflection of the dynamically allocated virtual particles. 42 e.
Implement the algorithm for the contribution to the boundary force from the virtual particles. f.
Implement a temperature boundary condition in the dynamic virtual particle allocation model. g.
Implement a bounce-forward algorithm for a solid boundary with arbitrary shape and orientation. h.
Revise and rewrite portions of the Velocity Verlet integration method for the position and momentum equations. i.
Develop and implement a Runge–Kutta integration algorithm for the entropy equation. j.
Implement the artificial incompressibility method for modeling liquid flows. k.
Implement a temperature power law for the shear viscosity, bulk viscosity and heat conductivity that appear in the momentum and energy SDPD equations. l.
Develop and implement algorithms for the evaluation of transport properties such as diffusion coefficient, shear viscosity and heat conductivity based on mean square displacement (MSD) and velocity autocorrelation function (VACF). m.
Develop analytical formulas of the self-diffusion coefficient based on the SDPD-fluid following Litvinov et al. (2009). 2.
Validate and verify the SDPD-DV implementation: a.
Validate the implementation of boundary particle and fluid particle loading, dynamic virtual particle allocation, periodic boundary particle allocation, fluid field properties sampling, and bounce forward reflection. Compare with analytical solutions for transient body-driven Poiseuille flow of water between 43 stationary infinite parallel plates of 10 -3 m height. Calculate the components of the forces attributed to the dynamically allocated virtual particles and compared with previous DPD investigations. b. Verify moving and no-slip boundary conditions achieved by the dynamic virtual particle method by comparisons of SDPD-DV simulations with transient Couette flow of water between stationary infinite parallel plates of 10 -3 m height. c. Verify the ability of the SDPD-DV to simulate curved 3D solid boundaries, and the evaluation of pressure by comparisons of SDPD-DV results of steady, low-Re incompressible flow over a cylinder of radius 0.02 m with results from FLUENT. d.
Perform simulations of Poiseuille flow between two infinite parallel plates at different temperature for verification of the non-isothermal SDPD-DV. e.
Perform simulations of Couette flow between two infinite parallel plates at different temperature for verification of the moving wall boundaries with constant temperature. 3.
Investigate mesoscopic flows using the isothermal and non-isothermal SDPD-DV implementations. a.
Perform simulations of H O(l) and N (g) at equilibrium states. Evaluate the self-diffusion coefficient, shear viscosity, translational temperature, thermal speed, and thermal fluctuations. Validate and verify by comparisons with analytical and experimental values. b. Examine the scale dependence in SDPD-DV and characterize the effects of particle mass, particle volume, smoothing length, and time step. 44 The presentation of this work is organized in the following manner. In Chapter 2, the overview of the SDPD-DV methodology is presented, as well as the mathematical and numerical aspect, as pertaining to its implementation with dynamic virtual particle and non-isothermal model, is presented in detail for each aforementioned code modification or addition. In Chapter 3, an extensive set of benchmark tests that cover the hydrodynamic and mesoscopic regimes is presented and used for verification, validation and error analysis of the SDPD-DV method. The first verification of SDPD-DV involves comparisons with analytical solutions for a body-force driven, transient, Poiseuille flow of water between parallel plates of 10 -3 m height. The second verification test involves the transient, Couette flow of water between parallel plates of 10 -3 m height. The third test used for verification involves the low-Reynolds number, incompressible flow over a cylinder of radius 0.02 m. An extensive set of SDPD-DV simulations of liquid water and gaseous nitrogen in mesoscopic periodic domains is also presented. In Chapter 4 the SDPD-DV code with entropy equation is applied to microscale non-isothermal case studies. The validation and verification includes simulation of equilibrium gaseous nitrogen in periodic domains and the evaluation of transport coefficients. Fluctuations in thermodynamic variables are evaluated and compared with analytical estimates. Results from SDPD-DV simulation of non-isothermal nitrogen Couette flow are verified with results from FLUENT. Conclusions and recommendations for future work are presented in Chapter 5. 45
2. SDPD-DV METHODOLOGY AND IMPLEMENTATION
We review in this chapter the basic elements of the SDPD model and present the discrete SDPD equations. We then present the major algorithmic features of the SDPD-DV implementation developed for simulation of mesoscopic fluids in wall-bounded domains. Algorithms presented include the particle loading indexing; the neighboring particle search; density and force evaluation on interior domains, solid boundary and periodic boundaries; pressure evaluation; integration of SDPD-DV equations; evaluation of particle transport properties; evaluation of sample-averaged particle properties. Material in this chapter appear in Gatsonis et al. (2013).
The SDPD derivation of Espanol and Revenga (2003) starts with the Navier-Stokes equations written in the material derivative form (i.e. following a fluid particle) with independent field (Eulerian) variables the density ( , ) t r , velocity ( , ) t v r , and entropy ( , ) S t r instead of the ( , ) E t r , , , D t tDt r r v , (1.57) , 1, , , 3 D t PDt t t t v r v vr r r , (1.58) , , , DS tT TDt t t r r r . (1.59) The derivation arrives first at a discrete SPH-type set of the deterministic Navier-Stokes equations for i r , i v and i S . The domain containing fluid of total mass F M and volume T V , following SPH, is discretized with a number of F N points each one representing a 46 thermodynamic closed system (equivalent to a material volume or a fluid particle). Each particle i has constant mass Fi F
Mm N , (2.1) and is described by independent variables position, velocity and entropy at time t ( ) ( ), ( ), ( ) , ( ) ( ), ( ), ( ) , ( ) i i i i i xi yi zi i t x t y t z t t v t v t v t S t r v . (2.2) The SDPD discretization is consistent with a set of F N thermodynamics systems (or Lagrangian fluid particles). Alternatively, the SDPD discretization can be considered as a set of F N grid nodes which are moving with the material (or particle) velocity. These two views are identical and can provide preferable viewpoints during analysis. The number density i d of the SDPD particle i follows SPH and is defined as a summation of neighboring particles as shown in Figure 3 Figure 3. Contribution of neighboring particles in particle i support domain to the number density evaluation. , i i jj d W r r h , (2.3) i j
47 where , ( , )
W r h is the interpolant (smoothing) function with a finite support h satisfying normalization condition , 1 W r h dr . (2.4) Several options are available for the interpolant function and a detailed list of rules to construct smoothing functions can be found on Liu and Liu (2003). The sum in Eq. (2.3) is extended to all the j particles that are within the finite support h of particle i , including the contribution from the particles i itself. The distance between particles i and j is , , ij ijj ij iji x y z r r r , (2.5) A dimensionally consistent particle volume i V is defined as the inverse of the particle number density i i V d . (2.6) The field density ( , ) t r at the particle position ( ) i t r of particle i is obtained as [ ( )] ( ) ( ) i i i i t t m d t r . (2.7) Similarly to the SPH approach, the discrete value at a point i r (or position of a particle) for any field variable ( , ) t r can be calculated by interpolation ( )( ) ( ) i j jji i jj W r rW r r r . (2.8) The deterministic SDPD equations consistent with the Naveir-Stokes Eq. (1.57)-(1.59)are: i i ddt r v , (2.9) 48 j ij iji i ij ij ij ij ij ijj j ji j i j i j P F Fd Pm Fdt d d d d d d v r v e e v , (2.10) ijii iji j i j FdST Tdt d d . (2.11) This deterministic dynamics is introduced into the GENERIC framework of Ortega (2005). This is done in order to introduce the stochastic part for the system dynamics. The introduction of thermal fluctuating terms leads to the GENERIC stochastic differential equations that contain (reversible and irreversible) deterministic terms, a term that relates the dissipative (irreversible) dynamics with stochastic terms, and a stochastic term (Eq. (38) in Espanol and Revenga, 2003) as B E Sdx L M k M dt dxx x x , (2.12) In the above dx is the stochastic term. The term L E x is the reversible part of the dynamics, and the second term
M S x is defined as the irreversible part. The matrices L is anti-symmetric, and M is symmetric and positive semi-definite. The following conditions must be satisfied by L and M
0, 0
S EL Mx x , (2.13)
0, 0
I E I SL Mx x x x . (2.14) The stochastic term dx satisfies T B dxdx k Mdt (2.15) which is an expression of the fluctuation-dissipation theorem . It is also required that
0, 0
E Idx dxx x , (2.16) 49 The term B k M x is from the stochastic interpretation of ˆIto process. Therefore, the matrix M is given by Ti j i jB Bij Ti j i jB B d d d dSk dt k dtM dS d dS dSk dt k dt v v vM v (2.17) The matrix L is defined by ijij ij L m
1L 1 (2.18) The derivatives of the energy and entropy with respect to the state variables result in
0, 01 ii ii
PdE Smx xT v (2.19) and the thermal noise is defined as ii dx ddS v (2.20) Eq. (2.9) (2.11) can be given in a matrix form as i i j i ji ii j jj B j Bi ij ij Bj j jii ii j i ji j jj B j B P d d d dS dd k dt S k dtk dtmS dSdS d dS dST dtk dt S k dt v v vr vvv L Mv vv (2.21) 50 The final result is a set of discrete particle equations that describe the deterministic and stochastic dynamics, referred to as the SDPD equations (Eq. (63) in Espanol and Revenga, 2003). We rewrite below the SDPD equations for the independent variables i , i v and i S in a format that is conducive to numerical implementation and allows also direct comparison with SPH and DPD. They are given by i i d dt r v , (1.62) d = d d i i C D R m t t v F F F , (1.63) i i Vi C R T dS E dt E dt E . (1.64) The terms appearing in the momentum equation (1.63) are broken into the conservative, dissipative and velocity random terms given by jiC ij ijj i j
PP Fd d F r , (2.22)
51 3 1 5 3 i j ijB BD ijj i j i ji ji j ijB B ij ij ijj i j i ji j
T T Fk kC C d dT TT T Fk kC C d dT T F ve e v , (2.23)
58 358 83 3 i j ijB iji j i jR ijj i j ijB iji j i j
T T Fk dT T d dT T Fk tr dT T d d WF eI W . (2.24) The terms in entropy equation (1.64) are divided into viscous, conductive and random terms. 51
22 2 i j j ijB B BV ijj i j i j i i ji jij i j ijBij ij ji j i j i j
T T T Fk k kE C C T T C d dT TF T T Fkd d m T T d d ve v , (2.25) ij ijBCd ij ijj ji j i i j
F FkE T Td d C d d , (2.26)
58 31 :2 58 83 3 + 4 i j ijB iji j i jR ij ijj i j ijB iji j i jijB i j ijj i j
T T Fk dT T d dE T T Fk tr dT T d dFk T T dVd d W e vI W . (2.27) A state equation is required to close the system ( , , ) eqi i i
E E N V S . (2.28) The state equation provides the particle temperature and pressure as eqi i ET S , (2.29) eqi
EP V . (2.30) The term ij i j ij e r r r in the SDPD equations is a unit vector, ij i j v v v , ij i j T T T and ij i j
F F W r r r r . (2.31) The terms and represent the shear and bulk viscosity; i C is the heat capacity at constant volume of particle i and is an extensive property, i V i C c m where V c is the specific heat capacity of fluid; the Boltzmann constant is B k and is the thermal conductivity of the fluid. 52 The velocity and entropy random terms contain I the identity matrix and ij d W the traceless symmetric part of a matrix of independents increments of the Wiener process ij d W .
12 3
Tij ij ij ij d d d tr d
IW W W W . (2.32) where the trace is defined as ij ij tr d d W W . (2.33) For B i k C the fluctuating terms reduce to zero and the set of Eqs.(1.62)-(1.64) reduce to a discrete SPH-form of the deterministic Navier-Stokes equations for i r , i v , and i S given by Eq. (2.9)-(2.11). Figure 4. General flow chart of SDPD-DV and post-processing.
Perform NNPSEvaluate FP DensityEvaluate FP Pressure and C p Evaluate FP Forces and S i Integration VV for r i and v i , RK for S i i =1, N F FPOutput Particle SamplesGenerate GridFPs Near PB Perform PBP MethodFPs Near SB Perform DVPAFPs Away from BoundariesInputs Generate Particle-Search Grid Load & IndexBPs & FPsFPs Near PB Perform PBP MethodFPs Near SB Perform DVPAFPs Away from Boundaries Penetrate PB Perform ReinjectionPenetrate SB PerformBFNo PenetrationPost-process for Field Properties Post-process for Particle & Transport PropertiesInitialize FPs Entropy In this section, we present the mathematical and computational aspects of the SDPD-DV developed in this work for simulation of wall-bounded domains. The general flow chart is presented in Figure 4. The physical domain is shown in Figure 5 and is characterized by arbitrary external and internal solid-wall boundaries, planar periodic inlets and outlets with equal areas, and the fluid region. A rectangular domain is also constructed to aid in the numerical implementation of the particle search and includes the physical domain as shown in Figure 5.
Figure 5. Physical wall-bounded domain with an exterior wall and an interior solid body showing the fluid particles (FPs), boundary particles (BPs), and virtual particles (VPs) used in the SDPD-DV. Inlets and outlets are considered periodic. The large rectangular cells are used for nearest neighbor particle search (NNPS). The boundary of the physical domain is discretized with a surface triangulation based on a surface length scale small enough to resolve the physical characteristics of the surface and much smaller than the smoothing length, h . A total of B N boundary particles ( BPs ) are placed on the vertices of the surface grid as shown in Figure 5. The BPs position and velocity are stored in the global particle list (0), (0), (0)( ) (0) 1,( ) (0) (0), (0), (0) k k kk k Bk k kx ky kz x y zt k Nt v v v r rv v (2.34) In addition, the normal vector entering the fluid domain at each BP is evaluated and later used in fluid particle-wall interaction. For a surface defined implicitly by ( , , ) 0
F x y z the normal vector is evaluated analytically by ( , , ) ( , , ) k k k k k k x y z F x y z n . For an arbitrary solid wall surface the local normal to each surface triangle surrounding a vertex is evaluated as n r r and the closest is assigned to the vertex as shown in Figure 6. Figure 6. Normal vector for a BP on a solid boundary.
The computational volume of the physical domain is populated with F ( 0) N t fluid particles ( FPs ) each with mass i m , corresponding to a total mass F M following Eq.(2.1). The FPs are 55 associated with a global particle list with index
B B F
1, 1 i N N N . The FPs are assigned upon initialization with position and a velocity as
B B F (0) (0), (0), (0) 1, 1(0) (0), (0), (0) i i i ii ix iy iz x y z i N N Nv v v rv (2.35) To simplify FP loading in cases of complex physical geometries, a mesh for the fluid region is generated with the number of tetrahedrons vertices to be equal to the total number of FPs in the domain.
The evaluation of the properties of a fluid particle i ( FP i ) require the identification of the nearest neighboring particles (NNP) in its support domain shown in Figure 7(c). A fluid particle j is considered as a nearest neighbor of i when it is located within the smoothing domain of particle i , and therefore ij r h , with determined by the smoothing function. An all-pair NNP search among F N particles in a domain requires ( ) F N operations per computational cycle. Numerous approaches have been developed for various types of particle simulations that require such a search. We implemented in this work, a NNPS approach that uses concepts from the linked-cell (Hockney and Eastwook, 1981), (Putz, 1998) and the Verlet-list algorithm (in't Veld, Plimpton and Grest, 2008). The approach uses a particle-search grid with size , , X Y Z
L L L as shown in Figure 7(c). This search-related grid contains C N rectangular cells generated with lengths , , max( ) CX CY CZ
L L L h . In case of a periodic (physical) boundary the search-grid face has to align with the physical boundary, as shown in Figure 5(a). The steps for the NNPS algorithm are summarized in Table 1. 56 Figure 7. Rectangular cells used for NNPS in the interior, near a solid and near a periodic boundary. Table 2. Nearest Neighboring Particle Search Algorithm used in SDPD-DV
S-1.
Update in the global particle list the current properties of each FP
1, 1
B B F i N N N . S-2.
For each cell C C N identify the BPs and
FPs contained in it. Generate C N cell particle lists each with ( ) FC N t
FPs and ( ) BC N t
BPs respectively. S-3.
For each cell particle list C C N loop through the ( ) FC N t
FPs and search all
FPs and all
BPs in the neighboring 27 cells. For each FP i create a list with ( ) i J t NN FPs and a list with the ( ) i K t NN BPs
This particle-search grid is used also in the implementation of periodic boundary conditions. For a FP residing in a cell with a periodic boundary face the NNPS generates the cells on the corresponding periodic boundary, as shown on Figure 7(c). All the FPs in FP i s support domain from the corresponding cells are indexed and become available for density and force evaluation. 57 The density evaluation procedure depends on the location of the FP . The three cases considered in this work include a FP away from a boundary (Figure 8(a)), a FP near a solid boundary (Figure 8(b)) and a FP near a periodic boundary (Figure 8(c)). (a) FP away from boundaries. (b) FP near a solid boundary. (c) FP near a periodic boundary. Figure 8. Support domain for a fluid particle (FP) for three cases considered in SDPD-DV. (a) FP Away from Boundaries
The particle number density (number of FPs in unit volume) at the location of the FP i is defined by the summation over the neighboring FPs in its support domain shown in Figure 8(a) and provides the thermodynamic volume i V , ( )
1( ) ( ) ( ) , , 1,( ) i L ti i j Fj i d t W t t h i NV t r r (2.36) where i L t is the number of particles in FP i ’s support domain. The summation includes the contribution from the FP i itself. Consequently, the mass density at the position i r is given by 58 ( ) [ ( )] ( ) ( ) ii i i i i mt t m d t V t r (2.37) This equation along with the normalization condition Eq. (2.4) satisfies the Lagrangian form of the continuity equation and therefore ensures that mass is conserved in a closed or in a periodic domain ( Espanol and Revenga, 2003). An alternative approach could be based on the continuity equation following SPH by Liu and Liu (2003). A common choice in SPH and SDPD for the interpolant is Lucy’s (1977) smoothing function r rW r h h h h (2.38) which, through Eq. (2.31) gives
315 1 .4 ijij rF h h (2.39) (b) FP Near Solid Boundaries: Dynamic Virtual Particle Allocation Method
When the support domain of a fluid particle falls outside a boundary (solid or periodic as shown in Figure 8(b) and Figure 8(c) respectively) the smoothing function approximation of Eq. (2.36) will result in error due to the with the unaccounted (truncated) part of the support domain. This boundary truncation error would lead to incorrect density evaluation and would affect the particle momentum and entropy evaluation. This issue has been addressed in SPH (Morris et al., 1997; Randles and Libersky, 1996; Liu and Liu, 2003; Li and Liu, 2007) and SDPD (Litvinov et al., 2008; Bian et al., 2012; Vazquez-Quesada, 2009). With the Lucy smoothing function Eq. (2.38) this occurs when iw r h , where iw r is the distance between the particle and the boundary as shown in Figure 8(b). It should be noted that due to the summation density approach followed in SDPD the density error does not affect the overall mass conservation but affects the total volume. 59 We describe the density evaluation method developed in this work to compensate for the boundary truncation error. For a FP near a solid wall we developed and implemented the dynamic virtual particle allocation (DVPA) method, by which the truncated portion of the support domain is dynamically filled with virtual particles (VPs). Approaches based on static and dynamic ghost particle allocation appeared in SPH (Morris et al., 1997; Randles and Libersky, 1996; Liu and Liu, 2003; Li and Liu, 2007) and SDPD (Litvinov et al., 2008; Bian et al., 2012). These virtual particles are included in the summation density of Eq. (2.36) following the algorithm described in Table 3. The VPs are generated as mirrored images of the FP i and its neighbors as shown in Figure 9(a). (a) Support domain for a FP near a solid boundary represented by BPs. (b)
Reflection of a FP near a solid boundary to create a virtual particle (VP).
Figure 9. Dynamic virtual particle allocation (DVPA) method used in SDPD-DV for density and force evaluation. A virtual particle (VP) is generated for a fluid particle (FP) near a solid boundary represented by boundary particles (BP).
60 In case of curved walls, we require that min h aR (2.40) where, the coefficient a and min R is the smallest radius of the curvature for the solid boundary in FP’s support domain. With such a constraint it is possible to consider the surface of the wall contained within the smoothing domain of a FP as quasi-flat and its normal vector as locally constant. The algorithm is summarized in Table 3. Table 3. Dynamic Virtual Particle Allocation and Density Evaluation Algorithm used in SDPD-DV.
S-1.
Loop through each FP i . S-2. Compute the particle number density i d using the FP j neighbors from Eq. (2.36) and FP j itself. S-3. Search the nearest-neighbor BP list to find the closest BP k . Compute ( ). iw i k k r r r n and compare with the average minimum distance FPs l between the FPs. A choice is min FPs l l provided by the surface triangulation. If min iw r l create a virtual particle VP i with its position VP r reflected across the normal to the surface plane, VP i iw k r r r n . The process is shown in Fig. 4(b). S-4. If VP i h r r add to particle number density i d the contribution due to the VP . S-5. Loop through each of the nearest-neighbors FP j and create a mirrored VP j following S-3. S-6. If VP j h r r then add to particle number density i d the contribution due to the VP j . 61 (c) FP Near Periodic Boundaries: Periodic Boundary Particle Allocation Method The density of a FP containing a periodic boundary within its support domain consists of contributions from the surrounding FPs and the FPs in the truncated part as shown in Figure 8(c). We implemented a periodic boundary particle allocation (PBPA) method, by which the truncated portion of the support domain is filled with copied particles from the periodic cells identified during the NNPS (Sec. 2.2.3). The density of the FP is then obtained using Eq. (2.36) where the ij r for a copied FP is given by ij Px ij Py ij Pzij x L y L z Lr (2.41) The values for , , Px Py Pz
L L L are determined by the location of the periodic boundary at the specified cell face. If there is only one periodic boundary at a cell face, as is the periodic outlet in Figure 5, , 0, 0
Px X Py Pz
L L L L . In a case of fully periodic boundaries, a corner cell has three periodic boundaries and , ,
Px X Py Y Pz Z
L L L L L L when copying the FP from the corresponding corner.
The full algorithm of SDPD, which include entropy equation, needs a closure by equation of states that provides pressure , i i i P P S and temperature , i i i T T S given by Eq. (1.65) and (1.66). Therefore, the internal energy is required as a function of i m , i S , and i V , in the form of , , i i i i E E m V S . (2.42) For ideal monatomic gases we consider the Sackur-Tetrode entropy equation given by Tetrode (1912) or Laurendeau (2005) 62 B V mES Nk N Nh , (2.43) where N is the number of molecules in the system, h is the Planck’s constant. For an SDPD fluid particle i which is considered as a thermodynamic subsystem, following Eq. (2.43), the entropy is given by i iii i B ii m EVS N k N hN , (2.44) where i N is the number of molecules in FP i , with volume i V . Then the internal energy of FP i is given by iii i i i i Bi i Sh NE N V S N kmV , (2.45) Therefore, the temperature and pressure of monatomic ideal gas are given by Espanol and Serrano (1999) ii i i i ii i B
ET E N V SS N k (2.46) i ii B ii i
E NP k TV V (2.47) And the heat capacity is given by i i B C N k (2.48) To initialize the entropy of FP i , we follow (2.44) by i B i Bi i N k mV kS Th N (2.49) 63 In Chapter 3, we consider also the isothermal SDPD implementation, where entropy equation is decoupled from the system, and temperature is therefore kept constant and initially assigned to the FPs. Closure in this case is obtained by an equation of state that provides ( , ) i i i
P P T . For an ideal gas flow we use i i i P RT (2.50) For a perfect gas with constant specific heats, the pressure is defined as (Batchelor, 1967) P e (2.51) For an incompressible flow, we follow Batchelor, (1967)
220 0 ii cP (2.52) where, c is an artificial sound speed and is the initial density. A similar law is given by Morris et al. (1997) and Liu and Liu, (2003), i i P c (2.53) The value for c should be low enough to avoid using very small time steps and high enough to be sufficiently close to the real fluid. As suggested by Morris et al. (1997) and Liu and Liu (2003), c can be estimated by
22 0 0 0 0 00 max , ,
V V F Lc L (2.54) where, is the relative density perturbation, V is a velocity scale, L is a length scale, is the kinematic viscosity, and o f is a body force per unit mass. Several alternative forms of Eq. (2.52) could be found as 64 P P (2.55) P B (2.56) P P b (2.57) where B and b are the coefficient chosen to keep the density fluctuation. In chapter 4, the full algorithm of SDPD-DV will be tested and validated. The force evaluation procedure depends on the location of the FP . As with density we considered three cases depicted in Figure 8. (a) FP Away from Boundaries For a FP with a position away from a boundary, as shown in Figure 8(a), the forces C F , D F and R F are evaluated directly from Eq. (2.22), (2.23) and (2.24). The summation is performed over all FPs in the support domain of the FP i provide by the NPPS. (b) FP Near Solid Boundaries: Dynamic Virtual Particle Allocation When a FP is in the proximity of a solid boundary ( iw r h ) a correction term is required also in the momentum Eq. (1.63) due to the presence of the wall and truncated domain. This correction term is based only on the local thermodynamic state of the fluid and the local geometry of the wall described by the BPs. For force evaluation, we use the DVPA method outlined for density in Sec. 2.2.4(b). and develop additional force terms to supplement the SDPD equations (2.22)-(2.24). 65 The thermodynamic properties of each VP are assigned to be identical to the properties of the corresponding FP that is mirrored, specifically _ j V j , _ j V j P P . (2.58) For each FP in the proximity of the wall the conservative C F in Eq. (2.22) becomes = j VPiC ij ij Cj N i j PP Fd d F r F (2.59) where the conservative force due to the VPs is _ _ _2 21: _ VP j VVP iC ij V ij Vj N i j V PP Fd d F r (2.60) In the above equations, N and V N are respectively the number of FP s and VPs contained within the support domain of a FP i . Generally, V N N , _ ij V ij F F and the identity sign is achieved only when the FP i is located exactly on the wall ( wall r ). Consistent with the above methodology we also introduce an additional term to the dissipative part of the dynamics to account for the truncated part of the domain. For each FP in the proximity of the wall the D F in Eq. (2.22) becomes
21: 21:
5= 1 3 1 5 + 3 i j ijB BD ijj N i j i ji ji j ij VPB B ij ij ij Dj N i j i ji j
T T Fk kC C d dT TT T Fk kC C d dT T F ve e v F (2.61) where the dissipative force due to the VPs in the truncated domain is 66 _21: _ _ _21:
5= 1 3 1 5 3 VV i j ijVP B BD ij Vj N i j i ji ji j ijB B ij V ij V ij Vj N i j i ji j T T Fk kC C d dT TT T Fk kC C d dT T F ve e v (2.62) For a given wall velocity W v the velocity of the VPs required to compute the velocity vector _ ij V v is defined as _ j V W j i j v v v (2.63) Substituting in the definition of the velocity vector _ ij V v it become _ _ ij V i j V i j W i j v v v v v v (2.64) The above formulation is sufficient to impose the non-slip condition at a solid wall boundary. For a planar wall and a FP i located on the wall _ _ , t t n nij V ij ij V ij e e e e , Eq. (2.64) enforces i W v v . For a wall without heat flux, we assign the VP temperature equal to the wall temperature given by _ j V w T T (2.65) If the temperature gradient along wall surface is not equal to zero, the w T would be the temperature of the nearest BP of FP j . Whereas if the heat flux through the surface of the wall is not zero, the VP temperature is assigned by _ j V w j T T T (2.66) (c) FP Near Periodic Boundaries
The force evaluation of a FP i containing a periodic boundary within its support domain follows the approach discussed in the density evaluation (Sec. 3.4.3). The FPs in the truncated part are 67 copied from the corresponding periodic cells. The force on the FP i is then obtained using Eq. (1.63) and Eq. (1.64) where the ij r for a copied particle is given by , , j Px j Py j Pzij i i i x L y L z Lx y z r (2.67) Once the particle neighbor list is constructed and the density, pressure and force is evaluated the integration of motion and momentum Eqs. (1.62)-(1.64) proceeds as indicated in Figure 4. The integration scheme used in this work is an implementation of the Velocity-Verlet scheme (Nikunen et al., 2003) for momentum equation, coupled with Runge-Kutta scheme for entropy equation, and summarized by the algorithm in Table 4. The choice for is based on required accuracy in the fractional change of the velocity estimate. The required t follows standard SPH conditions (Morris et al., 1997; Li and Liu, 2007) and once choice provides a posterior check, after the integration has proceeded, max ht v (2.68) (b) Near Solid Boundaries: Reflective Bounce Forward Condition In general the additional repulsive force Eq. (2.59) exerted on the FP i from the VPs is not enough to prevent FP i penetrating the solid wall. When such event occurs the FP i is reinserted in the fluid domain by a bounce forward reflection, this takes place after S-3 in the integration algorithm Table 4. The bounce forward reflection method is depicted in Fig. 5 and summarized by the algorithmic steps in Table 5 which are embedded in S-3 of Table 3. 68 (c) FP Near Periodic Boundaries In case that a FP penetrates a periodic boundary it is re-injected into the related periodic cell. This operation follows the bounce-forward method.
Table 4. Time integration in SDPD-DV. Velocity-Verlet is used for the particle position and momentum equation, and Runge-Kutta is used for the entropy equation.
S-1.
Calculate , n nC i i T F r , , , n n nD i i i T F r v , , n nR i i T F r and , , n n nVi i i i
E T r v , , n nCd i i E T r , , , n n nR i i i E T r v . S-2.
Update n ni i C D R t tm v v F F F , and calculate n Vi Cd Ri n n ni i i E dt E dt EdS T T T , then update n n ni i i T T f dS . S-3.
Update n n ni i i t r r v . S-4.
Calculate , n nC i i T F r , , n nR i i T F r , and , , n n nR i i i
E T r v . S-5.
Update n ni i C R tm v v F F , and n ni i T T . S-6.
Calculate , , n n nD i i i T F r v , , , n n nVi i i i E T r v and , n nCd i i E T r . S-7.
Update n ni i D tm v v F , and calculate n Vi Cd Ri n n ni i i E dt E dt EdS T T T , then update n n ni i i
T T f dS . S-8. If n ni ini v vv , update n ni i T T , loop over step 6. Calculate , , n n nD i i i T F r v , , , n n nVi i i i E T r v , , n nCd i i E T r and go to step S-2. 69 (a) Position reflection of a penetrating FP. (b)
Velocity reflection of a penetrating FP.
Figure 10. Reflective bounce-forward method in SDPD-DV applied in a case where a fluid particle (FP) penetrates a solid boundary represented by boundary particles (BP).
Table 5. Bounce-forward Method used in SDPD-DV in case of particle penetration.
S-1.
After the FPs position has been updated to the new time step n+1, check the sign of the scalar product ( n ni i P P n ), where n ni i P P is the vector between the FP n+1 position and the closest BP of n position, n is the wall normal unit. If n ni i P P n , particle penetrates, and proceed to S-2. S-2. If FP i penetrates the wall, compute the scalar product between the vector n ni i P P and the unit normal wall vector n obtaining the distance n np i i l P P n . S-3.
Compute the distance between the wall surface and FP i new position (inside the wall) ni P as p W l l r . S-4.
Reflect the FP i inside the domain and compute its new position as n niBF i l P P n . S-5.
Impose the velocity component of FP i normal to the wall to be opposite in sign. This is achieved by imposing the velocity of FP i after the reflection to be n n niBF i i v v v n . A property of a FP from the SDPD-DV simulation at a discrete time t k t is designated as ( ), ( ), ( ), ( ) ki i i i i
X t t S t T t r v etc. An SDPD-DV output (or sample) consists therefore of all particle properties , 1, ki F X i N . For steady SDPD-DV simulations we gather after reaching steady-state, m M independent SDPD-DV samples. For unsteady simulations, a number of M individual runs are performed in order to generate a sufficient number of independent samples. Transport coefficients such as diffusivity, shear viscosity are dynamical properties desired from SDPD simulations. The self-diffusion coefficient D can be evaluated for a system with dimensionality D d from the mean-square displacement (MSD) through the Generalized Einstein formula as shown by Allen and Tildesley (1997)
20 0 lim 2 i iD t tD d r r (2.69) where, is the delay time nK t , n is the time interval of each particle sample, and K the number of particle samples in delay time which is k k n t t t and k K k t t . The schema is shown in Figure 11. It can also be evaluated by the Green-Kubo formula through the velocity autocorrelation function (VACF) as introduced by Allen and Tildesley (1997) i iD D d t td v v (2.70) The term in the brackets in Eq. (2.69) and Eq. (2.70) denotes the time correlation function for time-dependent signals ( ), ( ) A t B t (Haile, 1997) 71
1( ) ( ) ( ) lim ( ) ( ) t ot t A t B t A t B t dtt (2.71) The function ( ) A t is sampled at t and ( ) B t after a delay time . The integral is evaluated over many time origins t shown in Figure 11. For a given number of time samples M and a number of particle samples K in delay times the time correlation can be approximated as a summation over M K number of available time origins (Haile, 1997)
1( ) ( ) ( )1
M Km t A t B tM K (2.72) The calculation can be improved by averaging over all particles in each sample,
10 0 0 01 1
1( ) ( ) ( ) ( )1 P NM K i im iP
A t B t A t B tNM K (2.73) Then the self-diffusion coefficient based on MSD is provided by P NM K i im iD P
D t td t NM K r r (2.74) And the self-diffusion coefficient based on VACF is given by P NM K i im iD P
D dt t td MN v v (2.75) The integral is approximated with a summation over all the discrete time in delay time shown by P NnK t M K i it t m iD P t t t tD d MN v v (2.76) In this work, we implemented the VAC algorithms following Haile (1997) and the MSD algorithm following Rapaport (1995). 72
Figure 11. Sampling used inn SDPD-DV for evaluation of the self-diffusion coefficient. The delay time is τ . The total number of samples is M , and K is the number of samples in delay time . The time origins are indexed from m= M-K+ (b) Shear Viscosity
Another transport coefficient we studied in this work is the bulk viscosity . An expression analogous to the Einstein diffusion Eq. (2.69) lead to the viscosity expression given by Allen and Tildesley (1987)
20 0 0 0 i yi xi i yi xix y i iB m v t r t m v t r tk TV , (2.77) where x y denotes a sum over the three pairs of distinct vector components ( xy , yz , and zx ) which improve the statistics. The alternative Green-Kubo form is given by Allen and Tildesley (1987) xy xyx yB V P t P dtk T , (2.78) where xy j xj yj xij yijj i j P m v v r fV (2.79) 73 The second term in xy P can be fulfilled with the force computation, and for Lennard-John potential the weighting function is ij ij ij ij f r r f r , which leads to xy yx P P . The numerical implementation follows section 2.2.8.(a) as shown in Figure 11. For analysis and visualization we construct also Eulerian (field) properties , , , , , t S t T t
V r r r . To obtain such properties we first generate a tetrahedral mesh over the fluid domain of interest. The vertices of the tetrahedron are designated as the nodes of the domain, each with an assigned index d d G and coordinates ( , , ) d d d d x y z r . An instantaneous fluid property associated with a node d at t k t , is obtained using the smoothing function approximation Eq. (2.8) shown in Figure 12, as ( )( , ) ( ) ( ) m kd i im m k m k id d d ii W r r XX t X d X W r r r (2.80) The summation is extended over all the FPs within a set distance d l from the vertex d r . For unsteady simulations, we obtain after reaching steady-state the sample-averaged field property at a node d at t k t , is given by ( , ) k m kd d dm X t X X r (2.81) For steady simulations, we obtain after reaching steady-state the sample-averaged field property at a node d is given by ( ) m kd d dm X X X r (2.82) 74 The choice of the Eulerian grid size and thus the interpolating distance d l must be compatible with the length scale of the phenomena under consideration. An excessively coarse (Eulerian) grid or interpolating distance could smooth out small scales of interest. (a) A typical post-processing grid. (b)
Sampling region near vertex of the grid.
Figure 12. Post-processing grid used for evaluation of field (Eulerian) instantaneous and time-averaged properties from SDPD-DV particle samples.
3. VERIFICATION, VALIDATION, AND ERROR OF THE ISOTHERMAL SDPD-DV
In this chapter we utilize the isothermal SDPD-DV code and perform verification and validation tests that cover the hydrodynamic and mesoscopic regimes. The first verification of SDPD-DV involves comparisons with analytical solutions for a body-force driven, transient, Poiseuille flow of water between parallel plates of 10 -3 m height. Physical insights on the force method in SDPD-DV are obtained by calculating the components of the forces attributed to the dynamically allocated virtual particles and compared with previous DPD investigations by Altenhoff et al. (2007) and Fedosov et al. (2008). The second verification test involves the transient, Couette flow of water between parallel plates of 10 -3 m height. The third test used for verification involves the low-Reynolds number, incompressible flow over a cylinder of radius 0.02 m. This benchmark test has been used in SPH by Morris et al. (1997), Vazquez-Quesada and Ellero (2012) and in DPD simulations by Kim and Phillips (2004). Our SDPD-DV results are compared with those obtained from ANSYS FLUENT (FLUENT . The final set of benchmark tests involves calculation of equilibrium states for liquids and gases in mesoscopic domains. This extended set of SDPD-DV simulations evaluates the translational temperature for liquid water and gaseous nitrogen, and the self-diffusion coefficient of liquid water. The SDPD-DV results are compared with analytical expressions introduced by Litvinov et al. (2009), Bird et al. (2007) and experiments Holz and Sacco (2000). The SDPD-DV liquid water simulations examine also the scale effects on the self-diffusion coefficient and the Schmidt number by varying the mass and size of the fluid particles (Vazquez-Quesada et al., 2009), (Litvinov et al., 2009). The material from this chapter can be found in Gatsonis et al. (2013). 76 The first test case involves an incompressible Poiseuille flow between two stationary infinite plane parallel plates as shown in Figure 13. The velocity profiles as predicted by SDPD-DV are compared to theoretical formulations. The SDPD-DV density is plotted along with the standard deviation. In addition, the boundary forces due to virtual particles are evaluated, and the distributions are plotted. This test case verifies the ability of the SDPD-DV to minimize density fluctuations and enforce the no-slip condition on solid boundaries.
The test involves an incompressible Poiseuille flow with density across two infinite parallel, stationary walls with separation height Z L as depicted in Figure 13. The fluid considered in the SDPD-DV simulations is H O(l) with -3
1, 000 kg m a and
10 kg m s . (a) 3D view (b) X-Z plane 2D view. Figure 13. SDPD-DV simulations of transient, body-force driven, planar Poiseuille flow. Physical domain showing the BPs and FPs. Table 6: Input parameters used in SDPD-DV simulations of transient planar Poiseuille flow, transient planar Couette flow, and flow over a cylinder.
Input Parameters Case Transient Poiseuille Transient Couette Flow Over Cylinder (m) X L -3 -3 (m) Y L -4 -4 (m) Z L -3 -3 m R N/A N/A 0.02 -2 (ms ) x f -4 N/A 1.5×10 -7 kg F M -7 -7 (K) T
300 300 300 -3 kg m a kg m s -3 -3 -3 -1 -1 kg m s
0 0 0 -1 -1
J kg K V c -1 -1 W m K -1 ms xw V N/A 1.25×10 -5 N/A F N B N C N
300 300 1323 (m) h -5 -5 -3 -2 (ms ) c -4 -4 -4 (s) t -4 -4 -2 min m l -6 -6 -4 m d l -4 -4 -3 M
50 50 100 The physical domain has
10 m X L , Y L ,
10 m Z L . The total mass of F M in the domain is represented by 9,000 fluid particles with a constant temperature of
300 K i T . The solid walls are represented by 8,282 boundary particles as shown in Figure 13. Periodic boundary conditions are imposed along the -axis x and -axis y . 78 Closure for pressure and density is obtained by using the artificial compressibility relation Eq. (2.52). Integration is carried out following Sec. 2.2.7 and a constant body force of
10 m/s f is imposed on each fluid particle. The DVPA method (Sec 2.2.4(b) and 2.2.6(b)) is applied for density and force evaluation. Over M simulations were performed to generate instantaneous samples used to derive the instantaneous particle properties. Fluid properties were then sampled on a structured grid with edge length
10 m d l . Input parameters in the SDPD-DV are listed in Table 6 (Gatsonis et al., 2013). The classic Poiseuille flow is due to an applied pressure gradient which can be replaced by a body force per unit mass f (ms -2 ) parallel to the x-axis . The flow starts from rest and develops a velocity profile given by Morris et al. (1997) and Papanastasiou et al. (1999)
22 233 20
Zx Z n Z Z zf fL nV z z L tz t nL Ln (3.1) The sample-averaged field density ( , ) d t r and velocity ( , ) d v t r profiles are plotted in Figure 14 along with the standard deviation. The SDPD-DV density ( ) d Z shown in Figure 14(a) has an error a a of less than 4% in the interior and less than 5% near the wall boundary. The field density averaged over the entire channel,
966 kg/m is used Eq. (3.1) to evaluate the analytical velocity profile. Comparisons between the analytical and SDPD-DV velocity profiles across the channel is plotted in Figure 14(b) and show them to be in excellent agreement for all times considered. The density and velocity results demonstrate the ability of our DVPA-based density and force evaluation method as well as integration algorithm implemented in SDPD-DV. 79 (a) Density distribution with standard deviation. (b) Velocity distribution. Figure 14. Sample-averaged fluid density ρ(Z d ) and velocity Vx ( Z d ,t) from SDPD-DV simulations of transient, body-force driven, planar Poiseuille flow. The domain has L X =10 -3 m, L Y =3×10 -4 m, L Z =10 -3 m, the fluid is H O(l) with ρ a =1,000 kg/m , T i =300 K, η =10 -3 kg·m -1 s -1 , and f x =10 -4 m/s . The analytical profiles are plotted for verification (Morris et al. 1997). To further investigate the DVPA method and its ability to enforce the no-slip boundary conditions we evaluate the boundary forces due to the virtual particles using Eq. (2.60) and Eq. (2.62). The distribution of the virtual particle force components are shown in Figure 15. As shown in Figure 15(a) and Figure 15(b), the normal components of conservative and dissipative VP forces are decreasing as the distance from solid wall increases. The negative and positive sign of the conservative normal virtual particle forces shows that they are perpendicular to the plates and pointing to the inner domain. Figure 15(c) and Figure 15(d) shows the tangential component of the conservative and dissipative virtual particle forces. Both contribute to the imposition of the non-slip condition. Our results show that the forces due to virtual particles in 80 SDPD-DV have the same qualitative behavior as the DPD boundary forces [Fig. 4, Altenhoff et al., 2007] and to the average wall forces [Fig. 4, Fedosov et al., 2008]. (a) Conservative normal component. (b) Dissipative normal component. (c) Conservative tangential component. (d) Dissipative tangential component.
Figure 15. Forces due to virtual particles from SDPD-DV simulations of transient, body-force driven, planar Poiseuille flow. The domain has L X =10 -3 m, L Y =3×10 -4 m, L Z =10 -3 m, the fluid is H O(l) with ρ a =1,000 kg/m , η =10 -3 kg·m -1 s -1 , and f x =10 -4 m/s . The test involves an incompressible flow across two infinite parallel walls with the top wall moving with a constant velocity as shown in Figure 16. The velocity profiles as predicted by SDPD-DV are compared to theoretical formulations. The SDPD-DV density is plotted along with the standard deviation. This benchmark test is used as further verification of the SDPD-DV and its ability to enforce no-slip and no penetration on a moving solid boundary.
The test involves an incompressible flow with density across two infinite parallel walls as depicted in Figure 13 with the top wall moving with a constant velocity xw V .The SDPD-DV simulation considers H O(l) with -3
1, 000 kgm a and
10 kg m s . The physical domain has
10 m X L , Y L , and
10 m Z L . The total mass of F M in the domain is represented by 9,000 fluid particles with a constant temperature of
300 K i T . The solid walls are represented by 8,282 boundary particles as shown in Figure 16 each with xw xw v V . Periodic boundary conditions are imposed along the -axis x and -axis y . Closure for pressure and density is obtained by using the artificial compressibility relation Eq. (2.52). Input parameters are listed in Table 6. 82 (a)
3D view (b)
X-Z plane 2D view.
Figure 16. SDPD-DV simulations of transient planar Couette flow. Physical domain showing the BPs and FPs.
The flow starts from rest and develops a velocity profile given by Morris et al. (1997) and Papanastasiou et al. (1999) nxw xwx n Z ZZ n nV V z tV zz t L LL n (3.2) The sample-averaged field density ( , ) d t r and velocity ( , ) d t V r profiles are plotted in Figure 17(a),(b) along with the standard deviation. The field density shown in Figure 17(a) has an error a a of less than 4% in the interior and less than 5% near the wall boundary. The analytical velocity profile from Eq. (3.2) using the average density
966 kg/m is plotted in Figure 17(b). A comparison between the analytical and SDPD-DV velocity profiles across the channel is plotted in Figure 17(b). The numerical results from SDPD-DV are in very good agreement with the analytical solution for all times considered. 83 (a) Density distribution with standard devian. (b) Velocity distribution.
Figure 17. Sample-averaged fluid density ρ(Z d )and velocity V x ( Z d, t ) from SDPD-DV simulations of transient, planar, Couette flow. The domain has L X =10 -3 m, L Y =3×10 -4 m, L Z =10 -3 m, the fluid is H O(l) with ρ a =1000 kg/m , T i =300 K, η =10 -3 kg·m -1 s -1 , and V xw =1.25x10 -5 m/s. The analytical profiles are plotted for verification (Morris, Fox and Zhu, 1997). The third benchmark test involves the SDPD-DV simulation of flow over a cylinder. This test appeared in several SPH (Morris et al., 1997; Vazquez-Quesada and Ellero, 2012) and DPD (Keaveny et al., 2005; Kim and Phillips, 2004) simulations. The sample-averaged, steady velocity and pressure fields from SDPD-DV and comparison with FLUENT are shown. The SDPD-DV flow field exhibits the anticipated behavior for the low-Reynolds number flow over a cylinder (Keaveny et al., 2005; Morris et al., 1997; Vazquez-Quesada and Ellero, 2012; Kim and Phillips, 2004). For further direct quantitative comparison the SDPD-DV and FLUENT, velocity and pressure fields are superimposed and plotted ( , 7.5 10 m, ) x y z . Additional verification 84 is obtained by comparing velocity and pressure field along contours. This test further verifies the ability of SDPD-DV to simulate curved solid boundaries and the implementation of the artificial compressibility method in SDPD-DV. The incompressible fluid of density is driven by a body-force x f over a cylinder of radius R . Periodic boundary conditions are assumed in the x, y and z direction as shown in Figure 18(a). The periodic boundary conditions can be realized as a flow in an infinite array of cylinders as depicted in Figure 18(b). For verification of the SDPD-DV results we also performed simulations using FLUENT over a domain shown in Figure 18(b). The SDPD-DV simulation considers H O(l) with -3
1, 000 kgm a and
10 kg m s driven by a body force x f . The physical domain has X L , Y L , Z L and the cylinder radius R . The total fluid mass of F M in the domain is represented by 17,632 fluid particles with a constant temperature of
300 K i T . The solid walls of the cylinder are represented by 7,440 boundary particles as shown in Figure 18(a). The initial fluid particle distribution is shown in Figure 18(a). The Reynolds number Re R V based on the cylinder radius and V in the average velocity field is Re 1 . Closure is obtained by using the artificial compressibility relation Eq. (2.52). Using Eq. (2.68) integration is carried out with t . Steady state was reached and M independent samples were used to evaluate fluid properties on a grid with d l . The input parameters for the SDPD-DV simulation are listed in Table 6. (a) (b) (c) Figure 18. Computational domains used in SDPD-DV and FLUENT simulations of steady, low-Reynolds number water flow with ρ a =1,000 kg/m , T i = 300 K, η =10 -3 kg·m -1 s -1 . The cylinder has R =0.02m and the body force per unit mass is f x = g x =1.5x10 -7 m/s . (a) SDPD-DV domain with L X =0.1 m, L Y =0.015 m, L Z =0.1 m showing the BPs on the cylinder and FPs in the domain. (b) FLUENT domain L X =0.3 m, L Y =0.015 m, L Z =0.3 m includes an array of cylinders to simulate the periodic flow. The insert shows the extend of the SDPD-DV physical domain. (c) Planes P1,P2,P3 and contours C1,C2,C3,C4 used for comparison of SDPD-DV and FLUENT results. In order to achieve periodic boundary conditions in ANSYS FLUENT, we placed 9 cylinders in a lattice as shown in Figure 18(b). The cylinders were placed apart so that there is no influence and simulation was performed with 246,132 cells. The physical domain for the FLUENT simulation has X L , Y L , Z L . The gravitational acceleration per unit mass x g in ANSYS FLUENT was equal to the SDPD-DV body force x f . For direct 86 comparison we sampled field properties on planes P1 ( y =0 m), P2 ( y =7.5×10 -3 m), P3 ( y =0.015 m) and along contours C1, C2, C3 and C4 shown in Figure 18(c). (a) SDPD-DV (b) FLUENT (a) SDPD-DV (b) FLUENT Figure 19. Sample-averaged V x (r) and P (r) on y = 0 m, y = 7.5×10 -3 m, y = 1.5×10 -2 m planes from SDPD-DV simulations. V x ( x,z ) and P (r) from FLUENT simulation. Steady, low-Reynolds number water flow with ρ a =1,000 kg/m , T i = 300 K, η =10 -3 kg·m -1 s -1 . The cylinder has R =0.02m and the body force per unit mass is f x = g x =1.5x10 -7 m/s . The overall flow field characteristics are shown in Figure 19. The sample-averaged, steady ( ) x d V r and ( ) d P r pressure fields from SPDP-BV and ( ) x V r and ( ) P r from FLUENT are shown on planes P1, P2, P3 parallel to x axis. The SDPD-DV flow field exhibits the anticipated behavior for the low-Reynolds flow over a cylinder (Keaveny et al., 2005; Morris et al., 1997; Vazquez-Quesada and Ellero, 2012; Kim and Phillips, 2004). The SDPD-DV properties are qualitatively and quantitatively similar to those obtained from FLUENT. For further direct quantitative comparison the SDPD-DV and FLUENT velocity and pressure fields are superimposed and plotted in Figure 20(a)-(b) on the plane ( , 7.5 10 m, ) x y z . The non-smooth character of the SDPD-DV pressure contours are due to fluctuations introduced by the artificial compressibility model applied to this incompressible flow. It is important to note also that velocities and pressures considered are very small and susceptible to numerical perturbations. Additional verification is obtained by comparing velocity and pressure field along contours C1, C2, C3, and C4 shown in Figure 20(c)-(d). The pressure comparison along the centerline C3 in Figure 20(d) shows the ability of our solid boundary density and force evaluation to accurately predict pressure in the ram and wake side of the cylinder where fluid particles are impinging the curved boundary. The overall very good agreement demonstrates the ability of our SDPD-DV implementation to simulate curved solid boundaries. 88 (a) V X ( r ) on y =7.5×10 -3 m plane. (b) P ( r ) on y =7.5×10 -3 m plane. (c) V X ( r ) along C1 and C2 . (d) P ( r ) along C3 and C4. Figure 20. Sample-averaged V x (r) and P (r) on y = 7.5×10 -3 m plane and along C1,C2,C3,C4 from SDPD-DV simulations. V x ( x,z ) and P (r) from FLUENT simulation. Steady, low-Reynolds number water flow with ρ a =1,000 kg/m , T i = 300 K, η =10 -3 kg·m -1 s -1 . The cylinder has R =0.02m and the body force per unit mass is f x = g x =1.5x10 -7 m/s . Z (m)0.00 0.02 0.04 0.06 0.08 0.10 V x - ( m / s ) X (m)0.00 0.02 0.04 0.06 0.08 0.10 P - ( P a ) -10-50510 SDPD-DV C3FLUENT C3SDPD-DV C4FLUENT C4 The last series of tests also serve for validation and verification since the SDPD-DV results for the self-diffusion coefficient are compared with data and analytical formulas. The tests also provide the opportunity to examine the scale-dependence of our SDPD-DV implementation and compare our results with those of Litvinov et al. (2009) and Vazquez-Quesada et al. (2009). These tests are intended as a demonstration of our SDPD-DV to simulate mesoscale flows at equilibrium states.
We performed SDPD-DV simulations of H O(l) with F M in the range
22 15 ,
1, 000 kg/m a ,
10 kg m s and
300 K i T . We performed also SDPD-DV simulations of N (g) with F M in the range
18 15 with a ,
10 kg m s and
300 K i T . In order to examine fluid particle scale effects we follow Vazquez-Quesada et al. (2009) and assume that the “size” of the fluid particle is given in terms of the SDPD variables as, / i i i i D V m (3.3) For comparison, we need also length scales for the real liquid and gas molecules. For a liquid with viscosity A , temperature A T and assuming that the fluid occupies a cubic lattice where the molecules simply touch each other, the appropriate scale is given by Bird et al. (2007) A A A A A
R V N m (3.4) where, A V is the molar volume and A N is the Avogadro number, A m is the mass of the molecule A. For a gas a relevant length scale can be considered the molecular diameter g D . The input 90 parameters and some derived variables from the SDPD-DV simulations are listed in Table 7 and Table 8. Table 7. Input and derived parameters in SDPD-DV simulations of mesoscale flows of liquid water at equilibrium states in rectangular domains with periodic boundaries.
Inputs, H O(l),
300 K T Derived SDPD-DV Case , ,
X Y Z L ( m) H O i mm FP N t (×10 -16 s) h (m) i i c ( )(2 ) i VR ii D (×m s -1 ) hii D (×m s -1 ) 1 5×10 -9 -9 -9 -9
2 8×10 -9 -9 -9 -10
3 1×10 -8 -9 -10 -10
4 5×10 -8 -8 -10 -10
5 1×10 -7 -8 -11 -11
6 1×10 -6 -7 -11 -11
7 1×10 -6 -7 -11 -11
8 1×10 -6 -6 -11 -11
9 1×10 -6 -6 -12 -12
10 1×10 -6 -7 -12 -12 Table 8. Input and derived parameters in SDPD-DV simulations of mesoscale flows of N (g) at equilibrium states in rectangular domains with periodic boundaries. Inputs N (g), a kg/m ,
300 K T Derived SDPD-DV Case , ,
X Y Z L ( m) N i mm FP N t (×10 -15 s) h (m) i i c ( ) i VD
1 1×10 -6 -7 -6 -7 -5 -6 Figure 21 shows the ( ) xi yi v v phase plots for the smallest and largest scales considered. The plots show that the equilibrium states sampled have particle velocities that are not scale free. To gain further insight, we use the fluid particle velocities and field mean velocities obtained from the SDPD-DV simulations to evaluate the average translational temperature for the entire system. Following the standard description for the thermal (or random) particle velocity ( ) ( ) ( , ) i i t t t C v V r (Gombosi, 1994)
22 2 21 F N yi yxi x zi zi iit iB F B v Vv V v Vm mT ck N k (3.5) Table 7 and Table 8 show the rms thermal speeds, i c evaluated from the SDPD-DV simulations. These thermal speeds are inversely related to the fluid particle masses considered and they are equal to those obtained with the equilibrium Maxwellian formula B i i k T m for 92 the SDPD fluid. Figure 22(a) shows the evolution of the translational temperature for Case 10. This behavior is typical for all cases considered and shows that the system quickly reaches equilibrium and the imposed FP temperature of
300 K i T . Figure 22(b) and Figure 22(c) depict the average translational temperature for all cases simulated in Table 7 and Table 8. It is clear that for the range of fluid particle masses and sizes considered the average translational temperature is scale-free and matches the imposed fluid particle temperature for both the H O(l) and N (g). We also evaluate the self-diffusion coefficient for the H O(l) simulations following the MSD Eq. (2.74) and VAC Eq. (2.76) methods. Various analytical formulas for the self-diffusivity are used for validation. The self-diffusion coefficient is given by Bird et al. (2007) by
B AAA A A k TD R (3.6) Assuming by analogy, that the SDPD liquid has a lattice scale i D from Eq. (3.3), then Eq. (3.6) can be expressed as B iii i k TD D (3.7) An additional expression for the self-diffusivity in a SDPD liquid is obtained following Litvinov et al. (2009) and using the smoothing function (2.38) used in this work. It provides h i B iii i h k TD m (3.8) 93 Figure 21. Phase plot v xi -v yi from SDPD-DV simulations of H O(l) with ρ a =1,000 kg/m , T i = 300 K, η =10 -3 kg·m -1 s -1 and N (g) with ρ a =1,184 kg/m , T i = 300 K, η =10 -5 kg·m -1 s -1 . Results show the scale effects of fluid particle size on velocity. The Schmidt number is then evaluated by Sc D (3.9) 94 (a) (b) (c) Figure 22. Average translational temperature from SDPD-DV simulations of H O(l) with ρ a =1,000 kg/m , T i = 300 K, η =10 -3 kg·m -1 s -1 and N (g) with ρ a =1,184 kg/m , T i = 300 K, η =10 -5 kg·m -1 s -1 . (a) Average translational temperature as a function of time (Case 1, H O(l)). (b) Average translational temperature and standard deviation as a function of m i and D i for H O(l). (c) Average translational temperature and standard deviation as a function of m i and D i for N (g). In Figure 23(a) we plot the self-diffusion coefficient for H O(l). The experimental value D is from Holz and Sacco (2009). With
26H O m and
3H O
1, 000 kg m then Eq. (3.4) provides
10H O R and the analytical value AA D is evaluated from Eq. (3.6). The values of ii D from the MSD Eq. (2.74) and VAC Eq. (2.76) are also plotted in Figure 23(a). For comparison we also evaluate ii D from Eq. (3.7) and hii D from Eq. (3.8) using SDPD-DV results for all the variables entering these expressions. Figure 23(a) shows that the self-diffusion coefficients obtained from MSD, VAC and the SDPD-liquid formulas are in good agreement been within a factor of about 3 of each other. The experimental 95 and analytical value of the self-diffusion coefficient is approached by SDPD-DV, as the mass of the SDPD particles approaches asymptotically the relevant physical scale. (a) (b) Figure 23. Self-diffusion coefficient and Schmidt number from SDPD-DV simulations of H O(l) , ρ a =1,000 kg/m , η =10 -3 kg·m -1 s -1 , T =300K for various sizes of the SDPD-DV fluid particles. Experimental and analytical estimates, as well as analytical estimates using SDPD-DV parameters are shown for validation and verification. It should be noted, that at the real molecule level the applicability of the SDPD model becomes tenuous. These results also show that for the fluid particle sizes considered in this work, the self-diffusion coefficient from SDPD-DV is not scale-free, a result that corroborates 96 earlier results of Litvinov et al. (2009) and Vazquez-Quesada et al. (2009). The Sc numbers from Eq. (3.9) using the various diffusion-coefficients are plotted in Figure 23(b). Similar to previous results, we show that the Sc number scales with the mass (or size) of the SDPD particles, identified in our work with either the SDPD-fluid size i D through Eq. (3.7) or the smoothing length h through Eq. (3.8). The results show also that with SDPD-DV we achieve Sc values close to the realistic ones. 97
4. VERIFICATION, VALIDATION AND ERROR ANALYSIS OF SDPD-DV
In this chapter, we perform validation and verification of SDPD-DV using the full-set of SDPD Eq. (1.62)-(1.64) as implemented in the SDPD-DV method discussed in Chapter 2. The validation and verification includes an extensive set of SDPD-DV simulations of gaseous nitrogen in mesoscopic periodic domains in equilibrium. The self-diffusion coefficient for N (g) and shear viscosity at equilibrium states are obtained through the mean-square displacement for the range of fluid particle masses (or sizes) considered. Additional verification involves SDPD-DV simulations of steady Couette N (g) flow between parallel plates. The top plate is moving at V xw =30m/s and separated by 10 -4 m from the bottom stationary plate. The top plate is assigned a constant temperature T =330K and bottom plate T =300K. The SDPD-DV field velocity and temperature profiles are compared with those obtained by FLUENT. We consider first N (g) systems in thermal equilibrium. The full algorithm of SDPD equations is considered with entropy provided by Eq. (1.64) and closed by Eq. (2.44)-(2.49). The simulations of N (g) were performed in rectangular domains with X Y Z
L L L in the range 0.25×10 -6 ~10×10 -6 m, with mass F M in the range
20 15 . The simulations examine the effects of fluid particle mass N / i m m and the smoothing length h . 98 Table 9: Input parameters in full-set SDPD-DV simulations of mesoscale flows of N (g) at equilibrium states in rectangular domains with periodic boundaries. Case Inputs N (g), a kg/m ,
300 K T , P =101942.4 XYZ L (×10 -6 m) N i mm FP N i V (×10 -8 m) h (×10 -8 m) t (×10 -15 s) (×10 -5 kg m -1 s -1 ) k (W m -1 K -1 ) 1 0.25 118 3375 1.67 3.7 1 1.79 0.026 2 0.25 118 3375 1.67 4.7 1 1.79 0.026 3 0.25 118 3375 1.67 5.9 1 1.79 0.026 4 0.25 118 3375 1.67 7.5 1 1.79 0.026 5 0.55 Table 10: Derived parameters in SDPD-DV simulations of mesoscale flows of N (g) at equilibrium states in rectangular domains with periodic boundaries. Case Derived SDPD-DV ( ) i V (×10 -8 m) i (kg/m ) T (K) P (Pa) i C D (×10 -9 m /s) (×10 -5 kg m -1 s -1 ) 1 1.66 1.193 299.2 102410 47.64 3.1 4 2 1.67 1.188 299.25 102030 47.64 4.1 2 3 1.67 1.185 299.25 101800 47.75 5.8 3 4 1.67 1.184 299.3 101750 47.64 8.5 2 5 3.667 1.184 299.41 101752 14.59 1.4 4.5 6 3.662 1.188 299.5 102155 14.63 2.4 1.5 7 3.665 1.185 299.43 101867 14.59 2.7 2 8 3.666 1.184 299.40 101775 14.56 4.2 1.5 9 6.678 1.178 299.41 101227 5.96 0.3 1 10 6.678 1.178 299.43 101228 5.96 0.7 1.5 11 6.668 1.183 299.42 101684 5.96 1.4 1 12 6.667 1.184 299.47 101729 5.96 2.9 1.5 13 16.7 1.177 299.41 101160 1.53 0.21 3 14 16.7 1.176 299.41 101052 1.52 0.3 3 15 16.7 1.176 299.42 101065 1.61 0.4 2 16 16.7 1.172 299.40 100752 1.52 0.65 2.5 17 66.9 1.173 299.4 100803 0.195 0.0365 1 18 66.7 1.171 299.4 100664 0.190 0.0259 1 19 66.7 1.181 299.4 101506 0.187 0.044 1 20 66.9 1.172 299.4 100670 0.187 0.1 3
100 For each N / i m m we considered four h ’s so that the resulting number of fluid particles within the support domain, i L , is 50, 100, 200, 360 upon initialization. Periodic boundary conditions are imposed on each side of the rectangular domain. Upon initialization we set a , ( 0) 300K i T t , ( 0) 1.79 10 kg/m s i t and ( 0) 0.026 W/m K i t . During the computation the transport coefficients appearing in the SDPD equations (1.62)-(1.64) are evaluated as functions of particle temperature i T t according to power law (White, 1974) ( )( ) nii
T tt T (4.1) nii
T tk t k T (4.2) where n is the power law coefficient of the order of 0.7, T the reference temperature, the heat conductivity of N (g) at T , and is the viscosity of N (g) at T . The heat capacity (J/K) is evaluated as Eq. (2.48). i i B C N k (2.48) In order to examine fluid particle scale effects we follow Vazquez-Quesada et al (2009) discussion and assume that the “size” of the fluid particle is given in terms of the SDPD variables as, / i i i i D V m (4.3) Input and some derived parameters from the SDPD-DV simulations are shown in Table 9 and Table 10. 101
Table 10 and Table 11 show that the average density , temperature T and pressure P exhibit minimal perturbations from the input value for the entire range of parameters considered. These rms thermal speeds, i c evaluated from the SDPD-DV simulations, are inversely related to the fluid particle masses considered and they are equal to those obtained with the equilibrium Maxwellian formula B i i k T m for the SDPD fluid.
Figure 24: Self-diffusion coefficient from SDPD-DV simulations of N (g) for values of h and m i / m N (g) used in Cases 1-20. Analytical estimates are from Eq.(4.4) and data by Holz and Sacco (2000). m i /m N (g) =118m i /m N (g) =1255m i /m N (g) =7544m i /m N (g) =1.2 10 m i /m N (g) =7.5 10 Eq. (4.4)Data (Bird et al., 2007)Smoothing Length h (m)10 -7 -6 S e l f- D i ff u s i v it y D ( m / s ) -11 -10 -9 -8 -7 -6 -5 Case 1 2 3 4Case 5 6 7 8 14Case 13 15 16Case 17 1819 20Case 9 10 11 12 m i /m N (g) =118m i /m N (g) =1255m i /m N (g) =7544m i /m N (g) =1.2 10 m i /m N (g) =7.5 10 Eq. (4.4)Data (Holz and Sacco, 2000))
102 We also evaluate the transport coefficient for the N (g) simulations following the MSD Eq. (2.74) for self-diffusion coefficient and Eq.(2.77) for viscosity. Validation and verification is obtained by comparison with experimental (Holz and Sacco, 2000) and analytical values. For a gas with molecular diameter A D , molecular mass A m and mass density A the self-diffusivity is given by Bird et al. (2007) as A B AAA A A m k TD D . (4.4)
The self-diffusion coefficient of SDPD-DV simulations from Eq. (2.74) for N (g) are plotted in Figure 24 for various smoothing lengths h and particle masses i m . For comparison, the values of AA D from Eq. (4.4) are evaluated and plotted in Figure 24 using N (g) molecular value. As it shown in Figure 24, the self-diffusivity of SDPD-DV is not scale free. The simulations show that D decreases with increasing mass ratio. This is because that the smaller the fluid particle is, the larger is its stochastic agitation. Figure 25 : Shear viscosity from SDPD-DV simulations of N (g) for value of m i / m N (g) used in Cases 1-20. Initial input viscosity η ( t =0) is plotted for verification. For a given mass ratio, increasing the h , increases D , since that the stochastic agitation is stronger in the bigger support domain. The values for calculated from the SDPD-DV m i /m N (g) S h ea r V i s c o s it y kg m - s - ) t
103 simulations using Eq. (2.77) are plotted in Figure 25. We plot for comparison the initial value ( 0) t which is also . It can be seen that shear viscosity is scale free and is not affected by the choice of particle mass or the smoothing length. Figure 26: Effects of time step on equilibrium temperature ̅ , translational temperature ̅ and pressure ̅ from SDPD-DV simulations of N (g). Results are for Case 9, 10, and 11. Figure 27: Effects of time step on equilibrium temperature ̅ , translational temperature ̅ and pressure ̅ from SDPD-DV simulations of N (g). Results are for Case 17, 18, and 19. t (s)10 -13 -12 -11 -10 T ( K ) m =45N m =90N m =180 Case 9Case 10Case 11 t (s)10 -13 -12 -11 -10 T r ( K ) m =45N m =90N m =180 t (s)10 -13 -12 -11 -10 P ( P a ) m =45N m =90N m =180 t T K P a P T K t s t s t s t (s)10 -13 -12 -11 -10 -9 -8 T ( K ) m =45N m =90N m =180 (s) t ( K ) T Case 17Case 18Case 19 t (s)10 -13 -12 -11 -10 -9 -8 T t r ( K ) m =45N m =90N m =180 (s) t ( K ) t T (s) t t (s)10 -13 -12 -11 -10 -9 -8 P X ( P a ) ( p ) P Table 11: Fluctuations in temperature, density, pressure and velocity, from SDPD-DV and analytical expressions. case Analytical Variance Variance in SDPD-DV T ( ) P V T ( ) P V
1 0.127 3.04×10 -6 -7 -6 -7 -6 -8 -6 -8 -2 -7 -2 -2 -8 -2
6 1.2×10 -2 -7 -2 -2 -8 -2
7 1.2×10 -2 -7 -2 -2 -9 -2
8 1.2×10 -2 -7 -2 -2 -9 -2
9 2×10 -3 -8
341 1.05×10 -2 -3 -8
392 1.05×10 -2
10 2×10 -3 -8
342 1.05×10 -2 -3 -8
329 1.04×10 -2
11 2×10 -3 -8
345 1.05×10 -2 -3 -8
359 1.05×10 -2
12 2×10 -3 -8
345 1.05×10 -2 -3 -9
345 1.06×10 -2
13 1.3×10 -4 -9 -4 -4 -8
159 6.92×10 -4
14 1.3×10 -4 -9 -4 -4 -9 -4
15 1.3×10 -4 -9 -4 -4 -8
524 6.73×10 -4
16 1.3×10 -4 -9 -4 -4 -8
500 6.69×10 -4
17 2.0×10 -6 -11 -5 -6 -7 -5
18 2.0×10 -6 -11 -5 -6 -7 -5
19 2.0×10 -6 -11 -5 -6 -8 -5
20 2.0×10 -6 -11 -5 -6 -8
416 1.05×10 -5
105 We examine also the effects of time step in the SDPD-DV simulations on the equilibrium characteristics. In Figure 26 and Figure 27, we plot the SDPD-DV temperature T averaged over the entire domain, the average translational temperature t T from (3.5), and the average SDPD-DV pressure P for Case 9, 10, 11 and Case 17, 18, and 19. These results show that time step can lead to a significant error depending on the FP mass and smoothing length. The results show also that for larger FPs the stochastic agitation is smaller. We compute also the variances in temperature, density, pressure and velocity for each case 1-20 according to Eq. (1.25), (1.22), (1.23) and (1.27). The results in Table 11 show that SDPD-DV are in very good agreement with the analytical values. In this section, we continue the verification of SDPD-DV by performing simulations of steady planar non-isothermal Couette flow.
The test involves an incompressible flow with density a across two infinite parallel walls as depicted in Figure 28(a) with imposed constant wall temperatures and the top wall is moving with a constant velocity xw V . The SDPD-DV simulation considers N (g) with density -3 a , initial viscosity ( 0) 1.79 10 kg m s t and initial heat conductivity -1 -1 ( 0) 0.026 W m K t . Closure for pressure and density is obtained by using the Eq.(2.47) and (2.46). Heat capacity V C is given by Eq.(2.48) with initial value ( 0) 2.42 10 J/K i C t . 106 The physical domain has X Z L and Y L as shown in Figure 28. Periodic boundary are imposed along the -axis x and -axis y . The fluid is represented by 10,890 FPs and total mass F M . The upper wall is located at Z and is assigned a velocity
30 ms xw V and temperature T . The temperature of the lower wall located at Z is set to T . The solid walls are represented by 10,560 BPs. Input conditions used in the SDPD-DV simulation are shown in Table 12. Table 12: Input parameters used in SDPD-DV non-equilibrium simulations of Couette flow.
Input Parameters Coutte (m) X L -4 (m) Y L -5 (m) Z L -4 -2 (ms ) x f N/A kg F M -13 Upper (K) T Lower (K) T -3 kg m a
0 kg m s t -5 -1 -1 ( 0) kg m s t -1 ( 0) J K i C t -14 -1 -1 t -1 ms xw V F N B N C N (m) h -6 (s) t -8
107 For verification we also performed a simulation using FLUENT in a domain shown in Figure 28 and input parameters shown in Table 12. The wall boundary conditions are assigned to both upper wall and lower wall with periodic boundary conditions to other boundaries. (a) SDPD-DV domain with L x =10 -4 m, L y =3×10 -5 m, L z =10 -4 m showing the BPs on the top and bottom and FPs in the domain. (b) FLUENT domain with L x =10 -4 m, L y =3×10 -5 m, L z =10 -4 m showing imposed wall boundary conditions and periodic boundary conditions. Figure 28: Computational domains used in SDPD-DV and FLUENT simulations of steady N (g) Couette flow with ρ a =1.184 kg/m , η =1.79×10 -5 kg·m -1 s -1 , κ =0.026 W·m -1 K -1 . The upper wall T = 330 K, lower wall T = 300 K and V xw =30m/s. The overall flow field characteristics are shown in Figure 29. We plot the SDPD-DV sample-averaged steady temperature T r and velocity x V r fields on planes X , X , and X . xw V T T W a ll B ound a r y W a ll B ound a r y P e r i od i c B ound a r y xw V T T (a) SDPD-DV (b)
FLUENT (c)
SDPD-DV (d)
FLUENT
Figure 29: Sample-averaged V x (r) and T (r) on y = 0 m, y = 1.5×10 -5 m, y = 3×10 -5 m planes from SDPD-DV and FLUENT simulations of steady state Couette flow. The domain has L X =10 -4 m, L Y =3×10 -5 m, L Z =10 -4 m, upper wall T = 330 K, lower wall T = 300 K, and V xw =30 m/s. The fluid is N (g) with ρ a =1.184 kg/m , η =1.79×10 -5 kg·m -1 s -1 and , κ =0.026 W·m -1 K -1 .
109 We also plot for comparison the steady-state T r and x V r from the FLUENT simulation. The velocity and temperature profiles are quantitatively similar to those obtained from FLUENT. For direct quantitative comparison, the sample-averaged field temperature ( ) d T r and velocity ( ) x d V r profiles are plotted in Figure 30. The SDPD-DV properties are in excellent agreement with FLUENT solutions. Figure 30: Sample-averaged fluid temperature T (r d ) and velocity V x (r d ) of steady-state Couette flow from SDPD-DV and FLUENT simulations. The domain has L X =10 -4 m, L Y =3×10 -5 m, L Z =10 -4 m, upper wall T = 330 K, lower wall T = 300 K and V xw =30 m/s. The fluid is N (g) with ρ a =1.184 kg/m , η =1.79×10 -5 kg·m -1 s -1 and , κ =0.026 W·m -1 K -1 . Z (m)0 2 4 6 8 10 T ( K ) V x ( m / s ) T (SDPD-DV) V x (SDPD-DV) T (FLUENT) V x (FLUENT)
5. SUMMARY, CONCLUSIONS AND RECOMMENDATIONS
This work presented the mathematical formulation and computational implementation of a smoothed dissipative particle dynamics method for wall-bounded domains (SDPD-DV) under development at the Computational Gas&Plasma Dynamics Lab (CGPL) at WPI. The SDPD-DV allows simulation of liquid and gaseous flows at mesoscopic and hydrodynamic scales. The SDPD-DV method and implementation utilizes fluid particles, boundary particles and dynamical allocated virtual particles. The boundary particles are placed on the physical boundary of the domain while fluid particles are used to discretize the fluid region of the domain. The algorithm for nearest neighbor particle search (NNPS) is based on a combination of the linked-cell and Verlet-list approaches. The NNPS utilizes large rectangular cells that organize the physical domain and increase computational efficiency. The dynamic virtual particle allocation (DVPA) method introduced in this work provides a formal way to minimize error in the density and force evaluation due to the truncated part of the support domain near a solid boundary. Models for the force components due to virtual particles are introduced near solid boundaries. A periodic boundary particle allocation method is used at periodic inlets and outlets. Integration of particle position and momentum equation is investigated with the implementation of the velocity-Verlet algorithm. The integration is supplemented by a bounce-forward algorithm to further prevent particle penetration. SDPD-DV outputs are sampled to obtain unsteady and steady particle properties. The self-diffusion coefficient is obtained by implementing the generalized Einstein (mean square displacement) and Green-Kubo (velocity 111 auto-correlation) methods. Field properties are obtained by sampling particle properties on a post-processing Eulerian grid using the smoothing function evaluation.
This work featured several modifications and revisions of existing algorithms in the SDPD-DV code as well as implementation of new algorithms in order to achieve a fully functional SDPD-DV code for wall-bounded mesoscopic computations. Specifically this work: 1.
Revised the periodic boundary conditions algorithm and the implementation of the periodic boundary cells searching list. 2.
Modified the evaluation of smoothing function in order to include the contribution of the self-density for each fluid particle. 3.
Modified the dynamic virtual particle allocation algorithm for the modeling of solid boundary in order to minimize the truncation error of density. 4.
Developed and implemented the boundary normal vector algorithm used in the reflection of the dynamically allocated virtual particles. 5.
Implemented the algorithm for the contribution to the boundary force from the virtual particles. 6.
Implemented a constant-wall temperature boundary condition in the dynamic virtual particle allocation model. 7.
Implemented a bounce-forward algorithm for a solid boundary with arbitrary shape and orientation. 8.
Revised the portions of the Velocity Verlet integration method for the position and momentum equations. 9.
Developed and implemented a Runge–Kutta integration algorithm for the entropy equation. 112 10.
Implemented the artificial incompressibility method for liquid flows. 11.
Implemented a temperature power law for the shear viscosity, bulk viscosity and heat conductivity that appear in the momentum and energy SDPD equations. 12.
Developed and implemented algorithms for the evaluation of transport properties such as diffusion coefficient, shear viscosity and heat conductivity based on mean square displacement (MSD) and velocity autocorrelation function (VACF). 13.
Developed analytical formulas of the self-diffusion coefficient based on the SDPD-fluid following Litvinov et al. (2009). Several benchmark tests were performed for validation and verification of the SDPD-DV method. The first verification of SDPD-DV involves the simulation of transient, planar Poiseuille body-force driven liquid water flow. The parallel plates are separated by a
10 m and H O(l) at
300 K T is driven by a body force per unit mass
10 ms x f . Numerical results show that density error is less than 4% in the interior and less than 5% near the walls. The velocity profiles are in very close agreement with theoretical ones. The second benchmark test involves transient, planar Couette water liquid flow at
300 K T . The parallel plates are separated by
10 m with the top plate moving at xw V . The numerical velocity profiles are in very good agreement with analytical solutions. These two benchmark tests verify major algorithmic parts of our SDPD-DV implementation and demonstrate the ability of SDPD-DV to enforce no-slip boundary conditions, avoid particle penetration of walls, and reduce the density error near boundaries. Further verification SDPD-DV is accomplished with a more challenging 3D benchmark test involving the body-force driven flow of H O(l) at
300 K T over a cylinder of radius R . The SDPD-DV simulations are performed in a domain with X L , 113 Y L , Z L and a body-force per unit mass of x f . FLUENT simulations are also performed in a lattice of cylinders with the flow driven by a x g . The SDPD-DV results are compared directly with FLUENT results and are in excellent agreement. This test further verifies the ability of SDPD-DV to simulate curved solid boundaries and the implementation of the artificial compressibility method in SDPD-DV. Additional SDPD-DV verification and validation features an extensive set of simulations used to obtain the equilibrium state of H O (l) and N (g) at
300 K T , and the self-diffusion coefficient of H O (l). In order to examine the scale-effects in SDPD-DV, the mass of the fluid particles for the H2O(l) SDPD-DV simulations is varied between 1.24 and 3.3×10 real molecular masses and their corresponding size is between 1.08 and 323 physical length scales. For the N (g) SDPD-DV simulations the mass of the fluid particles is varied between 6.37×10 and 6.37×10 real molecular masses and their corresponding size is between 2.2×10 and 2.2×10 physical length scales. The particle speeds are shown to depend on mass (or size) of the fluid particle. The average translational temperature from the SDPD-DV results is close to 300K and is scale-free for the range of particle masses (or sizes) considered. The self-diffusion coefficient for H O(l) is obtained using the MSD and VAC methods. Additional estimates are obtained from analytical expressions for the SDPD-fluid and values derived from the SDPD-DV simulations. Comparisons are also given with experimental data. The results show that MSD and VAC results are in very good agreement with the analytical SDPD estimates. The results also show that the self-diffusion coefficient from SDPD-DV is not scale-free and reaches the experimental value of H O(l) when the fluid particle mass asymptotically approaches the actual molecular size. The Schmidt numbers obtained from SDPD-DV are within the range expected for liquid water. 114 The full-set of SDPD equations implemented in SDPD-DV model is verified and validated with simulations in bounded and periodic domains that cover the hydrodynamic and mesoscopic regimes. Validation and verification is achieved with an extensive set of SDPD-DV simulations of gaseous nitrogen in mesoscopic periodic domains in equilibrium. The simulations of N (g) were performed in rectangular domains with in the range 0.25×10 -6 ~10×10 -6 m, with mass in the range 1.85×10 -20 ~1.184×10 -15 kg. The self-diffusion coefficient for N (g) at equilibrium states is obtained through the mean-square displacement for the range of fluid particle masses (or sizes) and smoothing length considered. The simulations show that both of the fluid particle mass and smoothing length affect the self–diffusion coefficient which is not scale free. This is because the smaller the fluid particle is, the larger is its stochastic agitation. For a given mass ratio, increasing the smoothing length, increases self-diffusivity, since the stochastic agitation is stronger in the bigger support domain. The shear viscosity for SDPD-DV simulation of N (g) is shown to be scale free and is not affected by the choice of particle mass or the smoothing length. The impact of fluid particle mass, smoothing length and timesteps are also presented. The results show also that for larger fluid particles the stochastic agitation is smaller. Additional verification involves SDPD-DV simulations of steady Couette N (g) flow. The top plate is moving at V xw =30m/s and separated by 10 -4 m from the bottom stationary plate. The top plate has a constant and the bottom plate . The SDPD-DV field velocity and temperature profiles are in excellent agreement with those obtained by FLUENT. This benchmark test further verifies the ability of our SDPD-DV to simulate non-isothermal flows and the implementation of the thermodynamics properties of virtual particles in SDPD-DV. 115 The limitations of the SDPD-DV method and code developed in the course of this research indicate several areas for future work. 1.
The dynamic virtual particle allocation method discussed in this work can lead to an imperfect representation of uneven boundary. In view of the need for wide applications with arbitrary geometries, there is need for modeling a solid boundary with dihedral or trihedral angle. With the implementation of the dynamic virtual particle allocation method, the reflection algorithm on flat wall is very mature. For curved boundary, higher resolution is required to avoid large errors from reflection distortion. However, for dihedral or trihedral angle, a new algorithm is needed to compensate for this distortion. 2.
Implementation of variable smoothing length would improve the accuracy of estimation with enough particles in smoothing length. Several algorithms with VSL can be found in SPH applications. Nevertheless, on a curved boundary the varying smoothing length would be larger than the radius that lead to incorrect representation. 3.
Development and implementation of a free surface boundary model with the DVPA method. A moveable boundary particle could be employed to form the free surface, as well as the virtual particle making up the truncation area. In turns, the model of boundary particle and fluid particle interaction will need to be added to simulate free surface. 4.
In case of more complicated boundary conditions, it is considerable to combine the dynamic virtual particle method and frozen ghost/virtual particle method. A wide application of SDPD-DV depends on the ability of handling complex flows. Considering the current SDPD-DV solver, future work could include: 116 1.
Development and implementation of a multiphase flow model. Particle interaction model need to be added to basic SDPD model for multiphase flow. 2.
Development and implementation of multispecies flow model. In the current solver, multispecies flow model capability can be developed by adding the multispecies particle interaction model. The force interpolation function of the fluid particle with another species fluid particle in support domain, the average coefficients of each force term are needed. 3.
Implementation of other forms of the continuity equation instead of summation form. 4.
Development and implementation of inflow and outflow boundary condition. In order to address pressure driven flows, the inlet and outlet boundary condition can be added by employing a buffer region or dynamic virtual particles. With a buffer region, the particle in the outflow can be re-injected into buffer region and will be assigned new properties. 5.
Implementation of heat flux boundary model. The heat flux boundary can be modeled by assigning different temperature to virtual particles. For example, with the same temperature as fluid particle, the virtual particle will ensure the zero heat flux on boundary. Assigning a constant temperature to virtual particle can model the heat bath. In addition to the algorithmic developments future work should include: 1.
Further study of the scale dependence of transport coefficients including the thermal conductivity. 2.
Implementation of a local diffusivity algorithm. This can be done either in sampling process or during the iterations. 117
REFERENCES
Allen, M.P. and Tildesley, D.J. (1987) Computer Simulation of Liquids, Oxford: Clarendon Press. Altenhoff, A.M., Walther, J.H. and Koumoutsakos, P. (2007) 'A Stochastic Boundary Forcing for Dissipative Particle Dynamics', Journal of Computational Physics, vol. 225, pp. 1125-1136. Amicarelli, A., Marongiu, J.-C., Leboeuf, F., Leduc, J. and Caro, J. (2011) 'SPH truncation error in estimating a 3D function', Computers & Fluids, vol. 44, pp. 279-296. Anurag Kumar, Y. (2009 ) 'From dissipative particle dynamics scales to physical scales: a coarse-graining study for water flow in microchannel', Microfluidics and Nanofluidics. Ata, R. and Soulaimani, A. (2005) 'A stabilized SPH method for inviscid shallow water flows', Int. J. Num. Method, pp. 139-159. Avalos, J.B. and Mackie, A.D. (1997) 'Dissipative Particle Dynamics with Energy Conservation', Eurphys. Lett., vol. 40, no. 2, pp. 141-146. Batchelor, G.K. (1967) An Introduction to Fluid Dynamics, Cambridge, UK: Cambridge Univ. Press. Benz, W., Hills, J.G. and Thielemann, F.K. (1989) 'Three-dimensional Hydrodynamical Simulations of Stellar Collisions', The Astrophysical Journal, vol. 342, pp. 986-998. Bian, X., Litvinov, S., Qian, R., Ellero, M. and Adams, N.A. (2012) 'Multiscale modeling of particle in suspension with SDPD', Physics of Fluids. Bird, R.B., Stewart, W.E. and Lightfoot, E.N. (2007) Transport Phenomena, New York: John Wiley & Sons. 118 Boek, E., Coveney, P., Lekkererker, H. and van der Schoot, P. (1997) 'Simulating the rheology of dense colloidal suspensions using Dissipative Particle Dynamics', Phys. Rev. E, vol. 55, pp. 3124-3133. Chatterjee, A. (2007) 'Modification to Lees–Edwards periodic boundary condition for dissipative particle dynamics simulation with high dissipation rates', Molecular Simulation, vol. 33, no. 15, pp. 1233-1236. Chaudhri, A. and Lukes, J.R. (2009) 'Multicomponent Energy Conserving Dissipative Particle Dynamics: A General Framework for Mesoscopic Heat Transfer Applications', Journal of Heat Transfer, vol. 131. de Groot, S.R. and Mazur, P. (1962) Non-equilibrium Thermodynamics , New York: Dover Publication Inc. Duong-Hong, D., Phan-Thien, N. and Fan, X. (2004) 'An implementation of no-slip boundary conditions in DPD', Comput Mech, vol. 35, pp. 24-29. Dzwinel, W., Yuen, D. and Boryczko, K. (2002) 'Mesoscopic dynamics of colloids simulated with dissipative particle dynamics and fluid particle model', J Mol Model. Ellero, M. and Adams, N.A. (2011) 'SPH simulations of flow around a periodic array of cylinders confined in a channel', Int. J. Numer. Meth. Engng, pp. 86: 1027-1040. Ellero, M., Kroger, M. and Hess, S. (2006) 'Multiscale Modeling of Viscoelastic Materials Containing Rigid Nonrotating Inclusions', Multiscale Modeling Simulation, vol. 5, no. 3, pp. 759-785. Ellero, M., Serrano, M. and Espanol, P. (2007) 'Incompressible smoothed particle hydrodynamics'. 119 Espanol, P. (1997) 'Dissipative Particle Dynamics with Energy Conservation ', Europhys. Lett., vol. 40, no. 6, pp. 631-636. Espanol, P. and Revenga, M. (2003) 'Smoothed Dissipative Particle Dynamics', Physical Review. Espanol, P. and Serrano, M. (1999) 'Thermodynamically Admissible Form for Discrete Hydrodynamics', Physical Review Letter, vol. 83, no. 22. Espanol, P. and Warren, P. (1995) 'Statistical Mechanics of Dissipative Particle Dynamics', vol. 30, no. 4. Fang, J., Parriaux, A., Rentschler, M. and Ancey, C. (2009) 'Improved SPH methods for simulating free surface flows of viscous fluids', Applied Numverical Mathematics, vol. 59, pp. 251-271. Fatehi, R. and Manzari, M.T. (2011) 'Error Estimation in smoothed particle hydrodynamics and a new scheme for second derivatives', Computers and Mathematics with Applications, vol. 61, pp. 482-498. Fedosov, D.A., Pivkin, I.V. and Karniadakis, G.E. (2008) 'Velocity limit in DPD simulations of wall-bounded flows', vol. 227. Filipovic, N., Ivanovic, M. and Kojic, M. (2009) 'A comparative numerical study between dissipative particle dynamics and smoothed particle hydrodynamics when applied to simple unsteady flows in microfluidics', Microfluid Nanofluid, pp. 7: 227-235. Flebbe, O., Munzel, S., Herold, H., Riffert, H. and Ruber, H. (1994) 'Smoothed Particle Hydrodynamics: Physical Viscosity and the Simulation of Accretion Disks', The Astrophysical Journal, vol. 431, pp. 754-760. FLUENT, A. (n.d) '13.2.5 Natural Convection and Buoyancy-Driven Flows'. 120 Frenkel, D. and Smit, B. (2002) 'Molecular Dynamics Simulation', in Frenkel, D. and Smit, B. Understanding Molecular Simulation from Algorithms to Applications, New York: Academic Press. Gatsonis, N.A., Potami, R. and Yang, J. (2013) 'Boundary Conditions in 3D Smooth Dissipative Particle Dynamics Simulations of Flows', under review in Journal of Computational Physics. Gombosi, T.I. (1994) Gaskinetic Theory, New York: Cambridge University Press. Haber, S., Filipovic, N., Kojic, M. and Tsuda, A. (2006) 'Dissipative particle dynamics simulation of flow generated by two rotating concentric cylinders: Boundary conditions', Physical Review. Hadjiconstantinou, N.G., Garcia, A.L., Bazant, M.Z. and He, G. (2003) 'Statistical error in particle simulations of hydrodynamic phenomena', Journal of Computational Physics, vol. 187. Haile, J.M. (1992). Molecular Dynamics Simulation Elementary Methods, New York: A Wiley-Interscience Publication. He, P. and Qiao, R. (2008) 'Self-consistent fluctuating hydrodynamics simulations of thermal transport in nanoparticle suspensions', Journal of Applied Physics. Hieber, S.E. and Koumoutsakos, P. (2008) 'An immersed boundary method for smoothed particle hydrodynamics of self-propelled swimmers', Journal of Computational Physics, vol. 227, pp. 8636-8654. Hockney, R.W. and Eastwook, J.W. (1981) Computing Simulations Using Particles, New York: McGraw-Hill. 121 Holz, S.R.H.M. and Sacco, A. (2000) 'Temperature-dependent Self-diffusion Coefficients of Water and Six Selected Molecular Liquids for Calibration in Accurate 1H NMR PFG Measurements', Chem. Chem. Phys., vol. 2. Hoogerbrugge and Koelman (1992) 'Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics', vol. 19(3). Hu, X. and Adams, N. (2006) 'A multi-phase SPH method for macroscopic and mesoscopic flows', Journal of Computational Physics. in't Veld, P.J., Plimpton, S.J. and Grest, G.S. (2008) 'Accurate and efficient methods for modeling colloidal mixtures in an explicit solvent using molecular dynamics', Comput. Phys. Comm. Keaveny, E., Pivkin, I., Maxey, M. and Karniadakis, G. (2005) 'A comparative study between dissipative particle dynamics and molecular dynamics for simple and complex geometry flows', The Journal Of Chemical Physics. Kim, J.M. and Phillips, R.J. (2004) 'Dissipative Particle Dynamics Simulation of Flow around Spheres and Cylinders at finite Reynolds Numbers', Chemical Eng. Science, vol. 59, pp. 4155-4168. Koumoutsakos, P. (2005) 'Multiscale Flow Simulation Using Particles', Annu. Rev. Fluid Mech., vol. 37. Kumar, A., Asako, Y., Abu-Nada, E., Krafczyk, M. and Faghri, M. (2009) 'From Dissipative Particle Dynamics Scales to Physical Scales: a Coarse-graining Study for Water Flow in Microchannel', Microfluids and Nanofluidics. Landau, L.D. and Lifshitz, E.M. (1980) Statistical Physics Part 1, 3rd edition, Oxford: Elsevier. 122 Landau, L.D. and Lifshitz, E.M. (1980) Statistical Physics, Part 2, Oxford: Pergamon Press. Laurendeau, N.M. (2005) Statistical Thermodynamics - Fundamentals and Applications, Cambridge University Press. Libersky, L. (1993) 'High Strain Lagrangian Hydrodynamics', J. Comp. Physics, pp. 67-75. Li, S. and Liu, W. (2007) Meshfree Particle Methods, New York: Springer Berlin Heidelberg. Litvinov, S., Ellero, M., Hu, X. and Adams, N. (2008) 'Smoothed dissipative particle dynamics model for polymer molecules in suspension', Physical Review. Litvinov, S., Ellero, M., Hu, X. and Adams, N.A. (2009) 'Self-diffusion Coefficient in smoothed dissipative particle dynamics', vol. 130. Litvinov, S., Ellero, M., Hu, X.Y. and Adams, N.A. (2010) 'A splitting scheme for highly dissipative smoothed particle dynamics'. Liu, G.R. and Liu, M.B. (2003) Smoothed Particle Hydrodynamics a meshfree particle method, Singapore: World Scientific Publishing Co. Pte. Ltd. Liu, M., Meakin, P. and Huang, H. (2007) 'Dissipative particle dynamics simulation of multiphase fluid flow in microchannels and microchannels networks', Physics of FLuids. Liu, M.B., Xie, W.P. and Liu., G.R. (2005) 'Modeling incompressible flows using a finite particle method', vol. 29. Lucy, L.B. (1977) 'A numerical approach to the testing of the fission hypothesis', Astron.. Mackie, A.D., Avalos, J.B. and Navas, V. (1999) 'Dissipative Particle Dynamics with Energy Conservation: Modeling of heat flow', PCCP. 123 Monaghan, J.J. (1994) 'Simulating Free Surface Flows with SPH', Journal of Compuational Physics. Monaghan, J.J. (2005) 'Smoothed particle hydrodynamics', Journal of Computational Physics, p. 399. Monaghan, J.J. (2009) 'A Turbulence Model for Smoothed Particle Hydrodynamics'. Morris, J.P., Fox, P.J. and Zhu, Y. (1997) 'Modeling Low Reynolds Number Incompressible Flows Using SPH', J. Comp. Physics, pp. 214-226. Nguyen, N.-T. and Wereley, T. (2002) Fundamentals and Applications of Microfluidics, Norwood: Artech House Inc. Nikunen, P., Karttunen, M. and Vattulainen, M. (2003) 'How would you integrate the equations of motion in Dissipative Particle Dynamics simulations', Computer Physics Communications, pp. 407-423. Novik, K. and Coveney, P. (1997) 'Using Dissipative Particle Dynamics to Model Binary Immiscible Fluids', International Journal of Modern Physics. Ortiz De Zarate, J.M. and Sengers, J.V. (2006) Hydrodynamics Fluctuations in Fluids and Fluid Mixtures, 1st edition, Oxford: Elsevier. Ottinger, H.C. (2005) Beyond Equilibrium Thermodynamics, New York: John Wiley & Sons Inc. Pan, W., Fedsov, D.A., Caswell, B. and Karniadakis, G.E. (2008) 'Hydrodynamics Interactions for Single Dissipative Particle Dynamics Particles and Their Clusters and Filaments', Phys. Rev. E. Papanastasiou, T.C., Georgiou, G.C. and Alexandrou, A.N. (1999) Viscous Fluid Flow, New York: CRC Press. 124 Pivkin , V. and Karniadakis, E. (2006) 'Coarse-graining limits in open and wall-bounded dissipative particle dynamics systems', THE JOURNAL OF CHEMICAL PHYSICS, vol. 124. Pivkin, I.V. and Karniadakis, G.E. (2005) 'A new method to impose no-slip boundary conditions in dissipative particle dynamics', Journal of Computational Physics, vol. 207, pp. 114-128. Pivkin, I. and Karniadakis, G. (2006) 'Controlling Density Fluctuations in Wall-Bounded Dissipative Particle Dynamics Systems', Physical Review Letters. Putz, K.A. (1998) 'Optimization techniques for parallel molecular dynamics using domain decomposition', Computer Physics Communications. Randles, P.W. and Libersky, L.D. (1996) 'Smoothed Particle Hydrodynamics some recent improvements and applications', Compu. Methods Appl. Mech. Engrg., vol. 139, pp. 375-408. Rapaport, D.C. (1995) The Art of Molecular Dynamics Simulation, Cambridge University Press. Revenga, M., Zuniga, I. and Espanol, P. (1999) 'Boundary Conditions in Dissipative Particle Dynamics', Computer Physics Communications, pp. 309-311. Ripoll, M. and Ernst, M.H. (2004) 'Model system for classical fluids out of equilibirum', Physical Review E, vol. 71. Ripoll, M., Ernst, M.H. and Espanol, P. (2001) 'Large Scale and Mesoscopic Hydrodynamics for Dissipative Particle Dynamics', Journal of Chemical Physics, vol. 115, no. 15. Ripoll, M., Espanol, P. and Ernst, M.H. (1998) 'Dissipative Particle Dynamics with Energy Conservation: Heat Conduction', International Journal of Modern Physics C, vol. 09, no. 08. 125 Serrano, M. (2006) 'Comparison between smoothed dissipative particle dynamics and Voronoi fluid particle model in a shear stationary flow', vol. 362. Serrano, M. and Espanol, P. (2001) 'Thermodynamically Consistent mesoscopic fluid particle model', Physical Review E. Smiatek, J., Allen, M.P. and Schmid, F. (2008) 'Tunable-slip boundaries for coarse-grained simulations of fluid flow', PHYSICAL JOURNAL E. Symeonidis, V. and Karniadakis, G.E. (2006) 'A family of time staggered schemes for integrating hybrid DPD models for polymers: Algorithms and Applications', Journal of Computational Physics. Symeonidis, V., Karniadakis, G.E. and Caswell, B. (2005) 'Dissipative Particle Dynamics Simulations of Polymer Chains: Scaling Laws and Shearing Response Compared to DNA Experiments', Phys. Rev. Lett. Symeonidis, V., Karniadakis, G.E. and Caswell, B. (2005) 'Simulation of λ-phage DNA in microchannels using dissipative particle dynamics', Bulletin of the Polish Academy of Sciences. Symeonidis, V., Karniadakis, G.E. and Caswell, B. (2006) 'Schmidt Number Effects in Dissipative Particle Dynamics Simulation of Polymers', Journal of Chem. Phys. Tetrode, H. (1912) 'The Chemical Constant of Gases and the Elementary Quantum of Action', Annalen der Physik, vol. 38, pp. 434-442. Thieulot, C., Janssen, L.P. and Espanol, P. (2005) 'Smoothed particle hydrodynamics model for phase separating fluid mixtures', PHYSICAL REVIEW E, vol. 72. Tiwari, A. and Abraham, J. (2006) 'Dissipative Particle Dynamics model for two phase flows', Physical Review. 126 Tiwari, A. and Abraham, J. (2006) 'Dissipative Particle Dynamics Simulations of Liquid Nanojet Breakup', Microfluid Nanofluid. Vázquez-Quesada, A. (2009) Smoothed Dissipative Particle Dynamics, the National University of Distance Education. Vazquez-Quesada, A. and Ellero, M. (2012) 'SPH simulations of a viscoelastic flow around a periodic array of cylinders confined in a channel', Journal of Non-Newtonian Fluid Mechanics, vol. 1, no. 8, pp. 167-168. Vazquez-Quesada, A., Ellero, M. and Espanol, P. (2009) 'Consistent scaling of thermal fluctuations in smoothed dissipative particle dynamics', vol. 130. Vazquez-Quesada, A., Ellero, M. and Espanol, P. (2009) 'Smoothed particle hydrodynamic model for viscoelastic fluids with thermal fluctuations', PHYSICAL REVIEW E, vol. 79. Vincenti, W.G. and Kruger, C.H. (1975) Introduction to Physical Gas Dynamics, Krieger Pub Co. Wang, L., Ge, W. and Li, J. (2006) 'A new wall boundary condition in particle methods', Computer Physics Communications, vol. 174, pp. 386-390. White, F.M. (1974) Viscous Fluid Flow, McGraw-Hill Inc. Whitworth, A.P., Bhattal, A.S., Turner, J.A. and Watkins, S.J. (1995) 'Estimating Density in Smoothed Particle Hydrodynamics', Aston. Astrophys., vol. 301, pp. 929-932. Willemsen, S.M., Hoefsloot, H.C.J. and Iedema, P.D. (2000) 'No-slip Boundary Condition in Dissipative Particle Dynamics', International Journal of Modern Physics C, vol. 11, no. 5, pp. 881-890. 127 Xu, Z. and Meakin, P. (2009) 'A phase-field approach to no-slip boundary conditions in dissipative particle dynamics and other particle models for fluid flow in geometrically complex confined systems', THE JOURNAL OF CHEMICAL PHYSICS, vol. 130. Yang, J., Gatsonis, N. and Potami, R. (2013) 'Multiscale Smooth Dissipative Particle Simulation of Non-Isothermal Flows', SIAM Conference on Computational Science and Engineering, Boston. Yang, J., Potami, R. and Gatsonis, N.A. (2011) 'Wall Boundary Conditions in 3D Smooth Dissipative Particle Dynamics Simulations', American Physical Society Division of Fluid Dynamics, Baltimore. Yang, J., Potami, R. and Gatsonis, N.A. (2012) 'Non-isothermal 3D SDPD Simulations', American Physical Society 65th Annual Meeting of the APS Division of Fluid Dynamics, San Diego. 128
APPENDIX A. SDPD Analytic Self-diffusion Coefficient
We follow the derivation procedure introduced by Litvinov et al. (2009) to form the 3D SDPD self-diffusion coefficient with Lucy function as the interpolation function. Review the SDPD equations, the hydrodynamic momentum equation in a Lagrangian description is given by de Groot and Mazur (1962) d Pdt v v v . (A.1) According to the derivation of self-diffusion coefficient for the SDPD fluid particles in Litvinov et al. (2009), one start from the following by neglecting the conservative forces i D Ri j i ji i ddt m m v F F . (A.2) Since the dissipative part is linear in velocity difference, the Eq (A.2) can be written as i i Ri ddt m v v F , (A.3) where, the coefficient iji j i i j ij ij Wm d d r r . (A.4) Assume the density and temperature are uniformly distributed, one can write Eq. (A.4) as following iji ji i ij ij Wm d r r , (A.5) where, i i i m d . For three dimensional cases, the kernel takes the form ijijij hW W fh h rr . (A.6) Replacing the summation in Eq. (A.5) with integration, the procedure is as following 129
22 20 00 0 iji j ij iji ii ij ij ii j i ji iii ii
Wr rm dV F F V F r dVr F r dr r drrF rW r drrW r i W r dr
Then the coefficient would be i W r dr (A.7) For Lucy function, we have
105 105 211 3 1 1 3 116 16 8 rr r r rW r dr dr dh h h hh h h h W r dr h (A.8) yield i h (A.9) The solution of Eq. (A.3) gives the following expression for the diffusion coefficient B k TD m (A.10) Substitute Eq. (A.9) to Eq.(A.10), we obtain B ii h k TD m (A.11) For incompressible flow, we set 130 (A.12) Then, Eq.(A.11) can be written as B ii h k TD m (A.13) 131 APPENDIX B. Boundary Method in 3D with Lucy’s Smoothing Function
In order to figure out a better way to correct truncation error, we used Vazquez-Quesada boundary condition method (2009a) and derived the formula for 3D using the Lucy function. The procedure is shown bellow. For 3D SDPD, the corrected density i d of the particle i is given by: i ij i ij r ij fluid k wall j fluid d W W W d n W r r r r r r (A.14) where denotes the missing volume of the particle i or the region belongs to the wall in the smoothing region of particle i . And n r is the distribution of wall particles. kk wall n r r r (A.15) If n r is equal to the corrected density i d , then i ij r ij fluid ij i ij fluid d W d n WW d d W r r r rr r r r (A.16) ijj fluid ii i ii W dd hd W rr r r v (A.17) where i i h d W r r r (A.18) Graph function i h
132 Integral equation (A.18)in the missing region arccos sin h d W d d drr W r r V r r (A.19) Set the origin at the particle i and i r rR H , where H is the smoothing length. Then arccos cosarccos cosarccos cosarccos sin sin sin sin h d d dR H RH W R105d d dR H RH 1 3R 1 R16 H1052 H d dRR 1 3R 1 R16 H105 d dRR 1 6 R 8R 38 cosarccos cosarccos cosarccos sin sin sin cos cos cos cos R105 d dR R 6 R 8R 3R8105 R 6 R 4R 3Rd8 3 5 3 7105 4 h 6 h 4h 3hd8 105 3 5 3 7 arccos cos cos cos cos cos
105 4 h 6 h 4h 3hd8 105 3 5 3 7105 4 h 6 h 4h 3hdr8 105 3r 5r 3r 7r105 4r h 3h 4h h8 105 6 r 10r 15r 14r105 for for (A.20) where i h h H For 2D case, we rewrite and correct the derivation 133 arccos cosarccos cosarccos cos cos cos cossin sin sincos cos cos h 1i 0 h h 1 322 0 h 2 4 5 6h 2 4 5 60 2 3 h 2 d dR R W R52 H d dR R 1 3R 1 RH10 1 h 3h 8h hd 10 2 2 5 2h10 2 210 arccos sin sin sinlncos cos cossin sin sincos cos cosarccos ln h4 54 265 3 02 2 2 3 23 22 25 arccos ln arccos ln
This can also be written as arcsin for ln ln for
18h 32h 16 h 1 h 3 6 h1 0 r Hh 6 36 h h 1 1 h0 r H (A.21) which is exactly the same as Vazquez-Quesada i h function. We plot the comparison the value of i h in 2D and 3D. 134 Figure 31: Function i h for Lucy kernel in 2-D and 3-D. Function h Given that hh h (A.22) In three dimension h h hh 1h h h H h
According to (A.20) for for
21 105h 315h 105hh 21h 0 r H16 16 16 16h 0 r H
Then the function h turns into for for (A.23) 135 The derivation is shown below: arccos cosarccos cosarccos cos sin cos sin cos sin cos sin co i i2 h 1 20 0 hh 1 2 250 hh 1 2 30 h h d Fd d dR F R H RH RH3152 d dR 1 R H RH RH4 H315 d dR 1 R R2H315 d2H r r r r r n arccos s cos cos cos
21 105h 315h 105h21h16 16 16 16 which is exactly equation (A.23)
Function i h Ψ After substituting the velocity field of the wall particles, obtaining the irreversible part as i wall i walli i iwall i i i i v V v VP a h b hd h d h Ψ (A.24) h is the equation (A.23), and i h Ψ is a tensor given by i ii i ii i h d F r r r rΨ r r r r r nr r r r (A.25) The tensor i ii i r r r rr r r r in equation (A.25)can be expressed as r rr r , and 136 sin cos sin sin cos sin cos cos sin sin cos sin sin sin cos sin sin cos cos sin cos sin cos r r r r δ δr r δ δ δ δ δ δδ δ δ δ δ δδ δ δ δ δ δ (A.26) where i j δ δ is the unit dyad. We express i h Ψ as i i i ii i i ii i h h h hh h h hh h δ δ δ δ δ δδ δ δ δ δ δ δ δδ δ δ δ Ψ Ψ Ψ ΨΨ Ψ Ψ ΨΨ Ψ (A.27) The integral form is arccos cosarccos cos sin cos sin cos i ii i ii i2 h 1 20 0 h2 h 1 230 0 h h d Fd d dR H RH RH F R315 d d dR R 1 R4 H r r r rΨ r r r r r nr r r r r rr rr rr r r rr r is only the function of and , so arccos cosarccos sin cos sin cos cos cos cos r rΨ r rr rr r (A.28) Integrate equation (A.28) by each dyad. For δ δ arccosarccos sin cos sin cos cos cos coscos sin cos cos cos cossin cos cos
315 1 h 2h hh d d4 H 60 4 5 6315 1 h 2h hd d4 H 60 4 5 6315 1 d4 H 2 δ δ Ψ arccos cos cos cos cos cos Where cos r
315 1 h 2h hh dr 1 r r4 H 60 4r 5r 6 r315 r h 2h hdr 1 r4H 60 4r 5r 6 r315 r h 2h h r h 2h hdr4H 60 4r 5r 6 r 60 4r 5r 6 r δ δ Ψ ln
315 r r h 2h h h 1 2h hdr4H 60 60 4r 5r 6 4 r 5r 6 r315 r r h 2h h h 1 2h hr4H 120 240 4 5r 6 4 2r 15r 24r1 1 2h h h315 120 240 5 12 84H ln
2h h15 24h h h 2h h h 2h hh120 240 4 5 12 8 15 24
Yield, one can get ln
315 1 h 3h 4h h hh h4H 240 24 16 15 24 4 δ δ
Ψ δ δ (A.29) For δ δ arccosarccos sin sin sin cos cos cos cos sin sin cos cos cos cos sin co
315 1 h 2h hh d d4 H 60 4 5 6315 1 h 2h hd d4 H 60 4 5 6315 14 H 2 δ δ Ψ arccos coss cos cos cos cos cos cos
42 4h 2 5 600 5 64 5 61 2 4 5 6h 41 2h ln
2h hr 5r 6 r315 1 h 3h 4h h h h4H 240 24 16 15 24 4
The integral result is ln
315 1 h 3h 4h h hh h4H 240 24 16 15 24 4 δ δ
Ψ δ δ (A.30) which gives i i h h δ δ δ δ Ψ Ψ
For δ δ arccosarccosarcco cos sin cos cos cos cos sin cos cos cos cos cos cos
315 1 h 2h hh d d4 H 60 4 5 6315 1 h 2h hd d4 H 60 4 5 6315 2 d4 H δ δ Ψ s cos cos cos Rearranging ln
315 1 5h 2h h hh h2H 240 16 5 12 4 δ δ
Ψ δ δ (A.31) 139 For δ δ arccos arccos sin sin cos sin cos cos cos cos sin cos sin cos cos cos cos
315 1 h 2h hh d d4 H 60 4 5 6315 1 h 2h hd d4 H 60 4 5 631 δ δ Ψ arccos sin cos cos cos cos for cos r ln ln
315 1 h 2h hh dr r2H 60 4r 5r 6 r315 r h 2h hdr2H 60 4r 5r 6 r315 r h 2h hr2H 240 4 5r 12r315 1 2h h h h2H 240 5 12 240 4 δ δ Ψ
2h hh 5 12
For δ δ i i h h 0 δ δ δ δ
Ψ Ψ (A.32) For δ δ arccos arccos sin cos cos sin cos cos cos cos cos sin cos cos cos cos
315 1 h 2h hh d d4 H 60 4 5 6315 1 h 2h hd d4 H 60 4 5 63154 H δ δ Ψ arccos sin cos cos cos cos
140 For δ δ i i h h 0 δ δ δ δ
Ψ Ψ (A.33) For δ δ arccos arccos sin cos sin sin cos cos cos cos sin sin cos cos cos cos
315 1 h 2h hh d d4 H 60 4 5 6315 1 h 2h hd d4 H 60 4 5 63154 H δ δ Ψ arccos sin cos cos cos cos For dyad δ δ i i h h 0 δ δ δ δ
Ψ Ψ (A.34) Substitute all components of i h Ψ into equation(A.27) ln + ln i i i i i ii i i i2 4 5 6 4 1 12 4 5 6 4 h h h h h hh h h h315 1 h 3h 4h h h h4H 240 24 16 15 24 4315 1 h 3h 4h h h h4H 240 24 16 15 24 4 δ δ δ δ δ δ δ δ δ δδ δ δ δ δ δ δ δ Ψ Ψ Ψ Ψ Ψ ΨΨ Ψ Ψ Ψ δ δ + ln
315 1 5h 2h h h h2H 240 16 5 12 4 δ δδ δ
Here we set nn δ δ , and i h Ψ can be written as i 1 i 2 i h h h Ψ 1 nn (A.35) where 141 ln
315 1 h 3h 4h h h h 0 r Hh 4H 240 24 16 15 24 40 r H (A.36) and ln
315 1 h 13h 8h 5h 3h h 0 r Hh 2H 480 48 32 15 48 80 r H (A.37) The functions i h , h , and h are plotted in Figure 32-Figure 34. Figure 32: Function i h for Lucy kernel in 2-D and 3-D. Figure 33: Function i h for Lucy kernel in 2-D and 3-D. Figure 34: Function i h for Lucy kernel in 2-D (Vazquez-Quesada formulation) and in 3-D. For the 2-D case r rh h hW rrF r r rr rh h hr rh h arccos cosarccos cosarccos cosarccos cos sin sin cos sin cos sin cos h 1 2i 0 hh 12 2 30 h3 h 1 22 24 0 hh 20 h 2 d d RH RH RH F R2 d dR R H F R60H2 d dR R 1 RH120 d dRH δ δ Ψ cosarccos cos cos cos cos cos R 1 R120 1 h h hd 1H 30 3 2 5 cos cos 1120 3 30 30 2 cos 5 3 cos2 cos 5cossin sin sin 1 sin sin sinln3 30 30 90 2 cos 5cos 3cos120 sin4 c hi h h h hh dH h hh h h hH h δ δ Ψ arccos4 5 52 3 0323 4 2 4 2 2 22 2 2 2 4 2 h h h hhh h h h h h h hhH h h h h h h
323 4 2 4 2 2 23 2 2 2 4 2 4 2
120 arccos 1 1 1 3 1ln3 90 4 15 20120 arccos 1 29 1 1 1 1 ln3 90 180 15 4 hh h h h h h h hH hh h h h h h h h hH h arccos 1 20 cosarccos 13 2 30 cos3 arccos 1 23 24 0 cosarccos 23 20 cos
2 cos cos 2 cos 60 2 cos 1120 cos 1 hi hh hh hh h h d d RH RH RH F Rd dR R H F RH d dR R RH d dR R RH δ δ Ψ
120 1 cos 30 3cos 2 cos 5cos120 cos 30 3 2 cos 5cos120 sin sin 1 sin sin ln30 90 3 2 cos 5co hi h h h hh dH h h hdH h h hH δ δ Ψ arccos0323 2 4 2 4 23 2 2 2 2 4 2 4 2 s1120 arccos 1 1 1 1 ln3 30 90 2 5120 arccos 1 1 1 1 1 1 ln3 30 90 90 2 5 h hh h h h h h hH hh h h h h h h h h hH h
120 arccos 1 1 1 1 1 ln3 45 90 5 4 h h h h h h h h hH h arccos 1 29 13 90 180120 01 1 1ln15 40 i h h h h h r Hh H h h h hh r H i h h h h h r Hh H h h h hh r H Figure 35: Function i h for Lucy kernel in 2-D (corrected) and 3-D.
1. Table used in the integral sin sintan 2 cos cos n n n dx ax Cax adx ax n dx forax a n ax n ax
23 24 35 4 26 5 dx x xC Cx xdx x Cx xdx x x Cx x xdx x x Cx x xdx x x x Cx x x xdx xx x x x Cx x x x xax dx C
2. Tensor r rr r
For 3-D The unit vector rr is given as sin cossin sincos rr Then the tensor r rr r is in the following expression sin cossin sin sin cos sin sin coscossin cos sin sin cos sin cos cos sin sin cos sin sin sin cos sinsin cos cos sin cos sin cos r rr r
This can be expressed in the sum of each i j δ δ sin cos sin sin cos sin cos cos sin sin cos sin sin sin cos sin sin cos cos sin cos sin cos i j i ji j r r r r δ δr r δ δ δ δ δ δδ δ δ δ δ δδ δ δ δ δ δ where and
0, 2
For the 2-D case 147 The unit vector rr is given as sincos rr Then the tensor r rr r becomes sin sin coscossin sin cos sin cos cos r rr r
This can be expressed as the sum of each dyad sin sin cos sin cos cos i j i ji j r r r r δ δr r δ δ δ δ δ δ δ δ
Where, and
0, 2