Dust Module

Warning: Beta feature

The dedicated dust dynamics module has not yet been exercised in a published science application with Quokka and should currently be treated as beta.

This module primarily implements two components: the dust transport term and the dust drag source term.

Equations for Gas-Dust System

$\begin{align} \frac{\partial \rho_{\mathrm{g}}}{\partial t} + \nabla \cdot (\rho_{\mathrm{g}} \mathbf{v}_{\mathrm{g}}) &= 0, \\ \frac{\partial (\rho_{\mathrm{g}} \mathbf{v}_{\mathrm{g}})}{\partial t} + \nabla \cdot (\rho_{\mathrm{g}} \mathbf{v}_{\mathrm{g}} \otimes \mathbf{v}_{\mathrm{g}} + P_{\mathrm{g}} \mathbf{I}) &= \sum_{n=1}^{N} \rho_{\mathrm{d},n} \frac{\mathbf{v}_{\mathrm{d},n} - \mathbf{v}_{\mathrm{g}}}{T_{\mathrm{s},n}}, \\ \frac{\partial E_{\mathrm{g}}}{\partial t} + \nabla \cdot \left[(E_{\mathrm{g}} + P_{\mathrm{g}}) \mathbf{v}_{\mathrm{g}}\right] &= \sum_{n=1}^{N} \rho_{\mathrm{d},n} \frac{\mathbf{v}_{\mathrm{d},n} - \mathbf{v}_{\mathrm{g}}}{T_{\mathrm{s},n}} \cdot \mathbf{v}_{\mathrm{g}} + \omega \sum_{n=1}^{N} \rho_{\mathrm{d},n} \frac{(\mathbf{v}_{\mathrm{d},n} - \mathbf{v}_{\mathrm{g}})^{2}}{T_{\mathrm{s},n}}, \\ \frac{\partial \rho_{\mathrm{d},n}}{\partial t} + \nabla \cdot (\rho_{\mathrm{d},n} \mathbf{v}_{\mathrm{d},n}) &= 0, \\ \frac{\partial (\rho_{\mathrm{d},n} \mathbf{v}_{\mathrm{d},n})}{\partial t} + \nabla \cdot (\rho_{\mathrm{d},n} \mathbf{v}_{\mathrm{d},n} \otimes \mathbf{v}_{\mathrm{d},n}) &= \rho_{\mathrm{d},n} \frac{\mathbf{v}_{\mathrm{g}} - \mathbf{v}_{\mathrm{d},n}}{T_{\mathrm{s},n}}, \end{align}$

where $\omega$ controls the level of frictional heating, with $\omega = 0$ turning it off and $\omega = 1$ depositing all dissipation into the gas.

Variable Storage

The dust cell-centred conserved variables ($\rho_{\mathrm{d}}$, $\rho_{\mathrm{d}}\vec{v_{\mathrm{d}}}$) are added to MultiFab.

Reconstruction and Riemann Solver

Dust reconstruction is performed together with gas using the same method. The Riemann Solver used is as follows:

In one dimension along the x-direction, given the left/right states $W_d^{L/R}$, one can provide the Riemann flux for conserved variables as follows. The density flux reads (Huang & Bai 2022):

$\begin{align*} F^{\text a}_x(\rho_d) = \begin{cases} \rho_d^L v_{d,x}^L & \text{if } v_{d,x}^L > 0, \, v_{d,x}^R \ge 0, \\ \rho_d^R v_{d,x}^R & \text{if } v_{d,x}^L \le 0, \, v_{d,x}^R < 0, \\ \rho_d^L v_{d,x}^L + \rho_d^R v_{d,x}^R & \text{if } v_{d,x}^L > 0, \, v_{d,x}^R < 0, \\ 0 & \text{else}. \end{cases} \end{align*}$

Similar expressions hold for the momentum flux for all directions.

This is implemented in src/dust/dustRiemannSolver.hpp and called in DustSystem::ComputeDustFluxes to compute the dust advection flux.

Time Integrator

A Strang-split method (Tedeschi-Prades et al. 2025) is used to integrate the dust-gas system. The Strang-split update can be expressed as:

$\mathbf{u}^{n+1} = \mathcal{D}_{\Delta t/2} \mathcal{H}_{\Delta t} \mathcal{D}_{\Delta t/2} \mathbf{u}^n$

where $\mathcal{D}$ is the dust-gas drag operator and $\mathcal{H}$ is the hydrodynamics operator (including both gas and dust transport). The hydrodynamics operator $\mathcal{H}$ is handled using the explicit RK2 scheme. The drag operator $\mathcal{D}$ is implemented in src/dust/DustDrag.hpp and called in QuokkaSimulation::addStrangSplitSourcesWithBuiltin via DustDrag::computeDustDrag.

Optional Picard iteration for dust–gas drag

Users may optionally enable Picard iteration for the dust–gas drag operator $\mathcal{D}$. When the stopping time depends on the gas or dust velocity, enabling iteration is required to maintain an implicit dust drag update. See Runtime parameters for details.

User-defined dust stopping time

For a given problem, users must define a problem-specific dust stopping time by implementing the DustDrag::ComputeReciprocalStoppingTime function (note that this function should return the reciprocal of the stopping time). An example can be found in the src/problems/DustDamping test.

Also, users can directly use the dust stopping time calculation helper DustDrag::ComputeReciprocalStoppingTimeKwok to compute the physical dust stopping time, following Kwok (1975) with an optional supersonic correction. The stopping time of dust $t_{\mathrm{s}}$ is given by:

$t_{\mathrm{s}} = \frac{\sqrt{\pi \gamma}}{2\sqrt{2}} \frac{a \rho_{\mathrm{gr}}}{\rho_{\mathrm{g}} c_{\mathrm{s}}} \times \begin{cases} \left( 1 + \dfrac{9\pi\gamma}{128} \left| \dfrac{\mathbf{v}_{\mathrm{d}} - \mathbf{v}_{\mathrm{g}}}{c_{\mathrm{s}}} \right|^2 \right)^{-1/2}, & \text{if supersonic correction is enabled,} \\[1.5em] 1, & \text{if supersonic correction is disabled.} \end{cases}$

When $\gamma=1$, this expression reduces exactly to the isothermal $t_s$. An example of its usage can be found in the src/problems/DustDampingCorrection test.

CFL Condition for Dust

For the dust-gas coupled system with N dust species, we use the following CFL condition:

$\Delta t_{\mathrm{CFL}} = C_{\mathrm{CFL}} \cdot \min_{\mathrm{cells}} \left( \frac{\Delta x}{\max\left( |v_{\mathrm{g}}| + c_s, \max_{n=1}^{N} |v_{\mathrm{d},n}|+c_s \right)} \right).$

Runtime Controls

The following input parameters tune the dust module and are documented in more detail in Runtime parameters:

enable_iter_stoptime – switch of iterative dust stopping time calculation.
omega – controls the level of frictional heating.
print_iteration_counts - switch to turn on/off printing of dust drag iteration counts for debugging.
dust.density_floor - the minimum dust density value allowed in the simulation.

Keyboard shortcuts

Quokka