A discontinuous Galerkin fast spectral method for the multi-species Boltzmann equation

TL;DR

This paper presents a fast spectral method combined with discontinuous Galerkin discretization for solving the multi-species Boltzmann equation, achieving high accuracy and efficiency in simulating gas mixtures with various collision kernels.

Contribution

It introduces a novel fast Fourier spectral method for the multi-species Boltzmann collision operator, coupled with DG discretization, improving computational efficiency and applicability to complex collision kernels.

Findings

Achieves spectral accuracy and reduced computational complexity.

Demonstrates efficiency and accuracy through numerical benchmarks.

Validates the method against DSMC simulations across various scenarios.

Abstract

We introduce a fast Fourier spectral method for the multi-species Boltzmann collision operator. The method retains the riveting properties of the single-species fast spectral method (Gamba et al. SIAM J. Sci. Comput., 39 pp. B658--B674 2017) including: (a) spectral accuracy, (b) reduced computational complexity compared to direct spectral method, (c) reduced memory requirement in the precomputation, and (d) applicability to general collision kernels. The fast collision algorithm is then coupled with discontinuous Galerkin discretization in the physical space (Jaiswal et al. J. Comp. Phys., 378 pp. 178--208 2019) to result in a highly accurate deterministic method (DGFS) for the full Boltzmann equation of gas mixtures. A series of numerical tests is performed to illustrate the efficiency and accuracy of the proposed method. Various benchmarks highlighting different collision kernels, different mass ratios, momentum transfer, heat transfer, and in particular the diffusive transport have been studied. The results are directly compared with the direct simulation Monte Carlo (DSMC) method.

Figures

Tables

Table 1. Table 1 : Spatially homogeneous Krook-Wu solution. L ∞ superscript 𝐿 L^{\infty} error ℰ ( i ) = ‖ ∂ t f e ​ x ​ a ​ c ​ t ( i ) − ∂ t f n ​ u ​ m ​ e ​ r ​ i ​ c ​ a ​ l ( i ) ‖ L ∞ , i = { 1 , 2 } formulae-sequence superscript ℰ 𝑖 subscript norm subscript 𝑡 superscript subscript 𝑓 𝑒 𝑥 𝑎 𝑐 𝑡 𝑖 subscript 𝑡 superscript subscript 𝑓 𝑛 𝑢 𝑚 𝑒 𝑟 𝑖 𝑐 𝑎 𝑙 𝑖 superscript 𝐿 𝑖 1 2 \mathcal{E}^{(i)}=\|\partial_{t}f_{exact}^{(i)}-\partial_{t}f_{numerical}^{(i)}\|_{L^{\infty}},\;i=\{1,2\} at t = 4 𝑡 4 t=4 for different mass ratios. N 𝑁 N is the number of discretization points in each velocity dimension and N ρ subscript 𝑁 𝜌 N_{\rho} is the number of Gauss quadrature points used in the radial direction. Number of quadrature points M 𝑀 M used on the sphere is fixed to 6 6 6 . A fixed velocity domain [ − 12 ,  12 ] 3 superscript 12 12 3 [-12,\;12]^{3} has been used for all cases.

$N$	$N_{\rho }$	$m_1/m_2=1$		$m_1/m_2=2$		$m_1/m_2=4$		$m_1/m_2=6$		$m_1/m_2=8$
		${ℰ}^{\left(1\right)}$	${ℰ}^{\left(2\right)}$	${ℰ}^{\left(1\right)}$	${ℰ}^{\left(2\right)}$	${ℰ}^{\left(1\right)}$	${ℰ}^{\left(2\right)}$	${ℰ}^{\left(1\right)}$	${ℰ}^{\left(2\right)}$	${ℰ}^{\left(1\right)}$	${ℰ}^{\left(2\right)}$
16	4	1.528e-03	1.528e-03	6.675e-03	4.444e-03	2.048e-03	3.633e-03	1.414e-03	4.941e-04	2.709e-04	2.481e-03
	8	2.114e-03	2.114e-03	7.795e-03	4.917e-03	2.253e-03	3.828e-03	1.425e-03	7.087e-04	2.715e-04	1.963e-03
	16	2.114e-03	2.114e-03	7.795e-03	4.917e-03	2.253e-03	3.828e-03	1.425e-03	7.086e-04	2.715e-04	1.963e-03
32	8	1.526e-04	1.526e-04	1.671e-03	2.237e-04	7.249e-03	9.770e-04	1.692e-02	1.701e-03	3.029e-02	2.928e-03
	16	1.873e-04	1.873e-04	1.729e-03	1.852e-04	8.018e-03	9.020e-04	1.935e-02	1.652e-03	3.338e-02	2.822e-03
	32	1.873e-04	1.873e-04	1.729e-03	1.852e-04	8.018e-03	9.020e-04	1.935e-02	1.652e-03	3.338e-02	2.822e-03
64	16	4.227e-08	4.227e-08	4.704e-06	4.749e-08	4.201e-04	3.460e-06	5.263e-03	4.056e-05	6.231e-02	4.053e-04
	32	4.227e-08	4.227e-08	4.754e-06	4.422e-08	4.043e-04	9.441e-07	5.153e-03	3.381e-05	6.186e-02	4.082e-04
	64	4.227e-08	4.227e-08	4.754e-06	4.422e-08	4.043e-04	9.441e-07	5.153e-03	3.381e-05	6.186e-02	4.082e-04

Table 2. Table 2 : Numerical parameters for normal shock wave [ 39 ] .

Parameter	Case NS-01	Case NS-02	Case NS-03	Case NS-04
Molecular mass: $m_1$ ($\times {10}^{27}kg$)	$6.63$	$6.63$	$6.63$	$6.63$
Molecular mass: $m_2$ ($\times {10}^{27}kg$)	$3.315$	$1.6575$	$3.315$	$3.315$
Mass Ratio: $m_2/m_1$	0.5	0.25	0.5	0.5
Mach number	1.5	1.5	1.5	3.0
Concentration: $n_{-}^{\left(2\right)}/n_{-}=n_{+}^{\left(2\right)}/n_{+}$	0.5	0.5	0.1	0.1
Non-dim physical space	$[-0.5,\mathrm{\hspace{0.17em}0.5}]$	$[-0.5,\mathrm{\hspace{0.17em}0.5}]$	$[-0.5,\mathrm{\hspace{0.17em}0.5}]$	$[-0.5,\mathrm{\hspace{0.17em}0.5}]$
Non-dim velocity space	${[-9,\mathrm{\hspace{0.17em}9}]}^3$	${[-15,\mathrm{\hspace{0.17em}15}]}^3$	${[-9,\mathrm{\hspace{0.17em}9}]}^3$	${[-15,\mathrm{\hspace{0.17em}15}]}^3$
$N^3$	${32}^3$	${64}^3$	${32}^3$	${48}^3$
$N_{\rho }$	$32$	$16$	$32$	$48$
$M$	12	12	12	12
Spatial elements	$8,\mathrm{\hspace{0.17em}16}$	16	16	16
DG order	3	3	3	3
Time step ($s\times {10}^8$)	$5.57,\mathrm{\hspace{0.17em}2.77}$	$1.64$	$2.77$	$1.64$
Viscosity index: ${\omega }_{ij}$	0.5	0.5	0.5	0.5
Scattering parameter: ${\alpha }_{ij}$	1	1	1	1
Ref. diameter: $d_{\text{ref},ij}$ ($\times {10}^{10}m$)	$2.17$	$2.17$	$2.17$	$2.17$
Ref. temperature: $T_{\text{ref},ij}$ ($K$)	273	273	273	273
Characteristic mass: $m_0$ ($\times {10}^{27}kg$)	$6.63$	$6.63$	$6.63$	$6.63$
Characteristic length: $H_0$ ($mm$)	30	30	30	30
Characteristic velocity: $u_0$ ($m/s$)	963.7	963.7	963.7	963.7
Characteristic temperature: $T_0$ ($K$)	223	223	223	223
Characteristic number density: $n_0$ ($m^{-3}$)	$2.889\times {10}^{21}$	$2.889\times {10}^{21}$	$2.889\times {10}^{21}$	$2.889\times {10}^{21}$
Upstream conditions (subscript -)
Velocity: $u_{-}$ ($m/s$)	1523.737	1669.171	1353.876	2707.753
Temperature: $T_{-}$ ($K$)	223	223	223	223
Mean free path: ${\lambda }_{-}={\left(\sqrt{2}\pi \left(n_{-}^{\left(1\right)}+n_{-}^{\left(2\right)}\right)d_{\text{ref},ij}^2\right)}^{-1}$ ($m$)	0.000827	0.000827	0.00148	0.00148
Number density: $n_{-}^{\left(1\right)}$ ($m^{-3}$)	$2.889\times {10}^{21}$	$2.889\times {10}^{21}$	$2.889\times {10}^{21}$	$2.889\times {10}^{21}$
Number density: $n_{-}^{\left(2\right)}$ ($m^{-3}$)	$2.889\times {10}^{21}$	$2.889\times {10}^{21}$	$3.209\times {10}^{20}$	$3.209\times {10}^{20}$
Downstream conditions (subscript +)
Velocity: $u_{+}$ ($m/s$)	888.847	973.683	789.761	902.584
Temperature: $T_{+}$ ($K$)	333.338	333.338	333.338	817.667
Number density: $n_{+}^{\left(1\right)}$ ($m^{-3}$)	$4.953\times {10}^{21}$	$4.953\times {10}^{21}$	$4.953\times {10}^{21}$	$8.669\times {10}^{21}$
Number density: $n_{+}^{\left(2\right)}$ ($m^{-3}$)	$4.953\times {10}^{21}$	$4.953\times {10}^{21}$	$5.502\times {10}^{20}$	$9.633\times {10}^{20}$
Initial conditions
Velocity: $u$ ($m/s$)	$u_{-}+\left(u_{+}-u_{-}\right)x/H_0$
Temperature: $T$ ($K$)	$T_{-}+\left(T_{+}-T_{-}\right)x/H_0$
Number density: $n^{\left(1\right)}$ ($m^{-3}$)	$n_{-}^{\left(1\right)}+\left(n_{+}^{\left(1\right)}-n_{-}^{\left(1\right)}\right)x/H_0$
Number density: $n^{\left(2\right)}$ ($m^{-3}$)	$n_{-}^{\left(2\right)}+\left(n_{+}^{\left(2\right)}-n_{-}^{\left(2\right)}\right)x/H_0$

Table 3. Table 3 : VHS and VSS model parameters for different mixture systems [ 7 ] .

