Development of Monte Carlo Methods in Hypersonic Aerodynamics

The advantage of Monte Carlo methods in the computational aerodynamics and application of these methods in rarefied fields are described in the present paper. The direct statistical simulation of aerodynamics processes with the solution of kinetic equations is considered. It is shown that the modern stage of the development of computational methods is impossible without the use of the complex approach (its physical nature, mathematical model, the theory of computational mathematics, and stochastic processes) to the development of algorithms. Main directions of the development of the direct simulation of Monte Carlo method in computational aerodynamics are discussed. Some calculation results by using Monte Carlo method are presented.


Introduction
The appearance of direct simulation Monte Carlo methods (DSMC) in various fields of applied mathematics is usually caused by the appearance of qualitatively new practical problems. The examples include the creation of nuclear weapons, space development, the study of atmospheric optics phenomena, and the study of physicochemical and turbulence processes. Some of efficient definition is as follows: The Monte Carlo methods are the methods designed for solving mathematical problems (for example, systems of algebraic, differential, or integral equations) based on the direct statistical simulation of physical, chemical, biological, economic, social, and other processes using the generation of random variables.
The first paper devoted to the Monte Carlo method was published in 1873 [1]. It described the experimental determination of π by a realization of the stochastic process of tossing a needle on a sheet of ruled paper. A striking example is the use of von Neumann's idea to simulate the neutron trajectories in the Los Alamos laboratory in 1940. Although the Monte Carlo methods require a large amount of computations, the absence of computers at that time did not discourage the researchers. As we known, the name of these methods comes from the capital of Monaco, which is famous for its Casino; indeed, the roulettes used in the casino are perfect tools for generating random numbers. The first paper usage of the Monte Carlo method that systematically expanded this method was published in 1949 [2]. In that paper, the Monte Carlo method was used to solve linear integral equations. In Russia, the Monte Carlo methods appeared after the Geneva International Conference on the Peaceful Uses of Atomic Energy [3].
The general scheme of the Monte Carlo method is based on the central limit theorem, which states that, the random quantity equal to the sum of a large number of random variables with the same expectation m and the same dispersions σ 2 . They have the normal distribution with the expectation N and the variance N σ 2 . Assume that we want to solve an equation or find the result of some process I. If we can construct the random variable ξ with the probability density p(x) such that the expectation of this variable is equal to the unknown solution M(ξ) = I, then we obtain a simple method for estimating the solution and its error: This implies the general properties of the Monte Carlo methods as below: 1. The absolute convergence to the solution with the rate 1/N. 2. An unfavorable dependence of the error ε on the number of trials: ε ≈ 1/√N (to reduce the error by an Development of Monte Carlo Methods in Hypersonic Aerodynamics order of magnitude, the number of trials must by increased by two orders of magnitude). 3. The main method of reducing the error is the variance reduction; in other words, this is a good choice of the probability density p(x) of the random variable ξ in accordance with the physical and mathematical formulation of the problem. 4. The error is independent of the dimensionality of the problem. 5. A simple structure of the computation algorithm (the computations needed to realize a proper random variable are repeated N times). 6. The structure of the random variable ξ can be generally based on a physical model of the process that does not require a formulation of the controlling equations as in regular methods; this fact is increasingly important for modern problems. We illustrate the main features of the Monte Carlo methods and the conditions under which these methods outperform the conventional finite difference methods or are inferior to them using the following example. Suppose that we want to evaluate the definite integral of a continuous function over the interval [a, b]: To evaluate this integral using the Monte Carlo method, we construct a random variable with the probability density p(x) such that is equal to I. Now, if we set ξ = f(x)/p(x) within the integration limits, then we have, by the central limit theorem, On the one hand, the evaluation of I by above formula can be interpreted as the solution of mathematically stated problem; on the other hand, it can be interpreted as a direct simulation of the area under the function of f(x). The evaluation of the one-dimensional integral I 1 by the Monte Carlo method corresponds to the computation of I using the rectangular rule with the step Δx ≈ 1/ √N and an error O(Δx). If f(x) is sufficiently "good", the integral I 1 in the one-dimensional case can be calculated accurate to O(Δx 2 ) using the trapezoid rule, accurate to O(Δx 3 ) using the parabolic rule, and to any desired accuracy without a considerable increase in the computational effort. In the multidimensional case, the difficulties in using schemes of a high order of accuracy increase; for this reason, they are rarely used for the calculation of n-dimensional integrals I n for n ≥ 3.
Let us compare the efficiency of the regular and statistical methods for the problem described above. Let n be the dimensionality of the problem, Y be the number of nodes on an axis, R = Y n be the total number of nodes for the regular methods, q be the order of accuracy, N be the number of statistical trials, and ν be the number of operations needed to process one node (to perform one statistical trial). Then, ε L = Y -q is the error of the regular methods, ε K = N -1/2 is the error of the statistical methods, L(ε) = ν⋅R = ν⋅ε -n/q is the number of operations when the problem is solved by a regular method, and K(ε) = ν⋅N = ν⋅ε −2 is the number of operations when the problem is solved by the Monte Carlo method. Then, in the case of an equal number of operations needed to obtain a solution with the same accuracy using each of the methods, we have n = 2q.
Therefore, for n ≥ 3 and q = 1 (first-order schemes), the Monte Carlo methods are preferable. For other classes of problems, the relation between the efficiency of the methods can be different.

