Solving the non-gray Callaway phonon Boltzmann transport equation allows a dual-relaxation-time approximation of separate normal and resistive scatterings and resolving the mode-dependent spectrum. The conventional iterative scheme (CIS) for deterministic solutions avoids a monolithic phase-space inversion, but its collision-source iteration can become prohibitively slow at a large characteristic length of a material. Existing synthetic acceleration schemes address either non-gray single-relaxation models or gray dual-relaxation models, leaving mode-resolved dual-relaxation transport without a dedicated acceleration framework. We develop a general synthetic iterative scheme (GSIS) for the stationary, linearized, non-gray Callaway equation, where synthetic approximations for the normal-process pseudo-temperature and phonon drift velocity are provided by exact energy and quasi-momentum balance laws closed with first-order Chapman-Enskog constitutive relations and non-equilibrium terms evaluated from the kinetic solution. The resistive-process pseudo-temperature is retrieved from the two quantities. Each iteration couples an upwind nodal discontinuous Galerkin kinetic sweep and a hybridizable discontinuous Galerkin solution of the synthetic equations to achieve high-order spatial discretization. A branch- and frequency-resolved Fourier analysis identifies the deterioration of CIS and shows that the GSIS contraction factor remains bounded away from unity for the considered graphene material. Asymptotic analysis indicates that GSIS reduces to a consistent discretization of a Guyer-Krumhansl-like equation and Fourier's law of heat conduction in the hydrodynamic and diffusive limits, respectively.

Modelling rarefied gas flow using the Boltzmann equation is vital in many areas. Due to the high dimensionality and coexistence of multiple characteristic scales, conventional solution strategies to this equation incur prohibitively high computational costs and are inadequate for rapid response in engineering design simulations. Based on proper generalised decomposition (PGD), we propose an a priori, asymptotic-preserving reduced-order method to solve the high-dimensional, parametrised Shakhov kinetic model equation. The method reduces the original problem to a few low-dimensional problems by formulating separated representations for the low-rank solution, thereby mitigating the curse of dimensionality. To capture the hydrodynamic asymptotics, we incorporated solutions of some synthetic equations into the PGD algorithm. This treatment allows the PGD solver to automatically reduce to a macroscopic solver for the Navier-Stokes equations, whose solution naturally exhibits low-rank structure. By treating the rarefaction parameter as an additional coordinate, a parametrised solution can be computed once and for all over the entire range of rarefaction, enabling fast multiple queries to any points in the parameter space. Numerical examples are presented to demonstrate the capability of the method to simulate rarefied gas flow with certain accuracy and a significant reduction in computational costs.

Understanding and controlling electron-phonon interactions is essential for optimizing thermal performance in microelectronic and thermoelectric devices, where strong non-equilibrium effects span multiple length and time scales. The coupled electron-phonon Boltzmann equations offer an insightful framework for modeling such transport, but their high dimensionality and stiffness pose significant computational challenges. Conventional iterative solvers typically suffer from excessive numerical dissipation and slow convergence in diffusive regimes. In this work, we develop a general synthetic iterative scheme (GSIS) that significantly accelerates the solution of coupled electron-phonon Boltzmann equations. The key innovation is the formulation of macroscopic synthetic equations that effectively capture diffusion-limit behavior, while precisely incorporating non-equilibrium corrections extracted from the kinetic level. During iterations, the kinetic solver provides high-order moment closures, while the macroscopic equations update the driving fields, ensuring rapid convergence. This strategic coupling across scales facilitates efficient global information exchange and robust multiscale resolution. Fourier analysis and numerical benchmarks show that GSIS can reduce computational time by up to three orders of magnitude compared to the conventional iterative scheme. Implemented with a high-order discontinuous Galerkin method, the GSIS is applicable to complex geometries and extendable to nonlinear and non-gray transport models.