Mixture	Ar-Kr	Ar-Kr
Collision kernel	VHS	VSS
Molecular mass: $m_1$ ($\times {10}^{27}kg$)	$66.3$	$66.3$
Molecular mass: $m_2$ ($\times {10}^{27}kg$)	$139.1$	$139.1$
Reference viscosity: ${\mu }_{\text{ref},1}$ ($\times {10}^5\text{Pa}\cdot s$)	$2.117$	$2.117$
Reference viscosity: ${\mu }_{\text{ref},2}$ ($\times {10}^5\text{Pa}\cdot s$)	$2.328$	$2.328$
Viscosity index: $({\omega }_{11},{\omega }_{22})$	$(0.81,\mathrm{\hspace{0.17em}0.8})$	$(0.81,\mathrm{\hspace{0.17em}0.8})$
Viscosity index: $({\omega }_{12},{\omega }_{21})$	$(0.805,\mathrm{\hspace{0.17em}0.805})$	$(0.805,\mathrm{\hspace{0.17em}0.805})$
Scattering parameter: $({\alpha }_{11},{\alpha }_{22})$	$(1,\mathrm{\hspace{0.17em}1})$	$(1.4,\mathrm{\hspace{0.17em}1.32})$
Scattering parameter: $({\alpha }_{12},{\alpha }_{21})$	$(1,\mathrm{\hspace{0.17em}1})$	$(1.36,\mathrm{\hspace{0.17em}1.36})$
Ref. diameter: $(d_{\text{ref},11},d_{\text{ref},22})$ ($\times {10}^{10}m$)	$(4.17,\mathrm{\hspace{0.17em}4.76})$	$(4.11,\mathrm{\hspace{0.17em}4.7})$
Ref. diameter: $(d_{\text{ref},11},d_{\text{ref},22})$ ($\times {10}^{10}m$)	$(4.465,\mathrm{\hspace{0.17em}4.465})$	$(4.405,\mathrm{\hspace{0.17em}4.405})$
Ref. temperature: $(T_{\text{ref},11},T_{\text{ref},22})$ ($K$)	$(273,\mathrm{\hspace{0.17em}273})$	$(273,\mathrm{\hspace{0.17em}273})$
Ref. temperature: $(T_{\text{ref},12},T_{\text{ref},21})$ ($K$)	$(273,\mathrm{\hspace{0.17em}273})$	$(273,\mathrm{\hspace{0.17em}273})$

Table 4. Table 4 : Numerical parameters for Fourier heat transfer. The molecular collision parameters for Ar-Kr system are provided in Table  3 .

Parameter	Case F-01	Case F-02	Case F-03	Case F-04	Case F-05	Case F-06
Mixture	Ar-Kr	Ar-Kr	Ar-Kr	Ar-Kr	Ar-Kr	Ar-Kr
Collision kernel	VHS	VHS	VHS	VHS	VHS	VHS
Non-dim physical space	$[0,\mathrm{\hspace{0.17em}1}]$	$[0,\mathrm{\hspace{0.17em}1}]$	$[0,\mathrm{\hspace{0.17em}1}]$	$[0,\mathrm{\hspace{0.17em}1}]$	$[0,\mathrm{\hspace{0.17em}1}]$	$[0,\mathrm{\hspace{0.17em}1}]$
Non-dim velocity space	${[-5,\mathrm{\hspace{0.17em}5}]}^3$	${[-5,\mathrm{\hspace{0.17em}5}]}^3$	${[-5,\mathrm{\hspace{0.17em}5}]}^3$	${[-9,\mathrm{\hspace{0.17em}9}]}^3$	${[-9,\mathrm{\hspace{0.17em}9}]}^3$	${[-9,\mathrm{\hspace{0.17em}9}]}^3$
$N^3$	${32}^3$	${32}^3$	${32}^3$	${64}^3$	${64}^3$	${64}^3$
$N_{\rho }$	$32$	$32$	$32$	$64$	$64$	$64$
$M$	12	12	12	12	12	12
Spatial elements	4	4	4	4	4	4
DG order	3	3	3	3	3	3
Time step (s)	$2\times {10}^{-8}$	$2\times {10}^{-8}$	$2\times {10}^{-8}$	$2\times {10}^{-8}$	$2\times {10}^{-8}$	$2\times {10}^{-8}$
Mass: $m_0$	$m_{\text{Ar}}=m_1$	$m_{\text{Ar}}=m_1$	$m_{\text{Ar}}=m_1$	$m_{\text{Ar}}=m_1$	$m_{\text{Ar}}=m_1$	$m_{\text{Ar}}=m_1$
Length: $H_0$ ($mm$)	1	1	1	1	1	1
Velocity: $u_0$ ($m/s$)	337.2	337.2	337.2	337.2	337.2	337.2
Temperature: $T_0$ ($K$)	273	273	273	273	273	273
Number density: $n_0$ ($m^{-3}$)	$1.680\times {10}^{21}$	$8.401\times {10}^{20}$	$1.680\times {10}^{20}$	$1.680\times {10}^{21}$	$8.401\times {10}^{20}$	$1.680\times {10}^{20}$
Left wall (purely diffuse) boundary conditions (subscript $l$)
Velocity: $u_l$ ($m/s$)	$0$	$0$	$0$	$0$	$0$	$0$
Temperature: $T_l$ ($K$)	263	263	263	223	223	223
Right wall (purely diffuse) boundary conditions (subscript $r$)
Velocity: $u_r$ ($m/s$)	$0$	$0$	$0$	$0$	$0$	$0$
Temperature: $T_r$ ($K$)	283	283	283	323	323	323
Initial conditions
Velocity: $u$ ($m/s$)	0	0	0	0	0	0
Temperature: $T$ ($K$)	273	273	273	273	273	273
Number density: $n^{\left(1\right)}$ ($m^{-3}$)	$1.680\times {10}^{21}$	$8.401\times {10}^{20}$	$1.680\times {10}^{20}$	$1.680\times {10}^{21}$	$8.401\times {10}^{20}$	$1.680\times {10}^{20}$
Number density: $n^{\left(2\right)}$ ($m^{-3}$)	$8.009\times {10}^{20}$	$4.004\times {10}^{20}$	$8.009\times {10}^{19}$	$8.009\times {10}^{20}$	$4.004\times {10}^{20}$	$8.009\times {10}^{19}$
Knudsen: $({\mathrm{Kn}}_{11},{\mathrm{Kn}}_{22})$	$(0.770,\mathrm{\hspace{0.17em}0.591})$	$(1.541,\mathrm{\hspace{0.17em}1.182})$	$(7.703,\mathrm{\hspace{0.17em}5.912})$	$(0.770,\mathrm{\hspace{0.17em}0.591})$	$(1.541,\mathrm{\hspace{0.17em}1.182})$	$(7.703,\mathrm{\hspace{0.17em}5.912})$
Knudsen: $({\mathrm{Kn}}_{12},{\mathrm{Kn}}_{21})$	$(0.782,\mathrm{\hspace{0.17em}0.540})$	$(1.564,\mathrm{\hspace{0.17em}1.080})$	$(7.820,\mathrm{\hspace{0.17em}5.399})$	$(0.782,\mathrm{\hspace{0.17em}0.540})$	$(1.564,\mathrm{\hspace{0.17em}1.080})$	$(7.820,\mathrm{\hspace{0.17em}5.399})$

Table 5. Table 5 : Numerical parameters for oscillatory Couette flow. The molecular collision parameters for Ar-Kr system are provided in Table  3 .

Parameter	Case OC-01	Case OC-02
Mixture	Ar-Kr	Ar-Kr
Collision kernel	VHS	VHS
Non-dim physical space	$[0,\mathrm{\hspace{0.17em}1}]$	$[0,\mathrm{\hspace{0.17em}1}]$
Non-dim velocity space	${[-5,\mathrm{\hspace{0.17em}5}]}^3$	${[-9,\mathrm{\hspace{0.17em}9}]}^3$
$N^3$	${24}^3$	${48}^3$
$N_{\rho }$	$24$	$48$
$M$	6	6
Spatial elements	4	4
DG order	3	3
Time step (s)	$2\times {10}^{-8}$	$2\times {10}^{-8}$
Characteristic mass: $m_0$	$m_{\text{Ar}}=m_1$	$m_{\text{Ar}}=m_1$
Characteristic length: $H_0$ ($mm$)	1	1
Characteristic velocity: $u_0$ ($m/s$)	337.2	337.2
Characteristic temperature: $T_0$ ($K$)	273	273
Characteristic number density: $n_0$ ($m^{-3}$)	$8.401\times {10}^{20}$	$8.401\times {10}^{20}$
Initial conditions
Velocity: $u$ ($m/s$)	0	0
Temperature: $T$ ($K$)	273	273
Number density: $n^{\left(1\right)}$ ($m^{-3}$)	$8.401\times {10}^{20}$	$8.401\times {10}^{20}$
Number density: $n^{\left(2\right)}$ ($m^{-3}$)	$4.004\times {10}^{20}$	$4.004\times {10}^{20}$
Knudsen number: $({\mathrm{Kn}}_{11},{\mathrm{Kn}}_{22})$	$(1.541,\mathrm{\hspace{0.17em}1.182})$	$(1.541,\mathrm{\hspace{0.17em}1.182})$
Knudsen number: $({\mathrm{Kn}}_{12},{\mathrm{Kn}}_{21})$	$(1.564,\mathrm{\hspace{0.17em}1.080})$	$(1.564,\mathrm{\hspace{0.17em}1.080})$
Left wall (purely diffuse) boundary conditions (subscript $l$)
Velocity: $u_l$ ($m/s$)	$(0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0})$	$(0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0})$
Temperature: $T_l$ ($K$)	273	273
Right wall (purely diffuse) boundary conditions (subscript $r$)
Velocity: $u_r$ ($m/s$)	$(0,\mathrm{\hspace{0.17em}50}sin\left(\zeta t\right),\mathrm{\hspace{0.17em}0})$	$(0,\mathrm{\hspace{0.17em}500}sin\left(\zeta t\right),\mathrm{\hspace{0.17em}0})$
Temperature: $T_r$ ($K$)	273	273
Period of oscillation: $\zeta$ ($s^{-1}$)	$2\pi /\left(5\times {10}^{-5}\right)$	$2\pi /\left(5\times {10}^{-5}\right)$
Velocity amplitude: $v_a$ ($m/s$)	$50$	$500$

Table 6. Table 6 : Numerical parameters for Couette flow. The molecular collision parameters for Ar-Kr system are provided in Table  3 .

