Solution of stationary boundary value problems for the Boltzmann equation by the Monte Carlo method
S. V. Rogasinsky
DOI: https://doi.org/10.1515/mcma.1999.5.3.263
Abstract:A new version of the direct Monte Carlo method for solving boundary value problems for the Boltzmann equation is presented. In contrast to the conventional approach, we do not solve the problem via stabilization in time; when evaluating functionals of the solution to Boltzmann equation, the random trajectories are stopped with probability one after a finite number of transitions. 1. Stationary boundary value problem for the Boltzmann equation The main issue of the present study is the construction and justification of a new Monte Carlo method for numerical solution to boundary value problems for the nonlinear Boltzmann equation in the stationary formulation. This class of problems is of practical interest because many physical processes are governed by stationary boundary value problems [1,2]· Let Gr and G7 are bounded domains in 3D with piecewise smooth boundaries Γ and 7, respectively. We assume for simplicity that G7 C Gr . We denote G = Gr \ G7. Let nr = nr(r) and n7 = n7(r) are external normal vectors to Gr and to G7 , respectively. The problem is formulated as follows: find a nonnegative function /(r, v), continuous in G χ Λ, satisfying in G the equation and the boundary conditions: /(r,v) = /r(r,v) if (nrv) < 0, for r € Γ, 264 S. V. Rogasinsky whilst for r 6 7, it satisfies the integral relation: , v) = fcy(v' -> v; r) (n7v')/(r, v')d v' if (i^v) > 0, where L·/ ' ' \ η ' Ί k(v , Vl-> v, Vl) σ(|ν vj, . Here σ(|ν' — νί|,Ω) is the differential cross-section of scattering of two particles, Ω is the angle characterizing the relative velocity vectors of the particles after the scattering. It is assumed that the nonnegative functions /r(r, v) and fc7(v' — » v;r) are given, and they are positive on the surfaces Γ and 7 , respectively, and -> v; r) dv = l r G 7 . In addition, the function /(r, v) should satisfy the condition + v)/(r,v)civdr < oo for an integer α > 1, where the integration is taken over the whole domain of velocity and spatial variables. We deal here with the problem of construction of the numerical Monte Carlo method, and therefore we assume that there exists a unique solution to the formulated problem. It will be convenient to give an equivalent formulation of the problem. To this end, we include the boundary conditions on 7 into the Boltzmann equation. Let us introduce the notation , ν! -> v, nv ην l ; "\ 0 , i f (nv)>0; (U > "\ 0 , if (nv) < 0; ir(r) is the generalized function (a simple layer) whose support is concentrated on the surface Γ , £7(r) is the simple layer on 7 , and g(r, v) = {-(nrv)-}/r(r,v)ir(r). Then the boundary value problem for the Boltzmann equation is reformulated as follows: find a nonnegative function /(r, v), continuous in G x jR and satisfying in r € G the equation Solution of stationary boundary value problems for the Boltzmann equation 265 ;/(r, v) = Sf [/, /] + ^(v' -> vlrH-Kv ΓΚΟΟ/ίΓ, v )dv , (1) and boundary conditions (nrv)-/(r,v) = (nrv)-/r(r,v) , r <E Γ, (2) (ηγν)/(Γ,ν) = 0. r G 7(3) The boundary conditions (2) on the surface Γ make it impossible to use the JV-particle Kolmogorov equation with a fixed number of particles to construct the direct Monte Carlo method, as in the homogeneous case [3]; in our case the number of particles varies and this should be properly taking into account. It is possible to do this approximately, by solving nonstationary boundary value problems and using the conventional splitting (over physical processes) technique and the stabilization method [1,4]. To this end, one solves the inhomogeneous η-particle equation with a constant number of particles in each time step Δ< of the splitting process; the boundary conditions are taken into account in the stage of spatial movement of particles [4]. This approach assumes that the stabilization method is applicable to solve the stationary boundary value problem for the Boltzmann equation (l)-(3). In this paper we do not use this conventional approach. 2. An auxiliary system of JV-particle equations In [6], we suggested a system of JV-particle equations which properly takes account of the change of particles caused by a flux of particles into the domain G, without affecting the probabilistic character of these equations. This is a crucial point in the use of the JV-particle equations for solving the nonlinear Boltzmann equations [5]. According to the approach presented in [3], we formulate a system of ΛΓ-particle equations which is a basis for constructing the direct Monte Carlo method for solving the problem (l)-(3). Let R„ = (ΓΙ, ..., rN), VN = (YI, ..., VN). We say that R„ lies in G C Λ (and write RN 6 G ) if rt· e G for all i = 1,..., N. Let us introduce the notation: Rl = { v : v e J?, v* < |v| < v* } is the domain of velocity variables. The constants v* and v* satisfy the condition 0 < v* < v*. Analogously we say that VN lies in JR* (and write VN e izj) if v; 6 R* for alH = 1,..., N. 266 S. V. Rogasinsky Denote by R^ the spatial coordinates of a system of N particles if there is at least one particle (indexed through i, l < i < N ,) such that rt e Γ . In this case we say that RN belongs to the boundary Γ. By RJ we denote the spatial coordiantes of a system of Ν particles if there is at least one particle (indexed through i, 1 < i < Ν ,) such that FJ G 7 . In this case we say that RN belongs to the boundarty 7. Let us define the indicator of the event that R^ belongs to the boundary 7 as follows: We define K(V'N -> VN|RN) as a function which determines a pair interaction in a system of Ν particles in the domain G. It satisfies the conditions K{V'N -> V„|R„) = K(V„ -> V'N\RN) , ->· VW|R^) = o, /f (v; -4 V^IRJ;) = o, if(vi -> V^RO Ξ ο . Let A(RW> Vw) = / tf (V„ -4 V^|Rw)dV; . (5) The function KC(V'N ->· V^IR^) is defined as a probability density function of velocities in the ΛΓ-particle system after an interaction, being in a state Rw , V^. This function has the same support as the function K(VN -> V„|R„) and they are related by K(V'H -> VN\RN] = l/7(RW) A(R„,V'K) KC(V'N -4 VW|RW) . The function ΚΊ(ν'Ν -4 V^IR^) describes the interaction of a system of Ν particles with the boundary 7 . It satisfies the conditions = 0, if Rjv^R^;