Numerical Resolution of a Three Temperature Plasma Model

Abstract

This paper is devoted to the numerical approximation of a three temperature plasma model: one for the ions, one for the electrons and one for the radiation (photons). A reformulation of the model is proposed that allows to build a convex combination-based scheme that unconditionally satisfies a maximum principle, at each sub-iteration of the non-linear iterative process. This yields a very robust scheme that can handle stiff source terms. In addition, the methodology is extended to include the contribution of a radiative flux (Rosseland diffusion approximation) and electronic and ionic conductivities (Spitzer–Härm diffusion approximation). Several numerical results are carried out to demonstrate the interest of the numerical approach.

Appendices

Appendix A: Stability of the State $T_i = T_e = T_r$.

Property 14

All states such that $T_i = T_e = T_r$ are equilibrium states (in the 0D case (no spatial variation) case with no source term).

Proof

We only consider here the simplified the 0D case (no spatial variation) (no spatial dependence) with no source terms. Setting

$$\begin{aligned} X= \begin{pmatrix} E_r \\ E_e \\ E_i \end{pmatrix}, \quad F(X) = \begin{pmatrix} c \sigma _P (\phi _e - \phi _r) &{} &{} \\ c \sigma _P (\phi _r - \phi _e) &{} + &{} \kappa (T_i - T_e) \\ &{} &{} \kappa (T_e - T_i) \end{pmatrix}, \end{aligned}$$

the model (1) writes under the following ordinary differential set of equations

$$\begin{aligned} {\dot{X}} = F(X). \end{aligned}$$

The unique fix point $X_0$ such that $F(X_0) = 0$ corresponds to the case of equal temperatures $T_r=T_e=T_i$. Now, considering the case of a small perturbation $\delta X$ around the equilibrium state $X_0$, one obtains

$$\begin{aligned} \dot{(X_0 + \delta X)} = \dot{\delta X} = F(X_0 + \delta X) \simeq F(X_0) + \partial _X F|_{X_0} \delta X = \partial _X F|_{X_0} \delta X. \end{aligned}$$

In this linear case, the following analytical expression of $\delta X$ is obtained

$$\begin{aligned} \delta X (t) = e^{\partial _X F|_{X_0} t} \delta X(0). \end{aligned}$$

The study of the eigenvalues of the matrix s $\partial _X F|_{X_0}$ shows that the matrix displays two negative eigenvalues and a third one which is zero. Therefore, $X_0$ is a stable equilibrium. $\square$

Appendix B: Practical Computation of the Terms $\beta _\alpha ^{k}$ and $\delta _{ie}^k$

Fig. 15

Representation of $\beta _\alpha ^{k}$ a function of $C_k$ using a constant reconstruction

Fig. 16

Representation of $\beta _\alpha ^{k}$ a function of $C_k$ using a linear reconstruction

The evaluations of the expressions (14) may be sensitive if the denominators become close to zero. However this difficulty is only numerical. Indeed, performing a Taylor expansion on $\beta _\alpha ^{k}$ in $T_\alpha ^k = T_\alpha ^n$ gives

$$\begin{aligned} \lim _{|T_\alpha ^k - T_\alpha ^n| \rightarrow 0} \beta _\alpha ^{k} = \frac{4 a (T_\alpha ^n)^3}{\rho C_{v,\alpha }}. \end{aligned}$$

(31)

Similarly, a Taylor expansion on $\delta _{ie}^k$ in $\phi _i^k = \phi _e^k$ leads to

$$\begin{aligned} \lim _{|T_i^k-T_e^k| \rightarrow 0} \delta _{ie}^k = \frac{1}{4 a (T_e^k)^3}. \end{aligned}$$

(32)

Coming back to the computation of $\beta _\alpha ^{k}$. The solution we consider consist in comparing the quantity $\displaystyle C_k = |T_\alpha ^k - T_\alpha ^n| / (T_\alpha ^k + T_\alpha ^n)$ with a threshold value $\varepsilon$. If $C_k$ is found larger than $\varepsilon$ then (14) is used otherwise one takes the expressions (31). In Fig. 15a, the value of $\beta _\alpha ^{k}$ is displayed as function of the value of $C_k$ for a large value of $\varepsilon$ (left) and a small value of $\varepsilon$ (right). Even for small $\varepsilon$ one observes a discontinuity between the limit value and the definition of $\beta _\alpha ^{k}$. In order to avoid this, one considers a linear reconstruction instead of a constant state. The limit value and the value in $\varepsilon$ are simply connected using a linear approximation. In Fig. 16a, it is observed that the discontinuity vanishes even using large value of $\varepsilon$.

The same procedure is applied for coefficient $\delta _{ie}^k$.

Remark

In practice, it has been observed that using this linear reconstruction procedure, the choice of the value of $\varepsilon$ does not have any impacts on the numerical results.

Appendix C: Temporal Discretisation of $\sigma _P$, $\kappa$ and $\sigma _R$

In the case of an explicit discretisation, they are fixed at time $t^n$ during the iterative process. In the implicit case, the quantity are computed at iteration $k+1$ while these coefficients are taken at step k. A third discretisation is studied, taking the half sum between the quantity at iteration k and time $t^n$. It is observed that taking small time steps all the discretisation give the same correct results. However, when using very large time step some differences appear. For example in Fig. 17, reconsidering the test case Marshak2A with large time step (only 100 time step and a coarse grid of 100 cells), it is observed that the explicit method is by far the less accurate while the implicit give the best accuracy. However, when looking at the computational time in Table 3, it is observed that the explicit method is the fastest while the implicit is the slowest.

Fig. 17

Temperature comparisons between different time discretisations of $\sigma _P$, $\kappa$ and $\sigma _R$ for the test case Marshak2A using the gradient conjugate method

Table 3 Computational time (in second) for different time discretisations

Remark

The results displayed in Fig. 17 are obtained using the conjugate gradient method. We mention that the results are very similar when using the Jacobi and Gauss–Seidel numerical schemes.