440 likes | 558 Views
Hybrid continuum-particle simulations for vacuum flows Sarantis Pantazis Physikalisch-Technische Bundesanstalt (PTB), Berlin School on Vacuum Gas Dynamics via Kinetic Theory Volos, 07-11 May 2012. Introduction.
E N D
Hybrid continuum-particle simulations for vacuum flows Sarantis Pantazis Physikalisch-Technische Bundesanstalt (PTB), Berlin School on Vacuum Gas Dynamics via Kinetic Theory Volos, 07-11 May 2012
Introduction In problems of practical interest, large variations of gas properties (such as pressure, density or temperature) frequently appear in the flow field. As a result, the local Knudsen number may have a range spanning several orders of magnitude, from well within thehydrodynamic regime up to the free molecular.
Introduction The Navier-Stokes equations, along with the equations of mass and energy conservation, are sufficient for the description of gas flow in a variety of situations. However, they are not valid for vacuum flows due to the violation of the continuum assumption and the breakdown of the Newton-Fourier-Fick constitutive laws. Therefore, it is physically inaccurate to use traditional CFD methods. This approach would lead to significant modelling errors for Kn>0.01 (or Kn>0.1 if appropriate BCs are used).
Introduction An approach based on the Kinetic Theory of Gases would be more suitable. The most commonly used methods involve • the solution of a system of kinetic equations, e.g. by the Discrete Velocity Method (DVM) • the evolution of a particle cloud, e.g. Molecular Dynamics (MD) or Direct Simulation Monte Carlo (DSMC) These methods are applicable for any value of the Knudsen number.
DSMC method The algorithm consists of the following steps: • Initial conditions • Free molecular flight and wall interactions • Simulation of molecular collisions • Indexation of particles • Sampling of macroscopic quantities Steps 2-5 are repeated at each time step. The iteration process continues until the statistical fluctuation is reduced to a certain level for steady state problems. The whole simulation is repeated multiple times for time dependent problems (ensemble averaging).
DSMC method However, significant difficulties appear when the DSMC method is applied in less rarefied gases. As the mean free path decreases: • the cell size must decrease (to avoid collisions with distant molecules) • the total number of particles must increase (to maintain a constant ratio of particles per cell) • the time step must decrease (to avoid “missing” collisions) • extended sampling must take place (to reduce noise) The factors stated above increase the computational effort significantly and practically prohibit the use of simple DSMC schemes for several practical applications.
Hybrid schemes The main concept behind hybrid schemes is to decompose the flow domain in different parts, according to the degree of non-equilibrium, and simulate each region with either traditional CFD or a kinetic approach. Our goal is to perform the domain decomposition and algorithm integration such that a reasonable compromise between accuracy and computational efficiency is reached. There exist several types of hybrid approaches. The most successful approaches involve the Navier-Stokes equations and the DSMC method. This type of simulations is very well known in the field of external, hypersonic, rarefied gas flows but may also be applied for vacuum flows.
CFD-DSMC hybrid The continuum/particle descriptions of the gas are performed as usual at the corresponding domains and results are combined at the interface, chosen based on the local conditions. The most important factors are: • Distribution function of incoming molecules. Either the Maxwell-Boltzmann or the Chapman-Enskog distribution function is chosen. • Correct determination of the interface location. The decomposition criterion must be a measure of the divergence from equilibrium at each position to ensure the accuracy of the physical description. At the same time, it is crucial for efficiency reasons to avoid using the DSMC method in places where it is not needed. • Proper exchange of interfacial information. This affects the properties of molecules entering the DSMC domain (distribution function, possible buffer zones, etc) and the macroscopic quantities for the CFD solver.
CFD-DSMC hybrid • Chapman-Enskog distribution • Description • Generation • Advantages and disadvantages • Interface location • Manual decomposition • Bird’s parameter • GLL Knudsen number • Thermal non-equilibrium • Chapman-Enskog criterion • Boyd criterion • Interface information exchange • Buffer Zones • Overlapping regions
C-E distribution Maxwell-Boltzmann distribution near equilibrium As we move away from equilibrium, the Chapman-Enskog expansion to 1st order gives [1] [1] A.L. Garcia, B.J. Alder, “Generation of the Chapman-Enskog distribution”, J. Com. Phys., 140, pp. 66-70 (1998).
C-E distribution Assume an envelope function g(c) Acceptance-Rejection algorithm: • If , we can simplify the acceptance relation to • Accurate and efficient choice[1]: A=1+30B, where [1] A.L. Garcia, B.J. Alder , “Generation of the Chapman-Enskog distribution”, J. Com. Phys., 140, pp. 66-70 (1998). Select ctry uniformly in [0, cmax] Select ftry uniformly in [0, g(ctry)] If ftry less or equal to f(ctry) accept
C-E distribution Advantages: • Generally considered to be more accurate than the Maxwell-Boltzmann distribution Disadvantages: • Influence on CPU time [2] • Possibly negative values [3] • Not significant improvement on results [3, 4] • [2] J.-S. Wu, Y.-Y. Lian, G. Cheng, R.P. Koomullil, K.-C. Tseng, “Development and verification of a coupled DSMC-NS scheme using unstructured mesh”, J.Comp.Phys., 219, pp. 579-607 (2006). • [3] P.V. Vashchenkov, A.N. Kudryavtsev, D.V. Khotyanovsky, M.S. Ivanov, “DSMC and Navier–Stokes study of backflow for nozzle plumes expanding into vacuum”, 24th Rarefied Gas Dynamics international Symposium, AIP conf. proc., 762, pp. 355–359 (2005). • [4] D.C. Wadsworth, D.A. Erwin, “Two-dimensional hybrid continuum/particle approach for rarefied flows”, 23rd Plasma Dynamics & Lasers Conference, AIAA-92-2975, pp. 1–7 (1992).
CFD-DSMC hybrid • Chapman-Enskog distribution • Description • Generation • Advantages and disadvantages • Interface location • Manual decomposition • Bird’s parameter • GLL Knudsen number • Thermal non-equilibrium • Chapman-Enskog criterion • Boyd criterion • Interface information exchange • Buffer Zones • Overlapping regions
Manual decomposition [6] F. La Torre, S. Kenjeres, J.-L. Moerel, C.R. Kleijn, “Hybrid simulations of rarefied supersonic gas flows in micro-nozzles”, Computers & Fluids 49, pp. 312-322 (2011). [5] O. Aktas, N.R. Aluru, “A combined continuum/DSMC technique for multiscale analysis of microfluidic filters”, J.Comp.Phys., 178, pp. 342-372 (2002).
Bird’s parameter Bird [7] proposed the following criterion for continuum breakdown in expanding flows: where ν is the collision rate, ρ is the density, D is the substantive derivative, s is the speed ratio. The threshold is 0.04. The second expression involving just the spatial derivative is for steady flows. The criterion was obtained from DSMC computations. [7] G.A. Bird, “Breakdown of translational and rotational equilibrium in gaseous expansions”, AIAA Journal, 8(11), pp. 1998-2003 (1970).
GLL Knudsen number As noted in Bird [8], the Gradient Length Local Knudsen number can be calculated using a reference length equal to the local gradient length of a quantity Q It can then serve as a decomposition criterion after taking the maximum value for several quantities [2]. A commonly chosen threshold is 0.02-0.05 [2,9]. [2] J.-S. Wu, Y.-Y. Lian, G. Cheng, R.P. Koomullil, K.-C. Tseng, “Development and verification of a coupled DSMC-NS scheme using unstructured mesh”, J.Comp.Phys., 219, pp. 579-607 (2006). • [8] G.A. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Oxford University Press, Oxford (1994). [9] W.L. Wang, I.D. Boyd, “Hybrid DSMC-CFD simulations of hypersonic flow over sharp and blunt bodies”, AIAA paper 3644 (2003)
Thermal non-equilibrium Thermal non-equilibrium [2] between translational and rotational temperatures (even vibrational for high temperatures) Serves as an additional criterion. A threshold value is around 0.03. [2] J.-S. Wu, Y.-Y. Lian, G. Cheng, R.P. Koomullil, K.-C. Tseng, “Development and verification of a coupled DSMC-NS scheme using unstructured mesh”, J.Comp.Phys., 219, pp. 579-607 (2006).
Chapman-Enskog criterion The C-E parameter may also be used as a breakdown criterion. Alternatively, Garcia et al. [10] proposed a simplified C-E parameter, such as The threshold in this case is 0.1-0.2. [10] A.L. Garcia, J.B. Bell, W.Y. Crutchfield, B.J. Alder, “Adaptive mesh and algorithm refinement using Direct Simulation Monte Carlo”, J.Comp.Phys., 154, pp. 134-155 (1999).
Boyd criterion The heat flux tensor is calculated by and the following relations must be valid in one-dimension [11] in order to approach continuum flow conditions. [11] I.D. Boyd, “Predicting breakdown of the continuum equations under rarefied flow conditions”, 23rd Rarefied Gas Dynamics international Symposium, AIP conf. proc., 663, pp. 899–906 (2003)
CFD-DSMC hybrid • Chapman-Enskog distribution • Description • Generation • Advantages and disadvantages • Interface location • Manual decomposition • Bird’s parameter • GLL Knudsen number • Thermal non-equilibrium • Chapman-Enskog criterion • Boyd criterion • Interface information exchange • Buffer Zones • Overlapping regions
Interface A layer of a few cells is often selected to act as an interface between the two methods. It is also called a buffer zone/overlap region/reservoir cells depending on the way it is used. A common characteristic between various implementations is that both methods are applied in that region. Also, the larger the size of the interface region, the smoother the transition should be. However, this practically means increasing the size of the DSMC domain and as a result the simulations become more demanding. A reasonable compromise must be found.
Interface DSMC region and buffer cells embedded within a continuum mesh [10]. Particles created in buffer zone, then move for one time step and finally get discarded if they remain in the buffer. [10] A.L. Garcia, J.B. Bell, W.Y. Crutchfield, B.J. Alder, “Adaptive mesh and algorithm refinement using Direct Simulation Monte Carlo”, J.Comp.Phys., 154, pp. 134-155 (1999).
Interface • A: DSMC region • B,C: interface • D: hydro • B: updated from DSMC • C: updated from NS • AB: BC for NS • CD: BC for DSMC [2] J.-S. Wu, Y.-Y. Lian, G. Cheng, R.P. Koomullil, K.-C. Tseng, “Development and verification of a coupled DSMC-NS scheme using unstructured mesh”, J.Comp.Phys., 219, pp. 579-607 (2006).
Interface Typical examples of the overlapping regions decomposition in the 1st and 15th iteration [2].
Solver characteristics Hybrid solvers may be categorized according to their algorithms and their characteristics: • Static/Dynamic • One-way/Two-way • State coupling/Flux coupling This distinction is important for both the accuracy and efficiency of the algorithms.
Static vs dynamic Static solver: The interface position is determined in the beginning and kept constant throughout the simulation of each region. However, the decomposition may be modified after each iteration step, when switching between the hydrodynamic and particle solvers. Static decomposition is frequently encountered in one-way simulations of steady state problems. Dynamic solver: The interface position changes during the simulation according to some criterion, such as the GLL Knudsen number. Dynamic decomposition is considered to be more favourable for time-dependent flows. It is more complex to apply but offers more flexibility.
One-way vs two-way One-way coupling (weak coupling): Information is not exchanged during a simulation. In order to communicate the necessary data, the simulation must be completed and data are then transfered to the other solver to serve as boundary conditions. Two-way coupling (strong coupling): The two algorithms run in a single simulation but one after the other at each time step. The necessary information is exchanged at the interface at each time step. One-way solvers imply aniterative procedure of simulations with alternating methods. This is very advantageous when only commercial (closed-source) packages are only available but renders time-dependent flows practically impossible.
State vs flux coupling State-based coupling (Schwarz technique): State properties are imposed at the interface as Dirichlet boundary conditions at the interface, according to the values obtained from the other solver. There are also a few works with Neumann/Robin boundary conditions in a similar fashion. State coupling suffers less from statistical noise [12] and is therefore more frequently used. Flux-based coupling: particles crossing the interface contribute to the fluxes through the corresponding faces. The fluxes may not be the same and may need to be modified to maintain conservation [13]. [12] H.S. Wijesinghe, N.G. Hadjiconstantinou, “Hybrid atomistic-continuum methods for multiscale hydrodynamics”, Int.J.Mult.Comp.Eng., 2(2), pp.189-202 (2008). [13] I.D. Boyd, T.R. Deschenes, “Hybrid particle-continuum numerical methods for aerospace applications”. Schematic of typical hybrid coupling techniques [13]
A simple implementation Solve the NS Eqs for the whole domain A simple implementation would be a static, one-way, state coupling solver Advantages: • Simplicity • Easily applicable with existing solvers • Numerical stability Disadvantages: • Iterative procedure (the process is repeated around 10 times) • Practically impossible for time-dependent flows Use the solution to find interface position Extract BCs for the rarefied domain at the interface Apply DSMC on the rarefied domain Use DSMC solution to determine BCs for the NS domain Solve the NS Eqs again for the corresponding part
A simple implementation Convergence history of density and total temperature [2]
Advanced schemes Determine interface position from initial conditions Dynamic, two-way, state-flux coupling solvers should be more flexible. Advantages: • More appropriate for time-dependent flows Disadvantages: • Programming difficulties • More sensitive to the interface treatment Solve the NS Eqs in their domain Exchange interface information Apply DSMC in the rarefied domain If time-dependent, repeat DSMC Exchange interface information Every few time steps re-adjust interface position. Advance time
Work in the literature Some recent papers are shown in chronological order and their characteristics are given next: [10] A.L. Garcia, J.B. Bell, W.Y. Crutchfield, B.J. Alder, “Adaptive mesh and algorithm refinement using Direct Simulation Monte Carlo”, J.Comp.Phys., 154, pp. 134-155 (1999). [16] Q.Sun, I.D. Boyd, G.V. Candler, “A hybrid continuum/particle approach for micro-scale gas flows”, 23rd Rarefied Gas Dynamics international Symposium, AIP conf. proc., 663, pp. 752–759 (2003). [5] O. Aktas, N.R. Aluru, “A combined continuum/DSMC technique for multiscale analysis of microfluidic filters”, J.Comp.Phys., 178, pp. 342-372 (2002). [3] P.V. Vashchenkov, A.N. Kudryavtsev, D.V. Khotyanovsky, M.S. Ivanov, “DSMC and Navier-Stokes study of backflow for nozzle plumes expanding into vacuum”, 24rd Rarefied Gas Dynamics international Symposium, AIP conf. proc., 762, pp. 355–360 (2005). [2] J.-S. Wu, Y.-Y. Lian, G. Cheng, R.P. Koomullil, K.-C. Tseng, “Development and verification of a coupled DSMC-NS scheme using unstructured mesh”, J.Comp.Phys., 219, pp. 579-607 (2006). [12] H.S. Wijesinghe, N.G. Hadjiconstantinou, “Hybrid atomistic-continuum methods for multiscale hydrodynamics”, Int.J.Mult.Comp.Eng., 2(2), pp.189-202 (2008). [14] A. Donev, J.B. Bell, A.L. Garcia, B.J. Alder, “A hybrid particle-continuum method for hydrodynamics of complex fluids”, Multiscale Model. Simul., 8(3), pp. 871-911 (2010). [6] F. La Torre, S. Kenjeres, J.-L. Moerel, C.R. Kleijn, “Hybrid simulations of rarefied supersonic gas flows in micro-nozzles”, Computers & Fluids, 49(1), pp. 312-322 (2011). [15] I.D. Boyd, “Hybrid particle-continuum methods for nonequilibrium gas and plasma flows”, 27th Rarefied Gas Dynamics international Symposium, AIP conf. proc., 1333, pp. 531-538 (2011).
Other improvements • Adaptive Mesh and Algorithm Refinement (AMAR) [10]: The DSMC domain is integrated as the last refinement level in a hydrodynamic solver. • Different time step for each method [10]. • Special treatments to ensure conservation principles [10]. • Variable time step [2]. [2] J.-S. Wu, Y.-Y. Lian, G. Cheng, R.P. Koomullil, K.-C. Tseng, “Development and verification of a coupled DSMC-NS scheme using unstructured mesh”, J.Comp.Phys., 219, pp. 579-607 (2006). [10] A.L. Garcia, J.B. Bell, W.Y. Crutchfield, B.J. Alder, “Adaptive mesh and algorithm refinement using Direct Simulation Monte Carlo”, J.Comp.Phys., 154, pp. 134-155 (1999).
Information Preservation The Information Preservation method [17] is a promising modification of the DSMC algorithm which is used to significantly reduce the noise in low Kn applications. Each particle is assigned an additional, “information” velocity, expressing the collective velocity of the real molecules that the particle represents. As a result, the sampling and collision processes slightly change. The main benefit lies in the fact that through this method, the sampling requirements are significantly reduced in comparison to classical DSMC. It has been seen [16] that its use in a hybrid solver results in less statistical noise. Therefore, it is expected that this approach as well as other noise reduction techniques may allow time-dependent simulations in the near future. [16] Q.Sun, I.D. Boyd, G.V. Candler, “A hybrid continuum/particle approach for micro-scale gas flows”, 23rd Rarefied Gas Dynamics international Symposium, AIP conf. proc., 663, pp. 752–759 (2003). [17] J. Fan, C. Shen, “Statistical simulation of low-speed rarefied gas flows”, J. Comp. Phys., 167, pp. 393-412 (2001).
Computational gain The benefits of a hybrid approach in terms of computational time varies. A CPU time reduction of 75-95% has been reported in [6], around 25% in [2] and 2 to 100 times in [5]. In other works [18] there is practically no computational gain or even an increase in CPU time instead [9]. However, the focus of the latter works was accuracy instead of efficiency. Clearly, the computational effort highly depends on the algorithm characteristics and implementation, as well as the nature of the problem itself. [2] J.-S. Wu, Y.-Y. Lian, G. Cheng, R.P. Koomullil, K.-C. Tseng, “Development and verification of a coupled DSMC-NS scheme using unstructured mesh”, J.Comp.Phys., 219, pp. 579-607 (2006). [5] O. Aktas, N.R. Aluru, “A combined continuum/DSMC technique for multiscale analysis of microfluidic filters”, J.Comp.Phys., 178, pp. 342-372 (2002). [6] F. La Torre, S. Kenjeres, J.-L. Moerel, C.R. Kleijn, Computers & Fluids, 49, pp. 312-322 (2011). [9] W.L. Wang, I.D. Boyd, “Hybrid DSMC-CFD simulations of hypersonic flow over sharp and blunt bodies”, AIAA paper 3644 (2003). [18] J.M. Burt, I.D. Boyd, “A multiscale particle approach for continuum/rarefied flow simulation”, AIAA paper 1184, (2008).
Restrictions: stability Hybrid algorithms inherit some of the advantages and disadvantages of both methods. One of the greatest problems is the possibility of numerical instability. Even though the DSMC method is unconditionally stable, traditional CFD methods are not. Therefore, fluctuations in the Monte Carlo field propagate through the interface and may cause divergence. It is therefore important to keep the statistical noise in the DSMC region as low as possible.
Restrictions: mesh Mesh quality for hydrodynamic simulations is a large field by itself. Cells must have good numerical properties (cell volume, determinant, skewness, etc.) and preferably a hexahedral shape. Even though the DSMC method has no restrictions regarding the cell shape, the cell size must be smaller than the mean free path. It is particularly important to satisfy the mesh quality and cell size restrictions everywhere, especially for dynamic solvers. This is also another factor which may cause divergence.
Restrictions: Time step It is well known that numerical simulations of hydrodynamic equations must be performed with a time step that satisfies the CFL criterion, in order to avoid divergence. The DSMC method on the other hand is unconditionally stable but may return unrealistic results if the time step is too high (more than the particle mean free time). The time step must satisfy simultaneously the time step criteria of both methods.
Restrictions: Time dependency In order to perform time-dependent DSMC simulations, an ensemble average procedure is usually applied. Many realizations of particles are performed and their sampling properties are averaged to extract the macroscopic quantities. In theory, one can also apply only one realization with many more particles to simulate a time-dependent phenomenon. This greatly simplifies the programming of the method but is also subject to hardware limitations. Therefore, it can only be used for very simple problems of low dimensionality. Also, the solver must be dynamic and two-way coupled in order to avoid the iterative solution procedure. Almost all work in the literature is steady state.
Conclusions • Hybrid schemes can significantly reduce the required computational effort. In fact, they are the only available choice for some flow configurations. • Care must be taken in the proper interface of the two regions through a layer of a few cells. • The characteristics of the algorithm play a very important role in the accuracy, efficiency and simplicity of the method. • Some restrictions regarding the mesh form and quality, as well as the time step must be taken into account • Time-dependent problems are still difficult to solve.