IMplicit-EXplicit schemes for BGK kinetic equations Sandra Pieraccini, Gabriella Puppo Politecnico di Torino Dipartimento di Matematica HYP06, Lyon, July 17-21 – p.1/22
Boltzmann equation Rarefied gas flow obeys Boltzmann equation which describes the evolution of the probability density f = f ( x, v, t ) of finding a particle at position x with microscopic velocity v at time t : ∂ t f + v · ∇ x f = Q ( f, f ) where Q ( f, f ) is the collision integral which accounts for particle interactions. HYP06, Lyon, July 17-21 – p.2/22
Collision integral Particles interact through collisions in which a particle located at the point ( x, v, t ) in phase space collides with a particle at ( x, v ∗ , t ) . The collision integral accounts for all possible interactions: � � ( f ( v ′ ) f ( v ′ Q ( f, f ) = σ ∗ ) − f ( v ) f ( v ∗ )) | ( v − v ∗ ) · n | dS dv ∗ ℜ 3 B − During collisions, particles exchange momentum, so their location immediately after impact will be given by ( x, v ′ , t ) and ( x, v ′ ∗ , t ) . The relation between v, v ∗ and v ′ , v ′ ∗ contains the physics of the interaction between particles. HYP06, Lyon, July 17-21 – p.3/22
Adding chemistry... The model described so far governs the motion of a monoatomic gas. To include more phenomena, the model becomes more complicated: For a polyatomic gas, the internal energy of the molecule must also be added as an independent variable. Modelling chemical reactions requires to consider distribution functions for all species involved and their reciprocal interactions in suitable collision integrals. Clearly the equations become more and more complicated, and their numerical solution extremely expensive. For this reason, simplified models for Boltzmann equation are introduced. HYP06, Lyon, July 17-21 – p.4/22
Motivation for BGK model The BGK model (Bhatnagar-Gross-Krook ’54) approximates Boltzmann equation for the evolution of a rarefied gas for small and moderate Knudsen numbers: mean free path Kn = characteristic length of the problem Lately, interest in this model has increased because: Several desirable properties have been shown to hold for the BGK model and its variants, such as BGK-ES, (Perthame et al. from 1989 on) HYP06, Lyon, July 17-21 – p.5/22
Motivation for BGK model The BGK model (Bhatnagar-Gross-Krook ’54) approximates Boltzmann equation for the evolution of a rarefied gas for small and moderate Knudsen numbers: mean free path Kn = characteristic length of the problem Lately, interest in this model has increased because: The BGK model has been extended to include more general fluids and can now be applied to the flow of a polytropic gas (Mieussens) and to mixtures of reacting gases (Aoki et al.) HYP06, Lyon, July 17-21 – p.5/22
Motivation for BGK model The BGK model (Bhatnagar-Gross-Krook ’54) approximates Boltzmann equation for the evolution of a rarefied gas for small and moderate Knudsen numbers: mean free path Kn = characteristic length of the problem Lately, interest in this model has increased because: New applications of kinetic models have appeared. For instance, fluid flow in nanostructures is described by Boltzmann equation, since it involves also moderate Knudsen numbers HYP06, Lyon, July 17-21 – p.5/22
Numerical schemes for the BGK model The development of numerical methods for the BGK model has started only recently. Yang, Huang ’95 This scheme is high order accurate in space, but only first order accurate in time HYP06, Lyon, July 17-21 – p.6/22
Numerical schemes for the BGK model The development of numerical methods for the BGK model has started only recently. Aoki, Kanba, Takata ’97 This is a second order scheme, designed for smooth solutions HYP06, Lyon, July 17-21 – p.6/22
Numerical schemes for the BGK model The development of numerical methods for the BGK model has started only recently. Mieussens, ’00 A second order scheme, where conservation is exactly enforced HYP06, Lyon, July 17-21 – p.6/22
Numerical schemes for the BGK model The development of numerical methods for the BGK model has started only recently. Andries, Bourgat, le Tallec, Perthame ’02 A stochastic Monte Carlo scheme HYP06, Lyon, July 17-21 – p.6/22
Numerical schemes for the BGK model The development of numerical methods for the BGK model has started only recently. Pieraccini, Puppo ’06 A non oscillatory high order scheme in space and time HYP06, Lyon, July 17-21 – p.6/22
BGK model The main variable is the probability density f that a particle be in the point x ∈ R d with velocity v ∈ R N at time t , thus f = f ( x, v, t ) . The evolution of f is given by: ∂f ∂t + v · ▽ x f = 1 τ ( f M − f ) , with initial condition f ( x, v, 0) = f 0 ( x, v ) ≥ 0 . Here τ is the collision time τ ≃ Kn , with τ > 0 and close to the hydrodynamic regime τ can be very small. HYP06, Lyon, July 17-21 – p.7/22
The Maxwellian f M is the local Maxwellian function, and it is built starting from the macroscopic moments of f : −|| v − u ( x, t ) || 2 � � ρ ( x, t ) f M ( x, v, t ) = (2 πRT ( x, t )) N/ 2 exp , 2 RT ( x, t ) where ρ and u are the gas macroscopic density and velocity and T is the temperature. They are computed from f as: ρ 1 � � � = f where � g � = R N g dv ρu v v 2 E E is total energy, and the temperature is: NT/ 2 = E − 1 / 2 ρu 2 . HYP06, Lyon, July 17-21 – p.8/22
The Maxwellian Thus the BGK equation ∂f ∂t + v · ▽ x f = 1 τ ( f M − f ) , describes the relaxation of f towards the local equilibrium Maxwellian f M . The local equilibrium is reached with a speed that is inversely proportional to τ . Thus the system is stiff for τ << 1 . HYP06, Lyon, July 17-21 – p.8/22
The Maxwellian Note that the Maxwellian f M has the same moments of f , namely: ρ 1 1 � � � � = f = f M ρu v v v 2 v 2 E HYP06, Lyon, July 17-21 – p.8/22
Conservation As in Boltzmann equation, the macroscopic moments of f are conserved: ∂ t � f � + ∇ x · � fv � = 0 , ∂ t � fv � + ∇ x · � v ⊗ vf � = 0 , � 1 � 1 2 � v � 2 f � 2 � v � 2 vf � ∂ t + ∇ x · = 0 . Moreover, for τ → 0 the macroscopic solution converges to the gas dynamic solution of Euler equations. HYP06, Lyon, July 17-21 – p.9/22
Conservation As in Boltzmann equation, the macroscopic moments of f are conserved: ∂ t � f � + ∇ x · � fv � = 0 , ∂ t � fv � + ∇ x · � v ⊗ vf � = 0 , � 1 � 1 2 � v � 2 f � 2 � v � 2 vf � ∂ t + ∇ x · = 0 . Thus a numerical scheme for the BGK model must be con- servative, and its numerical solution must converge to the Eu- ler solution as τ → 0 . HYP06, Lyon, July 17-21 – p.9/22
Entropy principle The BGK model satisfies an entropy principle: ∂ t � f log f � + ∇ x � vf log f � ≤ 0 , ∀ f ≥ 0 where equality holds if and only if f = f M . Thus the Maxwellian f M is the equilibrium solution of the system. The macroscopic entropy is: S = � f log f � Note that as τ → 0 , entropy is conserved on smooth solutions, as for Euler solutions. HYP06, Lyon, July 17-21 – p.10/22
Splitting schemes To circumvent the stiffness of the source term, the BGK equation can be integrated through operator splitting. HYP06, Lyon, July 17-21 – p.11/22
Splitting schemes First, the convective part of the equations is integrated explicitly in time: ∂f ∂t + v · ▽ x f = 0 In this fashion, space cells remain decoupled notwithstanding the term v · ▽ x . HYP06, Lyon, July 17-21 – p.11/22
Splitting schemes First, the convective part of the equations is integrated explicitly in time: ∂f ∂t + v · ▽ x f = 0 In this fashion, space cells remain decoupled notwithstanding the term v · ▽ x . Next, the relaxation part of the equations is integrated implicitly in time, to circumvent stability restrictions on the time step induced by the stiff term: ∂f ∂t = 1 τ ( f M − f ) Note however that this procedure is only first order accurate in time HYP06, Lyon, July 17-21 – p.11/22
IMEX RK schemes To achieve the advantages of operator splitting, coupled with high order accuracy in time, we use IMplicit EXplicit (IMEX) Runge -Kutta schemes. To set notation, we consider the following autonomous ODE: y ′ = f ( y ) + 1 εg ( y ) Let y n denote the numerical solution at time t n = n ∆ t . IMEX schemes are composed of an explicit and an implicit Runge- Kutta schemes, coupled together. HYP06, Lyon, July 17-21 – p.12/22
IMEX RK schemes The updated solution is computed as: ν ν b i f ( y ( i ) ) + ∆ t y n +1 = y n + ∆ t ˜ � � b i g ( y ( i ) ) ε i =1 i =1 where the stage values y ( i ) are given by: i − 1 i a il f ( y ( l ) ) + ∆ t y ( i ) = y n + ∆ t � � a il g ( y ( l ) ) ˜ ε l =1 l =1 HYP06, Lyon, July 17-21 – p.12/22
Recommend
More recommend
