Ideal and Resistive Magnetohydrodynamic Two-Dimensional Simulation of the Kelvin-Helmholtz Instability in the Context of Adaptive Multiresolution Analysis

i i “main” — 2017/8/14 — 15:37 — page 317 — #1 i Tema i Tendências em Matemática Aplicada e Computacional, 18, N. 2 (2017), 317-333 © 2017 Sociedade Brasileira de Matemática Aplicada e Computacional www.scielo.br/tema doi: 10.5540/tema.2017.018.02.0317 Ideal and Resistive Magnetohydrodynamic Two-Dimensional Simulation of the Kelvin-Helmholtz Instability in the Context of Adaptive Multiresolution Analysis A.K.F. GOMES1* , M.O. DOMINGUES2 and O. MENDES3 Received on November 28, 2016 / Accepted on May 10, 2017 ABSTRACT. This work is concerned with the numerical simulation of the Kelvin–Helmholtz instability using an ideal and resistive two-dimensional magnetohydrodynamics model in the context of an adaptive multiresolution approach. The Kelvin-Helmholtz instabilities are caused by a velocity shear and normally expected in a layer between two fluids with different velocities. Due to its complexity, this kind of problem is a well-known test for numerical schemes and it is important for the verification of the developed code. The aim of this paper is to verify the implemented numerical model with the well-known astrophysics FLASH code. Keywords: magnetohydrodynamics, Kelvin-Helmholtz instability, adaptive multiresolution analysis, numerical simulation, scientific computing. 1 INTRODUCTION The magnetohydrodynamics (MHD) theory describes the dynamics of a conducting fluid in presence of magnetic fields and constitutes an important tool to study the macroscopic behavior of plasmas [11]. In this context, the Kelvin-Helmholtz instability, which is commonly expected in boundary layers separating two fluids, is an important and a complex physical problem that can be studied with the MHD models, and should often occur in both astrophysical and space geophysical environments [6]. On the discretization of the MHD system, we use a finite volume (FV) method combined with an adaptive multiresolution (MR) approach to create a computational mesh refined where local structures are presented. The FV method is based on the integral form of conservation laws and guarantees the conservation of the model quantities. The MR for cell-averages was firstly introduced by Ami Harten [10]. The idea is to represent a set of data in different levels of resolution by using a wavelet transform. In the C ++ code named *Corresponding author: Anna Karina Fontes Gomes 1 Pós-graduação em Computação Aplicada, INPE, São José dos Campos, SP, Brasil. E-mail: 2 Laboratório Associado de Computação e Matemática Aplicada, INPE, São José dos Campos, SP, Brasil. E-mail: 3 Divisão de Geofı́sica Espacial, INPE, São José dos Campos, SP, Brasil. E-mail: i i i i i i “main” — 2017/8/14 — 15:37 — page 318 — #2 i 318 i MHD SIMULATION CARMEN [15] the MR algorithm is implemented for compressible Navier-Stokes and five more system of equations. The ideal MHD equations were added later to the CARMEN code [3, 8, 9], and it is employed herein. We use the FLASH code [7], developed by the Flash Center in University of Chicago, wellknown in astrophysics and space geophysics, to create a reference MHD solution to our results, since it is not possible to obtain a exact solution. The goal of this work is to verify the numerical results of CARMEN code for the Kelvin-Helmholtz instability problem by comparing them with the reference solution, which is obtained in a regular Cartesian mesh. The content is organized as follows. In Section 2, we briefly present both the MHD model and the MR approach we use to simulate the Kelvin-Helmholtz instabilities; in Section 3, the numerical methodology and implementation; in Section 4, the results and discussion. The final remarks are presented in Section 5. 2 THE MHD MODEL The ideal model describes the behavior of a perfectly conducting fluid under the influence of a magnetic field. By adding a resistive term to the MHD system, there is no magnetic flux conservation anymore, which can lead to a more diffusive behavior. The resistivity is associated to the parameter η, which comes from the Ohm’s law. For η = 0 it can trigger a different behavior to be studied in a plasma. It is important to note that when η → 0 the resistive MHD model becomes the ideal MHD equations, which describe the conservation of mass, energy, momentum and magnetic flux. In this work we consider η as a constant, but it is also possible to choose a scalar function. In this context, we introduce the resistive MHD equations. ∂ρ + ∇ · (ρu) = 0, ∂t    |B|2 ∂e +∇· e+ p+ u − B (u · B) = ∇ · [B × η(∇ × B)] , ∂t 2     2 |B| ∂ (ρu) + ∇ · ρut u + I p + − Bt B = 0, ∂t 2 (2.1b)   ∂B + ∇ · ut B − Bt u = −∇ × (η∇ × B), ∂t (2.1d) (2.1a) (2.1c) where ρ represents density, p the pressure, u = (u x , u y , u z ) the velocity vector, B = (Bx , B y , Bz ) the magnetic field vector, I the identity tensor of order 2, and γ the ratio of specific heats (γ > 1). The pressure is given by the constitutive law   |B|2 |u|2 − , p = (γ − 1) e − ρ 2 2 where e is the energy density. For the magnetic field, we have the Gauss’ law for magnetism ∇ · B = 0, which means there is no magnetic monopole in the solution of the MHD model. This Tend. Mat. Apl. Comput., 18, N. 2 (2017) i i i i i i “main” — 2017/8/14 — 15:37 — page 319 — #3 i GOMES, DOMINGUES and MENDES i 319 is an initial condition for the model. From the Faraday’s law ∇ × E = − ∂B ∂t (by applying the ∂(∇·B) divergence operator on the equation), we obtain ∂t = 0, i.e., this means there is no variation of the divergence of B over time. In numerical simulation, the divergence of B does not always vanish. Then it becomes necessary to implement a correction (or divergence cleaning) scheme so the solution will not lead to unphysical behavior or unwanted instabilities. In the next section, we present the numerical methodology for the simulation, including the divergence correction scheme. 3 NUMERICAL APPROACH To introduce the MHD simulation numerical methodology we firstly present the initial value problem for conservation laws of the form ∂U + ∇ · F(U) = S(U), ∂t U(x, y, t = 0) = U0(x, y), (3.1) (x, y) ∈ , (3.2) where U = (ρ , e, u x , u y , u z , Bx , B y , Bz ) is the vector of conservative variables, F = F(U) the flux tensor, S = S(U) the source term vector,  the domain and t the time. Using the definition of Equation (2.1), we have ⎛ ⎞ ρu   ⎞ ⎛ ⎜ ⎟ |B|2 0 ⎜ ⎟ u − B (u · B)⎟ ⎜ e+ p+ ⎜∇ · [B × η(∇ × B)]⎟ ⎜ ⎟ 2 ⎟ ⎜ ⎟ , S(U) = ⎜   F(U) = ⎜ ⎟. ⎜ ⎟ 2 ⎠ ⎝ 0 |B| ⎜ t t ⎟ ρu − B u + I p + B ⎜ ⎟ −∇ × (η∇ × B) ⎝ ⎠ 2 (3.3) ut B − Bt u In this section we describe the numerical methodology we use for the MHD simulation in the context of adaptive multiresolution approach. Divergence Cleaning approach. As discussed in the previous section, ∇ ·B = 0 is not satisfied numerically and it can lead to unphysical behavior in the numerical solution of the MHD model. We use the parabolic-hyperbolic correction [1], in which the errors are propagated and (...truncated)