The Monte Carlo Methods in Computational Aerodynamics
The dynamics of rarefied gas is described by the Boltzmann integro-differential kinetic equation for the single-particle distribution density like as Here, f = f(t, x, y, z, ξ x , ξ y , ξ z ) is the distribution den sity. f, f 1 , f′, f 1 ′, correspond to the molecules with the velocities ξ, ξ 1 and ξ', ξ 1 ', before and after the collisio n, g is the relative velocity of the molecules in binary collisions g = ξ-ξ 1 = ξ'−ξ 1 ', and b and ε are t he impact parameter and the azimuth angle for the coll ision.
The complex nonlinear structure of the collision integral and the large number of variables (seven in the general case) present severe difficulties for the analysis including the numerical analysis. The high dimension, the probabilistic nature of the kinetic processes, and complex molecular interaction models are the natural prerequisites for the application of the Monte Carlo methods. Historically, the numerical statistical methods in rarefied gas dynamics developed in three directions: the use of the Monte Carlo methods to evaluate the collision integrals in the regular finite difference schemes for solving the kinetic equations; the direct statistical simulation of physical phenomena, which is subdivided into two approaches: the simulation of trajectories of test particles by the Haviland method [4] and the simulation of the evolution of the ensemble of particles by the Bird method [5]; the construction of a stochastic process using the Ulam-Neumann procedure [6] corresponding to the solution of the kinetic equation.
The hierarchy of levels of the description of large molecular systems includes a wide range of approaches, and various descriptions of the molecular dynamics at different levels can be used for constructing efficient statistical simulation methods.
The most detailed level of description is a dynamical system. To describe a system consisting of large number N Universal Journal of Physics and Application 2(4): 213-220, 2014 215 of particles (a molecular gas is a system of this kind), one must specify the initial coordinates and velocity of each molecule r j , x j and the evolution equations of this system The solution of such a system is an unrealizable (cannot be solved in practice) problem even for a very rarefied gas. Indeed, at a height of 400 km (the most popular satellite orbits), one cubic centimeter contains 10 9 molecules. For this reason, a less detailed statistical description is used.
Following the Gibbs formalism, rather than consider a single system, an ensemble of systems in the 6N-dimensional Γ-space distributed according to the N-particle distribution function is considered. This function is interpreted as the probability of finding the system in the neighborhood dr 1 …dr N dξ 1 …dξ N of the point r 1 ,…, r N , ξ 1 , …, ξ N at the moment t:

