# A globally conservative finite element MHD code and its application to the study of compact torus formation, levitation and magnetic compression

The DELiTE (Differential Equations on Linear Triangular Elements) framework was developed for spatial discretisation of partial differential equations on an unstructured triangular grid in axisymmetric geometry. The framework is based on discrete differential operators in matrix form, which are derived using linear finite elements and mimic some of the properties of their continuous counterparts. A single-fluid two-temperature MHD code is implemented in this framework. The inherent properties of the operators are used in the code to ensure global conservation of energy, particle count, toroidal flux, and angular momentum. The code was applied to study a novel experiment in which a compact torus (CT), produced with a magnetized Marshall gun, is magnetically levitated off an insulating wall and then magnetically compressed through the action of currents in the levitation/compression coils located outside the wall. We present numerical models for CT formation, levitation, and magnetic compression, and comparisons between simulated and experimental diagnostics.

There are no comments yet.

## Authors

• 1 publication
• 1 publication
• ### Exposing and exploiting structure: optimal code generation for high-order finite element methods

Code generation based software platforms, such as Firedrake, have become...
11/07/2017 ∙ by Miklós Homolya, et al. ∙ 0

• ### Helicity-conservative finite element discretization for MHD systems

We construct finite element methods for the magnetohydrodynamics (MHD) s...
07/15/2020 ∙ by Kaibo Hu, et al. ∙ 0

• ### Geometrical discretisations for unfitted finite elements on explicit boundary representations

Unfitted (also known as embedded or immersed) finite element approximati...
09/09/2021 ∙ by Santiago Badia, et al. ∙ 0

• ### A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations

We develop general tools to construct fully-discrete, conservative numer...
06/26/2020 ∙ by Hendrik Ranocha, et al. ∙ 0

• ### On the spectrum of the finite element approximation of a three field formulation for linear elasticity

We continue the investigation on the spectrum of operators arising from ...
08/18/2021 ∙ by Linda Alzaben, et al. ∙ 0

• ### A Moving Mesh Method for Modelling Defects in Nematic Liquid Crystals

The properties of liquid crystals can be modelled using an order paramet...
10/05/2019 ∙ by Craig S. MacDonald, et al. ∙ 0

• ### A Large Scale Comparison of Tetrahedral and Hexahedral Elements for Finite Element Analysis

The Finite Element Method (FEM) is widely used to solve discrete Partial...
03/22/2019 ∙ by Teseo Schneider, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Ideally, a numerical discretisation method should, as well as accurately representing the continuous form of the mathematical equations that describe a particular physical system, reproduce the physical properties of the system being modelled. In practice, such properties are often expressed as conservation laws, and their maintenance can be just as important as standard numerical method assessment gauges like convergence, stability, accuracy, and range of applicability. Numerical solutions that contradict basic physical principles by, for example, destroying mass or energy are inherently unreliable when applied to novel physical regimes. On the other hand, many numerical methods do not have strong conservation properties, or conserve only some naturally conserved quantities, but can still regularly produce informative results that have a resemblance to real-world observations. However, as more complicated physical systems are modelled, physically incorrect numerical solutions can go unnoticed. This is especially true in simulating complex physical models, such as fluid dynamics or magnetohydrodynamics (MHD), where the maintenance of conservation laws puts more constraints on numerical schemes and, thus, helps to avoid spurious solutions. It is best to deal with the physical fidelity of the model at the numerical method design level [1], and try to replicate, in the discretised form, as many of the conservation laws inherent to the original physical system as possible.

Numerical methods with discrete conservation properties are well known in computational fluid dynamics (for example, [1, 2, 3, 4]), and in computational MHD [5, 6, 7]. In this paper we present a novel numerical scheme for a two-dimensional (axisymmetric) compressible MHD system, which is based on a continuous Galerkin linear finite element method on an unstructured triangular mesh and by construction has global (for the whole domain) conservation of mass, energy, toroidal flux and angular momentum. A novelty of the code is that all discrete spatial differential operators are represented as matrices, and the discretized MHD equations are obtained by simply replacing the original continuous differentiations with the corresponding matrix operators.

Note that by global conservation of a quantity in our numerical method, we imply that there is a discretised analogue of the continuity equation for that quantity, and, when integrated over the volume, its fluxes are completely cancelled in the interior of the domain, even though the explicit form of these fluxes are not always given. As shown in [8] and [1], global conservation for a method with local support ( local stencil) also implies local conservation. To enable the development of a numerical formulation with the aforementioned conservation properties, the discrete differential operators must obey a property equivalent, for a scalar field

