A large eddy lattice Boltzmann simulation of magnetohydrodynamic turbulence

TL;DR

This paper presents large eddy simulations of magnetohydrodynamic turbulence using an extended lattice Boltzmann model, demonstrating good agreement with direct numerical simulations for the Kelvin-Helmholtz instability.

Contribution

It introduces a novel LES approach for LB-MHD that incorporates non-commuting perturbations and physical constraints for subgrid modeling.

Findings

Simulations agree well with direct numerical results.

The method effectively captures magnetized Kelvin-Helmholtz instability.

Extension of LB-MHD to LES with non-commuting perturbations.

Abstract

Large eddy simulations (LES) of a lattice Boltzmann magnetohydrodynamic (LB-MHD) model are performed for the unstable magnetized Kelvin-Helmholtz jet instability. This algorithm is an extension of Ansumali et. al. (2004) to MHD in which one performs first an expansion in the filter width on the kinetic equations followed by the usual low Knudsen number expansion. These two perturbation operations do not commute. Closure is achieved by invoking the physical constraint that subgrid effects occur at transport time scales. The simulations are in very good agreement with direct numerical simulations.

Equations

$$\bar{\mathbf{u}}\left(\vec{r},\Delta \right)={\int }_{-\infty }^{\infty }\mathbf{u}\left(\vec{r}{}^{\prime }-\vec{r}\right)G\left(\vec{r}{}^{\prime },\Delta \right)d\vec{r}{}^{\prime }.$$

$${\partial }_t\text{vbox}\text{hrule}height=0.5pt\rho +\nabla \cdot \text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}=0,\nabla \cdot \overline{\mathbf{B}}=0$$

$$\begin{array}{r}\text{hfil}\begin{array}{rl}{\partial }_t\left(\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\right)& +\nabla \cdot \left(\frac{\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}}{\bar{\rho }}\right)=-\nabla \bar{p}+\nabla \cdot \left(\overline{\mathbf{B}}\overline{\mathbf{B}}\right)-\frac{1}{2}\nabla \left(\overline{\mathbf{B}}\cdot \overline{\mathbf{B}}\right)+\left(\xi +\frac{1}{3}\nu \right)\nabla \left(\nabla \cdot \text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\right)+\nu {\nabla }^2\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\\ & -\nabla \cdot \left\{\frac{6\nu }{6\nu +1}\frac{{\Delta }^2}{12\bar{\rho }}\left[\left({\partial }_{\beta }\left(\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\right)\right)\left({\partial }_{\beta }\left(\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\right)\right)-\frac{{\partial }_{\beta }\bar{p}}{\bar{p}}\left(\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\left({\partial }_{\beta }\left(\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\right)\right)+\left({\partial }_{\beta }\left(\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\right)\right)\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}-\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\frac{{\partial }_{\beta }\bar{p}}{\bar{p}}\right)\right]\right\}\\ & -\nabla \left\{\left(\frac{s_4}{4}+\frac{s_7}{20}-\frac{3s_8}{10}\right)\frac{{\Delta }^2}{12\bar{\rho }}\left[\left({\partial }_{\beta }\left(\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\right)\right)\cdot \left({\partial }_{\beta }\left(\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\right)\right)-\frac{{\partial }_{\beta }\bar{p}}{\bar{p}}\left(2\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\cdot \left({\partial }_{\beta }\left(\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\right)\right)-\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\cdot \text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\frac{{\partial }_{\beta }\bar{p}}{\bar{p}}\right)\right]\right\}\\ & -\frac{6\nu }{6\nu +1}\frac{{\Delta }^2}{12}\left\{\frac{1}{2}\nabla \left[\left({\partial }_{\beta }\overline{\mathbf{B}}\right)\cdot \left({\partial }_{\beta }\overline{\mathbf{B}}\right)\right]-\nabla \cdot \left[\left({\partial }_{\beta }\overline{\mathbf{B}}\right)\left({\partial }_{\beta }\overline{\mathbf{B}}\right)\right]\right\},\end{array}\end{array}$$

