A Partial Entropic Stabilization of Lattice Boltzmann MHD

TL;DR

This paper introduces a partial entropic stabilization method for lattice Boltzmann magnetohydrodynamics, improving simulation stability and accuracy in 2D MHD problems like Kelvin-Helmholtz instability.

Contribution

It extends the entropic lattice Boltzmann method to MHD by applying it only to the scalar distribution, enabling better handling of magnetic field reversals.

Findings

Achieved very good simulation results for MHD instabilities.

Benchmarking shows improved stability over previous models.

Successfully applied to complex 2D MHD phenomena.

Abstract

The entropic lattice Boltzmann algorithm of Karlin et. al. is partially extended to magnetohydrodynamics, based on the Dellar model of introducing a vector distribution for the magnetic field. This entropic ansatz is now applied only to the scalar particle distribution function so as to permit the many problems entailing magnetic field reversal. A 9-bit lattice is employed for both particle and magnetic distributions for our two dimensional simulations. The entropic ansatz is benchmarked against our earlier multiple relaxation lattice-Boltzmann model for the Kelvin-Helmholtz instability in a magnetized jet. Other two dimensional simulations are performed and compared to results determined by more standard direct algorithms: in particular the switch over between the Kelvin-Helmholtz/tearing mode instability of Chen et. al., and the generalized Orszag-Tang vortex model of Biskamp-Welter. Very good results are achieved.

Equations

$$\left({\partial }_t+{\partial }_{\gamma }{c}_{\gamma i}\right)f_i=\sum _j{X}_{ij}^{{}^{\prime }}\left(f_j^{\left(\mathrm{eq}\right)}-f_j\right)$$

$$\left({\partial }_t+{\partial }_{\gamma }{c}_{\gamma k}\right){\vec{g}}_k=Y^{{}^{\prime }}\left({\vec{g}}_k^{(\mathrm{eq})}-{\vec{g}}_k\right)$$

$$\displaystyle\begin{array}[]{ccccc}\sum_{i}f_{i}=\rho&,&\sum_{i}f_{i}\vec{c}_{i}=\rho\vec{u}&\mathrm{and}&\sum_{k}\vec{g}_{k}=\vec{B}\end{array}$$

$$f_i^{\left(\mathrm{eq}\right)}=w_i\rho \left[1+3\left({\vec{c}}_i\cdot \vec{u}\right)+\frac{9}{2}{\left({\vec{c}}_i\cdot \vec{u}\right)}^2-\frac{3}{2}{\vec{u}}^2\right]+\frac{9}{2}{w}_i\left[\frac{1}{2}{\vec{B}}^2{\vec{c}}_i^2-{\left(\vec{B}\cdot {\vec{c}}_i\right)}^2\right],i=0,..,8$$

$${\vec{g}}_k^{(\mathrm{eq})}=w_k\left[\vec{B}+3\left\{\left({\vec{c}}_k\cdot \vec{u}\right)\vec{B}-\left({\vec{c}}_k\cdot \vec{B}\right)\vec{u}\right\}\right],k=0,..,8$$

$$\begin{array}[]{ccc}M_{i}=\sum_{j=0}^{8}\mathrm{T}_{ij}f_{j}&,&\vec{N}_{k}=\sum_{q=0}^{8}\mathrm{T}_{kq}\vec{g}_{q}\end{array}$$

$$\mathrm{T}=\left(\begin{array}[]{c}\boldsymbol{1}\\ c_{x}\\ c_{y}\\ c_{x}c_{y}\\ c^{2}_{x}\\ c^{2}_{y}\\ c_{x}^{2}c_{y}\\ c_{x}c_{y}^{2}\\ c_{x}^{2}c_{y}^{2}\end{array}\right)=\left(\begin{array}[]{ccccccccc}1&1&1&1&1&1&1&1&1\\ 0&1&0&-1&0&1&-1&-1&1\\ 0&0&1&0&-1&1&1&-1&-1\\ 0&0&0&0&0&1&-1&1&-1\\ 0&1&0&1&0&1&1&1&1\\ 0&0&1&0&1&1&1&1&1\\ 0&0&0&0&0&1&1&-1&-1\\ 0&0&0&0&0&1&-1&-1&1\\ 0&0&0&0&0&1&1&1&1\end{array}\right)$$