Gas-radiation coupling critically influences hypersonic reentry flows, where extreme temperatures induce pronounced non-equilibrium gas and radiative heat transport. Accurate and efficient simulation of radiative gas dynamics is therefore indispensable for reliable design of thermal protection systems for atmospheric entry vehicles. In this study, a Boltzmann-type kinetic model for radiative gas flows is solved across a broad spectrum of flow and radiation transport regimes using the general synthetic iterative scheme (GSIS). The approach integrates an unstructured finite-volume discrete velocity method with a set of macroscopic synthetic equations. Within this framework, the kinetic model provides high-order closures for the constitutive relations in the synthetic equations. Simultaneously, the macroscopic synthetic equations drive the evolution of the mesoscopic kinetic system, significantly accelerating steady-state convergence in near-continuum regimes, as substantiated by linear Fourier stability analysis. Crucially, the algorithm is proven to be asymptotic-preserving, correctly recovering the continuum and optically thick limits—represented by the radiative Navier–Stokes–Fourier equations governing distinct translational, rotational, vibrational, and radiative temperatures—on coarse meshes independent of the mean free path. Numerical simulations of challenging benchmarks, including three-dimensional hypersonic flow over an Apollo reentry capsule, demonstrate that GSIS achieves orders-of-magnitude speedup over conventional iterative schemes in multiscale simulations of radiative gas flows while accurately capturing non-equilibrium effects and radiative heat transfer in hypersonic environments.

Experiments have demonstrated that the boiling heat transfer coefficient may be largely improved by engineering the heating surfaces using micro/nanostructures. While experiments are visually intuitive, the underlying mechanisms for this improvement remain unclear. In this paper, the boiling process on surfaces with various pillar-textured structures and wettability is studied using a two-dimensional lattice Boltzmann model. Our results demonstrate that wettability and surface structures significantly influence the bubble dynamics. The critical heat flux can be effectively increased through vortex generation structures such as pillar-groove structures, which disrupt the vapor film and then enhance the convective heat transfer. Our simulations also show that the width of grooves should exceed the characteristic length of the bubble radius, while not too large, can ensure a continuous bubble release and sufficient nucleation sites. Furthermore, a general relation based on the Arrhenius equation between the average heat flux (qave) and the standard deviation of the velocity field ( sigma ( u ) ) caused by vortexes, i.e., q ave = 6.8 x 10 - 4 exp ( - 3.9 x 10 - 3 / sigma ( u ) ) + 5.8 x 10 - 5 , is proposed. Results show that the convex corner of the pillar-groove structure can maximally enhance the strength of vortex convection within the fluid and maximumly increase the heat flux by 13.8%. Our work here provides a better understanding of bubble behaviors and boiling characteristics on surfaces with different wettability/pillar textures, which may facilitate the design of a strategy to improve the boiling heat transfer.

Numerical investigations and analyses are carried out on the interactions of low enthalpy hypersonic 30-55 degrees double wedge configuration, particularly focusing on steady cases at conditions similar to the experimental setup by Swantek & Austin [AIAA 2012-284], with Ma = 7 and h0 = 2.1MJ/kg. To achieve a steady solution, Reynolds numbers (Re) lower than those in the experiment are used. For increased accuracy, a third-order scheme WENO3 - PRM21,1 [Li et al., J. Sci. Comput., 88(3) (2021) 75-130] with improved resolution is employed. Meanwhile, three gas models, i.e., the perfect, equilibrium, and non-equilibrium gas models, are used to analyze the difference potentials that arise from the physical model. After validating the methods, grid convergence studies are first conducted at Ma = 7 and Re = 2.5 x 105/m, to determine the appropriate grid resolution for the main computations. Subsequently, comprehensive numerical studies are carried out on the steady interactions and their evolution at Ma = 7 and h0 = 2.1MJ/kg. Specifically: (a) The upper limits of Re are identified where the flows remain steady with the transmitted shock impinging on the aft wedge, and the corresponding interaction characteristics as well as differences in the three gas models are investigated qualitatively and quantitatively, e.g., the shock system, vortex structures, distributions such as pressure, Mach number, and specific heat ratio. Notably, a quasi-normal shock wave is observed within the slip line passage in the case of the perfect gas model. (b) The flow characteristics of the three models, including the interaction pattern, geometric features of triple points, impingements, and separation zone, are studied and compared for Re = 4, 3, and 2 x 104/m. Differences primarily emerge between the results of the perfect gas model and the real gas models. Specifically, a transmitted shock reflecting above the separation zone is observed in the case of the perfect gas model. The effect of the gas model on temperature and specific heat ratio distributions, as well as the heat transfer and pressure coefficients over the wedge surface are investigated. For an in-depth understanding, the shock polar method is applied for comparison with computational results, while a 1D flow model is proposed to explain the occurrence of the quasi-normal shock wave. Consequently, overall reasonable agreements are achieved. Finally, the effects of variations in Mach number and enthalpy are determined, by alternatively varying the two parameters around Ma = 7 and h0 = 2.1MJ/kg at Re = 4 x 104/m, focusing on alterations in interaction characteristics, thermodynamic properties, and aerodynamic performance.

