Effects of asymmetries in computations of forced vertical displacement events

Visco-resistive magnetohydrodynamic (MHD) computations with the NIMROD code (Sovinec C R et al 2004 J. Comput. Phys. 195 355) are applied to a model tokamak configuration that is subjected to induced vertical displacement. The modeling includes anisotropic thermal conduction within an evolving magnetic topology, and parameters separate the Alfvénic, resistive-wall, and plasma-resistive timescales. Contact with the wall leads to increasingly pervasive kink and tearing dynamics. The computed 3D evolution reproduces distinct thermal-quench and current-quench timescales, a positive bump in plasma current, and net horizontal forcing on the resistive wall. The MHD dynamo effect electric field, E f = − V ˜ × B ˜ , is analyzed for understanding the nonlinear effects of the fluctuations on the spreading of parallel current density and the resulting bump in plasma current. Forces on the resistive wall are consistent with Pustovitov’s analysis (Pustovitov V D 2015 Nucl. Fusion 55 113032); the plasma remains in approximate force-balance with the wall, so net force is accurately computed from integrating stress over the wall’s outer surface. Improvements to the modeling that are needed for predictive simulation of asymmetric vertical displacement events are discussed.


Introduction
Discharge-terminating disruptive events in tokamaks take many forms and involve a variety of plasma dynamics [1,2]. Whether it is a root cause or a consequence of other disruptive activity, vertical displacement poses a particularly significant risk, because it brings hot plasma in contact with surfaces that are not designed for extreme thermal loading. In addition, asymmetries that develop during vertical displacement events (VDEs) lead to horizontal forces on electrically conducting structures [3][4][5][6], and the forcing may rotate at rates that are comparable to mechanical harmonics [7]. The risks for ITER have motivated numerous experimental, analytical, and computational studies of disruptive dynamics and of means to mitigate their consequences. Here, we report on nonlinear visco-resistive magnetohydrodynamic (MHD) computations of vertical displacement in a model configuration and on the consequences of significant asymmetry that develops through contact with the wall. In these computations, vertical displacement occurs over the timescale of the resistive wall and is forced by externally imposed conditions on the magnetic field. Our analysis examines the force on the resistive wall during the current quench (CQ) and the spreading of parallel current density, which transiently raises the plasma current starting at the thermal quench (TQ).
To date, simulations of disruptions have been based on MHD modeling or on forms of reduced MHD [8][9][10][11][12][13][14][15][16][17]. The reasoning behind this choice is that macroscopic dynamics always arise during disruptions, and MHD provides a practical model for dynamics involving gross motions and changes in magnetic topology. Because disruptions in tokamak experiments involve electron and ion kinetics, radiation, neutral dynamics, and plasma-surface interaction, in addition to macroscopic dynamics, predictive simulation of disruptions in tokamaks will require far more comprehensive modeling [18] that is beyond present-day capabilities. Nonetheless, MHD modeling provides a basis for understanding disruptive dynamics and a conceptual framework from which to build comprehensive modeling. Previous axisymmetric studies evolve equilibrium relations quasi-statically to examine the generation of symmetric, open-field halo currents [11] and the influence of geometric details of surrounding structures [19]. Previous simulations of asymmetric VDEs by Strauss and coauthors examined the peaking of plasma current over the toroidal angle [15], the magnitude of horizontal forces as the ratio of external-kink time and resistive-wall time is varied [6], and the influence of changing boundary conditions on flow [20]. Recent simulations of a VDE in the NSTX experiment describe the destabilization of edge kink modes resulting from contact with solid surfaces [21].
Experimental studies of VDEs trigger them by turning parts of the control system off [22] or by programming the system to induce displacement [23]. The computations discussed here are similar to the latter approach. They start from a nominally up-down symmetric, double-null equilibrium, and one of the two divertor coils is effectively turned off at the start of the computation. Vertical displacement then proceeds on the resistive wall time , w t as eddy currents decay, which is the case whenever the wall slows the displacement. Running both axisymmetric and 3D nonlinear computations in this controlled scenario allows direct comparison, which helps us identify effects that result from asymmetry.
The remaining sections of this paper are organized as follows. Section 2 briefly describes the model and computational methods, including parameters and initial conditions. Section 3 discusses the computed results on displacement, TQ, CQ, and asymmetries. It also describes analysis of the spreading of parallel current density and forces on the wall. Section 4 provides a discussion of our findings and future efforts, and our conclusions are given in section 5.