$$\begin{array}[]{cccc}c_{x}=\left\{0,1,0,-1,0,1,-1,-1,1\right\}&,&c_{y}=\{0,0,1,0,-1,1,1,-1,-1\}&.\end{array}$$

$$\displaystyle\begin{split}\begin{array}[]{lll}{M}^{\left(\mathrm{eq}\right)}_{0}=M_{0}=\rho&{M}^{\left(\mathrm{eq}\right)}_{1}=M_{1}=\rho u_{x}&{M}^{\left(\mathrm{eq}\right)}_{2}=M_{2}=\rho u_{y}\end{array}\\ \begin{array}[]{ll}{M}^{\left(\mathrm{eq}\right)}_{3}=\rho u_{x}u_{y}-B_{x}B_{y}&{M}^{\left(\mathrm{eq}\right)}_{4}=\frac{1}{6}\left(6\rho u_{x}^{2}+2\rho-3\left(B^{2}_{x}-B^{2}_{y}\right)\right)\\ {M}^{\left(\mathrm{eq}\right)}_{5}=\frac{1}{6}\left(6\rho u_{y}^{2}+2\rho+3\left(B^{2}_{x}-B^{2}_{y}\right)\right)&{M}^{\left(\mathrm{eq}\right)}_{6}=\frac{1}{3}\rho u_{y}\\ {M}^{\left(\mathrm{eq}\right)}_{7}=\frac{1}{3}\rho u_{x}&{M}^{\left(\mathrm{eq}\right)}_{8}=\frac{1}{9}\rho\left(1+3u_{x}^{2}+3u_{y}^{2}\right)\\ \end{array}\end{split}$$

$$\begin{array}[]{lll}{N}^{\left(\mathrm{eq}\right)}_{\alpha 0}=N_{\alpha 0}=B_{\alpha}&{N}^{\left(\mathrm{eq}\right)}_{\alpha 1}=u_{x}B_{\alpha}-u_{\alpha}B_{x}&{N}^{\left(\mathrm{eq}\right)}_{\alpha 2}=u_{y}B_{\alpha}-u_{\alpha}B_{y}\\ {N}^{\left(\mathrm{eq}\right)}_{\alpha 3}=0&{N}^{\left(\mathrm{eq}\right)}_{\alpha 4}=\frac{B_{\alpha}}{3}&{N}^{\left(\mathrm{eq}\right)}_{\alpha 5}=\frac{B_{\alpha}}{3}\\ {N}^{\left(\mathrm{eq}\right)}_{\alpha 6}=\frac{1}{3}\left(u_{y}B_{\alpha}-u_{\alpha}B_{y}\right)&{N}^{\left(\mathrm{eq}\right)}_{\alpha 7}=\frac{1}{3}\left(u_{x}B_{\alpha}-u_{\alpha}B_{x}\right)&{N}^{\left(\mathrm{eq}\right)}_{\alpha 8}=\frac{B_{\alpha}}{9},\qquad\alpha=x,y\end{array}$$

$$f_i=k_i+s_i+h_i,i=\mathrm{0..8}$$

$$k_i=\sum _{j=0}^8\sum _{m=0}^2{\mathrm{T}}_{im}^{-1}{\mathrm{T}}_{mj}{f}_j,i=\mathrm{0..8}$$

$$s_i\in \left\{d,d\cup t,d\cup q,d\cup t\cup q\right\}$$

$$s_i=\sum _{j=0}^8\sum _{m=3}^5{\mathrm{T}}_{im}^{-1}{\mathrm{T}}_{mj}{f}_j,i=\mathrm{0..8.}$$

$$h_i=\sum _{j=0}^8\sum _{m=6}^8{\mathrm{T}}_{im}^{-1}{\mathrm{T}}_{mj}{f}_j,i=\mathrm{0..8.}$$

$$f_i^{{}^{\prime }}\equiv f_i\left(t+1\right)=f_i+2\beta \left(f_i^{\left(\mathrm{eq}\right)}-f_i\right)$$