Parameter	Case C-01
Mixture	Ar-Kr
Collision kernel	VSS
Non-dim physical space	$[0,\mathrm{\hspace{0.17em}1}]$
Non-dim velocity space	${[-7,\mathrm{\hspace{0.17em}7}]}^3$
$N^3$	${32}^3$
$N_{\rho }$	$32$
$M$	12
Spatial elements	4
DG order	3
Time step (s)	$2\times {10}^{-8}$
Characteristic mass: $m_0$	$m_{\text{Ar}}=m_1$
Characteristic length: $H_0$ ($mm$)	1
Characteristic velocity: $u_0$ ($m/s$)	337.2
Characteristic temperature: $T_0$ ($K$)	273
Characteristic number density: $n_0$ ($m^{-3}$)	$1.680\times {10}^{21}$
Initial conditions
Velocity: $u$ ($m/s$)	0
Temperature: $T$ ($K$)	273
Number density: $n^{\left(1\right)}$ ($m^{-3}$)	$1.680\times {10}^{21}$
Number density: $n^{\left(2\right)}$ ($m^{-3}$)	$8.009\times {10}^{20}$
Knudsen number: $({\mathrm{Kn}}_{11},{\mathrm{Kn}}_{22})$	$(0.793,\mathrm{\hspace{0.17em}0.606})$
Knudsen number: $({\mathrm{Kn}}_{12},{\mathrm{Kn}}_{21})$	$(0.803,\mathrm{\hspace{0.17em}0.555})$
Left wall (purely diffuse) boundary conditions (subscript $l$)
Velocity: $u_l$ ($m/s$)	$(0,-50,\mathrm{\hspace{0.17em}0})$
Temperature: $T_l$ ($K$)	273
Right wall (purely diffuse) boundary conditions (subscript $r$)
Velocity: $u_r$ ($m/s$)	$(0,+50,\mathrm{\hspace{0.17em}0})$
Temperature: $T_r$ ($K$)	273

Table 7. Table 7 : Numerical parameters for Ar-Ar self diffusion. The molecular collision parameters for Ar are provided in Table  3 .

Parameter	Case SD-01	Case SD-02
Mixture	Ar-Ar	Ar-Ar
Collision kernel	VSS	VSS
Non-dim physical space	$[0,\mathrm{\hspace{0.17em}1}]$	$[0,\mathrm{\hspace{0.17em}1}]$
Non-dim velocity space	${[-5.09,\mathrm{\hspace{0.17em}5.09}]}^3$	${[-5.09,\mathrm{\hspace{0.17em}5.09}]}^3$
$N^3$	${32}^3$	${32}^3$
$N_{\rho }$	$32$	$32$
$M$	12	12
Spatial elements	4	4
DG order	3	3
Time step (s)	$2\times {10}^{-8}$	$2\times {10}^{-8}$
Viscosity index: (${\omega }_{ij}$)	0.81	0.81
Scattering index: (${\alpha }_{ij}$)	1.4	1.4
Characteristic mass: $m_0$	$m_{\text{Ar}}=m_1$	$m_{\text{Ar}}=m_1$
Characteristic length: $H_0$ ($mm$)	1	1
Characteristic velocity: $u_0$ ($m/s$)	337.2	337.2
Characteristic temperature: $T_0$ ($K$)	273	273
Characteristic number density: $n_0$ ($m^{-3}$)	$1.680\times {10}^{21}$	$8.401\times {10}^{21}$
Initial conditions
Velocity: $u$ ($m/s$)	0	0
Temperature: $T$ ($K$)	273	273
Number density: $n^{\left(1\right)}$ ($m^{-3}$)	$1.680\times {10}^{21}$	$8.401\times {10}^{21}$
Number density: $n^{\left(2\right)}$ ($m^{-3}$)	$1.680\times {10}^{21}$	$8.401\times {10}^{21}$
Knudsen number: $({\mathrm{Kn}}_{11},{\mathrm{Kn}}_{22})$	$(0.793,\mathrm{\hspace{0.17em}0.793})$	$(0.159,\mathrm{\hspace{0.17em}0.159})$
Knudsen number: $({\mathrm{Kn}}_{12},{\mathrm{Kn}}_{21})$	$(0.793,\mathrm{\hspace{0.17em}0.793})$	$(0.159,\mathrm{\hspace{0.17em}0.159})$
Left boundary conditions (subscript $l$)
Ar1 enters: inlet boundary condition for Ar1
Velocity: $u_l$ ($m/s$)	$(0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0})$	$(0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0})$
Temperature: $T_l$ ($K$)	273	273
Number density: $n^{\left(1\right)}$ ($m^{-3}$)	$1.680\times {10}^{21}$	$8.401\times {10}^{21}$
Ar2 freely exits
Right wall (purely diffuse) boundary conditions (subscript $r$)
Ar2 enters: inlet boundary condition for Ar2
Velocity: $u_r$ ($m/s$)	$(0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0})$	$(0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0})$
Temperature: $T_r$ ($K$)	273	273
Number density: $n^{\left(2\right)}$ ($m^{-3}$)	$1.680\times {10}^{21}$	$8.401\times {10}^{21}$
Ar1 freely exits

Table 8. Table 8 : Numerical parameters for Ar-Kr mass diffusion. The molecular collision parameters for Ar-Kr system are provided in Table  3 .

Parameter	Case MD-01	Case MD-02
Mixture	Ar-Kr	Ar-Kr
Collision kernel	VSS	VSS
Non-dim physical space	$[0,\mathrm{\hspace{0.17em}1}]$	$[0,\mathrm{\hspace{0.17em}1}]$
Non-dim velocity space	${[-5.09,\mathrm{\hspace{0.17em}5.09}]}^3$	${[-5.09,\mathrm{\hspace{0.17em}5.09}]}^3$
$N^3$	${32}^3$	${32}^3$
$N_{\rho }$	$32$	$32$
$M$	12	12
Spatial elements	4	4
DG order	3	3
Time step (s)	$2\times {10}^{-8}$	$2\times {10}^{-8}$
Characteristic mass: $m_0$	$m_{\text{Ar}}=m_1$	$m_{\text{Ar}}=m_1$
Characteristic length: $H_0$ ($mm$)	1	1
Characteristic velocity: $u_0$ ($m/s$)	337.2	337.2
Characteristic temperature: $T_0$ ($K$)	273	273
Characteristic number density: $n_0$ ($m^{-3}$)	$1.680\times {10}^{21}$	$8.401\times {10}^{21}$
Initial conditions
Velocity: $u$ ($m/s$)	0	0
Temperature: $T$ ($K$)	273	273
Number density: $n^{\left(1\right)}$ ($m^{-3}$)	$1.680\times {10}^{21}$	$8.401\times {10}^{21}$
Number density: $n^{\left(2\right)}$ ($m^{-3}$)	$8.009\times {10}^{20}$	$4.004\times {10}^{21}$
Knudsen number: $({\mathrm{Kn}}_{11},{\mathrm{Kn}}_{22})$	$(0.793,\mathrm{\hspace{0.17em}0.606})$	$(0.159,\mathrm{\hspace{0.17em}0.121})$
Knudsen number: $({\mathrm{Kn}}_{12},{\mathrm{Kn}}_{21})$	$(0.803,\mathrm{\hspace{0.17em}0.555})$	$(0.161,\mathrm{\hspace{0.17em}0.111})$
Left boundary conditions (subscript $l$)
Ar enters: inlet boundary condition for Ar
Velocity: $u_l$ ($m/s$)	$(0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0})$	$(0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0})$
Temperature: $T_l$ ($K$)	273	273
Number density: $n^{\left(1\right)}$ ($m^{-3}$)	$1.680\times {10}^{21}$	$8.401\times {10}^{21}$
Kr freely exits
Right wall (purely diffuse) boundary conditions (subscript $r$)
Kr enters: inlet boundary condition for Kr
Velocity: $u_r$ ($m/s$)	$(0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0})$	$(0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0})$
Temperature: $T_r$ ($K$)	273	273
Number density: $n^{\left(2\right)}$ ($m^{-3}$)	$8.009\times {10}^{20}$	$4.004\times {10}^{21}$
Ar freely exits

Equations

$${\partial }_t{f}^{(i)}+v\cdot {\nabla }_x{f}^{(i)}=\sum _{j=1}^s{\mathcal{Q}}^{(ij)}(f^{(i)},f^{(j)}),t>0,x\in \Omega \subset {\mathbb{R}}^3,v\in {\mathbb{R}}^3,i=1,2,\dots ,s.$$

$${\mathcal{Q}}^{(ij)}(f^{(i)},f^{(j)})(v)={\int }_{{\mathbb{R}}^3}{\int }_{S^2}{\mathcal{B}}_{ij}(v-v_{\ast },\sigma )\left[f^{(i)}(v^{\prime })f^{(j)}(v_{\ast }^{\prime })-f^{(i)}(v)f^{(j)}(v_{\ast })\right]\mathrm{d}\sigma \mathrm{d}{v}_{\ast },$$

$$m_{i}v+m_j{v}_{\ast }=m_i{v}^{\prime }+m_j{v}_{\ast }^{\prime },m_i|v{|}^2+m_j|v_{\ast }{|}^2=m_i|v^{\prime }{|}^2+m_j|v_{\ast }^{\prime }{|}^2,$$