Modelling rarefied gas flow via the Boltzmann equation plays a vital role in many areas. Due to the high dimensionality of this kinetic equation and the coexistence of multiple characteristic scales in the transport processes, conventional solution strategies incur prohibitively high computational costs and are inadequate for rapid response for parametric analysis and optimisation loops in engineering design simulations. This paper proposes an a priori reduced-order method based on the proper generalised decomposition to solve the high-dimensional, parametrised Shakhov kinetic model equation. This method reduces the original problem into a few low-dimensional problem by formulating separated representations for the low-rank solution, as well as data and operators in the equation, thereby overcoming the curse of dimensionality. Furthermore, a general solution can be calculated once and for all in the whole range of the rarefaction parameter, enabling fast and multiple queries to a specific solution at any point in the parameter space. Numerical examples are presented to demonstrate the capability of the method to simulate rarefied gas flow with high accuracy and significant reduction in CPU time and memory requirements.

Modeling nano- and micro-scale heat conduction based on the phonon Boltzmann transport equation has gained increasing research interest due to the demand for better thermal performance of semiconductors. Nevertheless, the high dimensionality of the Boltzmann equation results in the so-called curse of dimensionality, presenting a bottleneck for efficient numerical solutions based on direct discretization. In practice, high-order numerical schemes such as discontinuous Galerkin finite element methods are preferable to reduce the degrees of freedom, thereby reducing the computational cost. However, when complex geometries emerge, cumbersome refinement is required to approximate the boundary of the computational domain if spatial meshes with straight-sided elements are employed, concealing the advantage of a high-order scheme. In this work, we extend the idea of the non-uniform rational B-splines enhanced finite element method. By embodying accurate geometric information, including parametric descriptions for curved boundaries and sampling information of rough surfaces reconstructed from scanning electron microscope images or by a random growth approach, into the faces of the elements adjacent to the physical boundary, the geometric inaccuracies and heavy refinement can be eliminated in a very coarse mesh. Strategies to define the polynomial basis and compute the integrals over the geometry-embodied elements are investigated. Numerical results, including heat conduction in a silicon ring, nano-porous media with circular pores, and a square domain with a rough boundary, show that to obtain solutions with the same order of accuracy, the discontinuous Galerkin method performed on accurate-geometry-embodied meshes can be 10-100 times faster than that implemented on straight-sided meshes. The efficiency of higher-order discretization methods is fully promoted, where fewer spatial elements combined with higher-order approximating polynomials are preferable to obtain solutions with high accuracy and reduced computational cost.

The formation of an artificial freezing curtain is severely constrained by groundwater seepage flow. Traditional numerical methods have limitations in addressing the interaction between seepage flow and freezing soil along a moving freezing front. An enthalpy-based lattice Boltzmann model was proposed to study the development of an artificial freezing curtain subjected to seepage flow. The proposed numerical model was validated against two analytical solutions. Finally, developing artificial freezing soil with a single freeze pipe and a row of three freeze pipes under seepage flow was discussed. For both cases, seepage flow restrained the formation and development of freezing soil. For a single pipe, the development of the freeze radius is restrained the most upstream, followed by midstream and downstream. For a row of three pipes, the closure position moved towards downstream compared with when there is no seepage flow. Compared with the initial ground temperature and pipe spacing, pipe diameter has less effect on the closure position. The closure time is approximately linearly related to the water-facing length, and the slope increases with seepage velocity.

Recent years have seen the emergence of new technologies that exploit nanoscale evaporation, ranging from nanoporous membranes for distillation to evaporative cooling in electronics. Despite the increasing depth of fundamental knowledge, there is still a lack of simulation tools capable of capturing the underlying non-equilibrium liquid-vapour phase changes that are critical to these and other such technologies. This work presents a molecular kinetic theory model capable of describing the entire flow field, i.e. the liquid and vapour phases and their interface, while striking a balance between accuracy and computational efficiency. In particular, unlike previous kinetic models based on the isothermal assumption, the proposed model can capture the temperature variations that occur during the evaporation process, yet does not require the computational resources of more complicated mean-field kinetic approaches. We assess the present kinetic model in three test cases: liquid-vapour equilibrium, evaporation into near-vacuum condition, and evaporation into vapour. The results agree well with benchmark solutions, while reducing the simulation time by almost two orders of magnitude on average in the cases studied. The results therefore suggest that this work is a stepping stone towards the development of an accurate and efficient computational approach to optimising the next generation of nanotechnologies based on nanoscale evaporation.