$$\text{to}{f}_i^{{}^{\prime }}=f_i-2\beta \Delta s_i-\beta \gamma \Delta h_i$$

$$S\left[f\right]=-\sum _i{f}_i\ln (\frac{f_i}{w_i}).$$

$$\sum _i\Delta h_i\ln (1+\frac{\left(1-\beta \gamma \right)\Delta h_i-\left(2\beta -1\right)\Delta s_i}{f_i^{\left(\mathrm{eq}\right)}})=0$$

$${\gamma }^{\ast }=\frac{1}{\beta }-\left(2-\frac{1}{\beta }\right)\frac{⟨\Delta s|\Delta h⟩}{⟨\Delta h|\Delta h⟩}$$

$$\text{where the inner product}⟨A|B⟩=\sum _i\frac{A_i{B}_i}{f_i^{\left(\mathrm{eq}\right)}}.$$

$$\vec{u}(t=0)=U_0{\text{sech}}^2(x)\hat{y},\vec{B}(t=0)=B_0\hat{y}$$

$$u_y(x,t=0)=-U_0\tanh (x),B_y(x,t=0)=B_0\tanh (x).$$

$$\vec{u}(x,y,t=0)=U_0\left[\sin (y+0.5)\hat{x}-\sin (x+1.4)\hat{y}\right]$$

$$\vec{B}(x,y,t=0)=B_0\left[\sin (y+4.1)\hat{x}-2\sin (2x+2.3)\hat{y}\right]$$

Full text

I Introduction

The lattice Boltzmann (LB) algorithm has proven to be an extremely interesting method for the solution of Navier-Stokes Succi (2001) flows because of its simplicity, extreme parallelizability and accuracy. Even though it is technically second order accurate it appears more comparable in accuracy to the pseudo-spectral computational methods. One of the major constraints on LB is that it is prone to numerical instability in certain parameter regimes: there is no inherent mechanism to enforce the LB distribution function to remain non-negative in time, particularly in strong turbulence simulations.

This has led to a generalization of the simple single-relaxation-time (SRT) LB collision operator to its multiple-relaxation-time (MRT) cousin Khirevich et al. (2015); d’Humières et al. (2002); Premnath and Abraham (2007). With these extra degrees of freedom one can achieve greater numerical stability, but the choice of these extra parameters is problem dependent, not known a prior, and can influence the fluid viscosity coefficient (and thus the associated Reynolds numbers). An alternate approach to achieving numerical stability is through an entropic principle and a discrete H-theorem Bösch et al. (2015a); Karlin et al. (2014); Chikatamarla and Karlin (2013); Boghosian et al. (2001, 2003, 2004); Keating et al. (2007); Vahala et al. (2009); Flint et al. (2016, 2015); Dellar (2001, 2003, 2009); Premnath et al. (2009); Bösch et al. (2015b); Pattison et al. (2008); Riley et al. (2008). In some respects, the entropic approach can be viewed as an optimal subset of MRT algorithms in which emphasis is placed on an algebraically determined entropy stabilizing parameter that is not directly dependent on the MRT collisional rates and which does not affect the fluid viscosity.

Here we partially extended these entropic ideas to magnetohydrodynamics (MHD). LB was first extended to MHD by Succi et. al. Succi et al. (1991) following on the heels of a cellular automata approach of Montgomery and Doolen Montgomery and Doolen (1987). An important breakthrough to a first principle LB-MHD model was by Dellar Dellar (2002) who introduced a vector kinetic equation for the magnetic distribution function. In conventional LB for Navier-Stokes, the zeroth moment of the scalar distribution function yields the density $\rho$, while the fluid velocity $\vec{u}$, is retrieved from the first moment. In the Dellar model for LB-MHD, the zeroth moment of the vector magnetic distribution function yields the magnetic field $\vec{B}$ itself. Not only does this permit moment closure at a lower level than in Navier-Stokes but it yields a consistent discrete approximation to $\divergence\vec{B}=0$ to machine round-off error Dellar (2002). However there is a significant difference between Navier-Stokes and MHD: in MHD there are many very important applications where there is a reversal in the magnetic field $\vec{B}$. Hence a first principle extension of an entropic principle, which relies on maximization of a concave function, cannot be applied to the magnetic field vector distribution function. Nevertheless, we shall apply the entropic stabilization scheme of Karlin et. al. Bösch et al. (2015a) to the scalar distribution function, and no such constraint on the vector distribution function. While some may insist that we are thus not forming any entropic stabilization as such, we will find that our partial entropic stabilization (as we shall call it) does permit numerically stable simulations at arbitrary small viscosity. This can not be achieved even in MRT LB-modeling, which only has static relaxation rates. We also shall show, from several different simulations, that there seems to be an increased stabilization in the LB-MHD algorithm due to this partial entropy constraint on the scalar particle distribution function. This can be attributed to the effect of the partial entropic parameter on the magnetic field because of the magnetic field coupling that exists within the velocity momentum equation/Navier-Stokes equation