$$\begin{array}{r}\text{hfil}\begin{array}{rl}{\partial }_t\overline{\mathbf{B}}& =\nabla \text{crossproduct}\left(\frac{\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\text{crossproduct}\overline{\mathbf{B}}}{\bar{\rho }}\right)+\eta {\nabla }^2\overline{\mathbf{B}}+\nabla \text{crossproduct}[\frac{{\Delta }^2}{12\bar{\rho }}\frac{6\eta }{6\eta +1}\{\left({\partial }_{\beta }\left(\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\right)\right)\text{crossproduct}\left({\partial }_{\beta }\overline{\mathbf{B}}\right)\\ & -\frac{{\partial }_{\beta }\bar{p}}{\bar{p}}\left(\left(\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\right)\text{crossproduct}\left({\partial }_{\beta }\overline{\mathbf{B}}\right)+\left({\partial }_{\beta }\left(\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\right)\right)\text{crossproduct}\overline{\mathbf{B}}-\frac{\left({\partial }_{\beta }\bar{p}\right)}{\bar{p}}\left(\text{vbox}\text{hrule}height=0.5pt\rho \mathbf{u}\right)\text{crossproduct}\overline{\mathbf{B}}\right)\}].\end{array}\end{array}$$

$$\nu =\frac{1}{3s_3}-\frac{1}{6}=\frac{1}{3s_4}-\frac{1}{6}$$

$$\xi =-\frac{1}{9}-\frac{1}{9s_4}-\frac{1}{15s_7}+\frac{2}{5s_8}$$

$$\eta =\frac{1}{3s_m}-\frac{1}{6}$$

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

$$\left({\partial }_t+{\partial }_{\gamma }{c}_{\gamma i}\right){\vec{g}}_i=s_m^{{}^{\prime }}\left({\vec{g}}_i^{(\mathrm{eq})}-{\vec{g}}_i\right),i=\mathrm{0...8}$$

$$f_i^{(\mathrm{eq})}=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}}_i^{(\mathrm{eq})}=w_i^{{}^{\prime }}\left[\vec{B}+3\left\{\left({\vec{c}}_i\cdot \vec{u}\right)\vec{B}-\left({\vec{c}}_i\cdot \vec{B}\right)\vec{u}\right\}\right],i=0,..,8$$

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

$$\displaystyle\mathrm{T}=\left(\begin{array}[]{c}\boldsymbol{1}\\ c_{x}\\ c_{y}\\ c_{x}c_{y}\\ c^{2}_{x}-c^{2}_{y}\\ 3c_{x}c^{2}_{y}-2c_{x}\\ 3c_{y}c^{2}_{x}-2c_{y}\\ 4\cdot\boldsymbol{1}-9\left(c^{2}_{x}+c^{2}_{y}-2c^{2}_{x}c^{2}_{y}\right)\\ 4\cdot\boldsymbol{1}-4\left(c^{2}_{x}+c^{2}_{y}\right)+3c^{2}_{x}c^{2}_{y}\end{array}\right)=\left(\begin{array}[]{rrrrrrrrr}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&-1&1&-1&0&0&0&0\\ 0&-2&0&2&0&1&-1&-1&1\\ 0&0&-2&0&2&1&1&-1&-1\\ 4&-5&-5&-5&-5&4&4&4&4\\ 4&0&0&0&0&-1&-1&-1&-1\end{array}\right)$$