and a vector field

, to

 ∫u∇⋅pdV+∫p⋅∇udV=∫up⋅dΓ (1.1)

so that the discrete forms of the differential product rule and divergence theorem are satisfied [9]. This is the essence of mimetic schemes, also known as support operator methods, to which this work is closely related [1, 10].

In section 2, we use the finite element method to derive various discrete differential matrix operators for a two-dimensional domain with axial symmetry in cylindrical coordinates, and discuss their properties. These operators satisfy the discrete form of equation 1.1, and constitute the base of the DELiTE (Differential Equations on Linear Triangular Elements) framework, which is implemented in MATLAB. In section 3, the discrete single-fluid two-temperature MHD equations are constructed within the framework. The formal process of equation development is detailed in [11]. Here, using inherent properties of the differential operators, we demonstrate global conservation of mass, toroidal flux, angular momentum, and energy. In section 4, an overview of the magnetic compression experiment is presented - more details can be found in [11, 12]. The magnetic compression experiment at General Fusion was a repetitive non-destructive test to study plasma physics applicable to magnetic target fusion compression. A compact torus (CT) is formed with a co-axial gun into a containment region with an hour-glass shaped inner flux conserver, and an insulating outer wall. External coil currents hold the CT off the outer wall (radial "levitation") and then rapidly compress it inwards. In section 5, we present the details of our numerical models for simulating CT formation, magnetic levitation, and compression within the DELiTE framework. The model includes coupling between the poloidal external vacuum field solution in the part of the domain representing the insulating wall, and the full MHD solutions in the remainder of the domain. In section 6, we present principal code input parameters and simulation results, and simulated diagnostics are compared with the experimentally measured counterparts. With inclusion of the insulating wall in the model, the effect of reduced plasma/wall interaction with an improved levitation/compression coil configuration, as observed in the experiment, is reproduced with MHD simulations. In the conclusion, section 7, we present a summary of principal code features and suggest possible further extensions.

## 2 Finite element discretisation and differential operators

In this section, we start with a brief description of the discretisation and finite element method, and introduce notations, then we derive several types of discrete differential operators used in the code and discuss their properties.

### 2.1 Discretisation and finite element method overview

The two-dimensional computational domain, with azimuthal symmetry in cylindrical coordinates, is represented by an unstructured triangular mesh. To develop the finite element discretisation, we drew inspiration from material presented in [13], which in turn, is based on material in [14], in which a finite element method is used to solve Laplace and Poisson equations in two dimensions. The freely-available mesh generator DISTMESH [15, 16, 17] was adapted to provide the computational grid. Nodes are located at triangle vertices. In the linear finite element method, any continuous function is approximated as a piecewise continuous function that is linear across each triangular element:

 u(r)≈U(r)=NeΣe=1Ue(r) (2.1)

where is the number of elements, and represents within element :

 Ue(r)=Ae+Ber+Cez (2.2)

Here, and are constants that are specific to element . These coefficients will be derived in section 2.2. Equivalent to equation 2.1, may be defined using a combination of basis functions as

 u(r)≈U(r)=NnΣn=1Unϕn(r) (2.3)

where is the number of nodes, represents the values of at node , and is the basis function associated with node .

The linear basis functions have the forms of pyramids, as depicted in figure 1, which indicates a portion of a computation mesh, with depictions of the basis functions for example nodes and . Each node is associated with a pyramid function , which has an elevation of one above the node, and falls linearly to zero at the immediately surrounding nodes. Each pyramid function has sides, where is the node-specific number of triangular elements surrounding node . Note that the pyramid shapes associated with basis functions for nodes on the domain boundary may have, in addition, one or more vertical sides. The remaining individual pyramid-side functions are defined by sections of planes that are tilted relative to the plane. Thus, each node is associated with tilted planes. In total, there are tilted planes defined in the solution domain, where . Each triangular element is associated with three tilted planes. The tilted planes are truncated to constitute pyramid sides by defining the functions representing the pyramid sides as

 ψen(r)=aen+benr+cenz (2.4)

with the additional truncating property that for all points located at which lie outside the triangular element associated with . The notation indicates that each pyramid side function is associated with a particular node , and with a particular triangular element . The coefficients , and are also associated with a particular node and element, and are such that at its associated node , and at the other two nodes in the triangular element . In summary, the basis functions , which define pyramid shapes with a peak elevation of one at node , can be expressed as

 ϕn(r)=enΣψen(r) (2.5)