In Section II, we present a moment-based representation for LB-MHD, while our partial entropic algorithm is outlined in Sec. III. In Sec. IV we present some two dimensional (2D) simulation results of our partial entropic LB-MHD algorithm: magnetic reconnection in the Kelvin-Helmholtz and the magnetic tearing instability as well as on the Biskamp-Welter profile for 2D MHD. We have concentrated on 2D MHD because of its much lower computational costs as compared to 3D LB-MHD. This is appropriate since 2D and 3D MHD there is a direct cascade of energy to small spatial scales - unlike 2D Navier-Stokes turbulence which exhibits an inverse cascade of energy to larger and larger spatial scales. Our partial entropic stabilization algorithm is readily extended to 3D.

II Moment Basis Representation for Multiple Relaxation Model for LB-MHD

There are quite a few MRT extensions Flint et al. (2016); Dellar (2009); Pattison et al. (2008); Dellar (2011) of the original SRT LB-MHD model of Dellar Dellar (2002). However, for simplicity, we shall work with only an SRT model for the vector magnetic field distribution $\vec{g}_{k}$, and an MRT model for the scalar distribution function $f_{i}$, where the subscripts denote the velocity streaming directions

[TABLE]

with the moments

[TABLE]

It is convenient to employ the summation convention only over the Greek indices which give the vector nature of the fields ($\gamma=1,2$ for 2D), while the Roman indices run over the corresponding (kinetic) lattice vectors $\vec{c}_{i},i=0..8$ for the 9-bit model in 2D (see Fig. 1). Summation over the Roman indices will always be made explicit. $X^{{}^{\prime}}_{ij}$ is the MRT collision operator for the evolution of $f_{i}$ while $Y^{{}^{\prime}}$ is the SRT for the evolution of $\vec{g}_{k}$. The MHD viscosity and resistivity transport coefficients are determined from these kinetic relaxation rates.

It is well known that the minimal LB representation of MHD equations on a square lattice is a 9-bit velocity streaming for $f_{i}$ and just 5-bit streaming for $\vec{g}_{k}$. This is because $\vec{u}$ is defined from the first moment of $f_{i}$ while $\vec{B}$ is defined as the zeroth moment of $\vec{g}_{k}$. It is convenient (and helpful for numerical stability) to employ the 9-bit streaming model for both kinetic equations. To recover the MHD equations in the Chapman-Enskog limit of the (discrete) kinetic equations, we take the well-known choice of relaxation distribution functions $f{{}_{i}}^{(eq)}$ and $\vec{g}_{k}^{\,(eq)}$

[TABLE]

The LB-MHD equations are typically solved by an operator-splitting method with time advancement coming from streaming and collisional relaxation. The excellent parallelization of the LB-MHD algorithm is now apparent: the discrete kinetic equations are solved using the streaming operation which is a simple shift of the data from one lattice point to another, while the collision step is a purely local operation with its evaluation requires only data from only that grid point. The streaming requires MPI only when the shift has data moving from one processor domain to another – and this can also be well parallelized. What makes LB so attractive is that the computationally difficult nonlinear convective derivatives $\vec{u}\cdot\nabla\vec{u}$, $\vec{u}\cdot\nabla\vec{B}$, $\vec{B}\cdot\nabla\vec{u}$ and $\vec{B}\cdot\nabla\vec{B}$ are replaced at the lattice kinetic level by simple linear advection and polynomial nonlinearities in the relaxation distributions.

