A new class of finite element variational multiscale turbulence models for incompressible magnetohydrodynamics
David Sondak, John N. Shadid, Assad A. Oberai, Roger P. Pawlowski, Eric C. Cyr, Tom M. Smith
AA new class of finite element variational multiscale turbulencemodels for incompressible magnetohydrodynamics
D. Sondak a, ∗ , J.N. Shadid b , A.A. Oberai c , R.P. Pawlowski b , E.C. Cyr b , T.M. Smith b a University of Wisconsin-Madison, Department of Mathematcis b Sandia National Laboratories, Computational Mathematics Department c Rensselaer Polytechnic Institute, Department of Mechanical, Aerospace, and Nuclear Engineering
Abstract
New large eddy simulation (LES) turbulence models for incompressible magnetohydro-dynamics (MHD) derived from the variational multiscale (VMS) formulation for finiteelement simulations are introduced. The new models include the variational multiscaleformulation, a residual-based eddy viscosity model, and a mixed model that combinesboth of these component models. Each model contains terms that are proportional to theresidual of the incompressible MHD equations and is therefore numerically consistent.Moreover, each model is also dynamic, in that its effect vanishes when this residual issmall. The new models are tested on the decaying MHD Taylor Green vortex at lowand high Reynolds numbers. The evaluation of the models is based on comparisons withavailable data from direct numerical simulations (DNS) of the time evolution of energiesas well as energy spectra at various discrete times. A numerical study, on a sequence ofmeshes, is presented that demonstrates that the large eddy simulation approaches theDNS solution for these quantities with spatial mesh refinement.
Keywords: turbulence models, magnetohydrodynamics, finite elements, variationalmultiscale formulation
1. Introduction
The study of fluid turbulence has a long and storied history that still lacks an ex-planation of the fundamental mechanisms driving this phenomenon [4, 43]. Despite this,the significance and ubiquity of turbulence drives us in trying to unlock its secrets. Anadditional layer of complexity is added when considering electrically conducting fluidssuch as liquid metals and plasmas and their interaction with a magnetic field. Thesetypes of fluids are even more common in the universe than non-electrically conductingfluids. They are found in stars, at the cores of planets and in the interstellar medium. Onearth, engineers try to exploit their properties to generate electricity and fusion energy.In most situations, these fluids are in a state of turbulence. In an important subset of ∗ Corresponding author
Email address: [email protected] (D. Sondak)
Preprint submitted to J. Comput. Phys. August 21, 2018 a r X i v : . [ phy s i c s . c o m p - ph ] D ec hese cases, liquid metals, and under certain conditions plasmas, can be described byincompressible magnetohydrodynamics (MHD) [15, 23, 24].The computational simulation of turbulence is extremely challenging due to its in-herent multiscale nature. A perfect numerical experiment (direct numerical simulation,DNS) would capture all motions down to the point at which energy is dissipated asheat by molecular dissipation. This is entirely infeasible for most problems and is likelyto remain so for the foreseeable future. The concept of large eddy simulation (LES)was introduced as a way of numerically resolving the largest-scales in the flow field whilemathematically modeling the effect of the small-scales on the large-scales; see [45] and [46]for an overview of LES for hydrodynamics. Several approaches for LES in MHD havebeen proposed. In [55, 25, 59, 38, 42, 33, 3, 26] many techniques are presented includingeddy viscosity models (EVMs), mean-field closures, and the Lagrangian-averaged MHDequations. A comprehensive review of LES models in astrophysics can be found in [47].Most simulations of MHD turbulence rely on spectral numerical methods (see [41]for background on spectral methods). Although very fast (exponential convergence) andrelatively straightforward to implement, these methods are not particularly amenable tothe simulation of flows in complex geometries. The finite element (FE) method is anattractive alternative that enables computational simulation in both simple and complexgeometries. In the context of fluid mechanics, and incompressible fluids in particular,there are a number of important issues that FE methods must address. These includesatisfaction of the incompressibility constraint, the inf-sup condition pertaining to con-ditions under which equal order elements for different fields may be used [30, 31], andthe control of spatial oscillations due to highly convected flow and unresolved internaland boundary layers. Among the most popular and widely implemented solutions tothese problems have been stabilized finite element methods [54, 20, 10, 12]. Around adecade after they were introduced, these methods were shown to derive from the generalframework of the variational multiscale (VMS) method [28, 32]. The VMS method wassubsequently used to develop LES turbulence models for hydrodynamics [27, 5].In spite of this work, LES models for MHD in the context of finite elements arerelatively unexplored. Recent work has focused on developing accurate and robust sta-bilized finite element methods for incompressible MHD [21, 11, 50, 2, 51]. We refer toclassical stabilized methods as incomplete turbulence models since they represent onlypart of the story coming from the VMS method. In the context of spectral formulationsVMS methods were used to derive new LES models for incompressible MHD in [53].However, the development of complete VMS based LES models for MHD that can beimplemented in a finite element context has not yet been explored. This is the centraltheme of this manuscript. In particular, the present work pursues two generalizationsto previous work on a finite element based VMS model for the large eddy simulation ofincompressible MHD flow: (1) the generalization of the spectral VMS turbulence modelsfor MHD that were developed in [53] and (2) the generalization of the VMS based finiteelement methods developed for incompressible MHD in [51].The previous finite element based VMS models for MHD discussed above neglectedsecond order correlations between unresolved scales (which we collectively refer to as“Reynolds” stresses) and included an incomplete version of the cross stresses. The presentwork fully accounts for the effects of the cross and Reynolds stresses in MHD throughthe full VMS model and a VMS-based eddy viscosity model. Our new models containseveral new terms which include two cross stress terms in the momentum equation that2re not present in classical stabilized formulations, two cross stress terms in the inductionequation that are not present in classical stabilized formulations, Reynolds stresses inboth the momentum and induction equations, and two VMS-based eddy viscosity modelsfor the Reynolds stresses in the momentum and induction equations.The layout of the remainder of the paper is as follows. Section 2 introduces the govern-ing equations both in the strong form and variational forms. Following this, section 2.3describes the finite element method and the general computational framework. Next,in section 3, a class of turbulence models is introduced for use with the finite elementmethod. A summary of numerical solution methods is presented in section 4. Section 5describes the test problem that was used to assess the model performance and discussesresults. Finally, section 6 reflects on the new models and considers improvements andother directions for the future.
2. Background and Numerical Method
The base model considered in this study is the incompressible resistive MHD systemthat is solved for a velocity field u ( x , t ) = [ u , u , u ] T and a magnetic induction B ( x , t ) = [ B , B , B ] T for x = [ x , x , x ] T in a domain Ω with boundary ∂ Ω. Weconsider an incompressible velocity field so that ∇ · u = 0. The momentum equation is ρ ∂ u ∂t + ρ ∇ · ( u ⊗ u ) = ∇ · T + J × B + f V . (1)In the momentum equation, (1), ρ is the fluid density, T is the fluid stress tensor, J × B isthe Lorentz force where J is the current density, and f V is a general body force. From thelow frequency Maxwell’s equations (specifically, Amp`ere’s law) we have J = ∇ × B /µ where µ is the magnetic permeability of free space. This allows the Lorentz force to berewritten as J × B = 1 µ ∇ · ( B ⊗ B ) − ∇ (cid:18) µ (cid:107) B (cid:107) (cid:19) . (2)Under the assumption of a Newtonian fluid, the stress tensor becomes T = − P I + Π (3)where P is the fluid pressure, I is the identity tensor, and the viscous part of the stresstensor is given by Π = − µ ( ∇ · u ) I + µ (cid:16) ∇ u + ( ∇ u ) T (cid:17) . (4)Using (2), (3), and (4) along with the divergence free constraint on the velocity fieldyields the final form of the momentum equation, ρ ∂ u ∂t + ρ ∇ · ( u ⊗ u ) − µ ∇ · ( B ⊗ B ) = ∇ · (cid:16) µ (cid:16) ∇ u + ( ∇ u ) T (cid:17)(cid:17) − ∇ P − ∇ (cid:18) µ (cid:107) B (cid:107) (cid:19) + f V . (5)3n equation for the magnetic induction is obtained by combining the low frequencyMaxwell’s equations (i.e. without a correction for the displacement current) with Ohm’slaw. This leads to ∂ B ∂t + ∇ · ( u ⊗ B − B ⊗ u ) = ∇ · (cid:16) λ (cid:16) ∇ B − ( ∇ B ) T (cid:17)(cid:17) (6)where λ is the magnetic diffusivity. The satisfaction of the solenoidal constraint on themagnetic induction, ∇ · B = 0, has garnered considerable attention in the numericalcommunity [56, 10]. We take the approach of introducing a scalar Lagrange multiplier ψ into the induction equation. This permits the enforcement of the solendoidal condition asa divergence free constraint on the magnetic field. Such an approach is common in finiteelement [10, 13] and finite volume methods [56, 16, 9]. Note that if a Fourier spectralmethod is used then the solenoidal constraints are satisfied identically and projectionsof ψ into the Fourier spectral space are also identically zero. Hence, in a pure Fourierspectral method the introduction of the scalar Lagrange multiplier is not necessary [53].The final form of the induction equation is ∂ B ∂t + ∇ · ( u ⊗ B − B ⊗ u ) = ∇ · (cid:16) λ (cid:16) ∇ B − ( ∇ B ) T (cid:17)(cid:17) − ∇ ψ + f I (7)where we have also included f I as a possible external electromagnetic force.We introduce some notation to simplify the subsequent exposition. First, denote thevector of solutions by U = [ u , P, B , ψ ] T . Then the nonlinear terms in the momentumequation are given by N V ( U ) = ρ u ⊗ u − µ B ⊗ B + 12 µ (cid:107) B (cid:107) I . (8)The viscous part of the stress tensor is F V ( U ) = µ (cid:16) ∇ u + ( ∇ u ) T (cid:17) . (9)Similarly, the nonlinear terms in the induction equation are given by N I ( U ) = u ⊗ B − B ⊗ u (10)and the magnetic diffusive flux is F I ( U ) = λ (cid:16) ∇ B − ( ∇ B ) T (cid:17) . (11)Putting this notation together, the final form of the incompressible MHD equations is ρ ∂ u ∂t + ∇ · (cid:16) N V ( U ) − F V ( U ) (cid:17) + ∇ P = f V (12) ∂ B ∂t + ∇ · (cid:16) N I ( U ) − F I ( U ) (cid:17) + ∇ ψ = f I (13) ∇ · u = ∇ · B = 0 . (14)4 .2. Variational Statement of the Governing Equations We use a finite element method to solve the MHD system ((12), (13), and (14)) ina domain Ω with either periodic boundary conditions for all variables or homogeneousDirichlet boundary conditions for the velocity and the induction fields. Both of theseboundary conditions imply that the space of trial and weighting functions are the same.We note that we have made this choice in order to simplify the notation. The methodsdescribed in this paper can easily be extended to more complex boundary conditions.Before discussing the numerical approach, we begin by stating the variational counterpartto the system described above:
Find U ∈ V s.t. ∀ W ∈ V A ( W , U ) = ( W , F ) (15)where W = [ w , q, v , s ] T is a vector of weighting functions and F = (cid:2) f V , , f I , (cid:3) T isa vector of forcing functions. The notation ( · , · ) represents an L inner product of twofunctions, ( u , v ) = (cid:90) u · v dΩ . (16)We must specify the function space V in which the solutions and weighting functionsreside. This function space consists of functions that are sufficiently smooth and satisfythe requisite boundary conditions and is given by V ≡ { W | W = [ w , q, v , s ] T ; w i , v i ∈ H (Ω) , q, s ∈ L (Ω) } . (17)In (17), L denotes the space of square-integrable functions while H is the spaceof square-integrable functions with square-integrable derivatives. The semilinear formin (15) is given by A ( W , U ) = A V ( W , U ) − ( ∇ q, u ) + A I ( W , U ) − ( ∇ s, B ) (18)( W , F ) = (cid:0) w , f V (cid:1) + (cid:0) v , f I (cid:1) (19)where A V ( W , U ) = (cid:18) w , ρ ∂ u ∂t (cid:19) − (cid:16) ∇ w , N V ( U ) (cid:17) − ( ∇ · w , P ) + (cid:0) ∇ w , F V ( U ) (cid:1) (20) A I ( W , U ) = (cid:18) v , ∂ B ∂t (cid:19) − (cid:16) ∇ v , N I ( U ) (cid:17) + ( v , ∇ ψ ) + (cid:0) ∇ v , F I ( U ) (cid:1) . (21)Note that no boundary terms appear in (20) and (21) as a result of the boundary con-ditions of the problem. We are now prepared to introduce the finite element method. The numerical solutionconsists of selecting specific functions from the function space V that provide reasonableapproximations to the actual functions that represent the solutions to the equations.5hat is U h ≈ U . Hence, we will select functions from a conforming, finite-dimensionalsubspace of V which we denote V h . Mathematically, V h ⊂ V . The parameter h indicatesthe finite dimensional space or a function that was chosen from the space V h . In practicalterms, h represents the mesh parameter and determines the size of the finite elementsbeing employed.We denote an element in physical space as K and an element in parametric space as K e . Coordinates in physical space are x = x i , i = 1 , , ξ = ξ i , i = 1 , ,
3. We consider an isoparamteric mapping, x ( ξ ) = (cid:88) a =1 x ea φ a ( ξ ) (22)where φ a ( ξ ) are trilinear hexahedral basis functions and x ea are the coordinates of theeight element nodes. The solution is given by U h ( ξ ) = (cid:88) a =1 U ea φ a ( ξ ) (23)and the discretized weighting functions are given by W h ( ξ ) = (cid:88) a =1 W ea φ a ( ξ ) . (24)The Galerkin statement is therefore: Find U h = (cid:2) u h , P h , B h , ψ h (cid:3) T ∈ V h s.t. ∀ W h = (cid:2) w h , q h , v h , s h (cid:3) ∈ V h A (cid:0) W h , U h (cid:1) = (cid:0) W h , F h (cid:1) (25)where A ( · , · ) is given by (18).A very serious pitfall of the Galerkin approximation is that the numerical solution isvery rarely able to capture all of the details of the flow field. This is particularly true forturbulent flows wherein multiple scales are active simultaneously. Very often, however,one is interested only in relatively large scales of the flow field. Considering the largescales in isolation is entirely inadequate and, at the very least, the effect of the smallscales on the large scales must be taken into account. This is the principle behind largeeddy simulation; that is, mathematical models (called turbulence models) are developedthat model the effects of the small scales on the large scales. This is the topic of the nextsection.
3. Turbulence Models
Several new turbulence models were developed and implemented in [53]. The modelswere implemented in a Fourier-spectral setting. There are several differences in the formthat the models take depending on the numerical method used. In particular, the modelsrequire more terms when using non-orthogonal basis functions as is done in the finiteelement setting. The first set of models that we present are derived from the variatonalmultiscale formulation under the assumption of scale separation.6 .1. Variational Multiscale Formulation
With the assumption of scale separation, we will model the effects below a certainlength scale and numerically resolve the motion above that length scale. Mathematicallywe say that the motions are represented by the finite dimensional functions U h . Tobe more precise, we say that the solutions U h are optimal when they are defined via aprojection operator P h . That is, U h = P h U ∈ V h . (26)The projection operator P h projects the analytical solution onto the finite dimensionalsubspace in an optimal manner. The user determines the meaning of optimal. Examplesinclude ensuring that the finite dimensional space contains functions that are nodallyexact or functions that minimize some type of error norm. The choice of the projectionoperator P h induces a fine-scale projection operator P (cid:48) which allows the introduction ofthe fine scale solution, U (cid:48) = P (cid:48) U ∈ V (cid:48) . (27)Thereafter we introduce a decomposition of the solution field into coarse-scales and fine-scales, U = U h + U (cid:48) . (28)This procedure is the starting point of the variational multiscale formulation as intro-duced in [28] and extended to large eddy simulation in [5]. The decomposition in (28) isalso used for the weighting functions so that W = W h + W (cid:48) . (29)These decompositions imply a function space decomposition of the form V = V h ⊕ V (cid:48) . (30)We note that the choice of projection operator will determine the nature of the functionspace decomposition. In [53], the projection operators were selected to be the high-pass and low-pass sharp cutoff filters. These projection operators result in orthogonalsubspaces which in turn possess convenient properties that can be exploited. For example,inner products between functions in V h and V (cid:48) are zero due to the orthogonality property.In this work, we leave the projection operators unspecified. A careful consideration ofprojection operators giving rise to orthogonal subspaces in a finite element context forfluid mechanics is presented in [10]. A thorough discussion of projection operators is alsoconsidered in [29].Introducing the decompositions (28) and (29) into (15) results in two coupled prob-lems, A (cid:16) W h , U h + U (cid:48) (cid:17) = (cid:0) W h , F (cid:1) , ∀ W h ∈ V h (31) A (cid:16) W (cid:48) , U h + U (cid:48) (cid:17) = (cid:16) W (cid:48) , F (cid:17) , ∀ W (cid:48) ∈ V (cid:48) . (32)7he goal is to solve the first of these problems for U h and the second for U (cid:48) . Up to thispoint these two problems are exact. The variational multiscale formulation proceeds byfinding an approximate solution to (32). The approximate solution to (32) is (see [5] forthe hydrodynamic case, [2] and references therein for MHD), U (cid:48) ≈ − τ R (cid:0) U h (cid:1) . (33)In (33), R (cid:0) U h (cid:1) is the residual of the coarse scales and is given by R (cid:0) U h (cid:1) = R V (cid:0) U h (cid:1) ∇ · u h R I (cid:0) U h (cid:1) ∇ · B h (34)where R V (cid:0) U h (cid:1) = ρ ∂ u h ∂t + ∇ · (cid:16) N V (cid:0) U h (cid:1) − F V (cid:0) U h (cid:1)(cid:17) + ∇ P h − f V (35)is the momentum residual and R I (cid:0) U h (cid:1) = ∂ B h ∂t + ∇ · (cid:16) N I (cid:0) U h (cid:1) − F I (cid:0) U h (cid:1)(cid:17) + ∇ ψ h − f I (36)is the induction residual. The parameter τ is the stabilization matrix and is an algebraicapproximation to the inverse differential operator on the fine scales. Moreover, thestabilization matrix is an intrinsic grid time-scale consisting of a combination of advectiveand diffusive time-scales. In this work, it is represented as a diagonal matrix operator τ = diag (cid:0) τ V , τ V , τ V , τ V c , τ I , τ I , τ I , τ I c (cid:1) . The diagonal components are given by τ V = 1 (cid:114)(cid:16) C V t ρ ∆ t (cid:17) + ρ u h · Gu h + ( C µ µ ) (cid:107) G (cid:107) + C B ρµ (cid:107) B h (cid:107) (cid:107) G (cid:107) (37) τ I = 1 (cid:114)(cid:16) C I t ∆ t (cid:17) + u h · Gu h + B h · GB h + C λ λ (cid:107) G (cid:107) (38) τ V c = 1 C V t tr ( G ) τ V (39) τ I c = 1 C I t tr ( G ) τ I (40)where the components of G are given by G ij = ∂ξ k ∂x i ∂ξ k ∂x j (41)and tr ( G ) = G ii . (Summation on i implied) (42)8he constants in these expressions are O (1) and do not change upon mesh refinementor coarsening. Ultimately, this fact does have implications for the influence of the stabi-lization parameter on the solution. The constants were selected to be C V t = 1, C I t = 1, C µ = 1, C B = 1, and C λ = 1. The VMS statement is: Find U h ∈ V h s.t. ∀ W h ∈ V h A V (cid:0) W h , U h (cid:1) − (cid:16) ∇ w h , N V C (cid:0) U h , U (cid:48) (cid:1) + N V R ( U (cid:48) ) (cid:17) − (cid:0) ∇ · w h , P (cid:48) (cid:1) + (cid:0) ∇ w h , F (cid:48) V ( U (cid:48) ) (cid:1) − (cid:0) ∇ q h , u (cid:48) (cid:1) + A I (cid:0) W h , U h (cid:1) − (cid:16) ∇ v h , N I C (cid:0) U h , U (cid:48) (cid:1) + N I R ( U (cid:48) ) (cid:17) + (cid:0) v h , ∇ ψ (cid:48) (cid:1) + (cid:0) ∇ v h , F (cid:48) I ( U (cid:48) ) (cid:1) − (cid:0) ∇ s h , B (cid:48) (cid:1) = (cid:0) w h , f V (cid:1) + (cid:0) v h , f I (cid:1) . (43)In (43) the cross stresses for the momentum and induction equations are respectivelygiven by N V C (cid:0) U h , U (cid:48) (cid:1) = ρ (cid:0) u h ⊗ u (cid:48) + u (cid:48) ⊗ u h (cid:1) − µ (cid:0) B h ⊗ B (cid:48) + B (cid:48) ⊗ B h (cid:1) + 1 µ B h · B (cid:48) (44) N I C (cid:0) U h , U (cid:48) (cid:1) = (cid:0) u h ⊗ B (cid:48) − B (cid:48) ⊗ u h (cid:1) + (cid:0) u (cid:48) ⊗ B h − B h ⊗ u (cid:48) (cid:1) (45)while the Reynolds stress contributions to the momentum and induction equations aregiven by N V R ( U (cid:48) ) = ρ ( u (cid:48) ⊗ u (cid:48) ) − µ ( B (cid:48) ⊗ B (cid:48) ) + 12 µ (cid:107) B (cid:48) (cid:107) (46) N I R ( U (cid:48) ) = ( u (cid:48) ⊗ B (cid:48) ) − ( B (cid:48) ⊗ u (cid:48) ) . (47)In (43) F (cid:48) V and F (cid:48) I are the viscous stress tensors arising from the subgrid scales where F (cid:48) V ( U (cid:48) ) = µ (cid:16) ∇ u (cid:48) + ( ∇ u (cid:48) ) T (cid:17) (48) F (cid:48) I ( U (cid:48) ) = λ (cid:16) ∇ B (cid:48) − ( ∇ B (cid:48) ) T (cid:17) . (49)It was demonstrated in [58] in the context of hydrodynamic turbulence that using the firstorder approximation to U (cid:48) , (33), results in a negligible contribution from the Reynoldsstresses (46). Following this example, we neglect the second order correlations givenby (46) and (47). Note that to date there is no precedent for dropping the term 12 µ (cid:107) B (cid:48) (cid:107) in (46). We ran simulations both with and without this term and found no discernableeffects on the results and therefore concluded that its contribution is negligible. Usinglinear finite elements the inner products that involve gradients of the subgrid scales canbe shown to be zero through an additional integration by parts. The final finite elementVMS statement is: Find U h ∈ V h s.t. ∀ W h ∈ V h , A V (cid:0) W h , U h (cid:1) − (cid:16) ∇ w h , N V C (cid:0) U h , U (cid:48) (cid:1)(cid:17) − (cid:0) ∇ · w h , P (cid:48) (cid:1) − (cid:0) ∇ q h , u (cid:48) (cid:1) + A I (cid:0) W h , U h (cid:1) − (cid:16) ∇ v h , N I C (cid:0) U h , U (cid:48) (cid:1)(cid:17) + (cid:0) v h , ∇ ψ (cid:48) (cid:1) − (cid:0) ∇ s h , B (cid:48) (cid:1) = (cid:0) w h , f V (cid:1) + (cid:0) v h , f I (cid:1) . (50)9f the fine scale space was designed to be truly orthogonal to the finite element space thenall inner products between functions in the fine scale and coarse scale spaces would vanishidentically. This is the situation that is found when using a Fourier-spectral method asin [53]. For a detailed, general discussion of this property the reader is referred to [10]and references therein. We briefly review some developments in eddy viscosity models for MHD. The generalform of the models presented in this work include any of the eddy viscosity modelsreviewed in this section. In particular, the eddy viscosity models reviewed here aredescribed in [55], [38], and [53].
In [53], the VMS-based turbulence modeling approach was compared to the classicalSmagorinsky LES model for MHD developed in [55] as well as an extension to this classicmodel that takes into account fundamental MHD physics [38]. The variational form of theMHD equations inclusive of an eddy viscosity model is:
Find U h ∈ V h s.t. ∀ W h ∈ V h A (cid:0) W h , U h (cid:1) + M (cid:0) W h , U h ; c , h (cid:1) = (cid:0) W h , F (cid:1) (51)where the model term M ( · , · ; · , · ) seeks to model the effects of the fine scales on the coarsescales. This model term generally depends on the coarse-scale solution, some modelparameters c , and the mesh parameter h . The classical Smagorinsky eddy viscositymodel for MHD was introduced in [55]. In a variational context the model term takesthe form, M (cid:0) W h , U h ; c , h (cid:1) = (cid:0) ∇ w h , ρν T ∇ s u h (cid:1) + (cid:0) ∇ v h , λ T ∇ a B h (cid:1) (52)where ∇ s u h = 12 (cid:16) ∇ u h + (cid:0) ∇ u h (cid:1) T (cid:17) ∇ a B h = 12 (cid:16) ∇ B h − (cid:0) ∇ B h (cid:1) T (cid:17) . The turbulent viscosity ν T and turbulent magnetic diffusivity λ T are given by ν T = C S V h S , S = √ ∇ s u h : ∇ s u h (53) λ T = C S I h (cid:12)(cid:12)(cid:12)(cid:12)(cid:114) µ ρ J h (cid:12)(cid:12)(cid:12)(cid:12) , J h = 1 µ ∇ × B h . (54)When the coefficients C S V and C s I are determined via a dynamic procedure [22, 39] themodels are referred to as the dynamic Smagorinsky eddy viscosity (DSEV) model.An extension to the model introduced in [55] was introduced in [38]. We refer to thisextension as the alignment-based dynamic Smagorinsky eddy viscosity (DSEVA) model.The model has the same form as (52) with the only difference being in the definition of10he turbulent diffusivities. These are now given by ν T = C S V h (cid:115) abs (cid:18) ∇ s u h : µ ρ ∇ s B h (cid:19) (55) λ T = C S I h sgn (cid:0) J h · ω h (cid:1) (cid:115) abs (cid:18)(cid:114) µ ρ J h · ω h (cid:19) . (56)In (55) and (56), sgn ( · ) gives the sign of its argument, abs ( · ) gives the absolute valueof its argument, and ω h is the fluid vorticity. The coefficients are still determined viathe classic dynamic procedure. The DSEVA model takes into account the alignment ofthe velocity field and magnetic induction which is a fundamental characteristic of MHDflows. A new type of eddy viscosity model was proposed in [40, 53]. This eddy viscositymodel was motivated via the VMS formulation and has the same form as (51) with themodel in the form of (52). The eddy diffusivities are given by ν T = λ T = Ch (cid:114) | u (cid:48) | + 1 µ ρ | B (cid:48) | (57)where C is a universal constant. The constant is given by C = 23 √ C / K hk h (cid:112) − β − / (58)where h is the mesh size, C K = 2 . β is a constant denoting the range of the numerical subgrid solutions, and k h is thecutoff wavenumber. This expression is derived from a more general expression in [53]which was determined by matching the numerical and exact dissipation rates. In (58)we have assumed that the energy spectrum exhibits a k − / power law behavior in theinertial range where k is the wavenumber. We note that the true inertial range behavior inMHD is still undetermined and that numerical evidence suggests a nonuniversal characterfor the MHD energy spectrum [36]. Recent theoretical studies have indicated that thespectrum in MHD turbulence is actually k − / [8]. However, a k − / energy spectrumleads to a universal expression for C that does not depend on the dissipation rate. Inlight of this, we choose C as given by (58) and acknowledge that modificiations to thederivation of C may be appropriate in the future.In a spectral setting the cutoff wavenumber is easily determined to be k h = π/h . Ina finite element setting we choose k h = π/ (2 h ) by reasoning that it requires four linearelements to resolve the smallest wavelength in the domain. Moreover, we choose β = 3 / β < / β = 3 / k h and 3 k h /
2. With these parameters and assumptions, we obtain C ≈ . u (cid:48) and B (cid:48) are given by (33) and are therefore based off of the residualof the coarse scales. Hence this model is automatically dynamic and it is conjecturedthat the constant C does not need to be adjusted dynamically. Because the expressions11or the eddy diffusivities involve formulae for the subgrid scales derived from the VMSformulation we refer to this eddy viscosity model as the residual-based eddy viscosity(RBEV) model. For large values of fluid and magnetic Reynolds numbers, Re and Rm , the second-order correlations ( N V R ( U (cid:48) ) and N I R ( U (cid:48) )) are important and we would expect that thepure VMS models given by (50) do not sufficiently account for the subgrid scales. Moti-vated by this, a mixed model for MHD was proposed in [53]. In the present formulation,this model takes the form, A V (cid:0) W h , U h (cid:1) − (cid:16) ∇ w h , N V C (cid:0) U h , U (cid:48) (cid:1)(cid:17) − (cid:0) ∇ · w h , P (cid:48) (cid:1) − (cid:0) ∇ q h , u (cid:48) (cid:1) + M V (cid:0) W h , U h ; c V , h (cid:1) + A I (cid:0) W h , U h (cid:1) − (cid:16) ∇ v h , N I C (cid:0) U h , U (cid:48) (cid:1)(cid:17) + (cid:0) v h , ∇ ψ (cid:48) (cid:1) − (cid:0) ∇ s h , B (cid:48) (cid:1) + M I (cid:0) W h , U h ; c I , h (cid:1) = (cid:0) w h , f V (cid:1) + (cid:0) v h , f I (cid:1) (59)where the additional model terms M V ( · , · ; · , · ) and M I ( · , · ; · , · ) are eddy viscosity modelsof the form (52). The eddy viscosities in the models could be given by any of (53)-(56)or (57). When the mixed model uses the RBEV model to account for the second ordercorrelations it is a fully residual-based model. For completeness, we present a general expression for the models. Note that theformulation contains parameters that may be adjusted depending on the desired model.Specific values of the parameters are provided in Table 1 that result in the particularmodels developed. A V (cid:0) W h , U h (cid:1) − b VVMS (cid:104)(cid:16) ∇ w h , N V C (cid:0) U h , U (cid:48) (cid:1)(cid:17) + (cid:0) ∇ · w h , P (cid:48) (cid:1)(cid:105) − c VVMS (cid:0) ∇ q h , u (cid:48) (cid:1) + b VEVM (cid:0) ∇ w h , ρν T ∇ s u h (cid:1) + A I (cid:0) W h , U h (cid:1) − b IVMS (cid:104)(cid:16) ∇ v h , N I C (cid:0) U h , U (cid:48) (cid:1)(cid:17) − (cid:0) v h , ∇ ψ (cid:48) (cid:1)(cid:105) − c IVMS (cid:0) ∇ s h , B (cid:48) (cid:1) + b IEVM (cid:0) ∇ v h , λ T ∇ a B h (cid:1) = (cid:0) w h , f V (cid:1) + (cid:0) v h , f I (cid:1) . (60)Note that (60) includes the DSEV and DSEVA models as possible choices by selectingthe eddy diffusivities according to (53)-(54) or (55)-(56) respectivly. Moreover, a fullyresidual-based model can be attained by selecting the eddy diffusivities according to (57).Finally, the parameter c V , IVMS deserves comment. This parameter is necessary dependingon the type of finite element basis functions that are being used. Classically, this termwas introduced to circumvent the inf-sup condition in mixed problems. That is, it allowsequal order elements for all of the fields in the problem. Note that this term comes forfree when using VMS-based methods. However, when using pure EVMs one needs to be12 able 1: A summary of the values of the parameters b i and c and corresponding models. Model b V , IVMS b V , IEVM c V , IVMS
VMS 1 0 1MM 1 1 1EVM 0 1 1 equal order elements0 any other timecareful and introduce the term only when appropriate as indicated in Table 1. For moredetails on the inf-sup condition and the addition of the stabilization term the reader isreferred to [30].
4. Summary of Fully-implicit AMG Preconditioned Newton-Krylov Solution
The incompressible restive MHD model, (12) - (14), supports a number of fast time-scales. These include elliptic effects due to the divergence constraints that support infi-nite speed pressure waves and enforce the solenoidal involution of the magnetic field, fasttime-scale parabolic effects due to momentum and magnetic diffusion, and Alfv´en wavepropagation. This system can exhibit behavior that is characterized by widely separatedtime-scales as well as overlapping time-scales. These characteristics make the robustand efficient iterative solution of these systems very challenging. In this context fully-implicit methods are an attractive choice that can often provide unconditionally-stabletime integration techniques. These methods can be designed with various types of sta-bility properties that allow robust integration of multiple-time-scale systems without therequirement to resolve the stiff modes of the system (which are not of interest since theydo not control the accuracy of time integration). In the computations presented in thisstudy, a fully-implicit multistep backward differentiation method of order three (BDF3)is employed to accurately and efficiently integrate the time evolution of the MHD system.More details on the fully-implicit integration and example order-of-accuracy verificationresults can be found in [50, 51].
The result of a fully-implicit solution technique is the development of very large-scale, coupled highly nonlinear algebraic system(s) that must be solved. Therefore, thesetechniques place a heavy burden on both the nonlinear and linear solvers and requirerobust, scalable, and efficient iterative solution methods. In this study Newton-basediterative nonlinear solvers ([17]) are employed to solve the challenging nonlinear systemsthat result from this application. These solvers can exhibit quadratic convergence ratesindependently of the problem size when sufficiently robust linear solvers are available.For the latter, we employ Krylov iterative techniques. A Newton-Krylov (NK) methodis an implementation of Newton’s method in which a Krylov iterative solution technique13s used to approximately solve the linear systems, J k s k + = − F k , that are generated ateach step of Newton’s method. For efficiency, an inexact Newton method ([17], [19]) isusually employed, whereby one approximately solves the linear systems generated in theNewton method by choosing a forcing term η k and stopping the Krylov iteration whenthe inexact Newton condition, (cid:107) F k + J k s k + (cid:107) ≤ η k + (cid:107) F k (cid:107) is satisfied. The particularKrylov method that is used in this study is a robust non-restarted GMRES methodthat is capable of iteratively converging to the solution of very large non-symmetriclinear systems such as developed from the discretized resistive MHD system provided asufficiently robust and scalable preconditioning method is available [50, 51]. As a scalable preconditioner for the discretized VMS resistive MHD system a fully-coupled algebraic multigrid method (AMG) is employed. In general AMG methods aresignificantly easier to implement and integrate with complex unstructured mesh simu-lation software than geometric multigrid methods [57, 48, 49]. AMG preconditionersassociate a graph with the matrix system being solved. Graph vertices correspond tomatrix rows for scalar PDEs, while for PDE systems it is natural to associate one vertexwith each nodal block of unknowns (e.g. velocities, pressure, magnetic field and the La-grange multiplier at a particular grid point). Note that since the VMS formulation allowsthe use of non-mixed FE interpolation (in this case for example all first order Lagrangeinterpolation on hexes) a simple block matrix (8x8) structure exists if all unknowns areordered consecutively at each FE vertex and a graph edge exists between vertex i and j if there is a nonzero in the block matrix which couples i ’s rows with j ’s columns or j ’s rows with i ’s columns. The specific AMG coarsening strategy that is employed isa parallel greedy graph aggregation technique which attempts to make an aggregate bytaking an unaggregated point and grouping it with all of its neighbors. Thus, it tendsto coarsen by a factor of three in each coordinate direction when applied to a standarddiscretization matrix obtained by a compact stencil on a regular mesh. In addition,a dynamic load-balancing package, Zoltan [18], is used to repartition coarsened opera-tors across processors. This generally improves parallel performance and gives betteraggregates on the next coarser level which is obtained by again applying the parallel ag-gregation technique. A more complete discussion of this strategy and the implementationfor this method can be found in [51]. The fully-coupled AMG preconditioner and theNewton-GMRES method described above have been demonstrated to provide a robustand scalable nonlinear/linear iterative solution strategy for a number of applications thatinclude MHD (see e.g. [50, 51, 37] and the references contained therein). To complete the brief description of the numerical method that is employed in thisstudy the convergence criteria that are used at each time-step to determine convergenceof the nonlinear/linear iterative methods must be stated. For the nonlinear solver twoconvergence criteria are employed. The first is a requirement of sufficient reduction in therelative nonlinear residual norm, (cid:107) F k (cid:107) / (cid:107) F o (cid:107) < − . In general, this requirement is setto be easily satisfied. The second convergence criterion is based on a sufficient decreaseof a weighted norm of the relative size of the Newton update vector. The latter criterion14equires that the correction, ∆ χ ki , for any grid point unknown, χ i , is small compared toits magnitude, (cid:12)(cid:12) χ ki (cid:12)(cid:12) , and is given by (cid:118)(cid:117)(cid:117)(cid:116) N u N u (cid:88) i =1 (cid:20) | ∆ χ i | ε r | χ i | + ε a (cid:21) < , where N u is the total number of unknowns, ε r is the relative error tolerance betweenthe variable correction and its magnitude, and ε a is the absolute error tolerance of thevariable correction, which sets the magnitude of components that are to be considered tobe numerically zero. In this study, the relative-error and absolute-error tolerance are 10 − and 10 − respectively. In practice this criteria controls convergence of the nonlinear solverto sufficient accuracy. For example, in the high Reynolds number 256 VMS simulation(see section 5.2) the maximum relative nonlinear norm decrease was O (10 − ). Finally ineach Newton sub-step the linear system is solved to a relative tolerance of η = 10 − withthe knowledge that the overall accuracy of the solution for each time step is controlledas described above for convergence of the Newton step.
5. Results and Discussion
The models are evaluated on the Taylor-Green vortex generalized to MHD as de-scribed in [44] and [36]. This problem begins with smooth initial conditions for thevelocity field and magnetic induction, u ( x , t = 0) = u sin ( x ) cos ( y ) cos ( z ) − cos ( x ) sin ( y ) cos ( z )0 . (61) B ( x , t = 0) = B cos ( x ) sin ( y ) sin ( z )sin ( x ) cos ( y ) sin ( z ) − x ) sin ( y ) cos ( z ) . (62)The domain is a periodic box of size [ − π, π ] . The energy of each field is determinedfrom their respective energy spectrum, K T , V , I = (cid:90) E T , V , I ( k ) d k (63)where K T , V , I represents the total energy, kinetic energy, or magnetic induction energy ata certain time and E T , V , I ( k ) is the total energy spectrum, the kinetic energy spectrum,or the magnetic energy spectrum at wavenumber k . We define the Reynolds number Re and magnetic Reynolds number Rm as Re = u rms L V ν (64) Rm = u rms L I λ . (65)15he characteristic length scales L V and L I are computed as L V , I = (cid:82) k E V , I ( k ) d k (cid:82) E V , I ( k ) d k . (66)The ratio of molecular viscosity to magnetic diffusivity, the magnetic Prandtl number( P r = Rm/Re ), is unity. The initial velocity and magnetic fields are scaled so that thetotal energy, which is the sum of the kinetic and magnetic energies, is initially equal to0 .
25. The initial state also has the kinetic and magnetic energies equal to each other.For details on properties of this flow field, the reader is referred to [36, 44, 14] andreferences therein. In particular, [36] discusses the nonuniversality of MHD by consideringmultiple Taylor-Green initial conditions and demonstrating that these different initialconditions lead to different energy cascades. In [14], the energy spectrum of the particularTaylor-Green vortex considered in this work is discussed as well as some implications ofnonuniversality for turbulence modeling in MHD. However, we also note the interestingconclusions in [1] on locality of MHD turbulence which is an essential ingredient for auniversal behvior of the energy cascade.
We first consider the case of Re = Rm = 1800 at the peak of dissipation. For thisproblem, the peak of dissipation occurs at t ∈ [3 . , . N VSUPG (cid:0) U h , U (cid:48) (cid:1) = ρ u (cid:48) ⊗ u h N ISUPG (cid:0) U h , U (cid:48) (cid:1) = (cid:0) u h ⊗ B (cid:48) − B (cid:48) ⊗ u h (cid:1) while the Reynolds stresses are neglected entirely. Figure 1 compares the total energyspectrum at t = 4 from the SUPG stabilization to the VMS and mixed models. Allenergy spectra have been averaged over neighboring modes to help remove the nodaloscillations that are typically present in energy spectra obtained from the Taylor Greenvortex. From these results it is clear that more than a straightforward SUPG stabilizationis required to sufficiently dissipate the unresolved highest wavenumber modes. Even thelow wavenumber modes are not accurately represented. This can be attributed to the factthat we are plotting the total energy and that the SUPG formulation does not accountfor cross coupling of subgrid fields. We note that with mesh refinement the SUPG resultswould improve because the effects of the subgrid fields would be diminished. Indeed, evenwith the VMS and mixed models there is still a slight pile-up of energy. Note, however,that the mixed model is able to dissipate the energy in the last modes despite the slightenergy pile-up. This is the desired behavior although we will have further comments onthe energy pile up in subsequent sections. 16 k -4 -3 -2 -1 E T ( k ) DNS: VMS : MM : SUPG : Figure 1: The total energy spectrum for the low Re case. Here the SUPG stabilized method is comparedwith the VMS model, the mixed model, and to a DNS simulation.
In Figure 2 we consider the time evolution of MHD energies for various LES meshresolutions ranging from N = 32 to N = 256. The LES results are compared to a spectralDNS at N = 512. It is observed that the models qualitatively converge monotonicallywith mesh refinement to the DNS results. t Plot LocationsUpper Left: K T ( t ) Upper Right: K I ( t ) Bottom: K V ( t ) Line Key
SpectralDNS: MM: MM: MM: MM: MM: MM: Figure 2: Time evolution of MHD energies from the mixed model (MM) LES MHD turbulence modelat various mesh resolutions. N = 32, N = 64, N = 128, and N = 256. The DNSresults have been low-pass filtered in order to compare with the FE LES solutions. Thishowever is not a straightforward comparison since the FE simulations clearly have amuch lower convergence rate (i.e. formally second-order) than the spectral method thatis exponentially convergent. It does however provide some qualitative assessment of therelative convergence of this class of models for this problem compared to the spectralDNS approximation. We first note that the solutions improve considerably with meshrefinement. K T DNS: VMS: MM: DNS: VMS: MM: t K T DNS: VMS: MM: t DNS: VMS: MM: Figure 3: Time evolution of total energy. The DNS results have been low-pass filtered such that thenumber of Fourier modes matches the number of finite elements in the LES model. K I DNS: VMS: MM: DNS: VMS: MM: t K I DNS: VMS: MM: t DNS: VMS: MM: Figure 4: Time evolution of magnetic energy. The DNS results have been low-pass filtered such that thenumber of Fourier modes matches the number of finite elements in the LES model. .020.040.060.080.100.12 K V DNS: VMS: MM: DNS: VMS: MM: t K V DNS: VMS: MM: t DNS: VMS: MM: Figure 5: Time evolution of kinetic energy. The DNS results have been low-pass filtered such that thenumber of Fourier modes matches the number of finite elements in the LES model.
Indeed, at N = 64 both models are in quite good agreement for times up to about t = 4, the peak of dissipation, with the mixed model performing marginally better. Weobserve, however, that both models have a tendency to overpredict the total energythrough the peak of dissipation for coarse discretizations. This is a result of the VMSapproach and the fact that we did not ensure the orthogonality of the subgrid scaleswith the resolved scales. In the spectral setting one can prove that the VMS method isglobally dissipative [58] (although local backscatter is possible). Such a demonstrationis facilitated by the orthogonality of two spaces involved. We note that for N = 128and beyond the over-prediction near the peak of dissipation is no longer present and theresults correspond to the DNS solution quite well.Next, we consider the total energy spectrum for relatively coarse FE discretizations.Figure 6 shows the total energy spectrum at t = 4 for the VMS model and mixed modelcompared to the DNS results. The LES results agree with the DNS results well for thelarge scales as expected. Both models exhibit a slight pile-up of energy but this energyis dissipated in the last mode or two as it should be. Further details of the solutionbehavior are uncovered in Figure 7 by considering the kinetic and magnetic energy spectraseparately, also at t = 4. The DNS data for the magnetic energy spectrum is matchedvery well. The models do not perform as well for coarse discretizations in replicating thekinetic energy. Indeed, we see that the observed pile-up of total energy in the small scalesstems directly from the kinetic energy. This implies that the models do not perform aswell for the momentum equation as for the induction equation. This lack of performanceis most likely due to highly active subgrid velocity scales. The good agreement in thetotal energy spectrum is mainly due to the fact that this flow field is dominated by themagnetic field. The energy in the velocity field is almost an order of magnitude less thanthat in the magnetic field. An interesting point is that the performance discrepancybetween the momentum and induction equation calls attention to the fact that we haveselected our eddy viscosities such that the turbulent magnetic Prandtl number is unity.It may be necessary to adjust the eddy viscosity portion of the model so that ν T (cid:54) = λ T .19 k -4 -3 -2 -1 E T ( k ) DNS: VMS : MM : Figure 6: The total energy spectrum for the VMS model and mixed model compared to the DNSsimulation at t = 4. k -4 -3 -2 -1 E I ( k ) DNS: VMS : MM : k E V ( k ) DNS: VMS : MM : Figure 7: Magnetic and kinetic energy spectra at t = 4. We can gain further insight into the nature and performance of the mixed model byinspecting the time evolution of the eddy viscosity. In Figure 8 we consider the timeevolution of the cell-averaged eddy viscosity at three different spatial resolutions. Weobserve that the eddy viscosity displays an automatic dynamic behavior. Initially, theeddy viscosity is zero because there are no unresolved scales in the initial conditionsand hence the residual is zero. However, as soon as other scales become active theeddy viscosity “turns on”. It ramps up to a peak value at the peak of dissipationwhich is physically consistent; that is, the eddy viscosity should work hardest aroundthe maximum of energy dissipation. We also note that with spatial mesh-refinement theeddy viscosities become less active. This is a direct result of the residual-dependenceof the eddy viscosities; the numerical residual becomes smaller as the mesh is refined20ecause the numerical solution is able to capture more scales of the flow field. t ν T , λ T
1e 4 Figure 8: Time evolution of the MHD eddy viscosity at three different resolutions. The inherent dynamicnature of the eddy viscosity is observed.
Although encouraging, we also know from Figure 6 that the eddy viscosity shouldbe slightly larger in order to dissipate the energy from the small scales more efficiently.Finally, we note that at the very beginning of the simulation the eddy viscosity exhibitsa small overshoot. This is the result of the finite element basis. The initial condition isnot contained in the finite element space. Thus, the residual is immediately nonzero. Ina Fourier spectral method of the Taylor-Green vortex in hydrodynamics this overshootis not present [40].Figures 9-11 show energy spectra at various resolutions at t = 4. k -7 -6 -5 -4 -3 -2 -1 E T ( k ) DNS: VMS : VMS : VMS : VMS : k DNS: MM : MM : MM : MM : Figure 9: Convergence of total MHD energy spectrum with mesh refinement at t = 4. In all cases, the LES results converge to the DNS results. At N = 128 and N = 25621he LES results are quite good. We note that for the VMS model there is a pile-up ofenergy in the last few modes of the kinetic energy spectrum at N = 256. However, thispile up is O (cid:0) − (cid:1) and is quite small. It is most likely due to the fact that the subgridscales at such high wavenumbers are very small and the models have some difficultyaccurately representing such small scales. The convergence in the case of the mixedmodel is seen to be much better with the high wavenumber modes being controlled muchmore effectively. The correspondence of the N = 256 MM discretization with the 512mode spectral DNS solution shows the remarkable ability of this model to reproduce thespectrum accurately. k -7 -6 -5 -4 -3 -2 -1 E I ( k ) DNS: VMS : VMS : VMS : VMS : k DNS: MM : MM : MM : MM : Figure 10: Convergence of magnetic energy spectrum with mesh refinement at t = 4. k -7 -6 -5 -4 -3 -2 -1 E V ( k ) DNS: VMS : VMS : VMS : VMS : k DNS: MM : MM : MM : MM : Figure 11: Convergence of kinetic energy spectrum with mesh refinement at t = 4. .2. High Reynolds Number Next we turn our attention to the same flow field at higher Reynolds numbers. TheReynolds numbers near the peak of dissipation are Re = Rm = 9000. Figure 12 presentsthe time evolution of energies obtained from DNS results and the mixed model at twodifferent resolutions. The results are quite good especially for N = 256 for which theLES solutions are only slightly overly dissipative for later times. K T ( t ) K I ( t ) t K V ( t ) Figure 12: Comparison of time evolution of total, magnetic, and kinetic energies between DNS results(brown circles) and mixed model results at resolutions of 128 (orange lambdas) and 256 (green clubs). We also consider, once again, plots of the total energy spectrum and compare thespectra from the models to DNS spectra from [36] in Figures 13 and 14. This time,however, the spectra have been averaged over certain time intervals which are indicatedin the respective figures. We have performed LES runs at various spatial resolutionsat N = 32 , , , times as many modes in the expansion compared to the number ofelements in the LES results. Given this difference in resolution, and the exponentialconvergence of the spectral method the correspondence is remarkable. Further, the abilityof the reasonably coarse LES type models to approximate the spectrum is encouragingin the context of controlling the computational costs of such simulations. Figure 13compares the performance of the LES models at N = 256 to the DNS results. Weobserve that both the VMS and mixed models perform quite well in the inertial rangeat both times. However, when the spectrum is averaged around the peak of dissipation(leftmost plot in Figure 13) we notice that the VMS model results in some energy pile-upin the smallest scales. We note that the energy pile-up is on the order of 10 − and istherefore rather small. Morever, in the last one or two modes, the VMS model beginsto disspate energy, albeit inadequately. On the other hand, the mixed model performsadmirably and matches the DNS spectrum nearly identically. The only difference betweenthe mixed model and VMS model is that the mixed model attempts to account for the23econd order correlations that the VMS model neglected. For high Reynolds numbersthese second order correlations are expected to be significant. One would expect thatnear the peak of dissipation these correlations are particularly active. This explains whythe spectrum from the VMS model exhibits some energy pile-up while the mixed modeldissipates energy at precisely the correct rate. k -5 -4 -3 -2 -1 h E T ( k ) i k − t ave = [ . , . ] DNS: VMS: MM: k -4 -3 -2 -1 k − t ave = [ . , . ] DNS: VMS: MM: Figure 13: Performance of VMS and mixed models for high Reynolds numbers at two different times.The total energy spectra have been averaged around the times indicated in the figure.
The rightmost plot in Figure 13 shows the spectra averaged around a later timethat is not near the peak of dissipation. We observe that the VMS model and mixedmodel both perform well but that the mixed model performs slightly better. The biggestdifference between this plot and the similar plot near the peak of dissipation is that theVMS model does not result in a pile-up of energy in the smallest scales. This makessense in light of the fact that the dissipation rate is not as intense and that the higherorder correlations are therefore less active. In addition, we note that the LES resultshave the correct inertial range scaling of k − . This is particularly interesting becausein the derivation of the mixed model constant we assumed that the energy spectrumwould follow a Kolmogorov scaling law. The good results are a further indication of thedynamic nature of the eddy viscosity portion of this model. That is, the model adapts asnecessary to the dynamics of the flow field. The exact value of the mixed model constantdoes not appear to be crucial.Finally, Figure 14 shows plots of the VMS and mixed models at various spatial resolu-tions in order to demonstrate qualitative convergence with mesh refinement. We observethat the 32 simulations are excessively dissipative as expected from linear finite elements.However, with subsequent refinements the spectrum matches the DNS better, especiallyin the inertial range. Each refinement shows improvement in agreement between the LESmodels and the DNS result. Clearly the correspondence of the spectrum for the mixedmodel is very impressive and indicates that this MHD LES model is performing verywell. 24 k -3 -2 -1 h E T ( k ) i VMS Model
DNS k Mixed Model
DNS Figure 14: Convergence of FEM VMS and mixed models demonstrated using 32 - 256 modes. Thespectra were averaged around t ∈ [6 , .
6. Conclusions
In this work, we have introduced a class of finite element turbulence models forincompressible magnetohydrodynamics. These models are derived from the variationalmultiscale formulation wherein the fine scales are assumed to be directly proportionalto the residual of the coarse scales. The general expression for the models is providedin (60) and specific model choices are summarized in Table 1. We briefly summarize thenew models:1. The variational multiscale formulation is given by (50) where the fine scales aredetermined from (33).2. Eddy viscosity models are given by (51) and (52). Different choices of the eddydiffusivites ν T and λ T result in different models. Dynamic Smagorinsky modelswould result from selecting either (53) and (54) for ν T and λ T respectively whilethe alignment-based dynamic Smagorinsky model would result from (55) and (56).We have introduced an alternative eddy viscosity model which we refer to as theresidual-based eddy viscosity model. In the RBEV model, the eddy diffusivitiesare proprotional to the unresolved velocity and magnetic field as determined viathe VMS method. In this way, the eddy viscosity model is fully residual based andtherefore results in a consistent numerical method. Moreover, we have assumed thatthe constant appearing in this model does not need to be determined dynamically.The RBEV model uses (57) for the eddy viscosities.3. The mixed model is given by (59). In this work we have used the VMS modelsalong with the RBEV model.The new models were tested on the Taylor-Green vortex generalized to MHD at twodifferent Reynolds numbers. Comparisons were made to existing DNS data in terms ofthe evolution of energies as well as energy spectra at various times. The mixed modelslightly outperformed the pure VMS model as expected. All models struggled slightly inreproducing the DNS energy spectra at the peak of dissipation except for the N = 256,high Reynolds number case in which the mixed model dissipated energy at the correct25ate and matched the DNS data very well. The models were able to eventually dissipatethe energy even in cases where slight energy pile-up was observed. It is interesting tonote that the models were able to capture the correct inertial range behavior for this flowfield even though it does not display a traditional Kolmogorov energy spectrum. This isin spite of the fact that the RBEV model consant was derived assuming a Kolmogorovspectrum. This behavior is not entirely surprising, however, considering that the eddyviscosity models are inherently dynamic and adjust automatically to the flow field.Certain numerical properties of the new models were explored in addition to physicalconsiderations. In particular, the nature of the eddy viscosities was considered as wellas convergence properties of the method. The eddy viscosities displayed their expecteddynamic behavior, peaking at the maximum of physical dissipation and decreasing there-after. With mesh refinement, the eddy viscosities decreased in strength because the roleof the subgrid scales was diminished. The numerical solution also improved with meshrefinement.Considering the encouraging performance of the new proposed models, nevertheless,there are still several potential areas of improvement. The RBEV constant could beimproved although it does not seem likely that a true universal constant for MHD can befound due to the nonuniversality of decaying MHD turbulence. This lack of universalitymay not be a significant hindrance to these models however due to their dynamic nature.Already, as shown in this work, the models are able to capture the correct inertial rangebehavior. The RBEV constant was derived in a spectral context and improvements tothis constant may involve considering alternative derivations in terms of finite elements.In general, the models appear to be quite sensitive to their constant parameters. Thisis especially true in the definition of the stabilization matrix found in (37) - (40). Aparameter sensitivity analysis may shed some light on the role that these constants playin affecting the solutions to different problems.The VMS method is not guaranteed to be globally dissipative; this feature depends onthe design of the finite element function spaces. The models developed in the present workdo not guarantee this global dissipativity and this fact is seen in the coarse discretizationplots of the time evolution of total MHD energy in Figure 3. Namely, just prior to thepeak of dissipation, the total MHD energy slightly exceeds the initial energy. A carefulconsideration of the finite element spaces as is done in [10] may help to overcome thisissue.There are several avenues left to explore in the future related to physical and numericalconsiderations. The optimal blending of parameters in the general mixed model presentedin (60) is still undetermined. For now, all of the models appearing in (60) are weightedequally. Some evidence exists that the eddy viscosity portion of the complete modelshould be weighted less than the VMS portion. This is based upon arguments concerningthe contribution of cross stresses and Reynolds stresses in hydrodynamic turbulence [35].One nice feature of the VMS-based models is the fact that they permit local backscat-ter. In MHD turbulence this feature is quite important as it permits the transfer ofenergy from the turbulent velocity fluctuations into the large scale magnetic field. Thisphenomenon has implications in models of the geodynamo. We have begun studies toquantify this phenomenon in relation to the turbulence models developed in this workboth in a spectral and finite element setting. Future work will also include the applica-tion of these models to complex geometries that include modeling liquid metals flows forfusion applications [52] and geodynamo simulations [34].26 cknowledgements This work was initiated with support from the Department of Energy (DOE) Of-fice of Science Graduate Fellowship (SCGF) under Contract No. DE-AC05-06OR23100.Support from NSF-DMS grant 1147523 is gratefully acknoweledged during preparationof this paper. Additionally the work of Shadid, Pawlowski, Cyr and Smith was par-tially supported by the DOE Office of Science Applied Mathematics Program at SandiaNational Laboratories under contract DE-AC04-94AL85000.27
1] H. Aluie and G. L. Eyink. Scale locality of magnetohydrodynamic turbulence.
Physical reviewletters , 104:81101, 2010.[2] Santiago Badia, Ramon Codina, and Ramon Planas. On an unconditionally convergent stabilized fi-nite element approximation of resistive magnetohydrodynamics.
Journal of Computational Physics ,234:399–416, 2013.[3] J. Baerenzung, H. Politano, Y. Ponty, and A. Pouquet. Spectral modeling of magnetohydrodynamicturbulent flows.
Physical Review E , 78:026310, 2008.[4] G. K. Batchelor.
The Theory of Homogeneous Turbulence . Cambridge: Cambridge UniversityPress, 1982.[5] Y. Bazilevs, V.M. Calo, J.A. Cottrell, T.J.R. Hughes, A. Reali, and G. Scovazzi. Variationalmultiscale residual-based turbulence modeling for large eddy simulation of incompressible flows.
Computer Methods in Applied Mechanics and Engineering , 197:173–201, 2007.[6] A. Beresnyak. Spectral slope and kolmogorov constant of mhd turbulence.
Physical Review Letters ,106:75001, 2011.[7] D. Biskamp.
Magnetohydrodynamic Turbulence . Cambridge: Cambridge University Press, 2003.[8] S. Boldyrev. Spectrum of magnetohydrodynamic turbulence.
Physical review letters , 96:115002,2006.[9] Luis Chac´on. A non-staggered, conservative, finite-volume scheme for 3d implicit extended magne-tohydrodynamics in curvilinear geometries.
Computer Physics Communications , 163(3):143–171,2004.[10] R. Codina. Stabilized finite element approximation of transient incompressible flows using orthog-onal subscales.
Computer Methods in Applied Mechanics and Engineering , 191:4295–4321, 2002.[11] R. Codina and N. Hern´andez-Silva. Stabilized finite element approximation of the stationarymagneto-hydrodynamics equations.
Computational Mechanics , 38:344–355, 2006.[12] R. Codina, J. Principe, O. Guasch, and S. Badia. Time dependent subscales in the stabilized finiteelement approximation of incompressible flow problems.
Computer Methods in Applied Mechanicsand Engineering , 196:2413–2430, 2007.[13] Ramon Codina and Noel Hern´andez. Approximation of the thermally coupled mhd problem usinga stabilized finite element method.
Journal of Computational Physics , 230(4):1281–1303, 2011.[14] V. Dallas and A. Alexakis. The origins of k − spectrum in the decaying taylor-green magnetohydro-dynamic turbulent flows. Phys. Rev. E , 88:053014, Nov 2013. doi: 10.1103/PhysRevE.88.053014.[15] P.A. Davidson.
An Introduction to Magnetohydrodynamics , volume 25. New York: CambridgeUniversity Press, 2001.[16] Andreas Dedner, Friedemann Kemm, Dietmar Kr¨oner, C-D Munz, Thomas Schnitzer, and MatthiasWesenberg. Hyperbolic divergence cleaning for the mhd equations.
Journal of ComputationalPhysics , 175(2):645–673, 2002.[17] John E Dennis Jr and Robert B Schnabel.
Numerical methods for unconstrained optimization andnonlinear equations , volume 16. Siam, 1996.[18] Karen Devine, Erik Boman, Robert Heaphy, Bruce Hendrickson, and Courtenay Vaughan. Zoltandata management services for parallel dynamic applications.
Computing in Science & Engineering ,4(2):90–96, 2002.[19] Stanley C Eisenstat and Homer F Walker. Choosing the forcing terms in an inexact newton method.
SIAM Journal on Scientific Computing , 17(1):16–32, 1996.[20] Leopoldo P Franca, Sergio L Frey, and Thomas JR Hughes. Stabilized finite element methods: I.application to the advective-diffusive model.
Computer Methods in Applied Mechanics and Engi-neering , 95(2):253–276, 1992.[21] J-F Gerbeau. A stabilized finite element method for the incompressible magnetohydrodynamicequations.
Numerische Mathematik , 87(1):83–111, 2000.[22] M. Germano, U. Piomelli, P. Moin, and W. H. Cabot. A dynamic subgrid-scale eddy viscositymodel.
Physics of Fluids A: Fluid Dynamics , 3:1760, 1991.[23] H. Goedbloed and S. Poedts.
Principles of Magnetohydrodynamics , volume 74. New York: Cam-bridge University Press, 2006.[24] H. Goedbloed, R. Keppens, and S. Poedts.
Advanced Magnetohydrodynamics , volume 1. New York:Cambridge University Press, 2010.[25] N.E. Haugen and A. Brandenburg. Hydrodynamic and hydromagnetic energy spectra from largeeddy simulations.
Physics of Fluids , 18:075106, 2006.[26] D.D. Holm. Lagrangian averages, averaged lagrangians, and the mean effects of fluctuations in fluiddynamics.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 12:518–530, 2002.[27] Thomas JR Hughes, Luca Mazzei, and Kenneth E Jansen. Large eddy simulation and the variational ultiscale method. Computing and Visualization in Science , 3(1-2):47–59, 2000.[28] T.J.R. Hughes. Multiscale phenomena: Green’s functions, the dirichlet-to-neumann formulation,subgrid scale models, bubbles and the origins of stabilized methods.
Computer Methods in AppliedMechanics and Engineering , 127:387–401, 1995.[29] T.J.R. Hughes and G. Sangalli. Variational multiscale analysis: The fine-scale green’s function, pro-jection, optimization, localization, and stabilized methods.
SIAM Journal on Numerical Analysis ,45:539–557, 2008.[30] T.J.R. Hughes, L. P. Franca, and M. Balestra. A new finite element formulation for computa-tional fluid dynamics: V. circumventing the babuˇska-brezzi condition: A stable petrov-galerkinformulation of the stokes problem accommodating equal-order interpolations.
Computer Methodsin Applied Mechanics and Engineering , 59:85–99, 1986.[31] T.J.R. Hughes, L.P. Franca, and G.M. Hulbert. A new finite element formulation for computationalfluid dynamics: Viii. the galerkin/least-squares method for advective-diffusive equations.
ComputerMethods in Applied Mechanics and Engineering , 73:173–189, 1989.[32] T.J.R. Hughes, G.R. Feij´oo, L. Mazzei, and J.B. Quincy. The variational multiscale method–aparadigm for computational mechanics.
Computer Methods in Applied Mechanics and Engineering ,166:3–24, 1998.[33] B. Knaepen and P. Moin. Large-eddy simulation of conductive flows at low magnetic reynoldsnumber.
Physics of Fluids , 16:1255, 2004.[34] Masaru Kono and Paul H Roberts. Recent geodynamo simulations and observations of the geo-magnetic field.
Reviews of Geophysics , 40(4):4–1, 2002.[35] R.H. Kraichnan. Eddy viscosity in two and three dimensions.
Journal of Atmospheric Science , 33(8):1521–1536, 1976.[36] E. Lee, M.E. Brachet, A. Pouquet, P.D. Mininni, and D. Rosenberg. Lack of universality in decayingmagnetohydrodynamic turbulence.
Physical Review E , 81:016318, 2010.[37] Paul T Lin, John N Shadid, Raymond S Tuminaro, Marzio Sala, Gary L Hennigan, and Roger PPawlowski. A parallel fully coupled algebraic multilevel preconditioner applied to multiphysicspde applications: Drift-diffusion, flow/transport/reaction, resistive mhd.
International Journal forNumerical Methods in Fluids , 64(10-12):1148–1179, 2010.[38] W.C.. M¨uller and D. Carati. Large-eddy simulation of magnetohydrodynamic turbulence.
ComputerPhysics Communications , 147:544–547, 2002.[39] A.A. Oberai and J. Wanderer. A dynamic approach for evaluating parameters in a numericalmethod.
International Journal for Numerical Methods in Engineering , 62:50–71, 2005.[40] A.A. Oberai, J. Lui, D. Sondak, and T.J.R. Hughes. A residual-based eddy viscosity model for theles of turbulent flows.
Comput. Methods Appl. Mech. Engrg. , 2014.[41] R. Peyret.
Spectral Methods for Incompressible Viscous Flow . New York: Springer-Verlag NewYork, Inc, 2010.[42] Y. Ponty, H. Politano, and J-F. P. Pinton. Simulation of induction at low magnetic prandtl number.
Physical Review Letters , 92:144503, 2004.[43] S. B. Pope.
Turbulent Flows . New York: Cambridge University Press, 2000.[44] A. Pouquet, E. Lee, M.E. Brachet, P.D. Mininni, and D. Rosenberg. The dynamics of unforcedturbulence at high reynolds number for taylor–green vortices generalized to mhd.
Geophysical andAstrophysical Fluid Dynamics , 104:115–134, 2010.[45] Robert S Rogallo and Parviz Moin. Numerical simulation of turbulent flows.
Annual Review ofFluid Mechanics , 16(1):99–137, 1984.[46] Pierre Sagaut.
Large eddy simulation for incompressible flows . Springer, 2002.[47] Wolfram Schmidt. Large eddy simulations in astrophysics. arXiv preprint arXiv:1404.2483 , 2014.[48] J.N. Shadid, RS Tuminaro, Karen Dragon Devine, Gary Lee Hennigan, and PT Lin. Performanceof fully coupled domain decomposition preconditioners for finite element transport/reaction simu-lations.
Journal of Computational Physics , 205(1):24–47, 2005.[49] J.N. Shadid, AG Salinger, RP Pawlowski, PT Lin, GL Hennigan, RS Tuminaro, and RB Lehoucq.Large-scale stabilized fe computational analysis of nonlinear steady-state transport/reaction sys-tems.
Computer methods in applied mechanics and engineering , 195(13):1846–1871, 2006.[50] J.N. Shadid, R.P. Pawlowski, J.W. Banks, L. Chacon, P.T. Lin, and R.S. Tuminaro. Towards ascalable fully-implicit fully-coupled resistive mhd formulation with stabilized fe methods.
Journalof Computational Physics , 229:7649–7671, 2010.[51] J.N. Shadid, R.P. Pawlowski, E.C. Cyr, R.S. Tuminaro, L. Chacon, and P.D. Weber. Scalableimplicit incompressible resistive MHD with stabilized fe and fully-coupled newton-krylov-amg. inpreparation, 2014.
52] Sergey Smolentsev, Ren´e Moreau, Leo B¨uhler, and Chiara Mistrangelo. Mhd thermofluid issues ofliquid-metal blankets: phenomena and advances.
Fusion Engineering and Design , 85(7):1196–1205,2010.[53] D. Sondak and A.A. Oberai. Large eddy simulation models for incompressible magnetohydrody-namics derived from the variational multiscale formulation.
Physics of Plasmas , 19:102308, 2012.[54] Tayfun E Tezduyar. Stabilized finite element formulations for incompressible flow computations.
Advances in applied mechanics , 28:1–44, 1991.[55] M.L. Theobald, P.A. Fox, and S. Sofia. A subgrid-scale resistivity for magnetohydrodynamics.
Physics of Plasmas , 1:3016, 1994.[56] G´abor T´oth. The · ¡ i¿ b¡/i¿= 0 constraint in shock-capturing magnetohydrodynamics codes. Journalof Computational Physics , 161(2):605–652, 2000.[57] RS Tuminaro, CH Tong, J.N. Shadid, KD Devine, and DM Day. On a multilevel preconditioningmodule for unstructured mesh krylov solvers: two-level schwarz.
Communications in numericalmethods in engineering , 18(6):383–389, 2002.[58] Z. Wang and A.A. Oberai. Spectral analysis of the dissipation of the residual-based variationalmultiscale method.
Computer Methods in Applied Mechanics and Engineering , 199:810–818, 2010.[59] A. Yoshizawa. Subgrid-scale modeling of magnetohydrodynamic turbulence.
Physical Society ofJapan , 60:9–12, 1991., 60:9–12, 1991.