Nonlinear MHD Halo Current Simulations

H.R. Strauss
New York University, New York, New York



A. Pletzer, W. Park, S. Jardin, J. Breslau
Princeton University Plasma Physics Laboratory, Princeton, New Jersey
Roberto Paccagnella
Istituto Gas Ionizzati del C.N.R., Padua, Italy

Abstract. Preliminary work is presented on halo current simulations in ITER. The first step is the study of VDE (vertical displacement event) instabilities. The growth rate is consistent with scaling inversely proportional to the resistive wall penetration time. The simulations have self consistent resistivity proportional to the ${-3/2}$ power of the temperature. Simulations have been done with temperature contrast between the plasma core and wall of $100$, to model the vacuum region between the core and resistive shell. Some 3D simulations are shown of disruptions competing with VDEs. The toroidal peaking factor can be as high as 3.

I. Introduction

Preliminary work is presented on halo current simulations in ITER. The first step is the study of VDE (vertical displacement event) instabilities. The growth rate is consistent with scaling inversely proportional to the resistive wall penetration time. The simulations have self consistent resistivity proportional to the ${-3/2}$ power of the temperature. Simulations have been done with resistivity contrast between the plasma core and wall of $1000$ times, to model the vacuum region between the core and resistive shell. Some 3D simulations are shown of disruptions competing with VDEs. Toroidal peaking factors are up to about 3.

The M3D (Multi-level 3D) project [1,2] carries out simulation studies of plasmas using multiple levels of physics, geometry, and grid models. In the present study is done with a resistive MHD model. M3D combines a two dimensional unstructured mesh with finite element discretization in poloidal planes [3], with a pseudo spectral representation in the toroidal direction. The unstructured mesh used in the calculations is shown in Fig.1. The part of the mesh adjacent to the outer wall (the ITER - FEAT first wall) was made using the ellipt2d package [4], which incorporates the Triangle mesh generation code.

The code includes a temperature equation, with thermal conduction along the magnetic field modelled by the artificial sound method [7]. The resistivity is proportional to $T^{-3/2}$, where $T$ is the temperature. The halo region between the plasma core and the wall is modelled as a cold resistive plasma. Simulations have been done with core temperature 100 times the halo temperature, for a resistivity contrast of $1000.$

The M3D code includes resistive wall boundary conditions, which match the solution inside the resistive wall to the exterior vacuum solution. The exterior problem is solved with a Green's function method, using A. Pletzer's GRIN code [5].

II. Resistive Wall Boundary Conditions

The plasma is bounded by an inner, thin, resistive wall. Surrounding this is an outer vacuum region, which can contain external current sources.

Represent the vacuum field as

$${\bf B}_v= \nabla \psi_v \times \nabla \phi + \nabla \lambda + I_0 \nabla \phi$$ (1)

Here $I_0$ is a constant which is equal to the constant part of $I$ in the plasma. Setting $\nabla\times {\bf B}_v= 0$ gives

$$\Delta^* \psi_v = R^2 \nabla {1\over R^2}\cdot \nabla \psi_v = 0$$ (2)

and

$$\nabla { \partial \psi_v \over \partial \phi } = 0$$

which can be satisfied if $\psi_v$ is a function only of $R, Z$ but not $\phi,$ i.e. only the $n = 0$ toroidal harmonic is included.

$$\psi_v = \overline{\psi_v}.$$

To satisfy $\nabla\cdot {\bf B}_v= 0,$

$$\nabla ^2 \lambda = 0$$ (3)

On the resistive wall boundary, integrating $\nabla\cdot {\bf B}$ across the thin shell gives

$$\hat{n}\cdot [[{\bf B}]] = 0$$

where the double bracket indicates a jump across the thin wall, and $\hat{n}$ is the outward normal from the plasma. In order that the normal component of the magnetic field be continuous, for the axially symmetric part of the solution,

$$\psi_v = \overline{\psi}_p.$$

For the toroidally varying part, we require

$${ \partial \lambda \over \partial n } = {\bf B}^p \cdot \hat{n}$$ (4)

where ${\bf B}^p$ is the magnetic field in the plasma. These boundary conditions determine the vacuum field.

The vacuum field is solved by the GRIN code. For an axially symmetric wall, the vacuum field is first Fourier expanded:

$$\lambda = \sum_{n=1}^{N} \lambda_n \exp{i n \phi}$$

From Green's identity one has an integral equation relating $\lambda_n$ to given ${\partial \lambda_n}/{\partial n}$ on the boundary contour [6]

$$\oint d\ell R \left\{ G_n \hspace{1mm} \frac{\partial \lam... ... n} - \frac{\partial G_n}{\partial n} \lambda_n \right\} = 0.$$

where $G_n(R, Z; R' Z')$ is the $n$th harmonic of the Green's function.

When discretized, this becomes a matrix equation relating the values of $\lambda_n$ and ${\partial \lambda_n}/{\partial n}$ at the dicrete mesh points of the boundary. Given a set of boundary points, $R_i, Z_i$ and the values $\psi_{vi}$ on the boundary, the $\phi$ indendendent, $n = 0$ part of the solution is

