# Dealing with spurious oscillations in particle tracking methods

I work on modelling high intensity discharge xenon-filled lamps. The model governing the discharge is quite complex and sadly includes fluid dynamics. After some time, I managed to implement a finite-difference particle-tracking lagrangian fluid dynamics for axisymmetric 1D discharge:

$$\frac{\partial}{\partial t} \left(\frac{1}{\rho}\right) = \frac{\partial(rv)}{\partial s} \\ \frac{\partial v}{\partial t} = - r \frac{\partial p}{\partial s} \\ \frac{\partial r}{\partial t} = v \\ \frac{\partial}{\partial t} \left( \mathcal{E} + \frac{v^2}{2}\right) = - \frac{\partial}{\partial s}(prv) - \frac{\partial(rW)}{\partial s} + \frac{f}{\rho} \\ W = -\lambda \frac{\partial T}{\partial r} \\ p = p(\rho, T) \\ \frac{\partial\mathcal{E}}{\partial T} = c_v(\rho, T)$$

Computational cells are organised so that lagrangian particle speeds and eulerian coordinates are defined on cell-interfaces and energy equation parameters (temperature, density, pressure, thermal capacity, thermal conductivity) are defined on cell centers (image below). Domain boundaries coincide with cell interfaces. .

Right-side boundary conditions are dirichlet for speed $$v(r)|_{r=R} = 0$$, keeping the scheme homogenous. For energy balance equation it gets a bit trickier as there is no boundary node to specify dirichlet condition at, so I intruduce fictious node in the depth of the wall material and interpolate boundary condition at cell interface:

$$\frac{\hat{T}_{last} + \hat{T}_{last + 1}}{2} = T_{boundary} \\->\\ \hat{T}_{last + 1} = 2T_{boundary} - \hat{T}_{last}$$

For computational efficiency reason I solve the equations separately. There's iterational process on energy equation; first I solve energy equation with implicit O(t^2)O(h^2) scheme, then I advance gas dynamics part on interpolated temperature profile with explicit scheme with time step several orders smaller to refine coefficients for energy equation. This is dictated by the extremely high cost to recalculate energy source function as there's radiative transfer over complex spectra involved.

Model works great for pulses with current of about 500A and duration greater than 50e-6 seconds, however once I go into kA ranges waves rebounding from the lamp wall tend to spuriously oscillate and destroy the model upon converging in the center ( ).

I tried some less radical methods of dealing with these oscillations like introducing artificial viscosity, however I found that my implementation of them greatly reduced stability due to loss of conservation.

With great dread I look forward to rewriting everything from scratch in terms of TVD scheme or with nested grids or with something equally time-consuming, and my question is: what are some less intrusive approaches to eliminating spurious oscillations that I overlooked?