Such an ensemble is described by the Liouville equation
From now on, the Liouville equation and all the subsequent kinetic equations following from the Bogolyubov chain including the last Boltzmann equation have a probabilistic nature. Although Liouville Eq. is simpler than dynamical system, it takes into account the collisions of N molecules and is very difficult to analyze. A less detailed description is achieved by roughening the description using s-particle distribution functions , which determine the probability to simultaneously find s particles independently of the state of the remaining (N-s) particles.
Following Bogolyubov's ideas, we obtain the chain of linked equations up to the single-particle distribution function F 1 = f(t, r, ξ) corresponding to the Boltzmann gas, which only takes into account the binary collisions: Following Boltzmann, we assume that the molecules are spherically symmetric and accept the molecular chaos hypothesis F 2 (t, r 1 , r 2 , ξ 1 , ξ 2 ) = F 1 (t, r 1 , ξ 1 )F 1 (t, r 2 , ξ 2 ). It is very interesting to consider particular case of Liouville's equation and of Bogolyubov's chain that describe a spatially homogeneous gas consisting of bounded number of particles and corresponding to two-particle collisions; in this case, on the final link of the chain, we obtain the Kac master equation [7] where ϕ 1 and ϕ 2 are the one-and two-particle distribution functions. In contrast to the Boltzmann equation, Kac master equation is linear, which will be used in the development and justification of efficient numerical direct statistical simulation schemes.
Returning to the Boltzmann equation, we easily obtain all the macroscopic parameters from the definition of the function f. For example, the number of molecules n in a unit volume of the gas is The mean velocity of the molecules, the strain tensor, and the energy flux are determined by the relations where c = ξ -V is the thermal velocity of the molecules. The mean energy of the heat motion of molecules is usually described in terms of the temperature Applying the Chapman-Enskog procedure to the Boltzmann equation, we obtain the hydrodynamical level of description. This sequentially yields the Euler, Navier-Stokes, Barnett, etc., equations:

Development of Monte Carlo Methods in Hypersonic Aerodynamics
Following the general logic of the presentation, we may assume that the dynamics of continua, being a particular case of the kinetic treatment of the gas motion, has some statistical features; this fact will be used below.

Construction of Statistical Simulation Methods
The investigation and justification of Monte Carlo methods, it is impossible to do without the kinetic equation that describes the phenomenon being simulated. It is important in order to be sure that the solution is correct and enable to use the results of many typical problems were first solved by the direct simulation methods.
The model equations, in contrast to the linearized equation, are not rigorously derived from the Boltzmann equation; moreover, they are much more nonlinear than the original equation. However, they can prove to be simpler in practical implementation.
From the practical point of view, the direct statistical simulation methods based on the Bird and Haviland approaches are naturally the most efficient, and their modifications have dominated the computational aerodynamics.
Presently, the leading place in rarefied gas dynamics is occupied by the Bird method; various modifications of this method developed by Russian researchers (see, e.g., [8][9][10][11][12][13][14][15][16][17][18][19]) improved the efficiency of the original method by several orders of magnitude. The idea of the method is to split the evolution of the system in a small interval of time into two physical processes: 1) Relaxation in accordance with the collision operator in the kinetic equation

2) Free molecular transfer
This is the well-known first-order splitting scheme with respect to Δt for any operator equation. In this case, this approach is attractive because it splits the dynamics of a very complex kinetic system into two clear physical processes. The distribution function is modeled by N particles that first collide in each cell between themselves with a given frequency during the time Δt and then move at the distance ξ j Δt during the time Δt.
The central role in the nonstationary statistical simulation method is played by the procedure used to count the collisions. A pair of particles is chosen for collision in accordance with the collision frequency independently of the distance between them in the given cell. The velocities of the particles after the collision are chosen according to the molecular interaction laws. Although the efficiency of the method depends on many parameters of the computation scheme (relaxation, splitting with respect to time, stabilization, time step, space grid, and so on), the main studies devoted to the improvement of the method focus on the improvement of the collision procedure and on reducing the statistical error because this is the main instrument that makes it possible to reduce the number of particles in the cells and thus decrease the computation time and the requirements for computer memory. For example, a modification of the collision procedure was proposed in [16] for a particular case of Maxwell's molecules. With this modification, the computation results are almost independent of the number of particles in a cell in the range from 40 to 6. (With the ordinary computation scheme, the number of particles in a cell must be about 30). In [9][10][11][12][13][14], a general method that is independent of the kind of molecules was proposed; in that method, the subsystem of particles in each cell is considered as an N-particle Kac model at the stage of collisions.
The simulation of a collision is reduced to a statistical realization of the evolution of model during the time Δt rather than to the realization of the Boltzmann equation. The collision time in the Kac model is calculated in accordance with collision statistics in the ideal gas following the Bernoulli scheme. This scheme makes it possible to use a considerably smaller number of particles in a cell and a finer grid. The analysis shows that the computation results are almost independent of the number of particles in a cell down to two particles. The point is that the Boltzmann equation requires the molecular chaos assumption to be satisfied; however, for the number of particles in a cell that can be processed by modern computers, this assumption is satisfied only with a systematic error. By contrast, Kac does not rely on this assumption; therefore, the collision is calculated as a Markov process. On the other hand, as N → ∞, the Kac model is completely equivalent to the spatially homogeneous Boltzmann equation. Thus, the approach developed by Belotserkovskii and Yanitskii provides a basis for constructing efficient numerical schemes for solving three-dimensional aerodynamic flow problems, and, on the other hand, it solves the important methodological problem of the equivalence of the numerical method and the solution of the kinetic equation.
A huge number of studies are devoted to the methods of the traditional use of statistical simulation. For this reason, we mainly consider aerodynamics problems in this paper. It has already been mentioned that the statistical methods are more efficient in practical problems of rarefied gas dynamics than the regular and semiregular methods. For the flow problems, which are the most important problems in aerodynamics, the statistical methods were successfully used for the calculation of aerodynamic characteristics of various (including very complex) bodies in free molecular and almost free molecular flows. The procedure, which was developed more than twenty years ago, is now implemented in standard computer programs and is widely used in many organizations. Applications in the case of smaller Knudsen Journal of Physics and Application 2(4): 213-220, 2014 217 numbers Kn present considerable computational difficulties due to the reduced mean free path of the molecules and, respectively, a finer step with respect to time and space; in the case of the direct statistical simulation; the number of particles that simulate the distribution function is also increased.