$$\displaystyle\mathrm{T_{m}}=\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}[]{rrrrrrrrr}1&\phantom{-}1&\phantom{-}1&1&1&\phantom{-}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}[]{ccc}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{array}[]{lll}M_{0}^{(eq)}=M_{0}=\rho&M_{3}^{(eq)}=\frac{\rho u_{x}\rho u_{y}}{\rho}-B_{x}B_{y}&M_{6}^{(eq)}=-\rho u_{y}\\ M_{1}^{(eq)}=M_{1}=\rho u_{x}&M_{4}^{(eq)}=\frac{\left(\rho u_{x}\right)^{2}-\left(\rho u_{y}\right)^{2}}{\rho}-B_{x}^{2}+B_{y}^{2}&M_{7}^{(eq)}=-3\frac{\left(\rho u_{x}\right)^{2}+\left(\rho u_{y}\right)^{2}}{\rho}\\ M_{2}^{(eq)}=M_{2}=\rho u_{y}&M_{5}^{(eq)}=-\rho u_{x}&M_{8}^{(eq)}=\frac{5}{3}\rho-3\frac{\left(\rho u_{x}\right)^{2}+\left(\rho u_{y}\right)^{2}}{\rho}\\ \end{array}$$

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

$$\overline{\left(XY\right)}=\text{vbox}\text{hrule}height=0.5ptX\text{vbox}\text{hrule}height=0.5ptY+\frac{{\Delta }^2}{12}\left({\partial }_{\beta }\text{vbox}\text{hrule}height=0.5ptX\right)\left({\partial }_{\beta }\text{vbox}\text{hrule}height=0.5ptY\right)+O\left({\Delta }^4\right)$$

$$\overline{\left(\frac{XY}{Z}\right)}=\frac{\text{vbox}\text{hrule}height=0.5ptX\text{vbox}\text{hrule}height=0.5ptY}{\text{vbox}\text{hrule}height=0.5ptZ}+\frac{{\Delta }^2}{12\text{vbox}\text{hrule}height=0.5ptZ}\left[\left({\partial }_{\beta }\text{vbox}\text{hrule}height=0.5ptX\right)\left({\partial }_{\beta }\text{vbox}\text{hrule}height=0.5ptY\right)-\frac{\left({\partial }_{\beta }\text{vbox}\text{hrule}height=0.5ptZ\right)}{\text{vbox}\text{hrule}height=0.5ptZ}\left(\text{vbox}\text{hrule}height=0.5ptX\left({\partial }_{\beta }\text{vbox}\text{hrule}height=0.5ptY\right)+\text{vbox}\text{hrule}height=0.5ptY\left({\partial }_{\beta }\text{vbox}\text{hrule}height=0.5ptX\right)-\frac{\text{vbox}\text{hrule}height=0.5ptX\text{vbox}\text{hrule}height=0.5ptY\left({\partial }_{\beta }\text{vbox}\text{hrule}height=0.5ptZ\right)}{\text{vbox}\text{hrule}height=0.5ptZ}\right)\right]+O\left({\Delta }^4\right)$$

$$\text{vbox}{\text{hrule}height=0.5ptM}_i^{(\mathrm{eq})}=M_i^{(\mathrm{eq})}\left(\text{vbox}{\text{hrule}height=0.5ptM}_0,\text{vbox}{\text{hrule}height=0.5ptM}_1,\text{vbox}{\text{hrule}height=0.5ptM}_2,\text{vbox}{\text{hrule}height=0.5ptN}_{x0},\text{vbox}{\text{hrule}height=0.5ptN}_{y0}\right)+{\Delta }^2\text{vbox}{\text{hrule}height=0.5ptM}_i^{(\Delta )}i=\mathrm{0....8}$$

