A Final Report for A Massively Parallel Hybrid Dusty-Gasdynamics and Kinetic Direct Simulation Monte Carlo Model for Planetary Applications NASA Grant NAGS-1 3404 from an unsolicited proposal to the Applied Information Systems Research Program for the period May 1,2003 to April 30,2004 Principal Investigator Michael R. Combi Department of Atmospheric, Oceanic and Space Sciences University of Michigan Space Research Building 2455 Hayward Street A m Arbor, MI 48 109 Tel: 1 (734) 764-7226 Fax: 1 (734) 647-3083 Email: [email protected] 1 1 Introduction In order to understand the global structure, dynamics, and physical and chemical processes occurring in the upper atmospheres, exospheres, and ionospheres of the Earth, the other planets, comets and planetary satellites and their interactions with their outer particles and fields environs, it is often necessary to address the fundamentally non-equilibrium aspects of the physical environment. These are regions where complex chemistry, energetics, and electromagnetic field influences are important. Traditional approaches are based largely on hydrodynamic or magnetohydrodynamic (MHD) formulations and are very important and highly useful. However, these methods often have limitations in rarefied physical regimes where the molecular collision rates and ion gyrofrequencies are small and where interactions with ionospheres and upper neutral atmospheres are important. At the University of Michigan we have an established base of experience and expertise in numerical simulations based on particle codes which address these physical regimes. The Principal Investigator, Dr. Michael Combi, has over 20 years of experience in the development of particle-kinetic and hybrid kinetichydrodynamics models and their direct use in data analysis. He has also worked in ground-based and space-based remote observational work and on spacecraft instrument teams. His research has involved studies of cometary atmospheres and ionospheres and their interaction with the solar wind, the neutral gas clouds escaping from Jupiter’s moon Io, the interaction of the atmospheres/ionospheres of Io and Europa with Jupiter’s corotating magnetosphere, as well as Earth’s ionosphere. Grants from the Applied Information Systems Research Program (AISRP) have supported the basic construction and first applications of three-dimensional, neutral and ion kinetic particle models, which have potential application to various space science problems. The purpose of this grant, which resulted from an unsolicited proposal for 1 year of support at $58,868, was to expand the original research and to develop the algorithms and operational code to enable seamless interaction between kinetic and hydrodynamic descriptions. This report describes our progress during the year. The contained in section 2 of this report will serve as the basis of a paper describing the method and its application to the cometary coma that will be continued under a research and analysis grant that supports various applications of theoretical comet models to understanding the inner comae of comets (grant NAGS- 13239 from the Planetary Atmospheres program). 2. Hybrid DSMC/CFD Approach Simulations of real gas flow often require numerical solutions with a wide variation of macroscopic parameters throughout computational domain. In some cases, the gas flow in a part of the domain could be in rarefied regime while another part may be in continuous regime. Rarefied flow conditions appear in different kinds of circumstances in the numerical study of cometary comae (e.g., the Knudsen layer at the surface, on the night side of nuclei, at larger distances from the nucleus and for weak comets in general). Even though the kinetic approach is generally appropriate physically to study gas flows even in the fluid limit (because the fluid equations themselves are approximations of the general kinetic formulation), it would make the numerical simulations of the flow very inefficient computationally. So, developing of efficient numerical approach that can be used for simulation of gas systems under continuum/rarefied condition is an important problem of modern numerical gas dynamics. An approach that utilizes DSMC and CFD methodologies in different subdomains is described in various papers [Bird, 1998; Ivanov et al., 2000; 2 Ivanov and Markelov, 2000; Hash and Hassan, 1996; Roy et al., 2002; Wang and Boyd, 2003; Wang et al., 20021. 2.1 Dusty-Gas Hydrodynamics The appropriate description for the physical state of any one of up to s species, p, in a multispecies dilute gas can be described by the Boltzmann equation, which is given as (Combi et al. 2004): where fp = fp (F,Z,t) and fd f, (T,c',t) are the full phase space velocity distribution functions for species p and q, cq is the velocity, r is the spatial coordinate, represents external forces, cTq is the magnitude of the relative velocity between particle p and q, 52 is the solid angle, the asterisks (*) indicate post-collision particles, the 1q subscript refers to the scattering target particles, and opqi s the total collision cross section between species p and q (which can in general be velocity dependent). In the collision integral on the right hand side of the equation, f* represents additions of particles, scattered into the region of velocity space in question, Le.: hetween F and C. + dc' - whereasfwithout the asterisk represents scattering out of that region. The Boltzmann equation makes neither assumptions nor any restrictions as to the form of the distribution functions. The description of the Euler equations for a dusty-gas cometary coma have been described in detail in many papers, including some by our own group at the University of Michigan (Combi et al. 1999, Korosmezey and Gombosi 1989). Various velocity moment expansions of the integrals over the distribution functions yield equations of conservation of mass, momentum and energy, which can be written for a dusty-gas neutral cometary coma as: -dP + v - (pu)= -'P, at 't du p-++(u*V)u+Vp=-F, at --@ + -(u1 - V)p+ -Y p (V - u) = -egd+ Qph- QIR y-1 c3t y-1 Y-1 respectively, and for any number of dust particle-size populations -d+P;v .(piui)=- 'Pi i = 1, ..., N, at st pi-a i + pi(vi V)V, = Fi i = ..., N df * respectively, where p is the gas mass density, u is the gas velocity, p is the gas pressure, Vi and pi are the velocity and mass density for dust particles of radius ai. The right hand sides 3 6P of all the equations contain the various source terms. The term - is the gas production 6t dP. source rate, and is the dust production source rate for particles of radius, ai. F is the gas- 2 dt dust drag force, which is related to the forces on the individual particle size populations as N F = -CF,, where the size dependent force is given by Fi = -3pi p C&, , where 4ai i-1 Pbi u-v si = k is the Boltzmann constant, Tis the gas temperature and m is the gas mean IjZkTlm’ molecular mass. The accommodation of gas via collisions with hotter dust yields the dust- gas heat exchange rate which is given as Qgd Here C, is the gas heat capacity at constant pressure and the rest of the coefficients can be defined under the assumption of diffusive reflection such that: Finally, pd is the bulk mass density of dust particles of radius ai; T is the gas temperature, Ti is the dust temperature assumed in equilibrium with solar radiation. Otherwise the other intermediate quantities are si, the relative Mach number between gas and dust, CD’, the dust- gas drag coefficient, TF, the recovery temperature, Ri’, the heat transfer function and Sti’. the Stanton number. Most of the dust-gas drag physics comes from the formulation by Finson and Probstein (1 968) with later corrections discussed by Wallis (1 982) and Kitamura (1 986). As described in the paper by Combi et al. (1999), we have successfully applied dusty- gas dynamics to a very collisionally thick comet, Hale-Bopp (1995 01). However, for most comets the gas becomes too rarefied both for productive comets at large distances from the nucleus, as well as for weak comets throughout most of the coma. 4 2.2 Continuum breakdown It is well known that there are two major limiting conditions of applicability of the computational fluid dynamics (CFD) approach. Fluid description of gas media is valid if velocity distribution function of gas molecules is close to Maxwellian and macroscopic parameters changes slowly on the scale of mean fi-ee path. The physical reasons [Boyd, 19951 for breakdown in expanding and compressing flows are different in several ways. In expanding flow, continuum breakdown generally occws when the collision rate becomes so low that thermal equilibrium cannot be maintained among the energy modes. In compressed flows, continuum breakdown occurs in general due to the very steep flow gradients in the shock front, and next to the body surface. A continuum breakdown parameter [ Wang and Boyd, 20021 is generally used as the criterion for switching between methods. In his pioneering work, Bird [ 19651 first advocated a semiempirical parameter for expanding flows. Another empirical parameter based on local flow gradients was later developed specifically for hypersonic flows. More recently, a new breakdown parameter based on the Chapman-Enskog perturbation expansion of the Boltzmann equation has been developed again for expanding flows. The breakdown of equilibrium was found to coincide with a certain value of the ratio of the logarithmic time derivative of density following the motion of the fluid to the collision frequency of the gas. This value had been proposed [lvunov et al., 1YYY; Boyd, 19%] as an empiricai breakdown criterion for use in engineering studies of systems which involve low-density expansion from continuum to highly rarefied conditions For 1D flows, V but v = L,S O A This form of P may be interpreted as the Lagrangian mean free path divided by the density scale length. Breakdown parameter (4) can be evaluated based on other macroscopic parameters of the gas system such as mean flow speed V or translational temperature T . In this case, the value of the breakdown parameter P is determined [ 11 11 as where P, , PT, P , are determined by substituting density number, temperature and velocity for p in the expression (4). The breakdown of the continuum assumption usually coincides with a degree of translational nonequilibrium for which anisotropic effects arise. Bird [Hush und Hussun, 19961 advocates a continuum breakdown parameter (3), with a limit value 5 P between 0.01-0.02 The Knudsen number is the key parameter in determination of the gas flow regime. It is commonly assumed [Boyd, 19951 that failure of the continuum approach occurs at Knudsen numbers around 0.01. Limiting value of Knudsen number 0.1 Kn is proposed in [Hush and Hussun, 19961. A local Knudsen number may be calculated [Boyd, 19951 throughout the flow field based on local flow properties and the length scale of the body. This parameter is termed the body-length local (BLL) Knudsen number: A Kn,,, = - D where h is the local mean free path, and D is a fundamental length scale of the problem. The limitation of the BLL definition is that a constant length scale, D, is employed. The local Knudsen number can better be formulated in terms of length scale determined by flow gradients. The gradient-length local (GLL) Knudsen number can be expressed [Boyd, 19951 in the following form where A again is the local mean free path, Q is a flow property (density or temperature), and 1 is some distance between two points in the flow field. This distance should be taken approximately along the line of the steepest gradients in the flow properties, streamlines or computational grid lines. In practical calculations, it is often more convenient to simplify expression (7) by evaluation [Wung and Boyd, 20021 dQldl as VQ without projecting it onto a preferential direction. Now the GLL Knudsen number for property Q can be written as A Kn, = -1VQl Q The Knudsen number of this form has a great physical meaning. When its value is much less than unity the flow can be regarded as locally slightly perturbed from equilibrium [ 111 1 that is a fundamental assumption of the fluid approach. A new parameter, Kn,, has also been max proposed [Wang and Boyd, 20021 that is modified from Boyd's Kn,, parameter. In the same way as it was done in expression (9,t he value of the local Knudsen number can be obtained based on different macroscopic parameters and define in the form Kn,,, = max(Kn,, Kn,, Kn, ) (9) In the rarefied gas dynamics literature it is now generally accepted that continuum breakdown is best predicted wherever the value of Kn,, exceeds 0.05. An entropy generation rate has been proposed [Chen et ul., 20031 as a continuum breakdown parameter since it has a stronger physical foundation than Knudsen numbers. The entropy generation is defined as s - r, dui 4, dT T dxj T2 axj gen Two non-dimensional parameters of entropy generation have been introduced for axisymmetric problems: where the former uses local and the latter uses global free stream conditions. The global parameter threshold is proposed to be 1, while the local parameter threshold is proposed to be 5 ~ 1 0 ~ . A Grad criterion has been proposed [Tallis and Mullinger, 19971 for determination of the Navier-Stokes (NS) solution. For a given known NS solution, moments of the Boltzmann equation, KG = (l,E,E - E,lE2k), need to be evaluated for the corresponding Grad Chapman-Enskeg distrihvtinn fiinction; giving the residual as where If the NS distribution is a good approximation of the distribution for the general Boltzmann equation, the value of the residualRNS will be close to zero, and the NS solution will be valid everywhere the residual is small. Similar criteria can be formulated to evaluate the discrepancy of distribution function fiom a Maxwellian which is useful for determination the area of validity of Euler equations. 2.3 The DSMCKFD Boundary and Interface Coupling kinetic and fluid approaches basically involves [Hush and Hassan, 19961 determination of the properties at the interface between the DSMC and CFD regions. The macroscopic properties are required by the conservation equations for the evaluation of the net fluxes and are necessary for the DSMC method for the specification of the distribution function used to initialize particles entering the DSMC region. At the end of each hybrid time step, particles enter from the fluid region, such that mass, momentum and energy fluxes are conserved along the interface. Several methods of maintaining the interface between these regions have been proposed. Marshak condition The most obvious way of coupling DSMC and CFD approaches is based on equating of half- fluxes at the interface. The half-fluxes can be defines as the flux of particles which velocities 7 comprising only of the normal velocity component. The DSMC half-fluxes are - 1 F”,, mivi AAT -(v*v) -2 On the fluid side, the half-fluxes are m m m FcFD= n j j $ @ f d u dv dq 1 Q= mivi m -(v v) * 2 The net flux is the sum of the incoming and outgoing half-fluxes. For the Marshak condition, the half-fluxes of the mass, momentum and energy of the particles crossing the interface during the DSMC step are measured. Theses half-fluxes are then added to the incoming half- fluxes from the fluid side and equated to the net flux. When assuming a Maxwellian distribution on the fluid side, this gives five equations and five unknowns, providing the macroscopic properties of density, velocity and temperature at the interface. Extrapolation offlow properties This approach of maintaining the interface between DSMC and CFD subdomains involves extrapolation of the macroscopic properties of density, velocity and temperature to the interface. The accuracy obtained is dependent on the order of extrapolation. The DSMC solution is allowed to evolve while accumulating samples for the determination of the macroscopic properties [Hash and Hassan, 19961. After a specified number of time steps, the solution is stopped, and the velocity and temperature are evaluated. At the first coupling, initial guesses for the velocity and temperature at the interface are taken from the full fluid solution. Subsequent initial guesses for the interface values are made using the values from the previous coupling. The two solutions are then smoothed across the interface region and new values for the interface are obtained. The solution for the fluid region is repeated, and the process iterated upon until the interface values converge. Because a good smoothing routine is required, a quadratic least-squares fit for the profiles across the interface is employed. High-order fits should be avoided because of the oscillations they can produce. Extrapolation of the net fluxes Instead of extrapolating the macroscopic parameters of the gas flow, net fluxes can be extrapolated across the interface. The main difficulty with this is that evaluation of the net flux constants must be very accurate in order to ensure that the calculated values conform 8 with the profile predicted by the DSMC method. In practice, the scatter in higher moments is so great that the method is essentially useless. Example of implementation of the method is presented in [Roy et al., 2002.1 where a cell-centered finite-volume approach is used for discretization of governing equations and second-order reconstructions of the interface fluxes is obtained via MUSCL extrapolation. Overlapping subdomains The approaches considered above are based on separation of the whole computational domain into a set of non-overlapping subdomains where an appropriate gas dynamics description is applied and where either macroscopic parameters or fluxes of the flow are used to couple the subdomains. Another approach, where kinetic and fluid subdomains overlap each other has been described in [Boyd et uZ., 20031. Figure. 1 Schematic of the interface between continuum and kinetic regions for the case of overlapping sub-domains Buffer cells and reservoir cells are used in the continuum domain as illustrated in Fig. 1. These buffer cells and reservoir cells are also treated as particle cells for particle movement and particle collisions. Then the interface becomes the internal cell edge for the DSMC. The reservoir cells are used to generate particles that can enter the DSMC domain, which avoids directly generating particles on the interface. The buffer cells, however, improve the quality of the particles that enter the DSMC. As a consequence, Chapman- Enskog or Maxwellian distribution functions can be assumed in the fluid region described by Navier-Stokes or Euler equations, respectively. Particles are generated in the buffer cells according to the distribution function based on the local macroscopic information. In the beginning of each time step, particles are generated in the reservoir cells according to the distribution function and at the end of the iteration all particles in the reservoir cells are removed. Particles in the buffer cells and reservoir cells are selected for collisions and move around using general DSMC procedures. The particles leaving particle cells (including buffer cells and reservoir cells) are removed. Application of the Hybrid DSMCEuler method to comet 21P/Borrelly In September 2001, NASA's Deep Space 1 spacecraft flew by comet 21P/Borrelly and obtained only the first set of close-up images and spectra of a cometary nucleus since the flyby of comet lP/Halley in 1986. Longer exposure images centered on the nucleus revealed very interesting, structured sets of dust jets emanating for various source regions on the dayside of the nucleus. Our first application for the hybrid DSMCEuler method is 1D- spherical gas coma relevant for the nucleus size and gas production rate of comet Borrelly. 9 We first solved for the flow using the Euler solver to get initial conditions for the hybrid DSMUEULER simulation, where the problem is solved in fluid approximation throughout whole domain. Even though the value of Knudsen number indicates that gas in some parts of the domain is under rarefied regime, this approach allows us efficiently generate initial conditions for hybrid simulation which would be “close” to converged solution. It also provides an interesting point of comparison to illustrate the weakness of the fluid approximation in some regions. To speed up convergence, simulations begin with a coarse mesh (i.e., not many computational cells and particles) and subsequently the mesh is refined and additional simulation particles added appropriately. To couple the DSMC and EULER subdomains, we used interpolation of macroscopic parameters on the boundary between subdomains as discussed above. At the start of the hybrid simulation, the EULER-solution was used for interpolation. Interpolated values were used to generate model particles that enter DSMC region and to calculate fluxes on the EULER side, where a standard MUSCL finite volume scheme is used. Model particles on the boundary were generated with a Maxwellian distribution function. The main difficulties in DSMUEULER hybrid simulation is due to statistical nature of DSMC method where a “large” sample in microscopic parameters is required to get acceptable noise level in the macroscopic flow field. In this work, different numbers of iterations were used on EULER and DSMC sides. While iterations on the EULER side are required only for convergence of the fluid solution, on DSMC side solution is iterated for decreasing of noise level also. Table. Physical parameters for Comet Borrelly for the ID spherically symmetrical case. Uniform Maxwellian at the surface Pure water coma Production rate of 102 9 s- 1 Zero bulk surface velocity Surface temperature = 180 K Model particles that cross boundary of computational domain are deleted. Euler domain Kn < 0.03 DSMC domain Kn >0.03 A Knudsen number in the form Kn, = -1VQI Q Figure 2 shows our first results using the new hybrid model compared with Euler hydrodynamics alone. The most obvious differences are (as should be expected) in the expanding coma. In the hybrid model the isotropic Maxwellian distribution is clearly not present, as evidenced by the divergence of the radial and surface-normal components of the temperature. As in the standard rocket nozzle problem, the radial temperature freezes out at a non-zero value. The additional effect of this is that the asymptotic radial expansion velocity is smaller in the hybrid model than in the Euler calculation. This is caused by the fact that the radial temperature (Le., the radial dispersion in velocities) is not transferred to directed radial motion. Another interesting result is subtler. In comparing a number of curves (but 10