Fitting of extended sub-grid scale models in compressible turbulent MHD

Fitting of extended sub-grid scale models in compressible turbulent MHD

Daniele Viganò    Carlos Palenzuela Departament de Física, Universitat de les Illes Balears, Institut d’Estudis Espacials de Catalunya and Institut d’Aplicacions Computacionals de Codi Comunitari (IAC3), Palma de Mallorca, Baleares E-07122, Spain
July 26, 2019

The use of Large Eddy Simulations in compressible, turbulent MHD is often limited to the implementation of purely dissipative sub-grid scale models. While such models work well in the hydrodynamic case due to the universality of the spectrum, they do not fully describe the complex dynamics of MHD, where the transfer of energy between internal, kinetic and magnetic energies at small scales are less trivial. For this reason, a sub-grid scale model based on the gradients of the fields entering in each non-linear term of the equations has already been proposed and studied in the literature, for the momentum and induction equations. Here, we consider the full set of compressible, ideal MHD equations, with an ideal gas equation of state, and we proposed a generalization of the gradient model, including the energy equation. We focus on the residuals coming from the whole set of equations, by filtering accurate high-resolution simulations of the turbulent Newtonian Kelvin-Helmholtz instability in a periodic box. We employ the same high-resolution shock capturing methods typically used in relativistic MHD, applicable in particular to neutron star mergers. The a-priori test, i.e. the fit between the sub-filter residuals and the model, allows us to confirm that the gradient model outperforms any other, in terms of accuracy of the fit and small deviations of the best-fit pre-coefficient from the expected value. Such results are validated for 2D and 3D, for a range of different problems, and are shown, for the first time, to hold also for the energy evolution equation. This paper is the first step, based on a solid theoretical and numerical basis, towards the near-future extension of the sub-grid scale gradient model to the relativistic MHD, and a future implementation in a full General Relativity LES.

I Introduction

In the context of computational fluid dynamics, direct numerical simulations (DNS) are employed to resolve the equations in all the dynamically relevant scales, including the dissipative ones. This means that the numerical grid size has to be smaller than the dissipative scale. Since the computational cost of DNSs scales as (Sagaut, Deck, and Terracol, 2006) (where is the Reynolds number), for many problems (including convection zones of stars, planetary atmospheres, accretion disks, and industrial applications such as gas turbines, steady turbulent flows…) it is unfeasible to have a small enough grid size and a large enough domain at the same time, thus representing all the scales of interest. Even though in the near future the increasing computation power will improve the achievable resolution, thus making DNSs feasible for more scenarios, currently there are many cases where we are still very far from it. In astrophysics, one of these cases, which is of particular interest due to the recent discoveries of gravitational waves, is the binary neutron star merger.

The two most widely used numerical approaches to overcome the lack of resolution in turbulent flows are the Reynolds Average Numerical Simulations and the Large Eddy Simulations (LESs). For both cases, the idea is to directly evolve the quantities on the resolved scales, i.e. from the integral scales (the largest ones, corresponding to dominion size or to the typical scale where energy is injected), down to the grid scales, and to describe the effect of the smaller ones through a suitable sub-grid scale (SGS) model.

While the purely hydrodynamic turbulence dynamics is relatively well understood due to the universality of the spectral properties, the dynamics becomes more problem-dependent when magnetic fields are included. Magnetic fields tend to break the usually assumed isotropy and self-similarity, and prevent one from defining universal laws, valid for different scenarios. Magneto-hydrodynamic (MHD) turbulence includes a non-linear dynamics, with transfer between magnetic and kinetic energy (dynamo mechanism, when the magnetic energy rises), and, at the same time, between different scales. Such transfers vary locally and in time, so that a statistical evaluation of them represents essential information. In general, the smallest scales are the most interesting ones in terms of interchange of kinetic and magnetic energy.

As a consequence, it is crucial to properly model the unresolved scale dynamics, by means of a SGS modeling, that is incorporated as additional extra terms in the MHD equations. Due to the variety of turbulent scenarios, LESs with SGS models represent an active field of inter-disciplinary research, including engineering, plasma physics and astrophysics (i.e., see for instance a review about the use of LESs in astrophysicsMiesch et al. (2015), which is also our main interest). A popular SGS model is the dissipative one, introduced by Smagorinsky half a century ago in the context of hydrodynamic direct cascadesSmagorinsky (1963) to mimic the fluid’s viscosity. Other SGS models have also been proposed relying mainly on dimensional analysis or physical assumptions like the dynamo mechanism.