Visco-resistive MHD modeling
The central objective of this study is to reproduce the major aspects of asymmetric disruption from vertical displacement. Modeling a TQ resulting self-consistently from asymmetric macroscopic instability requires thermal transport that is sensitive to changes in magnetic topology. The CQ occurs over a longer time-scale in experiments [1,[24][25][26], so numerical computations must also have multi-scale capabilities. In addition, net forcing on the wall only results when the wall is not an ideal conductor [27]. All of these properties are within the scope of time-dependent nonlinear non-ideal MHD models if the system of equations is solved with implicit numerical methods and is augmented with a resistive-wall and external magnetic response.
We separate the problem domain, which is axisymmetric, into the two subdomains shown in figure 1, where the resistive wall lies along their intersection. The visco-resistive MHD equations of the inner subdomain describe the evolution of particle density n n n , where total (electron plus ion) pressure is P nT is the adiabatic index. Relative to SI units, the Boltzmann constant is absorbed in T. The computations normalize the fields by using spatial dimensions and B | | of order unity, and by setting each of ion mass m, maximum-n, and 0 m to unity. The first term on the right side of equation (4) represents induction from the resistive-MHD electric field. We use the Spitzer dependence of The initial state of our idealized configuration has T varying by four orders of magnitude from the open-field region to the central value T 0 and resistivity varying by six orders of magnitude. Large resistivity at low temperature suppresses current density outside regions representing hot 1 More precisely, conditions at the magnetic axis of the equilibrium discussed below have R 0 =1.65 and B 0 =1.28, so the normalized τ A is 1.29. For convenience, discussion in the text considers τ A to be 1, however. plasma during the course of 3D dynamics. The second term on the right side of equation (4) is a numerical term that, together with a high-order representation, is used to keep magnetic divergence error small in the computations [28].
The closure relations represent anisotropic thermal conduction and anisotropic viscous stress where I is the identity tensor, and is the traceless rate of strain tensor. While the form of these relations is, technically, only suitable for collisional plasma with vanishingly small ion gyro-radii [29], it facilitates our present computations, together with the relatively modest anisotropy from the coefficient values 7.5 10 ,  n =´-|| Thermal conduction and viscous stress with these fixed coefficients are not accurate models of collisional transport in the edge plasma or of kinetic effects in the hot core. However, the unit vector b is along the 3D magnetic field, which makes the computed transport sensitive to the evolving magnetic topology, even with the modest anisotropy. The terms on the right side of equation (1) provide numerical smoothing of the particle density field. =´used in the computations are intended to have negligible impact on physically meaningful results. As a check on this, we find that repeating the linear computation of Section 3.1 with D n and D h reduced to 5 10 7 and 0, respectively, affects the growth rate by less than 0.5%. The outer subdomain is modeled by a simplified version of equation (4), where the electric field is E J h = with 0 h m having the fixed value of 100 so that resistive diffusion there is faster than other processes in the problem. This annular subdomain lies between the resistive wall and an outer conducting shell. The two subdomains are coupled by the thinwall approximation of where n is the unit normal along the surface of the wall, the parameter v w w w 0 h m º D is the ratio of the wall's magnetic diffusivity and thickness, and ΔB is the jump in magnetic field across the wall. For spatial scales of order unity, our v w -value of 10 3 implies that the time-scale for diffusion through the wall is 10 .
The initial central-η 0 value implies 10 w 3 t t h over the same spatial scale, hence .
We note that equations (1)-(4) describe the non-reduced form of MHD, where the toroidal component of magnetic field evolves according to equations (4) and (7), like the poloidal components. The outer boundary condition holds the total toroidal flux contained within the conducting shell constant, but magnetic flux may diffuse through the resistive wall, which supports poloidal and toroidal eddy currents that are subject to dissipation on the w t time-scale.
The conducting shell that surrounds the outer subdomain holds the normal component of B constant. The axisymmetric part of B n ·ˆalong that surface has contributions from external coils and from the initial plasma current. In addition, small error fields of order =respectively. Homogeneous Dirichlet conditions are set on the tangential components of V, and the inhomogeneous drift condition is the electric field within the resistive wall. This condition depends on time and space in the nonlinear computations, but its maximum value can be estimated from the initial equilibrium as v R I I I ,

@´-(ˆ· )
which is far smaller than the magnitude of flows that develop in the nonlinear computations, so this inhomogeneous condition is essentially equivalent to setting n V 0.
= · Mass and thermal energy escape the inner subdomain via diffusion and thermal conduction. Consistent with findings reported in [20], the concerns regarding the essentially impenetrable flow condition raised in [30] are not realized in these non-ideal computations, where magnetic flux is not frozen to cooled regions of plasma. However, we have found that results are sensitive to boundary conditions on T, due to its impact on resistivity in open-field regions [31]. More realistic modeling with boundary conditions inferred from sheath effects is under development, but it is beyond the scope of the present work and will be presented in a future publication.
The nonlinear computations start from the model equilibrium whose pressure and poloidal flux distribution are shown in figure 2(a). The equilibrium is the result of solving the Grad-Shafranov equation with the pressure and poloidal current profiles within the closed-flux region being P P P P a 4 1 and 8 where y is the normalized poloidal flux function that increases from 0 at the magnetic axis to 1 at the separatrix. The parameters P 1 10 , =´and I 2 e = produce the pressure and safety-factor profiles that are shown in figure 2(b). The equilibrium is determined numerically with the NIMEQ code [32], which has been modified to solve freeboundary computations. This up-down symmetric configuration has 15 axisymmetric coils outside the resistive wall at the positions shown in figure 1, and the two divertor coils are located at R 1.25 = and Z 1.45.

= 
The pressure distribution has a central β-value of 1.2%, and the particle density profile satisfies n n P P 0 0 1 5 = ( ) / so that edge-n is 1/10th of its central value.
The results section describes four initial-value computations. The first is a linear computation for the equilibrium of figure 2. The other three are nonlinear: a toroidally symmetric 2D computation with vertical forcing, a 3D computation with vertical forcing, and a 3D computation without vertical forcing. All are solved with the NIMROD code [28] using the implicit leapfrog time-advance that is described in [33] and the stabilization scheme from [34]. They use the same mesh over the R-Z plane, which has approximately 40 000 biquadratic elements in the inner subdomain and approximately 20 000 in the outer subdomain with some degree of concentration near the resistive wall in both. Vectors are expanded in cylindrical components with H 1 -elements, where the magnetic divergence constraint can only be satisfied approximately. However, even without considering the large uniform RB I e = f field, the error diffusion term in equation (4) keeps the divergence error smaller than one part in 1500 during the most violent phase of the 3D computations. Overall, results computed with the same mesh having bicubic basis functions are similar, providing evidence of spatial convergence over the R-Z plane, but running these computations through the entire CQ is computationally prohibitive at present. Toroidal variations in 3D NIMROD computations are represented by finite Fourier series. The toroidal resolution using wavenumbers n 0 21   in the 3D computations described here is based on experience with scanning resolution over a series of computations. Numerical quadratures over the toroidal angle use 64 evenly spaced points to preclude aliasing of quadratic nonlinearities. These computations include additional numerical damping of the n 20 = and n 21 = components at rates of 5 10 3 and 1 10 , 2 -respectively, to help avoid aliasing from high-order nonlinearities. NIMROD has been verified and benchmarked for many different models and applications (for example, 28,33,34). The resistive-wall implementation has been verified on cylindrical resistive-wall mode computations and is presently being benchmarked with the M3D-C1 code on a VDE instability using data from an NSTX discharge [35].
Growth-rates for the initial configuration are computed from the linearized versions of equations (1)-(4), where the equilibrium is a fixed distribution. In the nonlinear computations, the equilibrium density and temperature distributions become the initial conditions for n and T for the axisymmetric part of the solution. The axisymmetric part of the magnetic field is decomposed. The poloidal and toroidal fields from the plasma current become part of the initial conditions, and the externally generated toroidal field, i.e. the uniform RB I e = f part, is held constant. In the un-forced 3D computation, the poloidal field from all 15 external coils is also held constant. In the axisymmetric and forced 3D computations, field from the upper divertor coil, the unfilled rectangle in figure 1, is also made part of the initial condition, and field from the other 14 coils is held constant. This effectively turns the upper divertor coil off in these computations, leaving eddy current in the resistive wall in place of the upper divertor coil current. Also, the nonlinear computations are not driven from loop voltage or other external sources, and even without vertical displacement or asymmetric perturbation, this would lead to a gradual decay of plasma current and thermal energy from the start of the computations.

Results
For completeness, we first discuss the linear stability properties of the initial state without the forced VDE. We then describe general properties of the nonlinear evolution in the axisymmetric and 3D computations. Analysis of the current density and forces on the resistive wall follow.

Linear stability of the initial state
The equilibrium profile shown in figure 2 is linearly unstable to an external m 4, = n 1 = mode that has a growth rate of 4 10 .  where primes indicate d d , y is nontrivial for 1 y  in the equilibrium. Near the separatrix of this equilibrium, the above relation is dominated by the first term on the right. The fact that the mode is distributed over the poloidal angle indicates external-kink behavior, likely associated with this edge current density and its discontinuity across the separatrix. In our nonlinear computations, the relatively cool edge region just inside the separatrix is subject to resistive diffusion. However, as discussed in the next section, contact with the wall that results from vertical displacement tends to re-sharpen the edge.

Nonlinear evolution
The decay of eddy current in the nonlinear computations with forced vertical displacement transitions the configuration from being diverted to being limited within the first 1000 τ A . Over the first half of this period of time, the evolution of plasma current I p ( ) and thermal energy is nearly identical in the axisymmetric computation and the forced 3D computation. Figure 4 shows that discrepancies between the two computations develop over the second half of this period and increase thereafter, especially for the thermal energy. The spreading of the current-density distribution in the forced 3D computation, discussed later, weakens the attraction between the lower divertor coil and the plasma current, and figure 4(c) shows that the vertical motion is slowed as a result.
The asymmetric perturbations in the forced 3D computation achieve their maximum amplitude during the period t 500 1100 , A A   t t and as evident from the multiple local maxima in the fluctuation spectra shown in figure 5, saturation is not a simple process. The distortion of the plasma cross-section during this early saturation phase, shown by plasma pressure in figure 6, indicates that the m 3 = perturbation dominates prior to t 500 .
However, the m 2 = perturbation dominates by t 600 .
Plots of the n 1 = component of pressure (not shown) support the finding that the dominant poloidal wavenumber changes over time. The perturbations also alter the magnetic topology, producing magnetic islands near the edge of the distorting region of confinement and shrinking the volume of closed field lines. While the closed-flux region of the axisymmetric computation also decreases, the decrease in that case results from contact with the resistive wall, alone. This effect also occurs in the forced 3D computation, but chaotic scattering from the asymmetric perturbations is more significant, and anisotropic thermal conduction along the scattered field lines increases the rate of thermal energy loss. By t 1000 , there are no closed flux surfaces remaining in the forced 3D computation. The rate of thermal energy loss is largest at this time, when the highest-temperature region loses closed-flux confinement.
The 3D computation without vertical forcing also shows significant MHD activity. Resistive diffusion of the edge current again suppresses the m 4 = mode and excites the m 3 = external mode, but the core remains vertically centered ( figure 4(c)). The MHD activity is also less virulent. Comparing figures 5 and 7 shows that the initial saturation takes approximately twice as long without the wall contact. The m 2 = activity arises later in time, and figure 8 shows that while it again destroys all magnetic flux surfaces, there is recovery of a small central region of closed flux. This recovery does not occur in the computation with forced displacement. Noting from figure 4 that without wall contact, the plasma current continues with a relatively slow decay after thermal energy is diminished, the characterization of major versus minor disruption [24] generally describes the different behavior in our 3D computations with and without vertical forcing.
As a description of VDE evolution in the absence of asymmetric instabilities, the axisymmetric computation provides information on the source of free energy for MHD activity in the forced 3D computation during its vertical transient. With the condition of , contact with the wall scrapes-off the outer part of the equilibrium, while leaving the core relatively unchanged. Figure 9(a) shows the effect on the safety factor profile, and the most striking feature is the decrease in edge q-values over time. We surmise that the effect changes resonance conditions for different wavenumbers in the forced 3D computation. Moreover, the loss of  Figure 9(b) shows that the resulting edge-λ is nearly as large as the central λ-value, considerably larger than the initial edge-λ. The narrowness of this edge-current layer implies free energy for current-gradient-driven MHD modes, as also noted in [21], and enhances kink-type behavior, evidently followed by magnetic tearing, when asymmetries are allowed. Between the two 3D computations, this effect is only present in the one with vertical forcing, and it strengthens the MHD activity in that case. Contact with the wall also increases the edge pressure gradient in the axisymmetric computation. For q a 1, > ( ) this edge pressure gradient would tend to drive ballooning, but no distinct high-n activity appears in the forced 3D computation (figure 5). The loss of    A t = Each frame shows results from the 3D computation without vertical forcing at toroidal angle 0. f = edge thermal confinement from the low-n perturbations in the 3D computation precludes the development of a clear edge pressure pedestal, and ballooning does not arise.

Current bump and distribution of current density
The separation of the plasma current histories for the forced cases ( figure 4(a)) starts from the initial saturation of the asymmetric instabilities, which increases I p in the 3D case. Comparing figures 4 and 5(a) shows that the I p bump continues after t 1000 , which is when magnetic fluctuations and the rate of thermal energy loss are largest. This aspect is consistent with the discussion of the JET current spike in figure 29 of [25] and with the data for a major disruption in TFTR shown in figure 1 of [26]. Comparing the two frames in figure 10 shows that the parallel current density distribution, computed from the toroidally symmetric part of B, spreads while I p increases. The concentration of poloidal flux within the central-plasma region also decreases over this time, which can be observed from the decreased density of the equally spaced, poloidal-flux contours. The increase in I p simultaneous with the spreading of poloidal flux necessarily implies a reduction of inductance. This effect does not occur in the axisymmetric computation, where the distribution of poloidal flux within the plasma core remains relatively constant while the edge is removed by contact with the wall (figure 9(b)). Similar behavior has been noted in simulations of MHD activity excited by edge impurities for disruption mitigation [36,37].
Part of the change of poloidal flux can be attributed to resistive dissipation that is enhanced by the decreasing  temperature during the TQ. However, a stronger effect results from the correlation of magnetic and flow-velocity fluctuations, the MHD dynamo effect E V B , f º -á´ñ˜where 'fluctuation' refers to the toroidally asymmetric component of a field. Conceptually, this stems from mean-field theory of MHD for astrophysics [38] and was first appreciated for magnetic confinement in the context of magnetic relaxation in reversed-field pinches [39,40]. It was later applied to analyses of driven spheromaks [41,42], helicity injection in tokamaks [43] and spherical tokamaks [44][45][46], and the tokamak hybrid mode [47]. Averaging Faraday's law, t B E, ¶á ñ ¶ = -´á ñ and considering low frequencies leads to the Poynting theorem for the symmetric field, The term E J f á ñ · is included in the right side of equation (9) and is part of a nonlinear magnetic energy transport process that acts through MHD fluctuations [43,48]. Figures 11(a), (b) show this power density and the tor- which is when the I p bump begins in the forced 3D computation. To gauge the magnitude of E , f we also plot an approximation for the resistive electric field that acts on the toroidally averaged current density J . á ñ f Figure 11(c) shows the product of J á ñ f and resistivity which is computed with the averaged temperature. In the plot, locations where the current density values are below 1% of the maximum are set to 0 so that the image is not dominated by the outer, cold region where resistivity is orders of magnitude larger. Comparing figures 11(b), (c) shows that where the dynamo electric field is active, it is more than an order of magnitude larger than the resistive dissipation of J , á ñ f and approximately 100 times that of J h f on the magnetic axis at the start of the computation. Thus, the E f contribution generates net toroidal electric field, which alters the poloidal-flux and current-density distributions. The E J f á ñ · power density is mostly positive near the core and negative in the edge, which removes energy from B á ñ in the core and deposits it in the edge. This process spreads the current-density distribution. In addition, the fluctuationinduced loop voltage from E f f ·ˆredistributes poloidal flux; locations where RE f f ·ˆdecreases in the flux-normal direction push the poloidal flux distribution outward.
The consideration of correlated fluctuations helps describes the nonlinear effects of asymmetric instabilities on the toroidally averaged field, but it does not provide a full sense of their influence. Section 3.1 already noted that chaotic scattering extends over the entire plasma volume by the condition shown in figure 10(b). The local evaluation of λ at this time, figure 12(a), identifies considerable spatial structure that mean-field relaxation analysis does not. The spatial structure includes current sheets and regions of anti-parallel current, and these spatial structures extend beyond the region of closed poloidal flux contours shown in figure 10(b). This serves as a reminder that the poloidal flux function does not represent the magnetic topology of the 3D evolving field. The temperature distribution in the same poloidal plane, figure 12(b), also shows hot structures that extend beyond the region of closed flux contours. While the temperature distribution, and hence the electrical conductivity distribution, is not identical to the λ-distribution, a strong correlation is apparent from the figure.

Forces on the resistive wall
Our computation of net forces on the resistive wall follows the analysis given by Pustovitov in [27]. Two essential properties used in the analysis are that (1) together, the plasma and resistive wall form an electrically isolated system and (2) plasma inertia is negligible for the timescale of the evolution. The net force on the wall can then be computed from a surface integral of magnetic stresses over the outside of the wall. The thin-shell approximation of our model actually necessitates the use of stresses, because current per unit cross-section area is infinitely large. Here, we check the consistency of our where e ĵ is a Cartesian unit vector and the integral is broken into separate contributions over the inner and outer surfaces of the resistive wall. Noting that the net force acting on the plasma is F in - and, therefore, the implication from property 2 of F 0, in   we expect F F . j out j  The net horizontal force, which only results with toroidal asymmetries, peaks at a normalized amplitude of 0.08 at t 3000 A t = in the forced 3D computation ( figure 13(a)). To provide a sense of scale, if this force was generated from m 1, = n 1 = tilting of the plasma current, Noll's relation [3], would imply that the toroidally asymmetric displacement would only be 0.011. Nonetheless, we can check whether Pustovitov's property 2 holds for this relatively weak force. We find that the largest instantaneous realization of F F in out is the value of 3 10 2 that occurs briefly at t 600 .
A t @ Using the maximum over time of each integral, we find F F max max 5 10 . in out 3 @´-( ) ( ) These ratios are larger than what is discussed in [27] for cases of large forces in JET, but even here they are consistent with the analysis.
Separating the horizontal components of force, figure 13(b), shows that the orientation of the horizontal force vector changes over the simulated displacement event. This has been observed in previous MHD simulations of disruption [49], but effects outside the scope of MHD are likely needed to predict the rotation of forces.
Comparing figures 5(a) and 13(a) shows that the largest horizontal force on the wall occurs much later than the peak of the asymmetric magnetic fluctuations. In fact, the force  2 , here n is the unit outward normal along the wall. The maximum force density is larger at t 3000 A t = than at t 971 , A t = despite the smaller internal fluctuation energy. During the intervening time, the discharge has moved somewhat closer to the lower divertor, and the magnetic perturbations have more time to penetrate the resistive wall. The distribution over the toroidal angle has also changed such that there is less cancellation when integrating over the inboard and outboard sides.

Discussion
Visco-resistive MHD modeling greatly simplifies the entirety of physics that influences disruptions in tokamak experiments. Approximations that directly influence the computational results presented above are the lack of radiation modeling, the Dirichlet boundary conditions on temperature, and the lack of runawayelectron (RE) effects. The modeled TQ resulting from local anisotropic thermal transport in the presence of changing magnetic topology is fast in our forced 3D computation, relative to the CQ. However, a realistic kinetic model for parallel heat transport, or even a larger value of , c || would further separate the TQ and CQ timescales. Radiation from impurities is also expected to have an important role in the TQ, and as noted, it has not been modeled. The present computations rely on thermal conduction to the resistive wall, where the imposed Dirichlet condition maintains low temperature. Parallel conduction quickly cools open field-lines, leading to the thin halo region in our axisymmetric result. From [21] and our work in this area [31], we know that computed VDE evolution is sensitive to edge temperature, so more comprehensive modeling is needed. For example, boundary conditions for edge turbulence modeling have been developed using conditions at the magnetic pre-sheath entrance [50]. We are adapting this more realistic approach for our disruption computations [31], but because electrons are then largely insulated, it may also require a radiation model to expel electron heat that is conducted toward the surface. Finally, experimental CQs are often extended by the formation of RE beams [1,25,51]. We expect that practical computations in the near future will rely on reduced modeling of REs, such as the one developed in [52] and already applied in recent M3D simulations.
The possibility of kink-induced surface current, progressing in the direction opposite to that of the plasma current and traveling partly through conducting structures has been raised in [5]. While this wall-touching kink mode (WTKM) is predicted to be most problematic for the m 1, = n 1 = mode, which arises with sufficiently small q-values, the physics of reversed surface currents also arises with other external modes [53,54]. In Section 3.3, we note the existence of reversed parallel current density at the time of maximum-I p in the forced 3D computation, but the nonlinearly distorted current profile makes it difficult to relate the reversed parallel current density with any particular kink mode. However, at t 473 A t = when the m 3, = n 1 = mode dominates, the external kink distortion is clearer. The 0.15 l =isosurface in figure 15 is computed from the total current density at this time, and it lies along the outward bulge from the kink distortion. The existence of this helical current channel demonstrates that the physics of reversed surface currents is within the scope of resistive-MHD models using T-dependent resistivity to track an effective plasma surface.

Conclusions
This computational study of an asymmetric VDE in an idealized configuration considers a number of physical effects that are important for tokamak experiments. The 3D result reproduces distinct timescales for the TQ and CQ phases, despite the modest parallel heat conduction and absence of RE physics. The thermal energy collapse results from chaotic scattering of magnetic field lines associated with increasingly pervasive MHD modes and nonlinear coupling. In addition, the MHD activity generates nonlinear processes that redistribute the parallel current density. We argue that the effects on the symmetric components of B and J can be understood through analysis of the fluctuation-induced E f electric field, which is the MHD dynamo effect known from studies of current drive in RFPs, spheromaks, and other configurations. Because the fluctuations do not dissipate rapidly, the decrease in inductance and the resulting bump in plasma current that are associated with current-density spreading persist beyond the time when the rate of thermal energy loss is greatest. We also find that relative to the axisymmetric result, the spreading of current reduces the attractive force exerted by the current in the divertor coil. Thus, final termination of the plasma current through contact with the wall takes longer with the asymmetric distortions.
We have computed forces on the resistive wall by integrating the Maxwell stress over its surface. In the modeled case, the net horizontal force is not large relative to conditions of significant tilting of the toroidal current. Nonetheless, we find that the net force over the inner surface of the wall is much smaller than that over its outer surface. This is consistent with the analysis of [27], which argues that the plasma and wall must remain in approximate force balance on timescales of interest, so the net force can be accurately computed by integrating stress over the outer surface alone. Our force computation takes into account contributions from asymmetric conduction currents flowing from the plasma to the resistive wall, including any reversed currents of externalkink dynamics, demonstrated in section 4, which underlie the WTKM theory [5]. While the computed scenario does not develop a large m 1, = n 1 = distortion, our results support the possibility of reproducing WTKM physics. This work also supports the prospect of using Eulerianframe computation for simulating disruptions. We expect Lagrangian and other moving-frame approaches to be more efficient in cases where the plasma torus remains intact. However, instances of significant distortion and magnetic topology change, such as what our model case produces, would tangle Lagrangian meshes or at least remove many of the computational benefits of moving meshes. Nonetheless, 3D computations, such as the two presented here, are computationally intensive, and improved solver efficiencies and better use of recent computer hardware developments are needed to make these computations more practical. Predictive simulations of asymmetric VDEs are also expected to require detailed representations of external conductors, such as what has been developed for resistive-wall mode studies [55,56] and as discussed in [57], together with the physical model developments described in section 4.