In MRT-LB it is natural to perform the collisional relaxation in moment space (because of the local conservation of mass and momentum constraints) and the streaming in the distribution space $f_{i},\vec{g}_{k}$. There is a 1-1 map between these spaces. For the moment basis it is obvious to include the conservation moments (the zeroth and first moments of the $f_{i}$ and the zeroth moment of $\vec{g}_{k}$), while the remaining higher moments are somewhat arbitrary Bouzidi et al. (2001); Lallemand and Luo (2000). In particular, we consider the same constant $9\crossproduct 9$ $\mathrm{T}$- matrix that connects the scalar distributions $(f_{i},i=0..8)$ to their moments $(M_{i},i=0..8)$ as for the vector magnetic distributions $(\vec{g}_{k},k=0..8)$ with their moments $(\vec{N}_{k},k=0..8)$

[TABLE]

with

[TABLE]

where the Cartesian components of the corresponding $9$-dimensional lattice vectors are just

[TABLE]

For the scalar distributions, the 1st row of the $\mathrm{T}$-matrix is just the conservation of density while the 2nd and 3rd rows are just the conservation of momentum (2D). For the vector magnetic distributions the 1st row of the $\mathrm{T}$-matrix is the only collisional invariant.

With this moment basis, the MRT collisional relaxation rate tensor $X^{{}^{\prime}}_{ij}$ is diagonalized with the $T-$ matrix as a similarity transformation. It is convenient to denote this diagonal matrix with elements $X_{i}\delta_{ij}$. In the $D2Q9$ phase space, the relaxation rate $X_{j}$ is associated with the corresponding moment $M_{j}$, $j=0..8$. Similarly for the magnetic distributions in SRT, there is just a single collisional relaxation rate for each magnetic moment $\vec{N}_{k}$, and this will be denoted by $Y$.

In particular, the equilibrium moments can be written in terms of the conserved moments:

[TABLE]

[TABLE]

III Entropic Method and its Partial Extension to MHD

The Karlin group Bösch et al. (2015a); Karlin et al. (2014) introduces the entropic procedure for Navier-Stokes flows by separating the scalar lattice Boltzmann distribution into various moment-related groups. In particular,

[TABLE]

where the $k_{i}$ distributions correspond to those distributions with conserved moments, the $s_{i}$ distributions correspond to the stress/shear moments, and finally the $h_{i}$ distributions correspond to the remaining higher order moments. Thus for the $k_{i}$ distributions

[TABLE]

with the $m-$summation running from $m=0,1,2$ since there are 3 conserved moments. Similarly for $s_{i}$ and $h_{i}$. The $s_{i}$ distributions corresponding to the stress/shear moments will come from the set

[TABLE]

where $d$ is the deviatoric stress, $t$ is the trace of the stress tensor, and $q$ represents the third order moments. Here we choose, for simplicity, the moment contributions $s_{i}$ to be $d\cup t$ so that

[TABLE]

Moments 3, 4, and 5 are each second order moments in the $D2Q9$ model, and thus represent the second order quantities $d\cup t$. The moment contributions to $h_{i}$ are then all the remaining moments that do not contribute to either $k_{i}$ or $s_{i}$. Thus

[TABLE]

Karlin et. al. Bösch et al. (2015a); Karlin et al. (2014) now consider the entropy of the post-collisional state, and introduce a parameter $\gamma$ which yields an extremal to this entropy function. In MRT only some of the relaxation rates affect the transport coefficient under Chapman-Enskog expansions Pope (2000). The transport coefficient in Navier-Stokes simulations is first affected by the stress related distributions $(s_{i})$. The tunable parameter $\gamma$ is introduced to replace the relaxation rates for the higher order moment effects arising from the $(h_{i})$ distributions. In particular, one moves from the standard post-collisional distributions

[TABLE]

[TABLE]

where $\beta$ is related to the kinematic viscosity as $\nu=\frac{1}{6}\left(\frac{1}{\beta}-1\right)$ and $\Delta s_{i}=s_{i}-{s}^{\left(\mathrm{eq}\right)}_{i}$, $\Delta h_{i}=h_{i}-{h}^{\left(\mathrm{eq}\right)}_{i}$, while for the conserved moments $\Delta k_{i}=k_{i}-{k}^{\left(\mathrm{eq}\right)}_{i}=0$.

