High-order quadrature-based lattice Boltzmann models for the flow of ultrarelativistic rarefied gases

We present a systematic procedure for the construction of relativistic lattice Boltzmann models (R-SLB) appropriate for the simulation of flows of massless particles. Quadrature rules are used for the discretization of the momentum space in spherical coordinates. The models are optimized for one-dimensional flows. The applications considered in this paper are the Sod shock tube and the one-dimensional boost invariant expansion (Bjorken flow). Our models are tested against exact solutions in the inviscid and ballistic limits. At intermediate relaxation times (finite viscosity), we compare with the results obtained using the Boltzmann approach of multiparton scattering model for the Sod shock tube problem, as well as with a semianalytic solution for the nonideal Bjorken flow. In all cases our models give remarkably good results. We define a convergence test in order to find the quadrature order necessary to obtain convergence at a predefined accuracy. We show that, while in the hydrodynamic regime the number of velocities is comparable to that required for the more popular collision-streaming type models, as we go towards the ballistic regime, the size of the velocity set must be substantially increased in order to accurately reproduce the analytic profiles.

Figure 1

The flow structure in the inviscid limit of the Sod shock tube problem.

Figure 2

Density (a), pressure (b), temperature (c), fugacity (d) and velocity (e) profiles for decreasing values of the ratio $\left(\overline{\eta /s}\right)$, compared with the analytical results for the inviscid limit, presented in Sec. 4b1. The quadrature order is $Q_{\xi }=6$, the number of points is $Z=10000$, the time step was set to $5\times {10}^{-6}$ and the snapshot corresponds to $t=0.5$ (i.e., $100000$ iterations were performed). The initial conditions are $(n_{\text{L}},T_{\text{L}},P_{\text{L}})=(1,1,1)$ and $(n_{\text{R}},T_{\text{R}},P_{\text{R}})=(0.125,0.5,0.0625)$.

Figure 3

Profiles of $q$ (a) and $\mathrm{\Pi }$ (b) for decreasing values of $\left(\overline{\eta /s}\right)$. (c) The integrated values of $q$ and $\mathrm{\Pi }$ over small vicinities around the contact discontinuity and the shock front are shown in absolute value as functions of $\left(\overline{\eta /s}\right)$ in log-log scale (the integral of $\mathrm{\Pi }$ around the contact discontinuity is not shown, since $\mathrm{\Pi }$ vanishes in this region). The analytic results (4.43) and (4.44) are also represented for comparison, for values of the transport coefficients obtained with both the Chapman-Enskog (solid lines) and Grad (dashed lines) methods. The simulation parameters are presented in the caption of Fig. 2.

Figure 4

Profiles of the density (a), pressure (b), temperature (c), relative fugacity (d), and velocity (e) for the initial conditions in Eq. (4.34), for different orders of the polar quadrature. The profiles corresponding to $Q_{\xi }=200$ are overlapped with the analytic solutions given in Sec. 4c1. The profiles represent snapshots at $t=0.4$, obtained using a number of $Z=1000$ nodes with a time step $\delta t=5\times {10}^{-4}$. The initial state is described in the caption of Fig. 2.

Figure 5

Ballistic regime profiles of the heat flux (a) and pressure deviator (b), for different orders of the polar quadrature. The simulation parameters are described in the caption of Fig. 4.

Figure 6

Comparison between the profiles obtained with the LB model with $Q_{\xi }=500$ and $N_{\mathrm{\Omega }}=5$ and the results obtained using the BAMPS and vSHASTA methods, which were reported in Refs. [91, 93], for various values of $\left(\overline{\eta /s}\right)$. (a) Density $n$; (b) relative fugacity $\overline{\lambda }$; (c) pressure $P$; (d) velocity $\beta =u^z/u^0$; (e) pressure deviator $\mathrm{\Pi }$; (f) heat flux $q$. In (b), (e), and (f), the vSHASTA results are also included.

Figure 7

Comparison between the profiles of the pressure $P$, obtained with the LB model with $Q_{\xi }=500$ and $N_{\mathrm{\Omega }}=5$ (lines), and the results obtained using the BAMPS method reported in Ref. [92] (points), for various values of $\left(\overline{\eta /s}\right)$.

Figure 8

(a) The error $\epsilon$ (4.64) as a function of $N_{\mathrm{\Omega }}$ for various values of $\left(\overline{\eta /s}\right)$ ($Q_{\xi }=500$ is kept constant). The reference profiles considered here were obtained using $Q_{\xi }=500$ and $N_{\mathrm{\Omega }}=6$. (b) The minimum expansion order $N_{\mathrm{\Omega }}^{\text{conv}}$ required to reduce the error (4.64) of a model with $Q_{\xi }=500$ below the threshold ${\epsilon }_{\text{th}}=1\%$.

Figure 9