$$\displaystyle\left\{\begin{array}[]{l}\displaystyle v^{\prime}=\frac{v+v_{*}}{2}+\frac{(m_{i}-m_{j})}{2(m_{i}+m_{j})}(v-v_{*})+\frac{m_{j}}{(m_{i}+m_{j})}|v-v_{*}|\sigma,\\[12.0pt] \displaystyle v_{*}^{\prime}=\frac{v+v_{*}}{2}+\frac{(m_{i}-m_{j})}{2(m_{i}+m_{j})}(v-v_{*})-\frac{m_{i}}{(m_{i}+m_{j})}|v-v_{*}|\sigma,\end{array}\right.$$

$${\mathcal{B}}_{ij}=B_{ij}(|v-v_{\ast }|,\cos \chi ),\cos \chi =\frac{\sigma \cdot (v-v_{\ast })}{|v-v_{\ast }|},$$

$$B_{ij}(|v-v_{\ast }|,\cos \chi )=|v-v_{\ast }|{\Sigma }_{ij}(|v-v_{\ast }|,\chi ),$$

$${\Sigma }_{ij}(|v-v_{\ast }|,\chi )=\frac{b_{ij}}{\sin \chi }\left|\frac{\mathrm{d}{b}_{ij}}{\mathrm{d}\chi }\right|,$$

$$\chi =2{\cos }^{-1}\{(b_{ij}/d_{ij}{)}^{1/{\alpha }_{ij}}\},$$

$$d_{ij}=d_{\mathrm{ref},ij}\Bigg{[}\Bigg{(}\frac{2k_{B}T_{\mathrm{ref},ij}}{\mu_{ij}|v-v_{*}|^{2}}\Bigg{)}^{\omega_{ij}-0.5}\frac{1}{\Gamma(2.5-\omega_{ij})}\Bigg{]}^{1/2},$$

$$B_{ij}=b_{{\omega }_{ij},{\alpha }_{ij}}|v-v_{\ast }{|}^{2(1-{\omega }_{ij})}(1+\cos \chi {)}^{{\alpha }_{ij}-1},$$

$$b_{\omega_{ij},\,\alpha_{ij}}=\frac{d_{\mathrm{ref},ij}^{2}}{4}\Bigg{(}\frac{2k_{B}T_{\mathrm{ref},ij}}{\mu_{ij}}\Bigg{)}^{\omega_{ij}-0.5}\frac{1}{\Gamma(2.5-\omega_{ij})}\,\frac{\alpha_{ij}}{2^{\alpha_{ij}-1}}.$$

$$n^{(i)}={\int }_{{\mathbb{R}}^3}{f}^{(i)}\mathrm{d}v,{\rho }^{(i)}=m_i{n}^{(i)},u^{(i)}=\frac{1}{n^{(i)}}{\int }_{{\mathbb{R}}^3}v{f}^{(i)}\mathrm{d}v,T^{(i)}=\frac{m_i}{3n^{(i)}{k}_B}{\int }_{{\mathbb{R}}^3}(v-u^{(i)}{)}^2{f}^{(i)}\mathrm{d}v.$$

$$n=\sum _{i=1}^s{n}^{(i)},\rho =\sum _{i=1}^s{\rho }^{(i)},u=\frac{1}{\rho }\sum _{i=1}^s{\rho }^{(i)}{u}^{(i)}.$$

$$v_D^{(i)}=\frac{1}{n^{(i)}}{\int }_{{\mathbb{R}}^3}c{f}^{(i)}\mathrm{d}v=u^{(i)}-u,{\mathbb{P}}^{(i)}={\int }_{{\mathbb{R}}^3}{m}_{i}c\otimes c{f}^{(i)}\mathrm{d}v,q^{(i)}={\int }_{{\mathbb{R}}^3}\frac{1}{2}{m}_{i}c|c{|}^2{f}^{(i)}\mathrm{d}v,$$

$$\mathbb{P}=\sum _{i=1}^s{\mathbb{P}}^{(i)},q=\sum _{i=1}^s{q}^{(i)},p=n{k}_{B}T=\frac{1}{3}\text{tr}(\mathbb{P}).$$

$$\begin{array}{rl}& {\int }_{{\mathbb{R}}^3}{\mathcal{Q}}^{(ij)}(f^{(i)},f^{(j)})(v)\phi (v)\mathrm{d}v={\int }_{{\mathbb{R}}^3}{\int }_{{\mathbb{R}}^3}{\int }_{S^2}{\mathcal{B}}_{ij}(v-v_{\ast },\sigma )\left[f^{(i)}(v^{\prime })f^{(j)}(v_{\ast }^{\prime })-f^{(i)}(v)f^{(j)}(v_{\ast })\right]\\ & \cdot \frac{\phi (v)+\phi (v_{\ast })-\phi (v^{\prime })-\phi (v_{\ast }^{\prime })}{4}\mathrm{d}\sigma \mathrm{d}v\mathrm{d}{v}_{\ast },\\ & {\int }_{{\mathbb{R}}^3}{\mathcal{Q}}^{(ij)}(f^{(i)},f^{(j)})(v)\phi (v)\mathrm{d}v+{\int }_{{\mathbb{R}}^3}{\mathcal{Q}}^{(ji)}(f^{(j)},f^{(i)})(v)\varphi (v)\mathrm{d}v\\ =& {\int }_{{\mathbb{R}}^3}{\int }_{{\mathbb{R}}^3}{\int }_{S^2}{\mathcal{B}}_{ij}(v-v_{\ast },\sigma )\left[f^{(i)}(v^{\prime })f^{(j)}(v_{\ast }^{\prime })-f^{(i)}(v)f^{(j)}(v_{\ast })\right]\frac{\phi (v)+\varphi (v_{\ast })-\phi (v^{\prime })-\varphi (v_{\ast }^{\prime })}{2}\mathrm{d}\sigma \mathrm{d}v\mathrm{d}{v}_{\ast }.\end{array}$$

$$\begin{array}{rl}& {\int }_{{\mathbb{R}}^3}{\mathcal{Q}}^{(ij)}(f^{(i)},f^{(j)})\mathrm{d}v=0,\\ & {\int }_{{\mathbb{R}}^3}{\mathcal{Q}}^{(ij)}(f^{(i)},f^{(j)})m_{i}v\mathrm{d}v+{\int }_{{\mathbb{R}}^3}{\mathcal{Q}}^{(ji)}(f^{(j)},f^{(i)})m_{j}v\mathrm{d}v=0,\\ & {\int }_{{\mathbb{R}}^3}{\mathcal{Q}}^{(ij)}(f^{(i)},f^{(j)})m_i|v{|}^2\mathrm{d}v+{\int }_{{\mathbb{R}}^3}{\mathcal{Q}}^{(ji)}(f^{(j)},f^{(i)})m_j|v{|}^2\mathrm{d}v=0,\end{array}$$

$$\sum _{i,j=1}^s{\int }_{{\mathbb{R}}^3}{\mathcal{Q}}^{(ij)}(f^{(i)},f^{(j)})\ln f^{(i)}\mathrm{d}v\le 0.$$

$$\sum _{i=1}^s\left\{{\partial }_t{\int }_{{\mathbb{R}}^3}{f}^{(i)}\ln f^{(i)}\mathrm{d}v+{\nabla }_x\cdot {\int }_{{\mathbb{R}}^3}v{f}^{(i)}\ln f^{(i)}\mathrm{d}v\right\}\le 0,$$

$$\displaystyle f^{(i)}=\frac{n^{(i)}}{(2\pi R_{i}T)^{3/2}}\exp\Bigl{(}-\frac{(v-u)^{2}}{2R_{i}T}\Bigr{)}:=\mathcal{M}^{(i)},$$

$$\begin{array}{rl}& {\partial }_t{\int }_{{\mathbb{R}}^3}{f}^{(i)}\mathrm{d}v+{\nabla }_x\cdot {\int }_{{\mathbb{R}}^3}v{f}^{(i)}\mathrm{d}v=0,\\ & \sum _{i=1}^s\left\{{\partial }_t{\int }_{{\mathbb{R}}^3}{m}_{i}v{f}^{(i)}\mathrm{d}v+{\nabla }_x\cdot {\int }_{{\mathbb{R}}^3}{m}_{i}v\otimes v{f}^{(i)}\mathrm{d}v\right\}=0,\\ & \sum _{i=1}^s\left\{{\partial }_t{\int }_{{\mathbb{R}}^3}\frac{1}{2}{m}_i|v{|}^2{f}^{(i)}\mathrm{d}v+{\nabla }_x\cdot {\int }_{{\mathbb{R}}^3}\frac{1}{2}{m}_{i}v|v{|}^2{f}^{(i)}\mathrm{d}v\right\}=0,\end{array}$$

$$\begin{array}{rl}& {\partial }_t{n}^{(i)}+{\nabla }_x\cdot \left(n^{(i)}{u}^{(i)}\right)=0⟹{\partial }_t\rho +{\nabla }_x\cdot \left(\rho u\right)=0,\\ & {\partial }_t(\rho u)+{\nabla }_x\cdot (\rho u\otimes u+\mathbb{P})=0,\\ & {\partial }_{t}E+{\nabla }_x\cdot \left(Eu+\mathbb{P}u+q\right)=0,\end{array}$$

$$f^{(i)}(t,x,v)=\frac{n_{\mathrm{in}}^{(i)}}{(2\pi R_{i}T_{\mathrm{in}})^{3/2}}\exp\Big{(}-\frac{(v-u_{\mathrm{in}})^{2}}{2R_{i}T_{\mathrm{in}}}\Big{)},\quad x\in\partial\Omega,\quad v\cdot\hat{n}<0,$$

$$f^{(i)}(t,x,v)=(1-\alpha)f^{(i)}(t,x,v-2[(v-u_{w})\cdot\hat{n}]\hat{n})+\alpha\,n^{(i)}_{w}\exp\Big{(}-\frac{(v-u_{w})^{2}}{2R_{i}T_{w}}\Big{)},\quad x\in\partial\Omega,\quad(v-u_{w})\cdot\hat{n}<0,$$

$$n^{(i)}_{w}=-\frac{\int_{(v-u_{w})\cdot{\hat{n}}\geq 0}(v-u_{w})\cdot{\hat{n}}\,f^{(i)}\,\mathrm{d}{v}}{\int_{(v-u_{w})\cdot{\hat{n}}<0}(v-u_{w})\cdot{\hat{n}}\,\exp\Big{(}-\frac{(v-u_{w})^{2}}{2R_{i}T_{w}}\Big{)}\,\mathrm{d}{v}},$$

$$\hat{t}=\frac{t}{t_0},\hat{x}=\frac{x}{H_0},\hat{v}=\frac{v}{u_0},{\hat{m}}_i=\frac{m_i}{m_0},{\hat{f}}^{(i)}=\frac{f^{(i)}}{n_0/u_0^3},$$

$${\hat{\mathcal{B}}}_{ij}=\frac{{\mathcal{B}}_{ij}}{B_{0,ij}},$$

$$B_{0,ij}$$

$$\begin{array}{r}{\partial }_t{f}^{(i)}+v\cdot {\nabla }_x{f}^{(i)}=\sum _{j=1}^s\frac{n_0{H}_0}{u_0}{B}_{0,ij}{\int }_{{\mathbb{R}}^3}{\int }_{S^2}{\mathcal{B}}_{ij}\left[f^{(i)}(v^{\prime })f^{(j)}(v_{\ast }^{\prime })-f^{(i)}(v)f^{(j)}(v_{\ast })\right]\mathrm{d}\sigma \mathrm{d}{v}_{\ast }.\end{array}$$

$$\frac{u_0}{n_0{H}_0{B}_{0,ij}}=\frac{\frac{u_0}{n_0{B}_{0,ij}}}{H_0}={\mathrm{Kn}}_{ij}$$

Full text

††journal: Comput. Methods Appl. Mech. Engrg. [Submitted: 14th Dec. 2018]

\nobibliography

1 Introduction

The Boltzmann equation is an integro-differential equation describing the evolution of the distribution function in six-dimensional phase space. It governs the dilute gas behavior at the molecular level and its solution is required to accurately describe a wide range of non-continuum flow phenomena such as shocks, expansions into vacuum [3] as well as velocity and thermal slip at gas-solid interfaces [4, 5]. Most rarefied flows of technological interest involve gas mixtures with species diffusion playing a decisive role in turbulent, chemically reacting flows, and evaporation/condensation processes [6]. This paper focuses on the development and verification of a deterministic numerical solution to the full Boltzmann equation for gas mixtures.

The physics of Boltzmann equation is now most often simulated computationally using the direct simulation Monte Carlo (DSMC) method [7]. Based on the kinetic theory of gases, DSMC models the binary interactions between particles stochastically. The DSMC method can be rigorously derived as the Monte Carlo solution of the $N$-particle master kinetic equation [8]. Under the assumption that molecular interactions are Markov processes, in the limit of infinite number of particles $N\rightarrow\infty$, Wagner established the convergence of Bird’s DSMC method to the Boltzmann equation [9]. DSMC is widely used for simulating high-speed phenomena, whereas low-speed and unsteady flows are less tractable by stochastic simulations due to the inherent statistical noise.

To avoid the complexity of solving the full Boltzmann equation, many simplified multi-species kinetic models have been proposed and this is a very active research direction in the mathematical and engineering communities, see for instance some early works [10, 11, 12], and more recently [13, 14, 15, 16], and references therein. These simplified models perform better at low Knudsen numbers for flows in the slip and early transition regimes. Yet they often fail to capture the physics at high Knudsen numbers and for diffusion dominated flows at low Knudsen numbers (see [17, 18]). Consequently, in this work, rather than searching for a simple kinetic model to mimic some properties of the Boltzmann equation, we propose a deterministic evaluation of the full multi-species Boltzmann equation with an intention of correctly reproducing the mass, momentum, and energy transport in gas mixtures.

The main difficulty of numerically solving the full Boltzmann equation lies in its complicated collision term. Over the past years, the deterministic methods that approximate the Boltzmann collision operator/equation have undergone considerable development. This includes the discrete velocity methods, spectral methods, etc. The readers are referred to [19, 20] for a comprehensive review. In particular, the Fourier spectral method has been applied to solve the multi-species Boltzmann equation in the past. In [21], a spectral-Lagrangian Boltzmann solver was proposed for a multi-energy level gas for elastic/inelastic interactions with a Lagrangian based post-processing procedure to guarantee the conservation of macroscopic quantities. However, the method was implemented in a straightforward manner without any acceleration strategy and is therefore very expensive. In [22], a fast spectral method was introduced for the multi-species Boltzmann equation along with a strategy to treat large mass ratios. The method is based on the so-called Carleman representation which in its original form can only treat hard sphere molecules [23]. Extension to general collision kernels requires additional assumption on the kernel and parameter fitting/recalibration. This could be a reason that all the numerical tests were restricted to hard spheres in [22]. Recently, a fast Fourier spectral method for the single-species Boltzmann collision operator was introduced in [1]. The complexity for a single evaluation of the collision operator is reduced from $O(N^{6})$ (direct calculation) to $O(MN^{4}\log N)$, where $N$ is the number of discretization points in each velocity dimension, and $M\ll N^{2}$ is the number of discretization points on the sphere. Moreover, the method does not employ any assumptions or parametric fitting on the collision kernel, and is directly applicable for general molecular interactions. Based on [1], a discontinuous Galerkin fast spectral (DGFS) method was proposed in [2] for solving the full single-species Boltzmann equation. DGFS can produce high order spatially and temporally accurate solutions for low-speed and unsteady flows in micro-systems, and is amenable to excellent nearly-linear scaling characteristics on massively parallel architectures [24].

Along similar lines, we develop in this work the DGFS method for the multi-species Boltzmann equation and validate it on various benchmark tests including species diffusion which is of paramount importance in engineering applications. Specifically, we first generalize the method in [1] to derive a fast Fourier spectral method for the multi-species collision operator. The proposed method retains the riveting properties of the single-species fast spectral method including: (a) spectral accuracy in the velocity space, (b) reduced computational complexity compared to direct spectral method, (c) reduced memory requirement in the precomputation, and (d) applicability to general collision kernels. Next, we couple the fast collision algorithm with the discontinuous Galerkin discretization [2] in the physical space to result in a highly accurate deterministic method for the full Boltzmann equation of gas mixtures.

The rest of this paper is organized as follows. In section 2, we give an overview of the multi-species Boltzmann equation, the self/cross collision integrals, H-theorem, and the phenomenological collision kernels used in practical engineering applications. The nondimensionalization of the equation is performed in section 3. Section 4 introduces the fast Fourier spectral method for multi-species Boltzmann collision operator. The discontinuous Galerkin method for the full Boltzmann equation is described in section 5. Results of numerical experiments for the Krook-Wu solution, normal shock, Fourier flow, oscillatory Couette flow, Couette flow, and Fick’s diffusion are presented in section 6. Concluding remarks are given in section 7.

2 The multi-species Boltzmann equation

In this section, we give a brief description of the multi-species Boltzmann equation along with its basic mathematical properties.

Suppose we consider a many-particle system comprised of a mixture of $s$ species ($s\geq 2$). Each species is represented by a distribution function $f^{(i)}(t,x,v)$, where $t$ is time, $x$ is position, and $v$ is particle velocity. $f^{(i)}\,\mathrm{d}{x}\,\mathrm{d}{v}$ gives the number of particles of species $i$ to be found in an infinitesimal volume $\mathrm{d}{x}\,\mathrm{d}{v}$ centered at the point $(x,v)$ of the phase space. The multi-species Boltzmann equation describing the time evolution of $f^{(i)}$ is written as (cf. [25, 26])

[TABLE]

Here $\mathcal{Q}^{(ij)}$ is the collision operator that models binary collisions between species $i$ and $j$, and acts only in the velocity space:

[TABLE]

where $(v,v_{*})$ and $(v^{\prime},v^{\prime}_{*})$ denote the pre- and post- collision velocity pairs. During collisions, the momentum and energy are conserved:

[TABLE]

where $m_{i}$, $m_{j}$ denote the mass of particles of species $i$ and $j$ respectively. Hence one can parameterize $v^{\prime}$ and $v_{*}^{\prime}$ as follows

[TABLE]

with $\sigma$ being a vector varying on the unit sphere $S^{2}$. Finally $\mathcal{B}_{ij}=\mathcal{B}_{ji}(\geq 0)$ is the collision kernel characterizing the interaction mechanism between particles. It can be shown that

[TABLE]

where $\chi$ is the deviation angle between $v-v_{*}$ and $v^{\prime}-v^{\prime}_{*}$.

Given the interaction potential between particles, the specific form of $B_{ij}$ can be determined using the classical scattering theory:

[TABLE]

where $\Sigma_{ij}$ is the differential cross-section given by

[TABLE]

with $b_{ij}$ being the impact parameter. With a few exceptions, e.g. Hard Sphere (HS) model, the explicit form of $\Sigma_{ij}$ can be hard to obtain since $b_{ij}$ is related to $\chi$ implicitly. To avoid this complexity, phenomenological collision kernels are often used in practice with the aim to reproduce the correct transport coefficients. Koura et al. [27] introduced the so-called Variable Soft Sphere (VSS) model by assuming

[TABLE]

where $\alpha_{ij}$ is the scattering parameter, and $d_{ij}$ is the diameter borrowed from Bird’s Variable Hard Sphere (VHS) model (cf. eqn. (4.79) in [7]):

[TABLE]

with $\Gamma$ being the Gamma function, $\mu_{ij}=\frac{m_{i}m_{j}}{m_{i}+m_{j}}$ the reduced mass, $d_{\mathrm{ref},ij}$, $T_{\mathrm{ref},ij}$, and $\omega_{ij}$, respectively, the reference diameter, reference temperature, and viscosity index. Substituting the eqns. (9)-(11) into (8), one can obtain $B_{ij}$ as

[TABLE]

where $b_{\omega_{ij},\,\alpha_{ij}}$ is a constant given by

[TABLE]

In particular, the VHS kernel is obtained when $\alpha_{ij}=1$ and $0.5\leq\omega_{ij}\leq 1$ ($\omega_{ij}=1$: Maxwell molecules; $\omega_{ij}=0.5$: HS); and the VSS kernel is obtained when $1<\alpha_{ij}\leq 2$ and $0.5\leq\omega_{ij}\leq 1$.

Given the distribution function $f^{(i)}$, the number density, mass density, velocity, and temperature of species $i$ are defined as

[TABLE]

The total number density, mass density, and velocity are given by

[TABLE]

Further, the diffusion velocity, stress tensor, and heat flux vector of species $i$ are defined as

[TABLE]

where $c=v-u$ is the peculiar velocity. Finally, the total stress, heat flux, pressure, and temperature are given by

[TABLE]

It can be shown that the collision operator $\mathcal{Q}^{(ij)}$ satisfies the following weak forms:

[TABLE]

Using these weak forms, it is easy to derive

[TABLE]

and the well-known Boltzmann’s H-theorem

[TABLE]

(20) implies that the total entropy of the system decays with time:

[TABLE]

and the equality holds if and only if $f^{(i)}$ attains the local equilibrium

[TABLE]

where $R_{i}=k_{B}/m_{i}$ is the specific gas constant.

On the other hand, using (19), one can take the moments of eqn. (1) to obtain the following local conservation laws:

[TABLE]

which, using the previously defined macroscopic quantities, can be recast as

[TABLE]

where $E=3nk_{B}T/2+\rho u^{2}/2$ is the total energy. Note that this system is not closed. However, replacing $f^{(i)}$ by $\mathcal{M}^{(i)}$ in (24) yields a closed system, i.e., the compressible Euler equations. With more involved calculations (so-called Chapman-Enskog expansion), one can derive the Navier-Stokes equations. We omit the detail but mention that the heat flux term will contain the diffusion velocity $v_{D}^{(i)}$, a property unique to the mixtures (see for instance [26]).

For the eqn. (1), one can consider the in-flow equilibrium boundary condition:

[TABLE]

where $\hat{n}$ is the outward pointing normal at $x$, $n_{\mathrm{in}}^{(i)}$, $T_{\mathrm{in}}$ and $u_{\mathrm{in}}$ are the prescribed density, temperature and velocity. Another commonly used one is the Maxwell boundary condition:

[TABLE]

where $T_{w}$ and $u_{w}$ are the temperature and velocity of the wall, $n^{(i)}_{w}$ is determined from conservation of mass as

[TABLE]

and $\alpha$ is the accommodation coefficient, with $\alpha=1$ corresponds to purely diffusive boundary and $\alpha=0$ to purely reflective boundary.

3 Nondimensionalization

For easier manipulation, we perform a nondimensionalization of the eqn. (1). We first choose the characteristic length $H_{0}$, temperature $T_{0}$, number density $n_{0}$, and mass $m_{0}$, and then define the characteristic velocity $u_{0}=\sqrt{2k_{B}T_{0}/m_{0}}$ and time $t_{0}=H_{0}/u_{0}$. We rescale $t$, $x$, $v$, $m_{i}$, and $f^{(i)}$ as follows:

[TABLE]

and rescale the collision kernel as

[TABLE]

where

[TABLE]

Then the eqn. (1) becomes (dropping $\hat{~{}}$ for simplicity)

[TABLE]

The factor

[TABLE]

is the Knudsen number defined as the ratio of the mean free path and characteristic length scale, hence

[TABLE]

One can also define the “average” Knudsen number for each species $i$ as

[TABLE]

This is consistent with eqn. (4.76) in [7].

Therefore, the dimensionless Boltzmann equation for the VSS kernel (12) reads as

[TABLE]

with

[TABLE]

[TABLE]

Remark 1**.**

We adopt the VSS kernel in this paper for easy comparison with DSMC solutions. The fast algorithm for the collision operator does not rely on the specific form (37) (see Section 4).

In addition, we rescale the macroscopic quantities as

[TABLE]

then in rescaled variables (again dropping $\hat{~{}}$ for simplicity)

[TABLE]

and the Maxwellian (22) becomes

[TABLE]

Remark 2**.**

For the normal shock (see section 6.2.1), it is often convenient to define the so-called parallel ($T^{(i)}_{\parallel}$) and perpendicular ($T^{(i)}_{\perp}$) components of temperature as

[TABLE]

where subscripts $x$ and $y$ denote the first and second components of respective vector fields.

4 A fast Fourier spectral method for the multi-species Boltzmann collision operator

The main difficulty of numerically solving the multi-species Boltzmann equation (35) lies in the collision operator (36). In this section, we introduce a fast Fourier spectral method (in the velocity space) to approximate this operator. Discussion for the spatially inhomogeneous equation will be given in the next section.

We first perform a change of variables $v_{*}$ to $g=v-v_{*}$ in (36) to obtain

[TABLE]

where $\hat{g}$ is the unit vector along $g$ and

[TABLE]

Next we need to choose a finite computational domain $\mathcal{D}_{L}=[-L,L]^{3}$. This is based on the following criterion (similar discussion for the single-species case can be found in [28]).

Assume the support of functions $f^{(i)}$, $f^{(j)}$ can be approximated by a ball with radius $S$: $\text{Supp}(f^{(i)}(v),f^{(j)}(v))\subset\mathcal{B}_{S}$, then one has

$\text{Supp}(\mathcal{Q}^{(ij)}(f^{(i)},f^{(j)})(v))\subset\mathcal{B}_{\sqrt{1+m_{j}/m_{i}}S}$.

This is because if $|v|>\sqrt{1+m_{j}/m_{i}}S$, then $f^{(i)}(v)=0$; also $m_{i}|v^{\prime}|^{2}+m_{j}|v_{*}^{\prime}|^{2}\geq m_{i}|v|^{2}>(m_{i}+m_{j})S^{2}$, then either $|v^{\prime}|>S$ or $|v_{*}^{\prime}|>S$, so $f^{(i)}(v^{\prime})=0$ or $f^{(j)}(v_{*}^{\prime})=0$; either way $\mathcal{Q}^{(ij)}(f^{(i)},f^{(j)})(v)=0$. 2. 2.

It is enough to truncate $g$ to a ball $\mathcal{B}_{R}$ with $R=2S$:

[TABLE]

This is because if $2S<|g|=|v-v_{*}|\leq|v|+|v_{*}|$, then $|v|>S$ or $|v_{*}|>S$, so $f^{(i)}(v)=0$ or $f^{(j)}(v_{*})=0$; also $2S<|g|=|v-v_{*}|=|v^{\prime}-v_{*}^{\prime}|\leq|v^{\prime}|+|v_{*}^{\prime}|$, then $|v^{\prime}|>S$ or $|v_{*}^{\prime}|>S$, so $f^{(i)}(v^{\prime})=0$ or $f^{(j)}(v_{*}^{\prime})=0$; either way $\mathcal{Q}^{(ij)}(f^{(i)},f^{(j)})(v)=0$. 3. 3.

Since $|v|\leq\sqrt{1+m_{j}/m_{i}}S$ and $|g|\leq 2S$ in $\mathcal{Q}^{(ij)}(f^{(i)},f^{(j)})(v)$, we have

$|v_{*}|=|v-g|\leq|v|+|g|\leq(2+\sqrt{1+m_{j}/m_{i}})S$;

$|v^{\prime}|=\left|v-\frac{m_{j}}{m_{i}+m_{j}}g+\frac{m_{j}}{m_{i}+m_{j}}|g|\sigma\right|\leq|v|+\frac{2m_{j}}{m_{i}+m_{j}}|g|\leq(4m_{j}/(m_{i}+m_{j})+\sqrt{1+m_{j}/m_{i}})S$;

$|v_{*}^{\prime}|=\left|v-\frac{m_{j}}{m_{i}+m_{j}}g-\frac{m_{i}}{m_{i}+m_{j}}|g|\sigma\right|\leq|v|+|g|\leq(2+\sqrt{1+m_{j}/m_{i}})S$. 4. 4.

To avoid aliasing, need

[TABLE]

Remark 3**.**

From (45), it can be seen that the computational domain needs to be very large for large mass ratios $m_{j}/m_{i}\gg 1$. This is a common issue appearing in multi-species problems. Possible remedies include adaptive mesh in velocity space (cf. [29]), using an asymptotic model valid for large mass ratios (cf. [30]), or introducing independent velocity grid for each species wherein different collision types for every ($i,j$) pair are treated independently (cf. [31, 22]). In this paper, we only consider moderate mass ratios and postpone these studies to a future work.

Now we approximate $f^{(i)}$ (similarly for $f^{(j)}$) by a truncated Fourier series on $\mathcal{D}_{L}$:

[TABLE]

note here an abuse of notation: the summation over the 3D index $k$ means $-N/2\leq k_{i}\leq N/2-1$, where $k_{i}$ is each component of $k$. Upon substitution of $f^{(i)}$, $f^{(j)}$ into $\mathcal{Q}^{(ij)}(f^{(i)},f^{(j)})$ and a Galerkin projection to the same Fourier space, we obtain the $k$-th Fourier mode of the collision operator as

[TABLE]

with the weight

[TABLE]

Without special treatment, the summation (47) has to be evaluated directly, resulting in a computational cost of $O(N^{6})$. Furthermore, the weight $G^{(ij)}(l,m)$ needs to be precomputed and the storage requirement is $O(N^{6})$. This can quickly become a bottleneck even for moderate $N$. Motivated by our previous work for the single-species Boltzmann equation [1], we propose the following strategy to accelerate the direct summation as well as alleviate its memory bottleneck.

For the gain term (positive part) of $G^{(ij)}(l,m)$, we decompose it as

[TABLE]

where $\rho=|g|$ is the radial of $g$ and

[TABLE]

while for the loss term (negative part) of $G^{(ij)}(l,m)$,

[TABLE]

The idea is to precompute $F^{(ij)}(l+m,\rho,\sigma)$ and $G^{(ij)-}(m)$ up to a high accuracy, and approximate the integral in (48) on the fly using a quadrature rule:

[TABLE]

where for the radial direction, we use the Gauss-Legendre quadrature with $N_{\rho}=O(N)$ points (since the integral oscillates roughly on $O(N)$); for the integral over the sphere, we use the $M$-point spherical design quadrature [32, 33] (usually $M\ll N^{2}$).

Therefore, the gain term of the collision operator can be approximated as

[TABLE]

Written in the above form, we see that the inner sum is a convolution of two functions so that it can be evaluated efficiently in $O(N^{3}\log N)$ operations via the fast Fourier transform (FFT). Together with the outer sum, the total complexity of evaluating $\hat{\mathcal{Q}}^{(ij)+}_{k}$ (for all $k$) is $O(MN^{4}\log N)$ (recall the total number of quadrature points needed for $\rho$ and $\sigma$ is $O(MN)$).

On the other hand, the loss term of the collision operator can be written as

[TABLE]

which is readily a convolution, hence can be evaluated in $O(N^{3}\log N)$.

Putting both pieces together, we have obtained a fast algorithm of complexity $O(MN^{4}\log N)$ for evaluating the collision operator $\mathcal{Q}^{(ij)}(f^{(i)},f^{(j)})$, where $M\ll N^{2}$. In addition, the memory requirement to store the weight $F^{(ij)}(l+m,\rho,\sigma)$ and $G^{(ij)-}(m)$ is $O(MN^{4})$.

5 The discontinuous Galerkin method for the spatial discretization

The previously introduced fast spectral method allows us to compute the collision operator efficiently. To solve the full spatially inhomogeneous equation (35), we also need an accurate and efficient spatial and time discretization. Here we adopt the RKDG (Runge-Kutta discontinuous Galerkin) method [34] widely used for hyperbolic type equations. Since the transport term is linear in the Boltzmann equation, the application of DG method is straightforward. We give a brief description below for completeness.

We first decompose the physical domain $\Omega$ into $N_{e}$ variable-sized disjoint elements $D^{e}_{x}$:

[TABLE]

In each element $D_{x}^{e}$, we approximate the distribution function $f^{(i)}(t,x,v)$ for each species by a polynomial of order $N_{p}$:

[TABLE]

where $\phi_{l}^{e}(x)$ is the basis function supported in $D_{x}^{e}$, $K$ is the total number of terms in the local expansion, and $\mathcal{F}_{e,l}^{(i)}(t,v)$ is the elemental degree of freedom.

We form the residual by substituting the expansion (55) into the eqn. (35):

[TABLE]

where we used the quadratic property of the collision operator. We then require that the residual is orthogonal to all test functions. In the Galerkin formulation, the test function is the same as the basis function, thus

[TABLE]

Substituting (56) into (57) and applying the divergence theorem, we obtain

[TABLE]

where $\hat{n}^{e}$ is the local outward pointing normal and $F^{(i)}_{*}$ denotes the numerical flux. Specifically, the surface integral in the above equation is defined as follows

[TABLE]

with $\hat{n}^{e}_{E}$ and ${F}^{(i)}_{*,\,E}$ being the outward normal and numerical flux along the face $E$. In our implementation, we choose the upwind flux:

[TABLE]

where int and ext denote interior and exterior of the face $e$ respectively.

Finally, define the mass matrix $\mathcal{M}_{ml}$, stiffness matrix $\mathcal{S}_{ml}$, and the tensor $\mathcal{H}_{ml_{1}l_{2}}$ as

[TABLE]

then (58) can be written as

[TABLE]

for $1\leq m\leq K$, $1\leq i\leq s$.

(62) is the DG system we are going to solve in each element $D_{x}^{e}$ of the physical domain. The fast spectral method introduced in the previous section is used to evaluate the term $\mathcal{Q}^{(ij)}\left(\mathcal{F}^{(i)}_{e,\,l_{1}},\mathcal{F}^{(j)}_{e,\,l_{2}}\right)$. The second-order strong-stability-preserving (SSP) RK scheme [35] is applied for the time derivative.

5.1 Structure of $\mathcal{H}_{m\,l_{1}\,l_{2}}^{e}$: a spectral element approach

Needless to say, the main computational bottleneck when solving the system (62) lies in the term $\mathcal{Q}^{(ij)}\left(\mathcal{F}^{(i)}_{e,\,l_{1}},\mathcal{F}^{(j)}_{e,\,l_{2}}\right)$, whose complexity is $O(MN^{4}\log N)$ for given $i$, $j$, $l_{1}$, and $l_{2}$. For general polynomial basis (e.g., the modal DG basis), $\mathcal{H}_{m\,l_{1}\,l_{2}}^{e}$ is a full tensor, hence the total complexity to evaluate the collision part would be $O(s^{2}K^{2}MN^{4}\log N)$ (for all pairs of $(i,j)$ and $(l_{1},l_{2})$) inside each element $D_{x}^{e}$. This is still computationally demanding, even though we are equipped with the fast collision solver. Therefore, the sparsity of $\mathcal{H}_{m\,l_{1}\,l_{2}}^{e}$ would potentially save the computational cost since the collision operator only needs to be evaluated for $l_{1}$, $l_{2}$ such that $\mathcal{H}_{m\,l_{1}\,l_{2}}^{e}\neq 0$. It is known that in the spectral element method [36], if the nodal basis ([37]) is used and the interpolation points are chosen the same as the quadrature points, the mass matrix will become diagonal. Here to achieve better efficiency, we propose to use the same approach to treat the tensor $\mathcal{H}_{m\,l_{1}\,l_{2}}^{e}$. We present the 1D case for simplicity.

Suppose $D_{x}^{e}=[x^{e}_{l},x^{e}_{r}]$, with $x^{e}_{l}$ and $x^{e}_{r}$ being, respectively, the left and right ends of the element $D_{x}^{e}$. $h^{e}=|x^{e}_{r}-x^{e}_{l}|$ is the element size. The DG convention is to define an element in the standard interval $D^{(st)}=[-1,1]$ and map the standard element $D^{(st)}$ to the local element $D_{x}^{e}$ using an affine mapping

[TABLE]

Then

[TABLE]

where $\{\xi_{q},\,w_{q}\}_{q=1}^{N_{q}}$ are the quadrature points and weights.

Consider the Lagrange polynomials as basis functions, i.e.,

[TABLE]

where $\{\xi_{m}\}_{m=1}^{K}$ are the Gauss-Lobatto-Legendre (GLL) quadrature points. When $\{\xi_{q}\}_{q=1}^{N_{q}}$ are taken the same as $\{\xi_{m}\}_{m=1}^{K}$, $\phi_{m}^{(st)}(\xi_{q})=\delta_{mq}$, hence the mass matrix becomes diagonal. Similarly,

[TABLE]

For example, $N_{q}=K=3$ GLL quadrature yields

[TABLE]

Therefore, in this special case, the total complexity to evaluate the collision term $\mathcal{Q}^{(ij)}\left(\mathcal{F}^{(i)}_{e,\,l_{1}},\mathcal{F}^{(j)}_{e,\,l_{2}}\right)$ is reduced to $O(s^{2}KMN^{4}\log N)$ (for all pairs of $(i,j)$ and $(l_{1},l_{2})$ such that $\mathcal{H}_{m\,l_{1}\,l_{2}}^{e}\neq 0$). Of course, this improvement in efficiency comes with an accuracy loss which is quite complicated to analyze. Nevertheless, all numerical results presented in this paper are produced using the above described approach and the bulk properties such as density, temperature, etc. are found to be in good agreement with reference solutions (available finite difference solutions or DSMC solutions). A detailed study of numerical accuracy would be a subject of future work.

6 Numerical experiments

6.1 Spatially homogeneous case: Krook-Wu exact solution

For constant collision kernel, an exact solution to the spatially homogeneous multi-species Boltzmann equation can be constructed (see [38]). We use this solution to verify the accuracy of the proposed fast spectral method for approximating the collision operator. Considering a binary mixture, the equation simplifies to

[TABLE]

where $B_{ij}=B_{ji}:=\frac{\lambda_{ji}}{4\pi n^{(j)}}$ and $\lambda_{ij}$ is some positive constant. The exact solution is given by

[TABLE]

where

[TABLE]

Furthermore, the following condition needs to be satisfied

[TABLE]

For simplicity, we choose $n^{(1)}=n^{(2)}=1$, $\lambda_{11}=\lambda_{22}=1$, $\lambda_{12}=\lambda_{21}=1/2$ but vary the mass ratio $m_{1}/m_{2}$ in the following tests.

It is also helpful to take the derivative of eqn. (69), which yields

[TABLE]

where

[TABLE]

This allows us to check the accuracy of the collision solver without introducing time discretization error.

Figure 1 depicts the convergence behavior of the proposed fast algorithm with respect to $N$ for different mass ratios. Due to the isotropic nature of the solution, we observe that the errors remain relatively unaffected for different $M$ (number of quadrature points used on the sphere). On the other hand, the method exhibits a spectral convergence as $N$ (number of discretization points in each velocity dimension) increases. It is also clear that the accuracy deteriorates for large mass ratios (to keep the same level of accuracy, larger $N$ is needed). To understand the influence of $N_{\rho}$ (number of quadrature points in the radial direction), we list in Table 1 the errors of the method with respect to different $N_{\rho}$. It can be observed that the error is relatively unaffected upon reducing $N_{\rho}$ from $N$ to $N/2$.

Next we evolve the solution using the SSP-RK2 with time step $\Delta t=0.01$. Figure 2 illustrates the time evolution of the distribution function sliced along the velocity domain centerline, i.e., $f^{(i)}(:,N/2,N/2)$. It is observed that: a) the distribution function of the heavy particles becomes more skewed as the mass ratio increases; b) as time goes by, the distribution function tends toward the Maxwellian.