$$\begin{array}{rl}\text{vbox}{\text{hrule}height=0.5ptM}_3^{(\mathrm{eq})}& =\frac{\text{vbox}\text{hrule}height=0.5pt\rho u_x\text{vbox}\text{hrule}height=0.5pt\rho u_y}{\bar{\rho }}-\text{vbox}{\text{hrule}height=0.5ptB}_x\text{vbox}{\text{hrule}height=0.5ptB}_y\\ & +\frac{{\Delta }^2}{12\bar{\rho }}\left[\left({\partial }_{\beta }\text{vbox}\text{hrule}height=0.5pt\rho u_x\right)\left({\partial }_{\beta }\text{vbox}\text{hrule}height=0.5pt\rho u_y\right)-\frac{\left({\partial }_{\beta }\bar{\rho }\right)}{\bar{\rho }}\left(\text{vbox}\text{hrule}height=0.5pt\rho u_x\left({\partial }_{\beta }\text{vbox}\text{hrule}height=0.5pt\rho u_y\right)+\text{vbox}\text{hrule}height=0.5pt\rho u_y\left({\partial }_{\beta }\text{vbox}\text{hrule}height=0.5pt\rho u_x\right)-\frac{\text{vbox}\text{hrule}height=0.5pt\rho u_x\text{vbox}\text{hrule}height=0.5pt\rho u_y\left({\partial }_{\beta }\bar{\rho }\right)}{\bar{\rho }}\right)\right]\\ & -\frac{{\Delta }^2}{12}\left({\partial }_{\beta }\text{vbox}{\text{hrule}height=0.5ptB}_x\right)\left({\partial }_{\beta }\text{vbox}{\text{hrule}height=0.5ptB}_y\right)+\mathcal{O}\left({\Delta }^4\right).\end{array}$$

$$\begin{array}{rl}\text{vbox}{\text{hrule}height=0.5ptM}_3^{(\mathrm{eq})}& =\frac{\text{vbox}\text{hrule}height=0.5pt\rho u_x\text{vbox}\text{hrule}height=0.5pt\rho u_y}{\bar{\rho }}-\text{vbox}{\text{hrule}height=0.5ptB}_x\text{vbox}{\text{hrule}height=0.5ptB}_y+\frac{{\Delta }^2}{12\bar{\rho }}\left({\partial }_{\beta }\text{vbox}\text{hrule}height=0.5pt\rho u_x\right)\left({\partial }_{\beta }\text{vbox}\text{hrule}height=0.5pt\rho u_y\right)-\frac{{\Delta }^2}{12}\left({\partial }_{\beta }\text{vbox}{\text{hrule}height=0.5ptB}_x\right)\left({\partial }_{\beta }\text{vbox}{\text{hrule}height=0.5ptB}_y\right)+\mathcal{O}\left({\Delta }^4\right).\end{array}$$

Full text

1 Introduction

Recently [2] we derived a first-principles two dimensional (2D) magnetohydrodynamic (MHD) large eddy simulation (LES) model based on first filtering the lattice Boltzmann (LB) representation of MHD [3, 4, 5, 6, 7, 8, 9, 10] after which one applies the Chapman-Enskog limits to recover the final LES-MHD fluid equations. In essence, we extended to MHD the 2D Navier-Stokes (NS) LES-LB model of Ansumali et. al. [1] who exploited the non-commutativity of these two operations. (Of course, if one first applied the Chapman-Enskog limit to LB and then filtering one would land in the conventional quagmire of an LES closure problem.) A technical difficulty with the Ansumali et. al. model is that in 2D NS there is an inverse energy cascade to large spatial scales thereby rendering subgrid modeling non-essential. In 2D MHD, however, the energy cascades to small spatial scales as in 3D - and so makes it attractive to perform the LES-LB-MHD simulations in which there can be a substantial amount of excited subgrid modes. Here we present some preliminary LES-LB-MHD simulations of our model and compare the results with some direct numerical simulations (DNS). As Ansumali et. al. [1] did not perform any simulations on their LES-LB-NS model, these are the first such LES-LB-MHD simulations when one has first filtered the underlying LB representation, followed by the conventional small Knudsen number expansion.

The backbone of any LES [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 2] is the introduction of a spatial filter function to smooth out field fluctuations on the order of the filter width $\Delta$. Thus for the mean velocity

[TABLE]

