Electron thermal conductivity misbehaving

[trevilo]

The third order electron thermal conductivity model is misbehaving. Specifically, it is producing negative values for some combinations of inputs.

[trevilo]

The bad behavior was first identified in a test case that @shaering put together during development for PR #297 (see test/inputs/input.radDecay.ini, here). The pushed version uses

argon_minimal/third_order_thermal_conductivity = false

to avoid the problem but changing to

argon_minimal/third_order_thermal_conductivity = true

will lead to negative total thermal conductivity, which causes the temperature solve to fail.

The problem was traced to the electron contribution to the thermal conductivity becoming negative.

[trevilo]

By hacking test/test_argon_minimal.cpp, we can plot the electron thermal conductivity on a grid in (T, \alpha), where \alpha is the "ionization degree" as used in that function, which is n_e / (n_tot - n_e).

Doing so gives the attached plot, which clearly shows that there appear to be combinations where the third order model is singular (on a curve between approx 7900K and 8800K). For \alpha below the singularity, \kappa_e is negative.

kappa_e_badness.pdf

Assuming positive electron mole fraction (which is imposed in this test), negative \kappa_e is only possible if the factor (L11 - L12 * L12 / L22) that appears in the 3rd order model (see e.g., ArgonMixtureTransport::computeThirdOrderElectronThermalConductivity here) is negative, and when this factor is zero we get the singularity noted above.

So... what could cause this to go negative?

[trevilo]

Since the problem appears at low electron mole fraction, I decided to start the investigation by looking at the contribution of the electron-neutral collisions. This part involves e-Ar collision integrals that were curve fit by @dreamer2368 (see comments here). I haven't been able to locate the original fitting code, so instead I put together a python script to evaluate the fits against approximating the collision integral with quadrature (just trapezoid rule) using the BSR e-Ar momentum transfer cross section.

The script and the necessary input cross section are attached here for traceability. (The eAr.txt file is a python script, but github won't let me upload a .py file here).
eAr.txt
eAr_momentum.txt

The script generates the following results:
eAr_collision_integrals.pdf
eAr_determinant.pdf

The first plot shows the TPS fits compared to the quadrature results, which are not terrible, at least at first glance. But, computing the term that goes into \kappa_e (i.e., L11*L22 - L12*L12) from these fits leads to a very bad approximation of that term, as shown in the second plot. Specifically, compare the blue (quadrature) to orange (TPS fit)... the extra lines show what you get from either the "low T" or "high T" fits, because the original TPS fits use different forms for different temperature ranges. The oscillation is induced by the fact that the switch between forms occurs at different T for different collision integrals.

This is the source of the bad behavior in \kappa_e: marginal approximations of the the collision integrals leading to very bad behavior in this derived quantity.

[trevilo]

So, to fix the problem, we need to update these fits, which turns out to be very simple. Rather than using a high and low T form, I decided to try the following form:

Q = \sum_k c_k (\log(T))^k

for k= -1, 0, ..., 7. The least-squares fitting is implemented in the script attached to the previous comment, and it produces the following results:
eAr_collision_integrals_refit.pdf
eAr_determinant_refit.pdf

The collision integral fits are somewhat better than before, but much more importantly, the determinant result L11*L22 - L12*L12 is well-behaved. There is some non-negligible error at high T which we could improve. But, in practice, we never see T this high in the torch, so that is very low priority to me.