$$({ \partial \psi_v \over \partial n })_i = \sum_j K^0_{ij} \psi_{vj} + S_i.$$ (5)

The source term $S$ in (5) can be obtained from the applied external currents, or else using the ``virtual casing" method. In this method we first perform an ideal equilibrium calculation, with $\psi = 0$ on the boundary. The source term is found from

$$S = { \partial \psi \over \partial n }$$

where the right side is obtained from the ideal equilibrium.

The $n > 0$ part of the solution is

$$\lambda^n_i = \sum_j K^n_{ij} ({ \partial \lambda^n \over \partial n })_j$$

Now the magnetic field components in the plasma have to be matched using resistive evolution at the inner boundary, which is a thin resistive shell of thickness $\delta$ and resistivity $\eta_w.$ The boundary conditions are

$$\hat{n}\times [[{\bf E}]] = 0$$

where ${\bf E} = \eta {\bf J}.$ In the wall,

$${\bf J}= {1\over \delta}\hat{n}\times [[{\bf B}]].$$

Ohm's Law at the resistive wall is

$${\partial {\bf A} \over \partial t} = \nabla \Phi + {\eta_w \over \delta} \hat{n}\times [[{\bf B}]].$$

An equiilibrium, with ``virtual casing" source terms in the boundary conditions, is shown in Fig.2. The poloidal flux $\psi$ is shown in Fig.2(a) and the toroidal flux $I$ in Fig.2(b). The temperature is shown in Fig.3(a), and the toroidal current in Fig.3(b). A velocity perturbation is added to the equilibrium, with resistive boundary conditions, to obtain a VDE.

IV. VDE Simulations

The VDE instability growth rate is inversely proportional to the wall resisitive penetration time, or $\eta_w.$ This scaling is consistent with simulations, as will be shown below. To get the scaling it seems necessary to be in a regime in which

$$\tau_c > \tau_w > \tau_h.$$

Here $\tau_c = S \tau_A,$ and $\tau_h = (T_h/T_c)^{3/2} \tau_c,$ where $\tau_A$ is the Alfvéntime, $S = 10^4$ in the simulations, $T_h = 10^{-2} T_c,$ where $T_h$ and $T_c$ are halo and core temperatures, and $\tau_w = \delta_w / eta_w S \tau_A.$ We have chosen paramters in the regime

$$1 > 10^{-1} > {\tau_w \over \tau_c } > 10^{-2} > 10^{-3} = {\tau_h \over \tau_c }$$

It is not easy to get the VDE growth rate, even if the inequality is well satisfied. In Fig.4 are shown the growth rate as a function of time, for the case $\eta _w / \delta _w =$ (a) $0.01,$ and (b) $0.005.$ In Fig.4(a), the growth rate is large for small time as the equilibrium adjusts, and increases for large time as the VDE nonlinearly accelerates. For $35 \tau_A < t < 45 \tau_A$ there is a linear regime with $\gamma = 0.05.$ In Fig.4(b), the linear regime lasts longer, but the growth rates oscillates, about a mean of about $0.0025.$

For an order of magnitude variation in $\eta _w / \delta _w$, the growth rate of the VDE scales as

$$\gamma = 5.0 \eta_w / \delta_w.$$

The growth rate $\gamma$ as a function of $\eta _w / \delta _w$ is shown in Fig.5.

The nonlinear stage of the VDE is shown in Fig.6. The poloidal flux $\psi$ is shown in Fig.6(a) and the toroidal flux $I$ in Fig.6(b). The temperature is shown in Fig.7(a), and the toroidal current in Fig.7(b).

V. Disruption Simulations

In three dimensional simulations, disruptions can occur. In one scenario, a disruption causes a thermal quench, which in turn causes a current quench. This is accompanied by a VDE.

The following calculation is an example of this scenario. The time history of the toroidal peaking factor is shown in Fig.8(a), and the normalized peak current and peak temperature in Fig.8(b). The toroidal peaking factor almost reaches 3, but most of the time oscillates around 2. Fig.8(b) shows the current, plotted with a dashed line, which declines in value more slowly than the temperature, shown as a solid line. The temperature quench proceeds the current quench. The timescale is artificially rapid. In this simulation $S = 30,000.$ Fig.9 shows the poloidal flux $\psi$ at times $t = 71, 87, 126 \tau_A$. The early nonlinear stage of an internal kink occurs at time $t = 71 \tau_A$. The disruption occurs at the intermediate time $t = 87 \tau_A$. The VDE occurs at the later time, $t = 126 \tau_A$. The temperature is shown at the same times in Fig.10. In the intermediate time, the peak temperature is about $0.65$ of the initial peak temperature, and at the later time, it is about $0.1$ times the initial peak temperature. The toroidal flux is shown in Fig.11. There is significant intersection of the toroidal flux contours with the wall, indicating the flow of halo current. The halo current in this case is small. Fig.12 shows the electrostatic potential, which which is approximately the streamlines of the incompressible part of the flow. In Fig.12(a) the potential is typical of an internal kink, and is similar during the disruption phase of Fig.12(b). In the VDE phase of Fig.12(c) the flow has changed into a vertical flow into the divertor.