(a) The error $\epsilon$ (4.64) as a function of $Q_{\xi }$ for various values of $\left(\overline{\eta /s}\right)$. The reference profiles considered here were obtained using $Q_{\xi }=500$ and $N_{\mathrm{\Omega }}=5$. (b) The minimum quadrature order $Q_{\xi }^{\text{conv}}$ required to reduce the error (4.64) of a model with $N_{\mathrm{\Omega }}=\min (Q_{\xi }-1,5)$ below various values of the threshold ${\epsilon }_{\text{th}}$. The initial state for this analysis is given in Eq. (4.34).

Figure 10

The minimum quadrature order $Q_{\xi }^{\text{conv}}$ required to reduce the error (4.64) of a model with $N_{\mathrm{\Omega }}=\min (Q_{\xi }-1,5)$ below various values of the threshold ${\epsilon }_{\text{th}}$ for the initial state given in Eq. (4.61).

Figure 11

Evolution of $\tau n$ (a), $\tau P$ (b), and $\tau \mathrm{\Pi }$ (c) with respect to $\tau$. The $\tau$ factor is included in order to highlight the asymptotic limits (5.10). The numerical results corresponding to various $Q_{\xi }$ (dashed lines and points) are compared to the analytic solutions (solid lines) given in Eqs. (5.5) and (5.9), which are overlapped with the $Q_{\xi }=40$ curves.

Figure 12

Time evolution of (left) ${\tau }^{4/3}P$ and (right) ${\tau }^2\mathrm{\Pi }$ for ${\tilde{k}}_B{\tilde{T}}_0=0.3\mathrm{GeV}$ (top) and $0.6\mathrm{GeV}$ (bottom) for $4\pi \left(\overline{\eta /s}\right)\in \{0.1,0.5,1.0\}$ for partons ($\mu =0$). The numerical results (dotted lines and points) are compared with the first-order (dotted lines) and second-order (solid lines) hydrodynamics results obtained from Eqs. (5.18) and (5.19), respectively. The factors ${\tau }^{4/3}$ and ${\tau }^2$ multiplying $P$ and $\mathrm{\Pi }$ are inspired from the form of the first-order hydrodynamics solutions (5.18).

Figure 13

Evolution with respect to $\tau$ of (left) $T$ and (right) ratio ${\mathcal{P}}_L/{\mathcal{P}}_T$ of the longitudinal to the transverse pressure, for ${\tilde{k}}_B{\tilde{T}}_0=0.3\mathrm{GeV}$ (top) and $0.6\mathrm{GeV}$ (bottom), for various values of $\left(\overline{\eta /s}\right)$, under the assumption of vanishing chemical potential discussed in Sec. 5c. The numerical results (dotted lines and points) are validated using benchmark results (solid lines) as follows: for $\left(\overline{\eta /s}\right)=\infty$, the analytic formulas (5.9) and (5.14) are used; at $4\pi \left(\overline{\eta /s}\right)=0.3$, the numerical solution of the second-order hydrodynamics equations (5.19) is used; when $4\pi \left(\overline{\eta /s}\right)\in \{1,3,10\}$, the semianalytic solution reported in Ref. [101] is represented.

Figure 14

Time evolution of (left) $P$ and (right) $\mathrm{\Pi }$ for $4\pi \left(\overline{\eta /s}\right)\in \{0.1,0.5,1\}$ for initial temperatures ${\tilde{k}}_B{\tilde{T}}_0=0.3\mathrm{GeV}$ (top) and $0.6\mathrm{GeV}$ (bottom). The various curves show the numerical results obtained using LB for the cases of the ideal gas (dotted lines and points) and of the gas of partons (solid black lines), as well as the solutions of the first (dotted lines) and second (solid purple lines) hydrodynamics equations (5.24) and (5.25) for the ideal gas. The ideal and parton gas results seem to be overlapped for all tested values of $\left(\overline{\eta /s}\right)$ and ${\tilde{T}}_0$.

Figure 15

Evolution of the ratios (left) $T_{\text{ideal}\text{gas}}/T_{\text{partons}}$ and (right) $n_{\text{ideal}\text{gas}}/n_{\text{partons}}$ of the temperatures and densities obtained for the ideal ($P=nT$) and parton (${\mu }_{\mathrm{eq}}=0$) gases, for ${\tilde{k}}_B{\tilde{T}}_0=0.3\mathrm{GeV}$ (top) and $0.6\mathrm{GeV}$ (bottom) for $\left(\overline{\eta /s}\right)\in \{0.3,1,3,10,\infty \}$.

Figure 16

(a) Maximum $L_2$ norm (5.27) as a function of $Q_{\xi }$ for ${\tilde{k}}_B{\tilde{T}}_0=0.3\mathrm{GeV}$, under the assumption of a parton gas ($\mu =0)$. (b),(c) Dependence with respect to $4\pi \left(\overline{\eta /s}\right)$ of $Q_{\xi }^{\mathrm{conv}}$ for $L_2^{\text{th}}={10}^{-5}$ (b) and $L_2^{\text{th}}={10}^{-6}$ (c), obtained for ${\tilde{k}}_B{\tilde{T}}_0\in \{0.3\mathrm{GeV},0.6\mathrm{GeV}\}$, for both the parton and ideal gases. The $L_2$ norm is computed with respect to the reference model $\text{R-SLB}(5;100)$.