where the summation is over the pyramid side functions associated with node , each of which is non-zero only over its associated triangular element. The basis functions have the property that , where is the Kronecker delta, and are the coordinates of node . Noting that the volume of a pyramid is given by , where is the pyramid base area and is the pyramid height, this leads naturally, for any continuous function (including the piecewise linear approximation), to the property

 ∫ϕn(r)U(r)drdz=(Un+O(he))∫ϕn(r)drdz≈Unsn3 (2.6)

where , is the support area of node (area of the base of the pyramid function defined by ), pyramid height and is the typical element size. This identity is analogous to the integral property of the Dirac delta function. In deriving the property, it is assumed that the function is sufficiently smooth that it is approximately constant (to order ) in the support area of node . We neglect the term of order because our numerical scheme has overall the first order accuracy, as defined by the use of linear basis functions.

Another important property of the basis functions is related to partition of unity - the sum of all the basis functions in the domain, at any point in the domain ( at non-nodal locations, as well as at nodal locations), is equal to one. This property also hold for the pyramid side functions, .

In the following, the notation (or will be used to denote vectors of dimensions that contain node-associated quantities, while denotes vectors of element-associated quantities. The notation will be used to denote matrices of dimensions that operate on vectors of nodal quantities , to produce vectors of elemental quantities , while operates on vectors of elemental quantities to produce vectors of nodal quantities . The notations and will be used to denote square matrices with dimensions and respectively. Notations defining the various matrix dimensions are collected in table 1.

We introduce node-to-element averaging as

 e=13Me––––––––∗X–– (2.7)

Here, the connectivity matrix has dimensions . Each row of corresponds to a particular triangular element, and has only three non-zero entries. for the column indexes corresponding to the indexes of the nodes located at the vertices of the triangle with index The symbol represents regular matrix multiplication. The radial coordinates of element centroids are defined with , where contains the coordinates of the nodes. The superscript implies the transpose operation. Similarly, the axial centroid coordinates are defined as . The vector of nodal support areas , containing the areas of the bases of the pyramid functions associated with the nodes, is related to , the vector of elemental areas, as .

Volume integrals over the computational domain can be approximated in two ways, corresponding to nodal or elemental representations of the integrand function

 ∫u(r)dV ≈dV–––T∗U–– or ∫u(r)dV≈dVe–––––T∗Ue–––

where contains the volumes associated with each node, which are found by integrating the node-associated areas from to in the toroidal direction, and contains the elemental volumes. The factor of three in the former expression arises because each elemental area is shared by three nodes. Note that these two approximations do not give identical results for vectors related by equation 2.7, unless the original integrand function is constant. The symbol represents the Hadamard product, piecewise element-by element multiplication ( ), and the symbol represents Hadamard division, piecewise element-by element division ( ).

Defining , , and as the diagonal arrays constructed from , , , and , we define a volume-averaging operator

 Wn––––––––=R––––−1∗S––––−1∗Me––––––––T∗ˆS––––∗ˆR––––

that is used to map element-based quantities to node-based quantities, as

 =Wn––––––––∗Ue––– (2.8)

This operator satisfies the following identity:

 dV–––T∗(Q––∘(Wn––––––––∗Ue–––))=dVe–––––T∗(Qe–––∘Ue–––) (2.9)

where and are the discrete representations, defined at the element centroids, of the approximations to the continuous functions and , and is the discrete representation, defined at the nodes, of , and is related to by equation 2.7. A proof of this identity follows:

 dV–––T∗(Q––∘(Wn––––––––∗Ue–––)) =2π3(s–∘r–)T∗(Q––∘(R––––−1∗S––––−1∗Me––––––––T∗ˆS––––∗ˆR––––∗Ue–––)) =2π3Q––T∗(Me––––––––T∗ˆS––––∗ˆR––––∗Ue–––) =2π(ˆS––––∗ˆR––––∗Ue–––)T∗(13Me––––––––∗Q––) (transpose the scalar) =dVe–––––T∗(Qe–––∘Ue–––) (use equation 2.7)

Note that the matrix transpose relation

 (A∗B)T=BT∗AT (2.10)

for matrices and , is used to transpose the scalar on the right side of the equation in the second last step of the derivation above. In the particular case with , the identity becomes

 (2.11)

### 2.2 First order node-to-element differential operators

The node-to-element derivative operator matrices are defined such that

 U′er–––– =Dre––––––––∗U–– U′ez–––– =Dze––––––––––∗U––

Referring for example to the radial derivative operator, application of , a matrix with dimensions , to the vector , returns containing the values for the radial derivatives of inside the triangular elements. The node-to-element gradient and divergence operations, for cylindrical coordinates with azimuthal symmetry are

 ∇e––––––U–– =(Dre––––––––∗U––)^r+(Dze––––––––––∗U––)^z (2.12)

and

 ∇e––––––⋅P–– =((Dre––––––––∗(r–∘Pr–––))+(Dze––––––––––∗(r–∘Pz–––)))⊘re––

where and are the nodal representations of the and components the continuous vector field . A schematic of the node-to-element gradient operation mechanism is indicated in figure 2.

To derive the and operators we use the elemental representation of , equations 2.1 and 2.2. Equation 2.2 defines the values of the approximation at the nodes as

 Uε=Ae+Berε+Cezε (2.13)

where denotes the indexes of the nodes at the vertices of triangular element . In the following, the element-specific node indexes will be replaced with the indexes and for simplicity. From equation 2.2, the radial and axial first spatial derivatives of are constants across the element , and are given by , and . Expressions for and are found using equation 2.13:

 U1−U3 =Be(r1−r3)+Ce(z1−z3) U2−U3 =Be(r2−r3)+Ce(z2−z3)
 ⇒Be =∂Ue∂r=U1(z2−z3)+U2(z3−z1)+U3(z1−z2)2se Ce =∂Ue∂z=U1(r2−r3)+U2(r3−r1)+U3(r1−r2)−2se (2.14)

The triangle area is introduced here, assuming that vertices are number counterclockwise, noting that , where is the determinant of array , so that

 2se=(r1−r3)(z2−z3)−(r2−r3)(z1−z3)

The element-specific spatial derivatives can also be expressed in terms of the coefficients in the element-specific pyramid-side functions , and . As noted earlier, each pyramid-side function is specific to a particular node, and to a particular element, and is defined as , and has the property . This yields, for each element , the identity , or

 ⎛⎜⎝1r1z11r2z21r3z3⎞⎟⎠∗⎛⎜ ⎜⎝ae1ae2ae3be1be2be3ce1ce2ce3⎞⎟ ⎟⎠=⎛⎜⎝100010001⎞⎟⎠
 ⇒Ce (2.15)

Comparing equations 2.14 with equation 2.15, it is evident that

 ∂Ue∂r =U1be1+U2be2+U3be3=([]cccbe1be2be3)∗⎛⎜⎝U1U2U3⎞⎟⎠ ∂Ue∂z =U1ce1+U2ce2+U3ce3=([]cccce1ce2ce3)∗⎛⎜⎝U1U2U3⎞⎟⎠ (2.16)

Given the values of the approximation for the piecewise linear function on the triangle vertices, and the element-specific array evaluated using equation 2.15, equations 2.16 can be used to determine the values of the spatial derivatives of at the interior of each element. and are initially defined as sparse all-zero arrays. Reverting to node-specific notation, where are the indexes of the nodes at the vertices of element , the values and are inserted in row of with placements at the column indexes . Similarly, the values and , for each element, are used to assemble . The resulting derivative operators produce exact derivatives for nodal functions with linear dependence on and , and have first order accuracy () when applied to nonlinear functions. The operators introduced in the following subsections are all based on these node-to-element derivative operators, and so they all have the same accuracy.

Here, we will demonstrate a property of the node-to-element operators that will be used to demonstrate angular momentum conservation later in section 3.3.3. For any continuous scalar functions and , the identity (Stoke’s theorem) holds that

 ∫(∇f×∇g)⋅dS=∫∇×(f∇g)⋅dS=∫f∇g⋅dl (2.17)

With azimuthal symmetry, , where is an elemental area in the poloidal () plane. Hence, the discrete form of , using the node-to-element differential operators, is

 (∫(∇f×∇g)⋅dS)disc. =se––T∗((Dze––––––––––∗f––)∘(Dre––––––––∗g–)−(Dre––––––––∗f––)∘(Dze––––––––––∗g–)) =(Dze––––––––––∗f––)T∗ˆS––––∗(Dre––––––––∗g–)−(Dre––––––––∗f––)T∗ˆS––––∗(Dze––––––––––∗g–) =f––T∗Dze––––––––––T∗ˆS––––∗Dre––––––––∗g–−f––T∗Dre––––T––––––∗ˆS––––∗Dze––––––––––∗g– =f––T∗(Dze––––––––––T∗ˆS––––∗Dre––––––––−Dre––––––––T∗ˆS––––∗Dze––––––––––)∗g– =f––T∗(B––––T−–B––)∗g– =f––T∗A––––∗g– (2.18)

where . Setting and , where and are the vectors defining the nodal values of basis functions and respectively ( and ), equations 2.18 and 2.17 imply that

 Amn≡A(m,n)=ϕm–––T∗A––––∗ϕn–––=∫ϕm∇ϕn⋅dl (2.19)

We will show that each element of is zero at internal nodes, . There are four cases to consider: (1) and nodes and are not adjacent, (2) , (3) and are the indexes of adjacent internal nodes, and (4) and are the indexes of adjacent nodes located on the boundary of the computational domain. In case (1), is obviously zero. In case (2), . The first term scalar term here can be transposed, so that in case (2).

Figure 3 is a top-view of figure 1. Referring to equation 2.19 and to figure 3, the only finite contribution to in case (3) is from the contour integral along the boundary (depicted with dark blue and red dashed lines) of the diamond-shape representing the overlaping area of the basis functions and . It can be seen that along the part of the contour, and along the remaining parts, so that in case (3). In case (4), when and are the indexes of adjacent nodes located on the boundary of the computational domain, the integral would be finite along the light blue dashed line connecting nodes and (in case (4), the light blue dashed line would represent part of the computational domain boundary). Therefore, is finite only when and are the indexes of adjacent nodes located on the boundary of the computational domain, leading to the identity

 (Dze––––––––––T∗ˆS––––∗Dre––––––––−Dre––––––––T∗ˆS––––∗Dze––––––––––)int=0 (2.20)

### 2.3 First order element-to-node differential operators

The element-to-node derivative operator matrices are defined such that

 U′r––– =Drn––––––––––∗Ue––– U′z––– =Dzn––––––––––∗Ue–––

Referring for example to the radial derivative operator, application of the operator, which is a matrix with dimensions , to the vector , which contains the values of at the element centroids, returns the values for the radial derivatives of at the nodes. The element-to-node gradient and divergence operations, for cylindrical coordinates with azimuthal symmetry are

 ∇n––––––Ue––– =(Drn––––––––––∗Ue–––)^r+(Dzn––––––––––∗Ue–––)^z (2.21)

and

 ∇n––––––⋅Pe––– =((Drn––––––––––∗(re––∘Per–––))+(Dzn––––––––––∗(re––∘Pez–––)))⊘r– (2.22)

A schematic of the element-to-node gradient operation mechanism is indicated in figure 4. The element-to-node differential operators are derived in relation to the node-to-element operators by requiring the discrete form of equation 1.1 to hold:

 dV–––T∗[P––⋅(∇n––––––Ue–––)]+dVe–––––T∗[Ue–––∘(∇e––––––⋅P––)]=0 (2.23)

or

 dV–––T∗[U––∘(∇n––––––⋅Pe–––)]+dVe–––––T∗[Pe–––⋅(∇e––––––U––)]=0 (2.24)

In this derivation, it is assumed that the boundary term , there is no flux of the continuous vector field through the boundary. For arbitrary discrete nodal representations (equation 2.23) and (equation 2.24), a consequence of this assumption is that the element-to-node gradient operator produces valid results at the boundary nodes only if the original continuous function is zero at the boundary . Similarly, the element-to-node divergence operator produces valid results at the boundary nodes only if the original continuous function has no component perpendicular to the boundary . In the following, we will refer to these conditions as the natural boundary conditions. For the terms involving radial derivatives, equation 2.24 implies that

 dV–––T∗(U––∘((–Drn–––––∗(re––∘Per–––))⊘r–)) =−dVe–––––T∗(Per–––∘(Dre––––––––∗U––)) ⇒2π3s–T∗(U––∘(Drn––––––––––∗(re––∘Per–––))) =−2πse––T∗(re––∘Per–––∘(Dre––––––––∗U––)) ⇒13U––T∗(S––––∗(Drn––––––––––∗(re––∘Per–––))) =−(re––∘Per–––)T∗(ˆS––––∗Dre––––––––∗U––) ⇒13U––T∗(S––