Frequency-Domain GSIS for Rarefied Gas Flows

• Frequency-Domain GSIS is an advanced numerical framework that couples kinetic equations with high-order synthetic macroscopic corrections to simulate oscillatory rarefied gas flows efficiently in MEMS applications.
• The method employs a frequency-domain linearized Boltzmann–Shakhov formulation with iterative updates to achieve asymptotic preservation and super-convergence, drastically reducing computational iterations.
• GSIS delivers up to three orders of magnitude speed-up over conventional methods by tightly coupling mesoscopic kinetic updates with macroscopic moment corrections, enhancing simulation efficiency in complex microflows.

The frequency-domain General Synthetic Iterative Scheme (GSIS) is an advanced numerical framework for the efficient simulation of oscillatory rarefied gas flows, with direct relevance to microelectromechanical systems (MEMS) and applications requiring high-fidelity kinetic modeling. The formulation targets the periodic steady-state regime and achieves significant acceleration over conventional kinetic methods by tightly coupling the solution of the mesoscopic kinetic equation with high-order macroscopic (synthetic) equations. Through this iterative synergy, GSIS offers super-convergence, an asymptotic-preserving property, and up to three orders of magnitude speed-up for near-continuum regimes over traditional approaches (Li et al., 24 Jan 2026).

1. Frequency-Domain Kinetic Formulation

The central mathematical object of GSIS is the frequency-domain formulation of the linearized Boltzmann–Shakhov equation, which describes the response to small-amplitude wall oscillations $U_w(t) = \mathrm{Re}\{U_0 e^{i\omega t}\}$, with $U_0 \ll v_m$ (most probable molecular speed). By employing the time-harmonic ansatz

$$f(x, v, t) = \mathrm{Re}\left\{ \hat{f}(x, v) e^{i\omega t} \right\},$$

the kinetic equation reduces to

$$i\omega\,\hat{f} + v \cdot \nabla_x \hat{f} = \mathcal{C}[\hat{f}],$$

where $\mathcal{C}$ denotes the collision operator.

After nondimensionalization (characteristic length $L$, speed $v_m$, pressure $p_0$, viscosity $\mu$), the essential controlling parameters become the Strouhal number $S = \omega L / v_m$ and the rarefaction parameter $\delta_{rp} = p_0 L / (\mu v_m)$. Introducing $h$ via $\hat{f} = f_{eq}[1 + \xi h e^{i\omega t}]$ and $\xi = U_0 / v_m$, the frequency-domain linearized equation for the perturbation $h$ reads:

$$i S h + v \cdot \nabla h = \delta_{rp}\left[ \rho + 2 v \cdot u + (v^2 - \frac{3}{2})\tau + \frac{4}{15}(v^2 - \frac{5}{2})(v \cdot q) - h \right],$$

with $\rho, u, \tau, q$ representing the complex amplitudes of the perturbed density, velocity, temperature, and heat flux, respectively (Li et al., 24 Jan 2026).

2. High-Order Macroscopic Synthetic Equation Structure

Macroscopic variables are defined by velocity-space moments of $h$ (all quantities implicit in $e^{i\omega t}$):

• $\rho = \int h f_{eq} dv$
• $u = \int v h f_{eq} dv$
• $$\tau = \int \left(\frac{2}{3}v^2 - 1\right) h f_{eq} dv$$
• $$\Pi = \int 2\left(vv - \frac{v^2}{3} I\right) h f_{eq} dv$$
• $q = \int (v^2 - \frac{5}{2}) v h f_{eq} dv$

Multiplying the kinetic equation by relevant moment functions and integrating yields an unclosed macroscopic system:

$$i S \rho + \nabla \cdot u = 0 \ 2 i S u + \nabla \rho + \nabla \tau + \nabla \cdot \Pi = 0 \ \frac{3}{2} i S \tau + \nabla \cdot u + \nabla \cdot q = 0$$

The closure is constructed via decomposition:

$$\Pi = \Pi_{NSF} + HoT_\Pi,\quad q = q_{NSF} + HoT_q,$$

where the Navier–Stokes–Fourier (NSF) terms

$$\Pi_{NSF} = -\frac{1}{\delta_{rp}}\left[ \nabla u + (\nabla u)^T - \frac{2}{3}(\nabla \cdot u)I \right],\quad q_{NSF} = -\frac{15}{8 \delta_{rp}} \nabla \tau$$

are augmented by high-order terms (HoT) extracted from velocity moments of the kinetic residual at the half-iteration. The GSIS thus yields a closed macroscopic synthetic system, with the high-order corrections ensuring kinetic accuracy beyond the first-order constitutive laws (Li et al., 24 Jan 2026).

3. Iterative Coupled Algorithm

The GSIS iteration at outer step $n$ involves:

1. Kinetic half-step (CIS update):

Solve $$(iS + \delta_{rp}) h^{n+\frac{1}{2}} + v \cdot \nabla h^{n+\frac{1}{2}} = \delta_{rp}[\rho^n + 2 v \cdot u^n + (v^2 - \frac{3}{2})\tau^n + \frac{4}{15}(v^2 - \frac{5}{2})(v \cdot q^n)]$$.