Another branch of SGS models (including the multi-scale and multi-variational ones)Hughes et al. (2001) relies on self-similarity arguments, where basically the loss of information is estimated considering the average spatial correlations of the evolved fields at different scales, and extrapolating such information to infer the fields in the SGSs. These models are used especially in finite-elements schemes.

On the other hand, the gradient modelMüller and Carati (2002) is based on a mathematical formulation which relates the finite resolution to an effective filtering of the equations. By expanding in Taylor series the non-linear terms in any equation, one can approximate the SGS residuals with a functional form which only relies on mathematical assumptions. The model does not contain any arbitrary physical assumption, and just extrapolates to the SGS the trends seen at the smallest resolved scales.

In general, two steps can be taken in order to assess and compare the performance of any SGS model:

  • A-priori fitting. This test require the post-process filtering of the fields obtained with a relatively high-resolution simulation at a given time. The residuals coming from the filtering (i.e., the non-linear information contained between the filter size and the original resolution ) can be fit to a given SGS model, where the only free parameter is a pre-coefficient. The procedure can be easily repeated for different times, resolutions, filter sizes, initial conditions, in order to test the applicability of the tested SGS model. This procedure allows to identify the best-fitting functional form of the SGS model, but does not allow to evaluate the effect of the inclusion of the SGS in a simulation.

  • LES comparison, including different SGS models. Ideally, a low-resolution LES with a good SGS model should resemble a LES with a higher resolution. The LES comparison is a more demanding test, since it involves the dynamical feedback given by the implemented SGS terms. This test is much more informative and is crucial to evaluate the physical meaning of including a certain SGS tensor in a LES. However, the exploration of parameters is much more computationally expensive.

In this paper, we will perform the first step: an extensive a-priori assessment of the most popular explicit local SGS models present in the literature, leaving a battery of LES comparison tests for a forthcoming paper, which will obviously rely on the results found with the first approach. Works in this sense have been performed in the context of either incompressibleKessar, Balarac, and Plunian (2016) or compressible Grete et al. (2016, 2017a, 2017b); Grete (2017) MHD. These studies are among the most detailed in comparing different SGS models, employing either forced or decaying turbulence. Although they consider the momentum and induction MHD equations, the internal energy evolution is not included, since they impose either incompressibility or a polytropic/isothermal equation of state (EoS) as a closure Grete et al. (2017b). Here we extend these works by considering instead the compressible MHD equations with an ideal EoS.

We run several bounding box high-resolution simulations of the Kelvin-Helmholtz instability (KHI), which triggers the gradual development of turbulence. We do not include any forcing: we focus on dynamical, decaying turbulence, since we aim at applying them to simulations of astrophysical transient phenomena. We perform a-priori tests, where the filtered residuals are compared to what we would obtain with different SGS models applied to the information available at the filter resolution. We validate both the method and the SGS models, identifying the best-fit coefficients for a variety of initial data, resolutions and filters.

In § II we review and extend the formal description of a LES, seen as a consequence of separating scales by applying a filter to the solution. We describe in some detail the filtered full MHD equations in § III. In § IV we introduce different SGS models, with special emphasis on the gradient model. In § V we describe the platform and the numerical methods used in this work (in Appendix A we describe the validation of the code with a set of benchmark tests in 2D). In § VI we present the main highlights of the results coming from 3D KHI simulations, showing a variety of cases. Conclusions are drawn in § VII.

Ii Equivalent filtering in Large Eddy Simulations

In a LES, the smallest simulated scales (i.e., the ones corresponding to the grid size ) are not sufficient to capture entirely all the relevant dynamics, which might occur in the SGS scales. Since the SGS part of the solution is unknown by definition, it can only be either neglected, leading to the so-called implicit LES,gri (2007) where a low-pass filterSun and Domaradzki (2018) or the internal numerical dissipation (used for instance in solar physicsGhizaru, Charbonneau, and Smolarkiewicz (2010)) act as an intrinsic SGS model, or explicitly modeled through the information which is available, i.e. the one of the resolved scales. This explicit modeling, on which we focus, can be based either on mathematical or physical arguments.

