23P075-P
Fully-implicit MHD Simulation Using Nonconforming Vector
Finite Elements
非適合型ベクトル有限要素法による陰的 MHD シミュレーション
Yasuhiro Kagei and Shinji Tokuda1
影井康弘，徳田伸二

1Research Organization for Info. Sci. and Tech., 2-4, Shirakata, Tokai, Ibaraki, 319-1106, Japan
Japan Atomic Energy Agency, 2-166,Oaza-Obuchi-Aza-Omotedate, Rokkasho, Aomori, 039-3212, Japan
高度情報科学技術研究機構 319-1106 茨城県那珂郡東海村白方白根 2-4
1
原子力機構・青森研究開発センター 039-3212 青森県上北郡六ヶ所村大字尾駮字表舘 2-166
A novel vector finite element method is proposed to ensure that vector variables in the MHDequations satisfy the divergence-free and the curl-free constraints exactly in a general coordinatesystem. The method called ’Nonconforming Vector Finite Element Method’ were implemented ina single-fluid resistive MHD code and validations were performed.

1

Introduction

For long time simulations of MHD instabilities, it is essential to use a scheme satisfying thedivergence-free constraint on the magnetic field(∇ · b = 0). Another constraint is also important for the analysis of MHD modes close to themarginal stability. Fluid is nearly incompressible
and necessary to satisfy the condition ∇ · v ∼ 0
below the marginal stability. These constraints
are quite important because unphysical spurious
modes emerge and interact with physical modes
as violation of such a constraint. In this paper,
we propose a novel finite element method to ensure that vector variables in the MHD equations
satisfy the divergence-free and the curl-free constraints exactly in a general coordinate system.
The method, named ’Nonconforming Vector Finite Element Method’, consists of two ideas. The
first idea is that diﬀerent types of formulations(basis functions) are used to describe covariant
and contravariant components of a vector. A basis function of a covariant vector is determinedin accordance with the discrete ’curl’ operator,
and that of a contravariant vector is in accordance with the discrete ’divergence’ operator toensure both div-curl and curl-grad operators are
identically zero. Since the basis function is different among covariant and contravariant components, transform between them is not definedby the local metric tensor. Therefore, we introduce the second idea that the equation of thecovariant-contravariant metric transformation issubstituted into the weak form in which norm
conserving condition is imposed. We call thismethod ’nonconforming’ in the same manner as
in FEM theory. The proposed formulation wereimplemented in a single-fluid resistive MHD codeand validations were performed.

2Simulation Model

2.1

Nonconforming Vector Finite Element Formulation

We focus attention in this study on the cylindrical tokamak. Vector variables in the MHD
equations are represented as a combination ofradial basis functions. We adopt diﬀerent basis
functions for each component of a covariant vector (aµ for µ = s, θ, ϕ) and a contravariant vector(Aµ ),∑
as =as,j−1/2 cj−1/2 (s),(1)(aθaϕ)

j=∑( aaϕ,jjAs =)θ,j∑ej (s),

Asj ej (s),(2)(3)j()) ∑( θ
Aj−1/2Aθ
=cj−1/2 (s),AϕAϕj−1/2j(4)2.2 Fully Implicit Method
Explicit time stepping schemes are computationally expensive for the free boundary simulation in which functions have a fairly robust profile through the plasma surface. Fully implicit
approach is more appropriate, and the backwarddiﬀerentiation (BDF) algorithm is implemented
in this work for the linear MHD stability analysis.

3Ideal MHD Eigenvalue code
Nonconforming Vector FEM code

21.5

Vs

where ej (s) is the linear finite element andcj−1/2 (s) is the piecewise constant element. Sucha formulation makes the divergence-free and thecurl-free constraints satisfiable everywhere.

10.5

0
0

0.2

0.4

0.6

0.8

1sFig. 1: The radial profile of the m/n=2/1 Suydam mode.
The resonant surface, q=2, is located at s=0.5, where theSuydam index D =0.588. 500 radial grid points are used
and the amplitude is normalized by the norm ⟨v, v⟩1/2 .
The growth rate −5.632×10−3 for the ideal MHD parameter case is shown to be consistent with −5.858 × 10−3
obtained from the ideal MHD eigenvalue code.

Simulation Results
1.5e-11
N=500
N=200

√
i g ∇ · b / h b, b i1/2

1e-11

5e-12

0-5e-12-1e-11
00.2

0.4

0.6

0.8

1

250

300

s
3e-12h ∇ · b, ∇ · b i1/2 / h b, b i1/2

N=500
N=2002e-121e-120050

100150200t / τAFig. 2: The radial profile of the divergence of the magnetic field (upper fig.) and the temporal variation of theirnorm(lower fig.) for two diﬀerent cases of grid points(N ).
0.016
D = 0.588
D = 0.464
D = 0.343

0.014
0.012

√g ∇ · v / h v, v i1/2

Numerical simulations of MHD instabilities arecarried out to check the validity of our code. The
mode structure and growth rate are confirmed tobe consistent with that obtained from the linear
MHD eigenvalue solver using the conventional
FEM formulation in Ref. [1] as shown in Fig. 1.
The divergence constraint of the magnetic field,∇ · b = 0, is examined during the mode growth.
The measurement in Fig. 2 indicates that it issatisfied with an error of only ∼10−12 due to thediscretization of the spatial derivatives. Anotherdivergence constraint is examined for the velocity
field in Fig. 3. The incompressibility condition,∇ · v = 0, is violated around the resonant surface
position, but its amplitude tends towards zero as
Suydam index parameter is varied towards themarginal stability limit.
In the free boundary simulation by using thepseudo-vacuum model [2], the vacuum is modeled as a highly resistive, low density plasma, andthe density and resistivity profiles have a steep
gradient at the plasma surface. Fully-implicit
method allows the use of as large time steps as
Alfven time scale given by the core plasma parameters even when the simulation is run withthe pseudo-vacuum model. It is verified that thesimulation with the density ratio of 1/100 andthe resistivity ratio of 106 does not fail.

0.01
0.008
0.006
0.0040.0020References[1] R.Gruber and J.Rappaz: Finite Element Methods
in Linear Ideal Magnetohydrodynamics (SpringerVerlag, Berlin, 1985), Chap. 3, p.42-78.
[2] G.Kurita, M.Azumi, et al.: Nucl. Fusion 26(1986) 449.

-0.002
00.2

0.4

0.6

0.8

1sFig. 3: The radial profile of the divergence of the flow velocity. The incompressibility condition is violated aroundthe resonant surface position, s = 0.5, but the amplitudetends towards zero as Suydam index parameter(D) is varied towards 1/4, i.e., Suydam criterion.