A kinetic model is proposed for the nonequilibrium flow of dense gases composed of hard-sphere molecules, which significantly simplifies the collision integral of the Enskog equation using the relaxation-time approach. The model preserves the most important physical properties of high-density gas systems, including the Maxwellian at rest as the equilibrium solution and the equation of state for hard-sphere fluids; all the correct transport coefficients, namely, the shear viscosity, thermal conductivity, and bulk viscosity; and inhomogeneous density distribution in the presence of a solid boundary. The collision operator of the model contains a Shakhov model-like relaxation part and an excess part in low-order spatial derivatives of the macroscopic flow properties; this latter contribution is used to account for the effect arising from the finite size of gas molecules. The density inhomogeneity in the vicinity of a solid boundary in a confined flow is captured by a method based on the density-functional theory. Extensive benchmark tests are performed, including the normal shock structure and the Couette, Fourier, and Poiseuille flow at different reduced densities and Knudsen numbers, where the results are compared with the solutions from the Enskog equation and molecular dynamics simulations. It is shown that the proposed kinetic model provides a fairly accurate description of all these nonequilibrium dense gas flows. Finally, we apply our model to simulate forced wave propagation in a dense gas confined between two plates. The inhomogeneous density near the solid wall is found to enhance the oscillation amplitude, while the presence of bulk viscosity causes stronger attenuation of the sound wave. This shows the importance of a kinetic model to reproduce density inhomogeneity and correct transport coefficients, including bulk viscosity.

Highly rarefied gas flows through a rough channel of finite length with small bumps appended to its surfaces are investigated, by varying the accommodation coefficient α in Maxwell’s diffuse-specular boundary condition, the characteristic size and position of the bumps, and the channel length. First, we study the influence of the surface bumps and consider the rarefied gas flow in a unit channel with periodic boundary conditions to remove the end effect. It is found that the surface bumps have a significant impact on the flow permeability. When α is very small (i.e., nearly specular reflection of gas molecules at the channel surface), the apparent gas permeability is dramatically reduced, even in the presence of small bumps, to a value that is almost comparable to the one when fully diffuse gas-surface scattering is assumed. This impact can be taken into account through an effective accommodation coefficient, i.e., the permeability of the rough channel is taken equivalently as that of a smooth channel without bumps but having gas-surface scattering under the effective accommodation coefficient. Second, we study the end effect by connecting a smooth channel of length L_0 to two huge gas reservoirs. It is found that (i) the end correction length is large at small α . Consequently, the mass flow rate barely reduces with increasing L_0 rather than scales down by a factor of 1/L_0 as predicted by the classical Knudsen diffusion theory; and (ii) the end correction is related to the channel’s aspect ratio. Finally, based on the effective accommodation coefficient and end correction, we explain the exotic flow enhancement in graphene angstrom-scale channels observed by Geim’s research group (Keerthi et al, Nature 558:420–424, 2018).

The Enskog-Vlasov equation provides a consistent description of the microscopic molecular interactions for real fluids based on the kinetic and mean-field theories. The fluid flows in nano-channels are investigated by the Bhatnagar-Gross-Krook (BGK) type Enskog-Vlasov model, which simplifies the complicated Enskog-Vlasov collision operator and enables large-scale engineering design simulations. The density distributions of real fluids are found to exhibit inhomogeneities across the nano-channel, particularly at large densities, as a direct consequence of the inhomogeneous force distributions caused by the real fluid effects including the fluid molecules' volume exclusion and the long-range molecular attraction. In contrast to the Navier-Stokes equation with the slip boundary condition, which fails to describe nano-scale flows due to the coexistence of confinement, non-equilibrium, and real fluid effects, the Enskog-Vlasov-BGK model is found to capture these effects accurately as confirmed by the corresponding molecular dynamics simulations for low and moderate fluid densities.