Mathematically, the effect of finite resolution can be thought as equivalent to a low-pass spatial filter, where the filter size, , is given by the grid cell, . This separates, for a continuous field , the smoothed (resolved) part, from the sub-filter scale (SFS) fluctuations (or residuals), respectively:


Indicating the filter kernel with , the filtering operator over a field can be written as


An homogeneous isotropic low-pass filter is independent on the direction (i.e., ) and only smooths out fluctuations on length scales smaller than the filter size , leaving unchanged the variations of the solution at larger length scales. In addition, the filter operator is linear and commutes with spatial derivatives. Generically, it can be written for any dimension as


where is just the one-dimensional kernel function.

The simplest low-pass filter is the mean value in a cubic domain with size in each Cartesian direction , described by the normalized kernel


Despite the appealing simplicity of the box filter, which makes it very useful to perform numerical calculations, we will see below that it is not suitable for analytical calculations involving its derivatives, since they are not continuous. Therefore, at a formal level, it is more practical to introduce the normalized Gaussian kernel, which in the space domain can be written as


where defines the effective filtering width. Besides having the same zeroth and first moments, Gaussian and box filters have the same second moment if we set :


The filter is a useful mathematical tool to analyze not only the different scales on the solution, but also the structure of the evolution equations, allowing us to distinguish between the resolved or filtered fields and the unresolved ones on the SFS.

In order to see the effect of this scale-separation in a simple application, let us consider the well-known Burgers equation, a longstanding field of research in LES.D.  Love (1980) In the 1D case, given the velocity field , the equation reads:


The equation for the filtered fields can be obtained by applying the filtering operator to the equations. Since the filtering operator commutes with both spatial and time derivatives, we have


Notice that the filtered equation is equivalent to the original one except by a term proportional to which represents the SFS residuals, that is, the loss of small-scale information due to the filtering process related to the non linear terms. If in the region over which we average has a definite sign (which is usually the case except close to the shock), then and is negative-definite, so that represents the kinetic energy hidden in the SFS.

Regardless of the specific equations considered, these new terms always appear due to the non-commutativity of the filtering operator with the non-linear terms of the equations. In a LES, which effectively only evolves and cannot simulate the unresolved scales, the explicit expression of is a-priori unknown, and needs to be written as a function of the filtered fields in order to close the system of equations.

There are several ways of doing that, based mostly either on physical arguments or on an expected self-similarity of the solution. For instance, a popular, historical approach to fluid dynamics SGSsSmagorinsky (1963) is the setup of an artificial viscosity , namely


The viscosity parameter is often taken to be proportional to due to dimensional reasons and in order to ensure that this term vanishes in the continuum limit, thus guaranteeing numerical convergence to the continuum (or DNS) solution. The numerical value of the proportionality coefficient can be fixed by hand or estimated by means of dynamical procedures which assume self-similarity (thus similar to a multiscale/multivariational approach).Germano et al. (1991) Regardless of the value of , this functional form of is equivalent to the viscous Burger’s equation, where the viscosity is the numerical SGS model.

However, in many cases the physics involved does not consist only in dissipation, but can involve, for instance, inverse cascade, and scale-dependent transfers of energy between the fluid kinetic energy and the magnetic one. Therefore, in those cases any result coming from the LES will be biased by the a-priori choice of a given SGS physics.

If one wants to avoid introducing arbitrary physical assumptions, it is possible to compute these terms relying only on mathematical arguments by considering the analytical Taylor expansion of the SFS termsLeonard (1975); Müller and Carati (2002). The main hypothesis is that the discretization of the equations over a finite-size grid is approximated by a filtering operator acting on the equations, with a Gaussian kernel equivalent to a box filter with width , eq. (5). Its -dimensional Fourier transformed function is


The main idea is to compute an approximation of the inverse filtering operator based on gradient expansion of the filter kernel , that is, an approximation of the inverse Fourier transform of . The first step is to perform a Taylor expansion of the transformed function and its inverse in terms of the filter scale, that is,


Considering the expansions to the fields and , the application of the inverse Fourier transformation yields to an infinite series representation of the filter operator and its inverse in terms of gradient operators acting on the fields, namelyVlaykov et al. (2016)


