# 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?

Edit 1: added more information about energy balance equation, schema, and boundary conditions used.

• Have you considered that this could be due to your boundary condition at the right side of the domain? If so, what have you tried? – Spencer Bryngelson Jun 2 '19 at 1:12
• @SpencerBryngelson boundary conditions for gas dynamics part are dirichlet, namely fixed speed across boundary (0), and as far as I understood those require no additional work as they keep scheme homogenous, is my understanding incomplete? For the energy equation it is a bit tricky as energy nodes are offset by half step from velocity nodes so as to not deal with axial coordinate singularity, so to fix wall temperature I have to interpolate it with the help of fictitious node half step past the boundary, but my understanding is that the rest of the scheme is kept homogenous. – Dantragof Jun 2 '19 at 6:14
• I forgot to mention that 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. – Dantragof Jun 2 '19 at 6:38