6.2 Spatially inhomogeneous case

6.2.1 Normal shock with HS collision kernel

As a first example in the spatially inhomogeneous case, we consider the normal shock wave and compare our results with the finite difference solutions reported in [39]. Four cases are considered here whose numerical parameters are described in Table 2. The boundary conditions at upstream and downstream are the in-flow equilibrium boundary (see eqn. (25)). We solve the Boltzmann equation until the solution reaches a steady state. A convergence criterion of $(\|f^{n+1}-f^{n}\|_{L^{2}}/\|f^{n}\|_{L^{2}})$ /$(\|f^{2}-f^{1}\|_{L^{2}}/\|f^{1}\|_{L^{2}})<2\times 10^{-5}$ has been used, where $f^{n}$ denotes the distribution function at $n^{th}$ time step.

Figure 3 shows the bulk properties (number density, temperature, velocity, parallel/perpendicular temperature components) for Mach 1.5 normal shock with mass ratios $m_{2}/m_{1}=0.5$ and $m_{2}/m_{1}=0.25$. Based on these results, one can infer that DGFS recovers the normal shock reasonably well. In particular, from Figures 3(a), 3(b), we observe that a Mach 1.5 shock can be captured with just 8 elements within engineering $\pm 5\%$ accuracy. Note that the discontinuity in the flow profile is the characteristic of the DG method. The discontinuity expectedly vanishes upon refining the grid as in Figures 3(c), 3(d).