These expressions are absolutely convergent and formally accurate at all orders, since the Gaussian kernel is infinitely differentiable and with unbound support. In fact, it was found that these series converges for all canonical filters.

The unknown components of the SFS tensors usually have the form and . Applying the relations in Eq. (14) results in a series in terms of the individual fields, namely


which is accurate up to . The first order terms scale with the second moment of the filter, eq. (6). With this expansion is easy to model the SFS terms appearing in Eq. (8), namely


For a Cartesian grid with step , the pre-factor can be set to the equivalent value , according to eq. (6). In principle, corrections to these expression are expected due to the form factor (due to the kernel shape) and the contribution from higher orders. As a matter of fact, if the range of modeled SFSs, , becomes large, then the higher-order terms can be important. In other words, this approximated expression can help as long as we can capture most of the dynamics with our LES, leaving to the SGS the task of mimicking the small-scale contributions.

Regardless on these quite obvious caveats, notice that this prescription is very different from the viscous one given by eq. (9): one arises from physical considerations, while the other just from mathematical ones. After this simple example, we pass to the application of this formalism to the MHD equations.

Iii Filtering the Newtonian compressible MHD equations

iii.1 Compressible MHD equations

We consider the conservative formulation of the ideal Newtonian compressible MHD equations, consisting of the continuity, Euler, induction and energy evolution equations, respectively:


where the magnetized fluid is described by the total energy density defined as


the mass density , the fluid pressure , the internal energy density , and the velocity and magnetic field vectors, being the standard Kronecker delta. For simplicity, we have set the speed of light and the magnetic permeability . The conserved quantities of the equations , in their discretized versions, are numerically evolved. At each step of a simulation, the internal energy is recovered by inverting eq. (18).

The pressure, required to close the system of equations, is defined by an EoS, which in general can read . In this paper, we will present results with the ideal gas EoS, namely


where is the ideal gas adiabatic index. However, our formalism is general and can extend to any EoS.

iii.2 Filtered continuity equation and variable inversion

The filtered continuity equation is simply


where the Favre-filtered velocities, i.e., weighted by density(Miesch et al., 2015)


are the natural choice in compressible MHD, where the momenta are the evolved fields. Therefore, and are simply written as under the filtering operation.

We will use the same notation whenever we indicate an auxiliary field as a function of evolved filtered fields only. Specifically, in our case we can express the internal energy density as a function of the evolved filtered variables, eq. (18), defining


Note that there is a difference between what can be numerically evaluated with the resolved fields only, (i.e., combining the filtered evolved fields), and what cannot be, (i.e., filtering the local non-linear combinations of the variables):


which can be interpreted as the kinetic and magnetic energy density hidden in the SFSs. Note that, if such SFS contributions (usually positive-definite) are neglected in the variable inversion, this will cause a local artificial rise of the evaluated internal energy, .

Similarly, we can define the pressure as a function of all the other evolved variables . Its expression depends on the EoS, as we will explain below.

iii.3 Filtered momentum and induction equations

On the other hand, the filtered momentum and induction equations become, respectively


where , and indicate the tensors arising from the differences in the non linear products , and , respectively, when the filtering is considered:111In literature, different sign conventions are used. We use the same sign convention as Grete’s papers.Grete et al. (2015, 2016); Grete (2017) For instance, in Kessar’s workKessar, Balarac, and Plunian (2016), , are definite with an opposite sign, while has the same sign.


Notice that, by construction, and are symmetric tensors, while is anti-symmetric and its inclusion is equivalent to having a SFS electric field with components given by . Note that the SFS kinetic and magnetic energies defined in eq. (23) are related to and by


In the literature, sometimes the functional forms of the traces of and are defined separately from the ones of the off-diagonal components (related to the strain). This is useful as a tight analogy with the standard bulk and strain viscosity, but, since our formulation is general and does not rely on the underlying physics, we prefer to avoid this distinction.

The term associated to pressure is instead purely diagonal and has been only partially considered before. If isothermal EoS or incompressibility are assumed, pressure does not depend on the evolved field, and the only surviving term is the magnetic pressure, given by half the trace of . However, in the general case, the pressure depends non-linearly on the evolved fields, so that SFS residuals implicitly appear and is generally not zero. In the case of the ideal EoS, such term is proportional to eq. (23), so that:


Note that an ideal quasi-isothermal gas () can be taken as a shortcut to minimize the spurious increase of pressure.

iii.4 Filtered energy equation

While the momentum and induction equations have been explored in many different works, the energy evolution has usually been neglected (for incompressible and isothermal cases), or only partially considered. In particular, some authors have focused on the SFS kinetic and magnetic energy, studying in detail the inter-scale transfers between resolved and SFSsKessar, Balarac, and Plunian (2016); Vlaykov et al. (2016); Grete et al. (2017b). Those works aim at analyzing the energy transfer, but not to perform a detailed modeling of the SFS tensors appearing in the evolution equations. Hereafter, we propose two approaches. The first is a conservative one, similar to the previous equations, where only flux terms appear, and is apt for the implementations of SGS models in a LES. In the second one, we considered the already studied inter-scale energy transfers, appearing as source/sink terms.

Conservative approach (for SGS modeling in LES).
By following a similar approach as for the momentum and induction equations, we can write the energy evolution equation as


where the additional vectorial SFS terms on the right-hand side arise from the different non-linear terms in the original equation, as follows:


We distinguish three fundamental contributions. The first one is given by the SFS pressure advection term (), and the SFS convective term (). Combined, they represent the advection of the SFS enthalpy:


where we have exploited the identity to introduce the enthalpy density . In the case of an ideal EoS, this relation is reduced to


where the internal energy and pressure as functions of the filtered evolved quantities are given by eq. (22).

The second group of SFS terms in the right hand side of eq. (33) is given by the SFS advection of kinetic () and magnetic energy () densities. Note that, strictly speaking, the latter actually comes for one half from the magnetic energy advection, and for the other half from the magnetic pressure advection, similarly to . In order to keep the elegance of the enthalpy term described above, we prefer to keep the magnetic pressure inside this term.

The last term in the right hand side of eq. (33) is the cross-helicity term, which is partially related to the induction equation tensor, , but cannot be written only as a function of it (or a simple contraction with another filtered quantity).

The main disadvantage of this approach is the presence of triple products of three or four conserved fields in the residuals, which makes the calculation not straightforward. In this paper, we will not focus on this approach, which will be studied in depth in our generalized formulation to be published in a forthcoming paper.

Inter-scale approach (for analysis and fitting).
Alternatively, we can consider the evolution in time of the different components of the total energy density, . After some algebra, one can write the final expression


where .

This second approach is useful as an analysis tool, in the a-priori fitting tests, as in this paper.

iii.5 Final filtered equations

With these considerations, we can summarize the final set of equations to be evolved for the filtered fields are:






In the energy equation, we have reorganized the SFS source terms. Due to the (anti-)symmetric nature of the SFS tensors, we have that the integration by parts leaves a part of the terms as flux divergences, and another part which is the contraction of only the symmetric part of (i.e., the resolved fluid strain tensor ), and the anti-symmetric part of (i.e., the resolved current tensor ). The sink term is composed by the transfers of advected pressure, kinetic and magnetic energy from the resolved scales to the SFS, respectively:


In analogy with the viscous and resistive terms in non-ideal MHD (proportional to the strain and current tensors, respectively), at each point of the dominion represents the transfer of energy from the resolved scales into the SFSs (positive values) or viceversa (negative).

Note that, with this compact formulation, we need to give the prescriptions for four tensors, eqs. (48)-(51) and two scalars, related with pressure, i.e. eq. (54) and the term in eq. (52). The latter is automatically determined by and in the case of an ideal EoS. Once these six SFS tensors (having a total of 20 independent components in 3D) are prescribed, the sources and are as well, since they do not contain additional degrees of freedom.

Hereafter, to avoid confusion, we will indicate with an overline the SFS residuals tensors resulting from either the formal analysis or the post-process filtering of a numerical solution, and with a simple the analytical models of the SGS tensors that can be implemented in LES, which is the subject of the next section.

Iv SGS models