A thermodynamically consistent kinetic model is proposed for the non-equilibrium transport of confined van der Waals fluids, where the long-range molecular attraction is considered by a mean-field term in the transport equation, and the transport coefficients are tuned to match the experimental data. The equation of state of the van der Waals fluids can be obtained from an appropriate choice of the pair correlation function. By contrast, the modified Enskog theory predicts non-physical negative transport coefficients near the critical temperature and may not be able to recover the Boltzmann equation in the dilute limit. In addition, the shear viscosity and thermal conductivity are predicted more accurately by taking gas molecular attraction into account, while the softened Enskog formula for hard-sphere molecules performs better in predicting the bulk viscosity. The present kinetic model agrees with the Boltzmann model in the dilute limit and with the Navier–Stokes equations in the continuum limit, indicating its capability in modelling dilute-to-dense and continuum-to-non-equilibrium flows. The new model is examined thoroughly and validated by comparing it with the molecular dynamics simulation results. In contrast to the previous studies, our simulation results reveal the importance of molecular attraction even for high temperatures, which holds the molecules to the bulk while the hard-sphere model significantly overestimates the density near the wall. Because the long-range molecular attraction is considered appropriately in the present model, the velocity slip and temperature jump at the surface for the more realistic van der Waals fluids can be predicted accurately.

Sound wave propagation in rarefied flows of molecular gases confined in micro-channels is investigated numerically. We first validate the employed kinetic model against the experimental results and then systematically study the gas damping and surface force on the transducer as well as the resonance/anti-resonance in confined space. To quantify the impact of the finite relaxation rates of the translational and internal energies on wave propagation, we examine the roles of bulk viscosity and thermal conductivity in depth over a wide range of rarefactions and oscillation frequencies. It is found that the bulk viscosity only exerts influence on the pressure amplitude and its resonance frequency in the slip regime in high oscillations. In addition, the internal degree of freedom is frozen when the bulk viscosity of a molecular gas is large, resulting in the pressure amplitude of sound waves in the molecular gas being the same as in a monatomic gas. Meanwhile, the thermal conductivity has a limited influence on the pressure amplitude in all the simulated flows. In the case of the thermoacoustic wave, we prove that the Onsager–Casimir reciprocal relation also holds for molecular gases, i.e. the pressure deviation induced by the temperature variation is equal to the heat flux induced by the plate oscillation. Our findings enable an enhanced understanding of sound wave propagation in molecular gases, which may facilitate the design of nano-/micro-scale devices.

In rarefied gas flows, the spatial grid size could vary by several orders of magnitude in a single flow configuration (e.g., inside the Knudsen layer it is at the order of mean free path of gas molecules, while in the bulk region it is at a much larger hydrodynamic scale). Therefore, efficient implicit numerical method is urgently needed for time-dependent problems. However, the integro-differential nature of gas kinetic equations poses a grand challenge, as the gain part of the collision operator is non-invertible. Hence an iterative solver is required in each time step, which usu-ally takes a lot of iterations in the (near) continuum flow regime where the Knudsen number is small; worse still, the solution does not asymptotically preserve the fluid dynamic limit when the spatial cell size is not refined enough. Based on the general synthetic iteration scheme for steady-state solution of the Boltzmann equation, we pro-pose two numerical schemes to push the multiscale simulation of unsteady rarefied gas flows to a new boundary, that is, the numerical solution not only converges within dozens of iterations in each time step, but also asymptotically preserves the Navier-Stokes-Fourier limit in the continuum flow regime, when the spatial grid is coarse, and the time step is large (e.g., in simulating the extreme slow decay of two-dimensional Taylor vortex, the time step is even at the order of vortex decay time). The properties of fast convergence and asymptotic preserving of the proposed schemes are not only rigorously proven by the Fourier stability analysis for simplified gas kinetic models, but also demonstrated by several numerical examples for the gas kinetic models and the Boltzmann equation.

Abstract The Enskog-Vlasov equation provides a consistent description of the microscopic molecular interactions for real fluids based on the kinetic and mean-field theories. The present fluid flows in nano-channels are investigated by the Enskog-Vlasov-BGK model, which simplifies the complicated Enskog-Vlasov collision operator and enables large-scale engineering design simulations. The density distributions of real fluids are found to exhibit inhomogeneities across the nano-channel, particularly at large densities, as a direct consequence of the inhomogeneous force distributions caused by the real fluid effects including the fluid molecules' volume exclusion and the long-range molecular attraction. In contrast to the Navier-Stokes equation with the slip boundary condition, which fails to describe nano-scale flows due to the coexistence of confinement, non-equilibrium, and real fluid effects, the Enskog-Vlasov-BGK model is found to capture these effects accurately as confirmed by the corresponding molecular dynamics simulations for low and moderate fluid densities.