The Direct Simulation Method for Inviscid Ideal Gas Flows
The relationship between the kinetic and the continuous medium models can be illustrated using the direct statistical simulation (DSS) method for inviscid flows [21,22]. The DSS methods for flows in a rarefied gas (Bird's approach) are based on the splitting of the evolution of a system of particles into two physical processes on the time interval Δt. These processes are described by the equations They represent the change in the density of the distribution caused by the spatially homogeneous relaxation and the collision-free transfer, respectively. It is assumed that the time interval Δt is sufficiently small for these processes to be considered as independent; Δt is of the order of magnitude of the mean free path, and the size of the cell is of the order of magnitude of the local mean free path.
When the deviation from the equilibrium state is small, the stage of the collision (relaxation) simulation can be replaced by a procedure of generating particle velocities in accordance with the given distribution.
Inviscid flows can be modeled assuming that the result of relaxation in a cell is the establishment of an equilibrium distribution. The algorithm used to simulate the flows of an inviscid compressible ideal gas by the DSS method based on the use of the locally equilibrium distribution function f 0 consists of the following steps.

Simulation of the transfer of the particles in time
Δt with the boundary conditions taken into account. 6. Calculation of the macroparameters of the particles in the cells after the transfer. 7. Repetition of Steps (2)-(6) up to a certain moment in time. The initial values for the macroparameters are chosen using the conditions of the nonstationary problem or using an initial approximation for the solution of the stationary problem (Step 1). The DSS method can be modified for inviscid flows so as to reduce the computational effort for the generation of particle velocities in each cell. Namely, once a set of velocities corresponding to the unit temperature and the zero mean velocity is generated, it can be then used for various mean velocities and temperatures in each cell. The procedure of the random distribution of the particles over the cell (Step 2) at each time step makes it possible to avoid the correlations that can appear when the same set of thermal velocities is used. Since the precalculated set of velocities is fixed and the actual number of particles in the cells is variable, weight coefficients are used to ensure the validity of the conservation laws.
The number of particles in a cell that is necessary for the simulation is determined in such a way that the statistical error satisfies the error consistency condition δ ~ O(Δx, Δt). Here, N is the number of particles in the cell.
There is no need for tracing each particle during the entire calculation process because the distribution density is assumed to be given. The particles are traced individually only in the transfer process in order to count the macroparameters in the cells at the next time step.
The expressions for the components of the thermal velocity (Step 3) can be obtained by simulating the normally distributed random variable Here, α k are independent random numbers that are uniformly distributed on the interval (0, 1). In order to reproduce the mean velocity more accurately, it is reasonable to use the following symmetrized algorithm: the thermal velocities of the particles with the odd indexes are calculated, and the thermal velocities of the particles with the even indexes are set equal to the velocities of the corresponding odd particles with the opposite sign.
Since the number of the modeled particles is bounded, Development of Monte Carlo Methods in Hypersonic Aerodynamics the temperature of the generated set of velocities is slightly different from the temperature used to generate it. For this reason, a correction is required (Step 4). Let us write the actual temperature of the generated set: Here, R = k B /m and N is the number of particles in a cell. Multiply the velocities by √T/√N to make the temperature equal to the given one. Then, all the conservativeness conditions are satisfied: The coordinates of the particles at the step of the collision-free transfer (Step 5) are modified by the rule At this step, the particles are assigned the internal energies εl to take into account the internal degrees of freedom. When the particles are transferred, the boundary conditions must be satisfied. Concerning the rigid walls, it is assumed that the reflection of the particles is specular. It can be shown that this corresponds to the impermeability condition in the Euler equations. The conditions on the exterior boundaries through which the flow enters the domain and leaves it must correspond to the flow regime. Additional constraints can be taken into account, such as the condition that the fluxes of the mass, momentum, and energy are constant (the independency of these quantities of statistical fluctuations).
After the transfer step, the new values of the macroparameters in the cells are calculated with account for all the particles that arrived from the neighboring cells (Step 6): Under the local equilibrium condition, the internal degrees of freedom are taken into account by adding the term ε i = vk B T t /2, where T' is the temperature in the cell from which the particle arrived and v is the number of the internal degrees of freedom. In such a system, the number of particles, the mean velocity, and the temperature are the working macroscopic variables.
The time step ∆t is chosen from the condition that the particles belonging to the main group (i.e., the particles whose velocity differs from the macroscopic velocity (u, v, w) not greater than by the characteristic thermal velocity a) must be transferred not more than by a single cell: Here, a is the characteristic thermal velocity; e.g., a = √2RT=c √2/γ. It follows that CFL≤ 1.
The scale of the mean free path is not used in this model; therefore, the space size of the cells is chosen in the same way as for the continuous media.
There is no clear stability condition for this method: if CFL > 1, then the particles can pass several cells, the algorithm remains the same, but its accuracy is lower. The condition imposed on ∆t is an accuracy rather than a stability condition. The DSS method applied to inviscid flows is efficient in the regions of the flow in which the deviation from the equilibrium state is small. It reduces the computational cost of calculating the particle velocities after the collisions.

Results
The DSMC method is commonly used to simulate hypersonic rarefied flow problems, and the accuracy of the method depends directly on the accuracy of the gas-surface interaction model [23][24][25]. The Maxwell model is the most widely used and is based on classical thermodynamics in which it is assumed that molecules will either reflect diffusely from a surface with complete energy accommodation or will reflect specularly with no change in energy. An accommodation coefficient α τ is defined which specifies the fraction of molecules that will be scattered diffusely, with α τ = 0 giving complete specular reflection, and α τ = 1 giving complete diffuse reflection. In the (Cercignani-Lampic-Lord, CLL) model, the transformations of the normal and tangential components of velocity are assumed to be manually independent. The Lennard-Jones potential function reflected the fact that the attractive forces dominated at high distances, and the repulsive forces dominated at small distances. The calculation was carried out for aerospace vehicle which described in [26,27] with taking into account angles of attack from -90° to 90°. Parameters of the problem as follows: velocity relationship V ∞ /√2RT ∞ = 20, the ratio of specific heat is 5/3, temperature factors t w = T w /T ∞ = 0.001, 0.1, 1 (T w -wall temperature of vehicle, T ∞ -temperature of flow field, respectively). The total particles were using 5×10 6 particles. The tangential accommodation coefficient α τ = 0.9, 1 and the energy accommodation coefficient α n = 1. In figure 1 presented the results of the aerodynamic characteristics of aerospace vehicle on various accommodation coefficients. The effects of temperature factors are described in fig. 2. We can see that with the increment of temperature factors the values are increased.
The results of the aerodynamic characteristics of aerospace vehicle by using Monte Carlo method with different gas-surface interaction models are described in fig  3. In figure shows the coefficient of aerodynamic force C x  with application of different models (Maxwell, Lennard-Jones, CLL). The coefficient C x increase as the angle of attack increase. The multiple reflections have not been taken into account, since they are not significant for this body at the variation of the angle of attack. It is clear that the coefficient is sensitive to different models of interaction of molecules with surface. The method used in this paper is good for calculation aerodynamic characteristics of reentry vehicles.

Conclusions
The development of Monte Carlo method in the computational aerodynamics of rarefied gases is given, and application of these methods in unconventional fields is described. A short history of these methods is presented, and their advantages are discussed. Possible directions of the development of the statistical simulation methods are discussed. Results of the calculation of the aerodynamic characteristics of the aerospace vehicle at various degrees of angle of attack with various gas-surface interaction models are presented.