The explicit expression of the SGS tensors represents the closure of the system of equations in a LES. By definition, in a LES it is impossible to know the exact numerical values or the functional forms of the SFS tensors and sources above, and . The SGS models considered here are local, i.e. depend only on the derivatives of the fields at the smallest resolved scales: we do not consider the multi-scale (or self-similar) models, where a second filter is performed in order to evaluate the SGS model. We also neglect the structure models, based on different possible combinations of mutual contractions of the four first-derivative tensors (fluid strain, magnetic strain, current, vorticity), maintaining the parity needed by SGS models in the equations.(Speziale, 1991; Kosović, 1997; Lu and Rutland, 2016) In general, they do not provide a good fit (a-priori tests) and are also known to be numerically less stable when implemented in a LES. Below we list the ones we explicitly compare in this work.

iv.1 Gradient model (extended formulation)

One can consider the analytical Taylor expansion of the SFS terms appearing in the MHD equations, under the hypothesis of having a filter with a Gaussian kernel. To apply the presented derivation self-consistently to the compressible SFS stresses, one has to take into account the compressibility effects from the mass-weighted filtering operator (i.e., the Favre filtering), which is given by


Using the previous relations one can write


One can solve for , derivate this expression and use the recurrence relation to obtain, at linear order in ,


so that the following expression, function of instead of , holds


The correction, is zero in the incompressible case. Using the rules described before, the expansions for the SFS terms (48)-(50) give, at leading order in :


This formulation, with as in eq. (6), is known as the gradient modelMüller and Carati (2002) (sometimes called generically non-linear(Grete et al., 2016)) for MHD.

If we consider the dependencies of the total (fluid plus magnetic) pressure for an ideal gas with index , as given by eqs. (23) and (32), then eqs. (51)-(52) are expressed by:


Besides the well-known magnetic pressure contribution , there is an additional term, related to the fluid pressure , proportional to times the SFS kinetic and magnetic energy. Note again that for an isothermal EoS (), the fluid contribution vanishes and the only term left is the one related to the magnetic pressure, in agreement with previous formalism.Grete et al. (2016) The proper modeling of includes a much more complicated expression. Since in a LES we will use only conservative terms (i.e., no source terms will be modeled), we will avoid to model and study this term in detail.

While the induction and momentum equations obtained had already been explored in depth,Miesch et al. (2015); Grete (2017) the pressure- and enthalpy-related terms in the momentum and energy equations are a novelty, as far as we know.

The gradient model, despite its relative simplicity, strongly differs from the linear models, explained below and largely employed in literature, because there is no a-priori dissipative nature. Moreover, since they reproduce the dynamics seen as the smallest resolved scales where the inter-energy transfers are more effective, they allow a strong coupling between velocity and magnetic fields.

iv.2 Dissipative model (Smagorinsky)

The basic purely HD model (Smagorinsky model)Smagorinsky (1963) accounts for the viscosity at small scales, introducing an effective viscosity tensor given by


where is a factor to be determined, and the proportionality to is taken for dimensional reasons. For the isotropic HD turbulent cases, where the universal laws apply, the value of numerically in spectral methods, being around (Lu and Rutland, 2016; Germano et al., 1991), depending also on the problem. The coefficient can be determined with a dynamical procedure, instead of being fixed, as done for a channel flow.Germano et al. (1991); Oberai and Wanderer (2005)

A subsequent workTheobald, Fox, and Sofia (1994) extended the Newtonian Smagorinsky model to the MHD, introducing the analogous turbulent magnetic diffusivity in the induction equation, proportional to the current density tensor. We will test then the following tensors:


which are the ones commonly used in the literature, properly adapted (factors with introduced on dimensional basis) and extended with the magnetic tensor , that has been introduced only by analogy to the fluid strain, but do not corresponding to any dissipative mechanism. In literature, different normalizations are present, taking trace and off-diagonal components separately, or considering different dependences with . However, they do not show in general a substantial, systematic improvement,Grete et al. (2017a) and we only test the expression above.

Note also some authorsMüller and Carati (2002); Kessar, Balarac, and Plunian (2016) consider a mixed model, given by a superposition of the Smagorinsky model with the gradient model, basically in order to improve the numerical stability, compared to the gradient model alone.

By definition, the contractions of eq. (67) with the fluid strain tensor and of eq. (69) with the current tensors provide positive values of (as long as the magnetic term is not dominant) and respectively: this shows their dissipative nature.

iv.3 Cross-helicity model

This model(Müller and Carati, 2002) considers the dissipation of cross-helicity, which is given by two contributions, proportional to the fluid and magnetic strain rate tensors, , and to the product of current and vorticity , so that