Figure 4 shows the bulk properties (number density, temperature, and velocity, parallel/perpendicular temperature components) for Mach 1.5 and Mach 3 normal shock for mass ratio $m_{2}/m_{1}=0.5$ at low concentration $n_{-}^{(2)}/n_{-}=0.1$. Again, we observe a fair agreement with the reference solutions.

6.2.2 Solver configurations

In the sections that follow, we consider the standard benchmark cases of Fourier heat transfer, oscillatory Couette flow, Couette flow, and Fick’s diffusion problem at different Knudsen numbers for different collision kernels including VHS and VSS kernels. The results are compared with those obtained from DSMC with equivalent molecular collision models.

All the cases, unless otherwise noted, employ Argon-Krypton mixture. The collision model parameters are tabulated in Table 3 (as provided in [7]). The reference diameters are selected so as to maintain the reference viscosity (cf. eqn. (4.62) in [7]). Note that the viscosity index ($\omega_{ij}$) and scattering index ($\alpha_{ij}$) are empirical parameters, which are calibrated against experiments so that DSMC simulations reproduce experimental observations. The values of these parameters need to be recalibrated for different temperature ranges and different molecules. It is worth noting that there are hundreds of works on recalibration of transport coefficients. In the present DGFS formulation, no recalibration is needed (since the fast collision solver works for general collision kernels), i.e., one can directly use the HS/VHS/VSS model parameters from DSMC literature.