In general, the filtering results in the standard closure problem. Previous LB-LES-NS modeling [28, 27] have first considered the Chapman-Enskog expansion followed by filter and thus have concentrated on the Smagorinsky closure for the subgrid stresses. It has been pointed out [29] that in the conventional NS-LES closure, the subgrid stresses are assumed to be in equilibrium with the filtered strain. However, in LB-LES-NS the stresses relax towards the filtered strain at a rate dictated by the current eddy-viscosity, thereby permitting some spatio-temporal memory effects that are absent when applying LES directly to the continuum equations. In essence this gives an edge to any LB - approach.

In the Ansumali approach, however, one first performs a perturbation expansion in the filter width $\Delta$ followed by the standard LB Chapman-Enskog expansion in the Knudsen number $\mathrm{Kn}$. These two perturbation expansions do not commute. Closure is now achieved by making the physically plausible assumption that eddy transport effects occur at the transport time scale and this results in the scaling $\Delta\simeq\order{\mathrm{Kn}^{1/2}}$. One still retains the LB effects of spatial -temporal memory as noted earlier [28, 27] .

2 LES-LB-MHD Model

For completeness we briefly review the essentials of our LES-LB-MHD model [2] that yields the following closed set of filtered MHD equations (without further approximations) for the filtered density $\left(\bar{\rho}\right)$, momentum \left(\,\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.99997pt\rho\mathbf{u}\kern 0.0pt}}}\right) and the magnetic field $(\overline{\mathbf{B}})$ in multiple relaxation time (MRT)

[TABLE]

where $s_{3}\ldots s_{8}$ are relaxation rates, and in this isothermal model, the pressure is directly related to the density $p=\rho c_{s}^{2}=\frac{\rho}{3}$, in lattice units ($c_{s}$ is the sound speed). The transport coefficients (shear viscosity $\nu$, bulk viscosity $\xi$ and resistivity $\eta$) are determined from the LB-MRT for the particle distribution function (the $s^{\prime}s$) and the single magnetic distribution function relaxation rate $s_{m}$:

[TABLE]

[TABLE]

[TABLE]

We now summarize our computational LB-LES-MHD model that underlies Eqs. (2). For 2D MHD, we consider an LB model with 9-bit lattice

[TABLE]

with the moments $\sum_{i}f_{i}=\rho$, $\sum_{i}f_{i}\vec{c}_{i}=\rho\vec{u}$, and $\sum_{k}\vec{g}_{k}=\vec{B}$. Here the summation convention is employed on the vector nature of the fields (using Greek indices) while for Roman indices, correspond to the corresponding lattice vectors for the kinetic velocities $\vec{c}_{i}$, there is no implicitly implied summation. The lattice is just the axes and diagonals of a square (along with the rest particle $i=0$). $s^{{}^{\prime}}_{ij}$ are the MRT collisional relaxation rate tensor for the $f_{i}$ while the SRT $s^{{}^{\prime}}_{m}$ is the collisional relaxation rate for $\vec{g}_{i}$. These kinetic relaxation rates determine the MHD viscosity and resistivity transport coefficients. (Of course, more sophisticated LB models can be formed by MRT on the $\vec{g}_{k}$ equations, but for this first reported LB-LES-MHD simulation we will restrict ourselves to the simpler SRT model)

A convenient choice of the relaxation distribution functions, will under Chapman-Enskog, yield the MHD equations

[TABLE]

where the $w^{\prime}s$ are appropriate lattice weights. In the operator-splitting solution method of collide-stream, it is most convenient to perform the collision step in moment-space (because of collisional invariants of the zeroth and 1st moment of $f_{i}$, and the zeroth moment of $\vec{g}_{i}$.), while the streaming is optimally done in the $(f_{i},\vec{g}_{i})$-space. Moment space $(M_{i},\vec{N}_{i})$ is defined by

[TABLE]

with the 1-1 constant transformation matrices, $\mathrm{T}$ and $\mathrm{T}_{\mathrm{m}}$ given by

[TABLE]

The $x$ and $y$ components of the $9$-dimensional lattice vectors are

[TABLE]

In terms of conserved moments, we can write

[TABLE]

2.1 Filtering LB