1. Moment update: Compute $$\{ \rho^{n+\frac{1}{2}}, u^{n+\frac{1}{2}}, \tau^{n+\frac{1}{2}}, \Pi^{n+\frac{1}{2}}, q^{n+\frac{1}{2}} \}$$ via velocity-space quadrature.
2. High-order term evaluation: Use above formulas involving the kinetic residual at $n+\frac{1}{2}$.
3. Synthetic system solve: Solve the high-order closed macroscopic system for $[\rho^{n+1}, u^{n+1}, \tau^{n+1}]$.
4. Distribution correction: $$h^{n+1} = h^{n+\frac{1}{2}} + [\rho^{n+1} - \rho^{n+\frac{1}{2}}] + 2 v \cdot [u^{n+1} - u^{n+\frac{1}{2}}] + (v^2 - 3/2)[\tau^{n+1} - \tau^{n+\frac{1}{2}}]$$.

The iteration proceeds until the macroscopic residual norm drops below a user-specified tolerance, e.g., $10^{-5}$ (Li et al., 24 Jan 2026).

4. Asymptotic-Preserving Properties via Chapman–Enskog Analysis

A Chapman–Enskog expansion

$$h = h_0 + \delta_{rp}^{-1} h_1 + O(\delta_{rp}^{-2})$$

yields $$h_0 = \rho + 2 v \cdot u + (v^2 - 3/2)\tau + (4/15)(v^2 - 5/2)(v \cdot q)$$ as the leading term and $h_1 = -v \cdot \nabla h_0$ at first order. Substitution demonstrates $HoT_\Pi, HoT_q = O(\delta_{rp}^{-2})$, reducing the high-order corrections to negligible levels as $\delta_{rp} \to \infty$. As a result, $\Pi = \Pi_{NSF} + O(Kn^2)$ and $q = q_{NSF} + O(Kn^2)$, ensuring that the GSIS is asymptotic preserving (AP). On coarse spatial meshes, the system correctly limits to the Navier–Stokes–Fourier equations in the continuum regime, without requiring mesh resolution proportional to the Knudsen number (Li et al., 24 Jan 2026).

5. Stability, Super-Convergence, and Error Decay

A Fourier stability analysis considers iteration error modes,

$$Y^{n+1}(x,v) = h^{n+1} - h^n = e^n \bar{Y}(v) e^{i \theta \cdot x}, \quad \Phi^{n+1}(x) = M^{n+1} - M^n = e^{n+1}\alpha e^{i\theta \cdot x},$$

where $|\theta| = 1$. For the Conventional Iterative Scheme (CIS), the decay factor $e$ approaches $1-O(1/\delta_{rp}^2)$ in the near-continuum regime, producing "false convergence"—the final error scales as $\delta_{rp}^2 \epsilon$ even for small residual $\epsilon$. For GSIS, the error decay factor satisfies $\rho(G) = O(1/\delta_{rp}^2)$ uniformly in $S$, ensuring error reduction by several orders of magnitude per iteration—so-called super-convergence (Li et al., 24 Jan 2026).

Empirical stability data show that while CIS iteration counts diverge ($\sim 10^3-10^4$) for large $\delta_{rp}$, GSIS consistently converges in $O(10)$ iterations across all rarefaction levels and Strouhal numbers.

6. Numerical Performance and Benchmark Results

The GSIS's computational efficacy was validated on two benchmark problems:

Geometry	$N_{cell}$	$\delta_{rp}$	$S$	CIS iters	GSIS iters	CIS CPU(s)	GSIS CPU(s)
Shear-driven flow between eccentric cylinders	$6\!-\!12\times10^3$	1	1.0	64	24	20	14
		10	1.0	414	25	134	14
		100	0.1	10,047	29	5,349	25
		1,000	0.001	37,570	27	23,787	32
Squeeze-film damping (oscillating cantilever)	$2.8\times10^4$	1	1.0	58	27	67	66
		10	1.0	1,091	29	1,209	67
		100	1.0	10,725	30	11,765	80
		1,000	0.001	—	27	—	60

For all tested regimes, GSIS retains $20$–$30$ iteration convergence throughout the rarefaction range, while CIS iteration counts increase by orders of magnitude. CPU times reflect up to a $10^3$-fold speed-up for GSIS as $\delta_{rp}$ increases (near-continuum regime) (Li et al., 24 Jan 2026).

7. Summary and Importance

The frequency-domain GSIS provides a deterministic kinetic solver coupled with a high-order macroscopic (synthetic) equation, utilizing high-order stress/heat-flux extraction from a velocity-space half-step. Its principal attributes are:

• Super-convergence in iteration counts ($O(1)$ as Kn → 0)
• Asymptotic preservation (Navier–Stokes–Fourier recovery on coarse grids)
• Robust stability across all frequencies and rarefaction
• Computational acceleration up to three orders of magnitude over conventional split transport–collision iterations

This scheme enables tractable simulation of oscillatory rarefied gas flow phenomena in complex geometries and high-dimensional spaces relevant in MEMS and microfluidics, with practical impact on the computational modeling of thermally or mechanically forced microflows (Li et al., 24 Jan 2026).