SPARTA [17] has been employed for carrying out DSMC simulations in the present work. It implements the DSMC method as proposed by Bird [7]. The solver has been benchmarked [17] and widely used for studying hypersonic, subsonic and thermal gas flow problems [40, 41, 42, 43, 44]. In this work, cell size less than $\lambda/3$ has been ensured in all test cases. A minimum of 30 DSMC simulator particles per species per cell are used in conjunction with the no-time collision (NTC) algorithm. Each steady-state simulation has been averaged for a minimum 100,000 steps so as to minimize the statistical noise.

More specifically, for all the cases except oscillatory Couette flow, DSMC-SPARTA simulations employ 500 cells, $>$100 particles per cell, a time step of $2\times 10^{-9}$ sec, 1 million unsteady time steps, and 100 million steady time steps. These DSMC parameters have been in part taken from [2] where the authors investigated the single-species rarefied gas flow problems. The parameters have been selected partially to minimize the statistical fluctuations and linear time-stepping errors inherent to DSMC simulations. We, however, note that these parameters are very conservative from a numerical simulation perspective.

6.2.3 Fourier heat transfer of Argon-Krypton mixture using VHS collision kernel

In the current test case, we consider the effect of temperature gradient on the solution. The coordinates are chosen such that the walls are parallel to the $y$ direction and $x$ is the direction perpendicular to the walls. The geometry as well as boundary conditions are shown in Figure 5. We consider six cases for a range of temperature gradients and rarefaction levels. The numerical parameters for these six cases are given in Table 4.