The phonon Boltzmann transport equation with dual relaxation times is often used to describe the heat conduction in semiconductor materials when the classical Fourier's law is no longer valid. For practical engineering designs, accurate and efficient numerical methods are highly demanded to solve the equation. At a large Knudsen number (i.e., the ratio of the phonon mean free path to a characteristic system length), steady-state solutions can be obtained via the conventional iterative scheme (CIS) within a few iterations. However, when the Knudsen number becomes small, i.e., when the phonon transport occurs in the diffusive or hydrodynamic regime, thousands of iterations are required to obtain converged results. In this work, a general synthetic iterative scheme (GSIS) is proposed to tackle the inefficiency of CIS. The key ingredient of the GSIS is that a set of macroscopic synthetic equations, which is exactly derived from the Boltzmann transport equation, is simultaneously solved with the kinetic equation to obtain the temperature and heat flux. During the iteration, the macroscopic quantities are used to evaluate the equilibria in the scattering terms of the kinetic equation, thus guiding the evolution of the phonon distribution function, while the distribution function, in turn, provides closures to the synthetic equations. The Fourier stability analysis is conducted to reveal the superiority of the GSIS over the CIS in terms of fast convergence in periodic systems. It is shown that the convergence rate of the GSIS can always be maintained under 0.2 so that only two iterations are required to reduce the iterative error by one order of magnitude. Numerical results in wall-bounded systems are presented to demonstrate further the efficiency of GSIS, where the CPU time is reduced by up to three orders of magnitude, especially in both the diffusive and hydrodynamic regimes where the Knudsen number is small.

The temperature jump problem in rarefied molecular (diatomic and polyatomic) gases is investigated based on a one-dimensional heat conduction problem. The gas dynamics is described by a kinetic model, which is capable of recovering the general temperature and thermal relaxation processes predicted by the Wang–Chang Uhlenbeck equation. Analytical formulations for the temperature jump coefficient subject to the classical Maxwell gas–surface interaction are derived via the Chapman–Enskog expansion. Numerically, the temperature jump coefficient and the Knudsen layer function are calculated by matching the kinetic solution to the Navier–Stokes prediction in the Knudsen layer. Results show that the temperature jump highly depends on the thermal relaxation processes: the values of the temperature jump coefficient and the Knudsen layer function are determined by the relative quantity of the translational thermal conductivity to the internal thermal conductivity; and a minimum temperature jump coefficient emerges when the translational Eucken factor is 4/3 times of the internal one. Due to the exclusion of the Knudsen layer effect, the analytical estimation of the temperature jump coefficient may possess large errors. A new formulation, which is a function of the internal degree of freedom, the Eucken factors, and the accommodation coefficient, is proposed based on the numerical results.

The general synthetic iterative scheme (GSIS) is extended to find the steady-state solution of the nonlinear gas kinetic equation, resolving the long-standing problems of slow convergence and requirement of ultra-fine grids in near-continuum flows. The key ingredient of GSIS is the tight coupling of gas kinetic and macroscopic synthetic equations, where the constitutive relations explicitly contain Newton's law of shear stress and Fourier's law of heat conduction. The higher-order constitutive relations describing rarefaction effects are calculated from the velocity distribution function; however, their constructions are simpler than our previous work (Su et al., 2020 [28]) for linearized gas kinetic equations. On the other hand, solutions of macroscopic synthetic equations are used to accelerate the evolution of gas kinetic equation at the next iteration step. A rigorous linear Fourier stability analysis of the present schemes in periodic system shows that the error decay rate of GSIS can be smaller than 0.5, which means that the deviation to steady-state solution can be reduced by 3 orders of magnitude in 10 iterations. Other important advantages of the GSIS are: (i) it does not rely on the specific form of Boltzmann collision operator, and (ii) it can be solved by sophisticated techniques in computational fluid dynamics, making it amenable to large scale engineering applications. In this paper, the efficiency and accuracy of GSIS are demonstrated by a number of canonical test cases in rarefied gas dynamics, covering different flow regimes.