Skip to content

Neutral-Neutral

Neutral-Neutral Collisions

Neutral-neutral collisions in GUERNICA are currently modeled with a BGK operator. This operator is enabled by including neutral_neutral under operators in [components]. In the current implementation, it relaxes the neutral distribution function toward an equilibrium distribution \(f_M\) according to

\[ \left.\frac{\partial f_n}{\partial t}\right|_{\mathrm{nn}} = -\nu \left(f_n - f_M\right), \]

where \(f_n\) is the current neutral distribution function, \(f_M\) is the BGK equilibrium distribution, and \(\nu\) is the neutral-neutral collision frequency.

The important point is that \(f_M\) is not formed by simply evaluating the continuum Maxwellian using the local neutral density, flow velocity, and temperature. In the continuum, that would be the natural choice, since a Maxwellian with the same mass, momentum, and energy is the correct BGK equilibrium. However, in GUERNICA the velocity space is discretized with the Discrete Velocity Method (DVM), so the moments are computed as discrete sums over the velocity grid. If one naively evaluates the continuum Maxwellian on that discrete grid, the resulting discrete distribution will generally not reproduce the exact discrete mass, momentum, and energy of the current state. As a result, the BGK operator would no longer preserve those quantities at the discrete level.

To avoid that problem, GUERNICA computes \(f_M\) so that it explicitly matches the discrete conserved moments used by the solver. Rather than defining \(f_M\) directly from \((\rho,u,T)\), the code writes it in exponential form as

\[ f_M(v_x,v_y,v_z) = \exp\left( a_0 + a_{1x} v_x + a_{1y} v_y + a_{1z} v_z + a_2 |v|^2 \right), \]

where the five coefficients \((a_0,a_{1x},a_{1y},a_{1z},a_2)\) are chosen so that the discrete moments of \(f_M\) match the target neutral moments:

  • density \(\rho\),
  • momentum components \(\rho u_x\), \(\rho u_y\), \(\rho u_z\),
  • and energy \(E\).

So in practice, the code solves the nonlinear system

\[ \sum_{\mathbf v} w_{\mathbf v} f_M(\mathbf v) = \rho, \]
\[ \sum_{\mathbf v} w_{\mathbf v} v_x f_M(\mathbf v) = \rho u_x, \qquad \sum_{\mathbf v} w_{\mathbf v} v_y f_M(\mathbf v) = \rho u_y, \qquad \sum_{\mathbf v} w_{\mathbf v} v_z f_M(\mathbf v) = \rho u_z, \]
\[ \sum_{\mathbf v} w_{\mathbf v} |\mathbf v|^2 f_M(\mathbf v) = E, \]

where the sums are taken over the discrete velocity grid with the DVM quadrature weights. This is the sense in which the BGK equilibrium in GUERNICA is a discrete Maxwellian: it has the exponential Maxwellian form, but its coefficients are determined so that the discrete conserved quantities are matched exactly rather than only in the continuum limit.

The solver computes these coefficients in two stages. First, it builds an initial guess from the corresponding analytic Maxwellian using the local neutral moments. In the code this gives

\[ a_0 = \log\left(\frac{\rho}{(2\pi T)^{3/2}}\right) - \frac{|u|^2}{2T}, \qquad a_{1x} = \frac{u_x}{T}, \qquad a_{1y} = \frac{u_y}{T}, \qquad a_{1z} = \frac{u_z}{T}, \qquad a_2 = -\frac{1}{2T}, \]

with \(T\) inferred from the local energy and flow speed. This serves only as an initial guess. The actual equilibrium used by the BGK operator is then obtained by solving the full discrete moment-matching system.

That nonlinear system is solved in the code with Newton’s method. At each update, the BGK operator:

  1. computes the discrete moment residuals between the current trial \(f_M\) and the target moments,
  2. assembles the corresponding \(5\times5\) Jacobian,
  3. solves for the Newton correction,
  4. updates the coefficient vector \((a_0,a_{1x},a_{1y},a_{1z},a_2)\),
  5. and repeats until convergence.

This is why the BGK operator in GUERNICA preserves the intended moments much more faithfully than a naive pointwise Maxwellian evaluation on the DVM grid would. The equilibrium is still Maxwellian in form, but it is specifically the Maxwellian whose discrete mass, momentum, and energy agree with the current neutral state.

To turn neutral-neutral collisions on in the input file, include the operator in [components] and add a [neutral_neutral] section. A minimal setup looks like this:

[components]
species = n
operators = neutral_neutral

[neutral_neutral]
model = bgk
nu = 1.0

Here:

  • model = bgk selects the BGK neutral-neutral collision model.
  • nu sets the relaxation frequency controlling how quickly \(f_n\) is driven toward the discrete equilibrium \(f_M\). Larger nu gives stronger relaxation.