Figure 6 shows the variation of normalized temperature along the domain length for different initial mixture densities: a–b) $\Delta T=20$ (Case F-01, F-02, F-03), and c–d) $\Delta T=100$ (Case F-05, F-06, F-07). The results are compared against DSMC. We note minor ($1-2\%$) discrepancy between DGFS and DSMC for Krypton in the bulk-region away from the walls. Note however that the amount of predicted temperature jump is consistent between DSMC and DGFS for both species.

6.2.4 Oscillatory Couette flow of Argon-Krypton mixture using VHS collision kernel

In the current test case, we consider the effect of transient momentum transport for verifying the temporal accuracy of the DGFS. The schematic remains the same as in the previous test case. The left wall is at rest, and the right wall moves with a velocity of $u=(0,v_{a}\sin{(\zeta t)},0)\;m/s$, where $v_{a}$ is the amplitude of oscillation. The simulation parameters are given in Table 6. The present case is run for two different wall velocities: a) $v_{a}=50\;m/s$, and b) $v_{a}=500\;m/s$. Argon-Krypton mixture with VHS collision model is taken as the working gas. Specifically for DSMC simulations, the domain is discretized into 50 cells with 100000 particles per cell (PPC). For $v_{a}=50\;m/s$ case, a time step of $2\times 10^{-10}$ sec is employed. For $v_{a}=500\;m/s$ case, a time step of $2\times 10^{-11}$ sec is employed. The results are averaged for every 1000 (Navg) time steps. These DSMC simulation parameters have been taken from [2]. Note that such low DSMC time steps are particularly needed for obtaining time accurate results since the time stepping is inherently linear in traditional DSMC method [7].

Figure 7 illustrates the results for the oscillatory Couette flow along the domain length for different $v_{a}$. Ignoring the statistical noise, we observe a good agreement between DGFS and DSMC. Note in particular that for both species, the amount of slip at the left wall are different – which is in accordance with the conservation principles. Moreover, the amount of slip is consistent between DSMC and DGFS.

6.2.5 Couette flow of Argon-Krypton mixture using VSS collision kernel

Phenomenological scattering models are designed and calibrated (against experiments) so as to recover the correct transport properties. VSS model, in particular, recovers two transport properties: a) viscosity and b) diffusion [7]. Couette flow serves as a test case for reproducing the correct viscosity coefficient (the test case for reproducing the correct diffusion coefficient is provided in the later sections). In the current test case, the schematic remains the same as in the previous test case. The left and right parallel walls move with a velocity of $u_{w}=(0,\mp 50,0)\;m/s$. The simulation parameters are given in Table 6. Argon-Krypton mixture with VSS collision kernel is taken as the working gas.

Figure 8 illustrates the velocity and temperature along the domain length for both species. Ignoring the statistical noise, we observe an excellent agreement between DGFS and DSMC.

The viscosity $\mu^{(i)}$ can be recovered from the 1-D Couette flow simulations using the relation between shear-stress and velocity-gradient [45, 44]:

[TABLE]

For consistency, we use eqn. (74) for both DSMC and DGFS. For computing the derivative in (74), we use centered finite difference for DSMC, and the polynomial derivative for DGFS. Figure 9 illustrates the variation of viscosity along the domain for both species. It is observed that: (a) the viscosity is lower for the heavier (Kr) species since the mixture contains $\sim 32\%$ Kr, and $\sim 68\%$ Ar; and (b) both DSMC and DGFS match well within the expected statistical scatter inherent to DSMC simulations. Note that, in the present simulation, we use DG scheme with $K=3$ which implies that the underlying polynomial is quadratic. Hence all the bulk properties including velocity should be a quadratic polynomial. Since the viscosity (74) contains the derivative of the velocity, the overall reconstructed viscosity should be linear, as we observe in Figure 9(a). Upon increasing $K$, we recover the smooth high order polynomial for viscosity as illustrated in Figure 9(b).

From a computation viewpoint, DSMC-SPARTA simulations with 500 cells, 900 particles per cell, a time step of 2e-9 sec, 1 million unsteady time steps, and 20 million steady time steps, on 24 (Intel-Xeon Gold) processors took 26086.45 sec. The parameters have been selected to minimize the statistical fluctuations and linear time-stepping errors inherent to DSMC simulations. On the other hand, DGFS simulations on a single (Titan X Pascal) GPU with 4 elements, $K=3$, $N^{3}=32^{3}$, $M=12$ took 6020.19 sec to achieve $(\|f^{n+1}-f^{n}\|_{L^{2}}/\|f^{n}\|_{L^{2}})/(\|f^{2}-f^{1}\|_{L^{2}}/\|f^{1}\|_{L^{2}})<5\times 10^{-6}$. Note that these are representative simulation times for indicating the computational efforts required in DGFS and DSMC for 1-D simulations. Our experience shows that even heavily tuned codes can be further improved. A detailed comparison between CPU and GPU performance is subject of future study.

6.2.6 Self diffusion of Argon-Argon mixture using VSS collision kernel

In the current test case, we consider the effect of diffusive transport. The schematic remains the same as in the previous test case. Argon-Argon mixture with VSS collision kernel is taken as the working gas. To differentiate between two types of Argon, we tag the molecules as $\text{Ar}_{1}$ and $\text{Ar}_{2}$. At the left boundary, $\text{Ar}_{1}$ enters and exits at the right boundary. At the right boundary, $\text{Ar}_{2}$ enters and exits at the left boundary. The molecules enter the domain with zero mean velocity. The simulation parameters are provided in Table 7.

Figure 10(a) shows the variation of concentration $n^{(i)}/n$ along the domain. Since the species-1 enters from the left boundary and exits at right, we observe a drop in species-1 concentration as we move towards the right boundary. Conversely for species-2, since the species-2 enters from the right boundary and exits at left, we observe a drop in species-2 concentration as we move towards the left boundary. The non-linearity of the concentration profile and the associated slip at the boundaries can be explained through the low mixture density, and the fact that the flow is in slip regime. It is also worth noting that throughout the domain at any given $x$ location, the sum of the concentrations of two species is unity which asserts that the numerical formulation is conservative.

Figure 10(b) shows the variation of diffusion velocity along the domain. Since the species-1 enters from the left boundary and exits at right, we observe a low net diffusion speed (magnitude of the diffusion velocity) for the first species and a high diffusion speed for the second species. Conversely at the right boundary, since the species-2 enters from the right boundary and exits at left, we observe a low diffusion speed for the second species and high diffusion speed for the second species.

Figure 11 illustrates the temperature profile along the domain, where we observe a drop in temperatures of the two species. Based upon these results, it can be inferred that DGFS can resolve the strong gradients in temperature and diffusion velocity with just 4 elements and $K=3$ within engineering accuracy.

For this test case, the self-diffusion coefficient is given as (cf. eqn. (12.18) in [7])

[TABLE]

While writing this equation, it is also assumed that the coefficient of thermal diffusion is low, and therefore the effect of temperature gradient is negligible (see eqn. (8.4.7) in [46]). Note that this equation is an approximation to the diffusion equation, derived from leading order Chapman expansion (see section 8.4 in [46]), and therefore, strictly speaking, the values computed from this equation might not be fully accurate especially for the rarefied flows, since the higher order terms have not been accounted for. It is worth noting that the full diffusion equation based on moment of the distribution function is highly non-trivial from a computation perspective, and is therefore rarely implemented in DSMC codes. For consistency, we use (75) for both DSMC and DGFS. For computing the derivatives in (75), we use centered finite difference for DSMC, and the polynomial derivative for DGFS. In particular, for DSMC simulations, we used 500 cells, 2000 particles per cell, a time step of $1\times 10^{-8}$ sec, and averaged the results for 1 million time steps to minimize the statistical scatter in diffusion coefficients.

Figure 12 illustrates the variation of self-diffusion coefficient along the domain as a function of scattering parameter $\alpha_{ij}$. It is observed that: (a) the diffusion coefficient increases with increase in $\alpha_{ij}$ in accordance with the VSS model (cf. eqn. (3.75) in [7]), (b) both DSMC and DGFS match well within the expected statistical scatter inherent to DSMC simulations, and (c) with increase in number density from Case 01 to Case 02, the diffusion coefficient decreases in accordance with (75). Note that, in the present simulation, we use DG scheme with $K=3$ which implies that the underlying polynomial is quadratic. Hence all the bulk properties including number density should be a quadratic polynomial. Recall that the diffusion (75) contains the derivative of the number density, and hence the overall reconstructed diffusion coefficient should be linear, which is what we observe in Figures 12(a),12(b). Upon increasing $K$, we recover the smooth high-order polynomial for diffusion coefficient.

6.2.7 Mass diffusion of Argon-Krypton mixture using VSS collision kernel

In the current test case, we consider the effect of mass diffusion. The conditions remain the same as in previous case, except that Argon-Krypton mixture with VSS collision model is taken as the working gas. More specifically, Argon enters the left boundary and exits at the right boundary; and Krypton enters through the right and exits at left. The molecules enter the domain with zero mean velocity. We consider two cases with different initial number density. The numerical parameters for both the cases are given in Table 8.

Figure 13(a) shows the variation of concentration profile for the two species. We observe that the concentration of Argon remains greater than Krypton throughout the domain, except for a small portion near the right boundary. This can be directly inferred from the mass/momentum conservation principle i.e., the heavier species diffuses slower and the lighter species diffuses faster. Therefore, after a sufficiently long time, the concentration of lighter species will be greater than that of heavier species in the major part of the domain. As in the self-diffusion case, the sum of the concentrations of both species is unity throughout the domain at any given $x$ location. The effect of the momentum conservation is more pronounced in the Figure 13(b) where we observe a higher diffusion speed for the lighter species and a lower diffusion speed for the heavier species.

7 Conclusions

A fast spectral method for the multi-species Boltzmann collision operator has been proposed in this work. The method is designed to handle the cross-molecular interactions between dissimilar species with moderate mass ratios. In particular, it is applicable to general collision kernels which allows us to directly compare our results against the well-known stochastic DSMC solutions. The fast collision algorithm in the velocity space was then coupled with the discontinuous Galerkin discretization in the physical space to yield highly accurate numerical solutions for the full spatially inhomogeneous Boltzmann equation. The DG-type formulation employed in the present work has advantage of having high order accuracy at the element-level, and its element-local compact nature (and that of our collision algorithm) enables effective parallelization on massively parallel architectures.

To validate our solver, extensive numerical tests were performed, including the spatially homogeneous Krook-Wu solution for Maxwell molecules where the exact solution is known; normal shock wave for HS where finite difference solutions are available for comparison, and Fourier, oscillatory Couette, Couette, self diffusion, mass diffusion problems where different collision kernels (VHS and VSS), Knudsen numbers, and mass ratios were considered and the results were compared well with DSMC solutions.

Acknowledgement

This work was supported by NSF grant DMS-1620250 and NSF CAREER grant DMS-1654152.