In order to maximize the entropy $S\left[f\right]$

[TABLE]

one now writes the entropy in terms of the post-collisional state and the $\gamma$ parameter. The critical point of the entropy Bösch et al. (2015a); Karlin et al. (2014) determines the tunable parameter $\gamma$ from

[TABLE]

This is a rather computationally expensive root-finding procedure having to be done at every point of the grid and at every time step. Karlin et. al Bösch et al. (2015a); Karlin et al. (2014) noted that if one invokes the simple small argument expansion $\log(1+x)=x+...$ one can then readily determine the entropic factor algebraically. The parameter determined algebraically is denoted by $\gamma^{*}$ :

[TABLE]

[TABLE]

On substituting $\gamma^{*}$ back into the new post-collisional state (Eq. 18) a maximal entropy state has been determined for Navier-Stokes flows. The Karlin group successfully tested this approximation for the tunable parameter $\gamma^{*}(\vec{x},t)$ in various simulations of 2D and 3D Navier-Stokes Bösch et al. (2015a); Karlin et al. (2014). One thus sees that this emtropic algorithm is a subset of MRT - but it has a dynamic entropic parameter determined at every lattice point and every time step algebraically for entropic stabilization as opposed to the static relaxation times for typical MRT models.

Clearly, this analysis does not simply carry over to LB-MHD with possible non-positive vector magnetic distributions. Hence we make the ansatz for our partial entropic algorithm that the entropic parameter in LB-MHD is still determined by Eq. (21) for the corresponding LB-MHD $\Delta h$ and $\Delta s$. The validity of our ansatz will now be tested against various 2D MHD simulations.

