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 , 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  and , 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
so that the discrete forms of the differential product rule and divergence theorem are satisfied . 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 . 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 , which in turn, is based on material in , 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:
where is the number of elements, and represents within element :
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
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
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
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
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
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
that is used to map element-based quantities to node-based quantities, as
This operator satisfies the following identity:
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:
|(transpose the scalar)|
|(use equation 2.7)|
Note that the matrix transpose relation
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.2 First order node-to-element differential operators
The node-to-element derivative operator matrices are defined such that
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
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.
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:
The triangle area is introduced here, assuming that vertices are number counterclockwise, noting that , where is the determinant of array , so that
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
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.
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
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
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
2.3 First order element-to-node differential operators
The element-to-node derivative operator matrices are defined such that
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
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:
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