In applying filtering to the LB Eqs. (6) and (7), only the nonlinear terms in the relaxation distributions, Eqs. (8) and (9), require further attention. On applying perturbations in the filter width $\Delta$ we immediately see that

[TABLE]

and

[TABLE]

for arbitrary fields $X$, $Y$, and $Z$. Moreover since collisions are performed in moment space, we need first to transform from $f^{(\mathrm{eq})},\vec{g}^{(\mathrm{eq})}$ to $M^{(\mathrm{eq})},\vec{N}^{(\mathrm{eq})}$ and then apply filtering in terms of the filtered collisional invariants \hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.99997ptM\kern 0.0pt}}}_{0},\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.99997ptM\kern 0.0pt}}}_{1},\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.99997ptM\kern 0.0pt}}}_{2},\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.99997ptN\kern 0.0pt}}}_{x0},\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.99997ptN\kern 0.0pt}}}_{y0}

[TABLE]

where the $O(\Delta^{2})$ term arises from the nonlinearities. In particular for the \hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.99997ptM\kern 0.0pt}}}_{3}^{(\mathrm{eq})}-term :

[TABLE]

Similarly for the other filtered equilibrium moments.

3 LES-LB-MHD Simulation

The filtered LB equations are now solved, with streaming performed in distribution space and collisions in moment space. As this is the first simulation on the LB-filtered-LES approach we have made a significant number of simplifications. We first restrict the evolution of the filtered scalar distribution function to an SRT collision operator. In this case the relaxation rates $s_{i}$ are all equal so that the 3rd term in Eq. (2b) is automatically zero. Moreover since nearly all LB simulations are quasi-incompressible at the fluid level, we neglect (filtered) density gradients in the moment representation of the collision operator. Thus, for example, we approximate \hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.99997ptM\kern 0.0pt}}}_{3}^{(\mathrm{eq})} by

[TABLE]

Also, since the last term in Eq. (2c) is dependent on the filtered density gradient its effects at the filtered MHD level will not be significant when we code the filtered LB system.

It should be noted that as in regular LB-MHD, the filtered $\nabla\cdot\overline{\mathbf{B}}=0$ is maintained to machine accuracy.

There is a little subtlety in that not all the spatial derivatives in the filtered collision moments can be determined from local perturbed moments [21, 4]. This limitation is thought to arise from the low D2Q9 lattice. It is expected that on a D3Q27 lattice the linearly independent set of derivatives can be represented by the now larger number of local perturbed moments.

While we solve the filtered LB equations, resulting in the filtered LES MHD Eqs. (2), there is some similarity in our final MHD model with that of the ”tensor diffusivity” model of Müller-Carati [12]. However it must be stressed that we are performing a first principles derivation of the eddy transport coefficients from a kinetic (LB) model while Muller-Carati propose an ad hoc scheme of minimizing the error between two filters at each time step in their determination of their model’s transport coefficients.

We will now evolve in time our filtered LB equations and consider the magnetized Kelvin-Helmholtz instability in a sufficiently weak magnetic field so that the 2D velocity jet is not stabilized [30]. The initial jet velocity profile is $U_{y}=U_{0}\sech^{2}{\left(\frac{2\pi}{L}4x\right)}$. The corresponding vorticity is shown in Fig. 1. The initial Reynolds number is chosen to be Re = $U_{0}L/\nu=50$k = const., with $U_{0}=4.88\times 10^{-2}$ and $B_{0}=0.005U_{0}$. The viscosity and resistivity on a grid of $1024^{2}$ are $\nu=\eta=10^{-3}$ and scale with the grid to maintain a constant Re and a constant magnetic Reynolds number $U_{0}L/\eta$. The initial perturbation to the fields are: $U_{y}=0.01U_{0}\sin{\left(\frac{2\pi}{L}4x\right)}$, $B_{y}=0.01B_{0}\sin{\left(\frac{2\pi}{L}4x\right)}$, $U_{x}=0.01U_{0}\sin{\left(\frac{2\pi}{L}4y\right)}$, and $B_{x}=0.01B_{0}\sin{\left(\frac{2\pi}{L}4y\right)}$. Note that initially $\nabla\cdot\vec{B}=0=\nabla\cdot\vec{U}$.