Summarizing, our partial entropic LB-MHD algorithm consists of the following steps (c.f., Karlin et. al. Bösch et al. (2015a):

Compute the conserved moments ($\rho$,$\mathbf{u}$,$\mathbf{B}$) (Eq. 7, 10, 11) 2. 2.

Evaluate the equilibria $\left(f_{i}^{(\mathrm{eq})}\left(\rho,\mathbf{u},\mathbf{B}\right),\vec{g}_{k}^{(\mathrm{eq})}\left(\rho,\mathbf{u},\mathbf{B}\right)\right)$ (Eq. 5) 3. 3.

Compute $s$ and $s^{(\mathrm{eq})}$ (Eq. 13, 14) 4. 4.

Compute $\Delta s_{i}=s_{i}-{s}^{\left(\mathrm{eq}\right)}_{i}$ 5. 5.

Compute $\Delta h_{i}=h_{i}-{h}^{\left(\mathrm{eq}\right)}_{i}=f_{i}-{f}^{\left(\mathrm{eq}\right)}_{i}-\Delta s_{i}$ 6. 6.

Evaluate $\gamma^{*}$ (Eq. 21) 7. 7.

Relax (Collide): $f_{i}^{{}^{\prime}}$ (Eq. 18), and corresponding $\vec{g}_{k}^{{}^{\prime}}$.

Standard LB-MHD is recovered for entropy parameter: $\gamma(\vec{x},t)=const.=2$. As mentioned earlier, there is no attempt made to find a corresponding maximal entropy state for the magnetic distribution function since the magnetic field in most problems of interest undergoes field reversal. (e.g., in magnetic field reconnection..). However the effect of working with the maximal entropy state for the particle distribution function will have direct effects on the evolution of the magnetic field due to the coupling of the $\vec{B}$-field in the relaxation distribution function ${f}^{\left(\mathrm{eq}\right)}$ as well as the coupling of the fluid velocity $\vec{u}$ in ${\vec{g}}^{\left(\mathrm{eq}\right)}$.

IV Partially Entropic LB-MHD Simulations

We first have benchmarked our partially entropic LB-MHD code against our earlier (totally non-entropic) MRT LB-MHD simulations of a Kelvin-Helmholtz jet instability in a magnetic field Flint et al. (2016). Here we show the physics recovered by the variations in the partially entropic parameter $\gamma^{*}$ and its variations away from the MRT value of $\gamma^{*}(\vec{x},t)\equiv 2.0$ for sufficiently weak axial $\vec{B}$ that the jet is unstable. Some runs were then performed to examine the increased numerical stability in the parameter regime of the mean velocity $\vec{u}$ and magnetic field $\vec{B}$ due to the partially entropic algorithm. Following this we consider the interplay between Kelvin-Helmholtz instability and the tearing mode instability and qualitatively compare our results to that of Chen et. al, Chen et al. (1997). Finally we qualitatively compare our simulations with the Biskamp-Welter profile.

IV.1 Magnetized Kelvin-Hemholtz Jet Instability

We now consider the partially entropic-LB-MHD algorithm for the breakdown of a Kelvin-Helmholtz jet in a weak magnetic field. In our simulations, the initial parameters are so chosen that there is a direct cascade of energy to small spatial scales (indicating the existence of a magnetic field) but the magnetic field is sufficient weak so as not to stabilize the jet, Fig. 2

[TABLE]

The evolution of the vorticity, $\omega$, the current, $j$ as well as the entropic stabilization parameter $\gamma^{*}$ for this 2D jet is plotted in Fig. 3. With the (dimensionless) choice of $B_{0}=0.005U_{0}$, the jet breaks into a Kelvin-Helmholtz vortex street ($t\leq 266k$). There is then further symmetry breaking as the vortex street is broken up, leading to vortex-vortex reconnection (ala 2D Navier-Stokes turbulence), as well as the generation of small scales eddies (characteristic of 2D MHD) for $t>266k$. One notices that the partial entropy parameter $\gamma^{*}$ in fig. 3 deviates from the ordinary LB-MHD value of $\gamma^{*}(\vec{x},t)=2$ wherever there are significant number of small eddies. These are regions of steep gradients and it is in these regions where the partial entropic stabilization of the simulation occurs. It is important to note that this partial entropy stabilization is occurring from local information at each lattice site. This is reminiscent of LB where gradients can be computed from local moments of perturbed distributions: e.g., in large eddy simulation modelings in the Smagorinsky model, the mean velocity gradients are determined from simple local moments. For stronger magnetic fields, the jet will be stabilized and is of little interest for our partial entropic-LB-MHD model, Flint et al. (2016). A spectral plot of the total energy of the Kelvin-Helmholtz simulation at $t=500$k is presented in fig. 4 with a slope of $k^{-\frac{5}{3}}$. This spectral plot corresponds to the final timestep in Fig. 3.

IV.1.1 Stability Improvements from the Entropic Algorithm

Some numerical stability boundaries were investigated between ordinary LB-MHD and our partial entropic-LB-MHD with the $\gamma^{*}$ parameter on a grid of $1024^{2}$ for the Kelvin-Helmholtz jet. We found that the partial entropic-LB-MHD algorithm permitted a maximal mean velocity $\vec{U}_{0,max}$ to be increased by a factor of $2$ in a purely Navier-Stokes turbulence simulation (i.e., no $\vec{B}$-field) while the velocity maximum could be increased by a factor of $8$ when there was a strong stabilizing $\vec{B}$-field. As regards the magnetic field (at fixed $\vec{U}_{0}$), the partial entropic-LB-MHD algorithm permitted an increase by a factor of $2$ in $\vec{B}_{0}$. In the partial entropic-LB-MHD algorithm the viscosity could become arbitrary small, while ordinary LB-MHD the minimum stable viscosity was $10^{-5}$ when $\vec{B}_{0}=0$, and $10^{-2}$ when there was a strong stabilizing $\vec{B}_{0}$. No substantial stability limits were found on the achievable minimum resistivity.

It should be stressed that the computational overhead of computing this entropic parameter $\gamma^{*}$ is quite small, primarily because it is determined algebraically from local information only.

IV.2 Chen et. al. Profile

Chen et. al. Chen et al. (1997) has considered the linear and nonlinear evolution of Kelvin-Helmholtz (velocity shear) vs. the tearing mode (magnetic shear) instabilities in 2D compressible MHD. Their closure includes an evolution equation for the enthalphy as well as various resistivity profiles using standard CFD techniques. Their initial profiles are

[TABLE]

Thus our comparisons can only be qualitative, and we only consider the Chen et. al. Chen et al. (1997) simulations when they keep their resistivity constant. Typically, when the velocity is below the Alfven speed, it stabilizes the tearing mode and so reduces the reconnection rate. However, if the velocity is above the Alfven speed the Kelvin-Helmholtz instability sets in. In our first partial entropic LB-MHD simulation, we consider super-Alfven velocity shear flow and the Kelvin-Helmholtz induced magnetic islands due to reconnection in Fig. 5. In Fig. 5a we show the simulation results of case 13 in Chen et. al. for the magnetic field lines and compare them to those arising from our partial entropic LB-MHD model for resistivity $\eta=0.001$, Fig. 5b.

For the case of no initial shear, large magnetic islands are formed. A corresponding snapshot is given of the magnetic field lines from the case 5 simulation in Chen et. al., Fig. 6a, and from our entropic LB-MHD model, Fig. 6b. In Fig. 7 we show the partial entropic LB-MHD evolution of the magnetic field lines for this initial zero velocity shear flow parameter set of Fig. 6. It seems for the case considered here, the enthalpy equation in Chen Chen et al. (1997) does not play a significant role.

IV.3 Biskamp-Welter Profile

We now consider the model of Biskamp and Welter Biskamp and Welter (1989) for decaying 2D MHD turbulence, using their initial profiles

[TABLE]

(These are a generalization of the canonical Orszag-Tang vortex). A snapshot of the current lines are shown in Fig. 8 and compared with those from the Biskamp-Welter simulation.

In Fig. 9 we plot the corresponding 2D entropy parameter $\gamma^{*}(x,y)$ at this time snapshot. The lattice points at which $\gamma^{*}(x,y)\neq 2$ correspond to points where there are effects of in our partial entropic LB-MHD algorithm. The energy dissipation rate for this Biskamp-Welter case is shown in fig. 10. This can be compared with figure 20 in Biskamp and Welter (1989).

V Conclusion

We have extended the Karlin Bösch et al. (2015a); Karlin et al. (2014) entropic Navier-Stokes algorithm to LB-MHD and tested the ensuing model on 3 different problems: velocity shear flows exhibiting Kelvin-Helmholtz and/or tearing instability, a generalized Orszag-Tang vortex and magnetized jet instability. We considered the D2Q9 model for both the particle and vector magnetic field distributions. The partial entropy algorithm is applied only to the particle distributions while in using a vector distribution for the magnetic field one must allow for magnetic field reversals. Hence we do not have a fully entropic LB-MHD model. The algorithm clearly extends immediately to 3D, but because of the much greater computational costs we have restricted our simulations to 2D for we can still capture turbulence effects of the generation of small scale motions since in 2D MHD energy cascades to small scales. We have found good agreement with the CFD simulations of Chen et. al. and Biskamp and Welter. The partial entropic algorithm permits much larger ranges of velocity and magnetic field amplitudes than could be found in standard LB-MHD algorithms. This greater numerical stability is achieved at a quite small increase in computational costs since Karlin et. al. have determined a simple algebraic approximation to the full entropic parameter. This approximation is then carried over as an ansatz for our 2D LB-MHD model. Moreover the extreme parallelization of this partial entropic LB-MHD algorithm is retained since this algebraic entropic parameter $\gamma^{*}$ is determined purely from local information at each lattice site. The accuracy of the under-resolved Navier-Stokes simulations of Bösch et. al. Bösch et al. (2015b) portend that this new (partial) entropy method could be a possible subgrid model in itself. In some sense, this is the spirit behind our pushing the magnitude of $U_{0}$ and $B_{0}$. We are not trying to claim rigorous error bounds on various equilibria representations. This partial entropic LB-MHD algorithm is a subset of MRT models in which there is now a dynamical relaxation rate determined for quasi-stabilization of the fluid flow by a well-defined procedure as opposed to the standard static MRT relaxation rates.

VI Acknowledgments

This work was supported by an NSF and AFOSR grant. The computations were performed on DoD Supercomputer Topaz.