In Fig, 2 we compare the evolution of vorticity in time from DNS on a $2048^{2}$ grid with that determined from our LES-LB-MHD model on a $1024^{2}$ grid. The DNS simulations are determined by solving the direct unfiltered LB Eqs. (6) and (7). For constant Reynolds number simulations at different grid sizes, the kinematic viscosity is adjusted appropriately. Thus on halving the spatial grid, a DNS time step of $2t_{0}$ corresponds to time step $t_{0}$ in LES-LB-MHD.

At relatively early times the jet profile width slightly widens while within the vorticity layers the Kelvin-Helmholtz instability will break these layers into the familiar vortex street (Fig. 2). Since we have chosen a weak magnetic field insufficient to stabilize the jet, the vorticity streets break apart with like vortex-vortex reconnection (Fig. 2). There is very good agreement between DNS and LES-LB-MHD with filter width $\Delta=2$ (in lattice units) on a grid $L/2$.

Finally, we show the corresponding vorticity (Fig. 3), total energy spectrum (Fig. 4), and current (Fig. 5) plots at t = 780k for simulations on $1024^{2}$ grids and their counterparts on $2048^{2}$ grids at time t = 1.56M. We consider 4 cases: (a) DNS on $1024^{2}$, (b) filtered LB-LES-MHD on $1024^{2}$ grid and small filter width, $\Delta=1$, (c) DNS on $2048^{2}$ and (d) filtered LES-LB-MHD on a $1024^{2}$ grid but with filter width $\Delta=2$. The effect of the filter width $\Delta$ in our LB-LES-MHD model on the evolution of the vorticity is evident when comparing Fig. 3b to Fig. 3d - both in location and strength of the main vortices as well as in the fine grained small scale vorticity. As the filter width increases to $\Delta=2$ (fig. 3d), there is stronger agreement now with the DNS (fig. 3c) on $L^{2}$ grid with our LES-LB-MHD filtered model on $(L/2)^{2}$ grid. This shows that the subgrid terms are now influencing larger scales with some accuracy. The spectral plots (fig. 4) are somewhat similar in all simulations with a very localized Kolmogorov energy spectrum. Presumably this is because the turbulence is limited and relatively weak. There appears to be good agreement in both the vorticity and current between DNS and LES-LB-MHD with $\Delta=2$ on half the grid.

4 Conclusion

Here we have presented some preliminary 2D filtered SRT LB-MHD simulation results based on an extension of ideas of Ansumali et. al. [1] that leads to a self-consistent LES-LB closure scheme based solely on expansions in the filter width $\Delta$ and invoking the constraint that any eddy transport effects can only occur on the transport time scales. We find very good agreement between DNS and our LES-LB-MHD models. This warrants further investigation of other filters used in LES, as well as in dynamic subgridding commonly used in LES of Navier-Stokes turbulence. Finally, an exploration of the effects of MRT on this LES algorithm should be quite interesting as a somewhat unexpected term related to the gradient of a pressure appears in the subgrid viscosity. This term reveals that higher-order moments (not stress related) can have a first order effect on the subgrid viscosity when MRT is employed. Given that this subgrid pressure term relies on the existence of higher order moments, it suggests that the extra parameters in lattice Boltzmann (ie. the distribution velocities/moments) are introducing new physics naturally absent from LES in computational fluid dynamics. It would be very interesting to see whether this new term enhances the LES accuracy or increases stability at even higher Reynold’s flow. Further study could include how this term effects other, well-established LES approaches in computational fluid dynamics. These ideas are under consideration.

5 Acknowledgments

This work was partially supported by an AFOSR and NSF grant. The computations were performed on Department of Defense supercomputers.