The implementation of a resistive-wall extension to the non-linear MHD-code
JOREK via a coupling to the vacuum-field code STARWALL is presented along with
first applications and benchmark results. Also, non-linear saturation in the
presence of a resistive wall is demonstrated. After completion of the ongoing
verification process, this code extension will allow to perform non-linear
simulations of MHD instabilities in the presence of three-dimensional resistive
walls with holes for limited and X-point plasmas.Comment: Contribution for "Theory Of Fusion Plasmas, Joint Varenna - Lausanne
  International Workshop, Villa Monastero, Varenna, Italy (27.-31.8.2012)",
  accepted for publication in Journal of Physics Conference Serie

Chu M S

E Nardon

E Strumberger

G T A Huysmans

Hertout P

Hirsch C

Huysmans G

Hölzl M

I Chapman

Igochine V

K Lackner

M Hölzl

Merkel P

P Merkel

R McAdams

S Günter

Strumberger E

Journal of Physics : Conference Series

English

arXiv

Crossref

Journal of Physics Conference Series

Coupling JOREK and STARWALL Codes for Non-linear Resistive-wall Simulations

Hölzl, M.

Merkel, P.

Huysmans, G.

Nardon, E.

Strumberger, E.

McAdams, R.

Chapman, I.

Günter, S.

Lackner, K.

MPG.PuRe

Coupling JOREK and STARWALL Codes forNon-linear Resistive-wall SimulationsM Hölzl1, P Merkel1, G T A Huysmans2, E Nardon3,E Strumberger1, R McAdams4,5, I Chapman5, S Günter1, K Lackner11Max-Planck-Institute for Plasmaphysics, EURATOM Association, Boltzmannstr. 2, 85748Garching, Germany2ITER Organisation, Route de Vinon sur Verdon, St-Paul-lez-Durance, France3CEA, IRFM, CEA Cadarache, F-13108 St Paul-lez-Durance, France4York Plasma Institute, University of York, York, YO10 5DD, UK5Euratom/CCFE Fusion Association, Culham Science Centre, Abingdon, OX14 3DB, UKE-mail: mhoelzl@ipp.mpg.deAbstract. The implementation of a resistive-wall extension to the non-linear MHD-codeJOREK via a coupling to the vacuum-field code STARWALL is presented along with firstapplications and benchmark results. Also, non-linear saturation in the presence of a resistivewall is demonstrated. After completion of the ongoing verification process, this code extensionwill allow to perform non-linear simulations of MHD instabilities in the presence of three-dimensional resistive walls with holes for limited and X-point plasmas.1. IntroductionPlasma instabilities growing on fast time-scales such as vertical displacement events, disruptions,external kink modes, or edge localized modes are connected with a time-dependent magneticperturbation outside the plasma. The mirror currents induced in conducting structures aroundthe plasma1 by the varying field perturbations act back onto the instabilities. Thus, linear andnon-linear dynamics of the plasma can be strongly influenced by the wall. An external kinkmode driven by strong edge current densities, for instance, is converted into a resistive wallmode (RWM) in the presence of a conducting wall close to the plasma edge and grows on thetimescale of resistive wall-current decay as described in great detail in Refs. [1, 2].We describe the status of a resistive wall extension for the non-linear MHD code JOREK [3, 4]as well as first benchmark results. The implementation is done by coupling JOREK to a modifiedversion of the STARWALL code [5] which determines the magnetic field in the vacuum regionin the presence of three-dimensional walls with holes. The “response” of conducting structuresto magnetic field perturbations is represented by response-matrices computed by STARWALL.These matrices are used in JOREK to keep track of wall-currents and calculate the tangentialmagnetic field required for the modified boundary conditions as discussed in Section 3.2.The article provides some information about the STARWALL code in Section 2. The couplingequations are derived in Section 3. First benchmarks of the resistive wall extension are presentedin Section 4. Conclusions and a brief outlook are given in Section 5.1 For simplicity, all conducting structures are denoted wall in the following.2. The STARWALL CodeSTARWALL solves the vacuum magnetic field equation outside the JOREK computationaldomain in presence of a three-dimensional conducting wall with holes as a Neumann-likeproblem. The continuity of the magnetic field component normal to the boundary of the JOREKcomputational domain (called interface) is used as a boundary condition. The wall is representedby infinitely thin triangles which allows to approximate realistic wall structures very well usingan effective resistance ηw/dw where ηw and dw denote the resistivity of the wall material2 andits thickness, respectively. Wall currents are assumed constant within each triangle such thatthey can be described by current potentials Yk at the triangle nodes [5]3.In case of an ideally conducting wall, solving the vacuum field equation results in an algebraicexpression for the magnetic field component tangential to the interface4 in terms of the normalcomponent. Wall currents don’t need to be considered explicitly as they are instantaneouslygiven by the normal field at the interface. Using the ideal wall response-matrix M̂ id calculatedby STARWALL, this can be written in the following way5:Btan =∑ibi Btan,i =∑ibi∑jM̂ idi,j Ψj . (1)Here, bi denotes a JOREK basis function consisting of a 1D Bezier basis function along thepoloidal direction on the interface and Fourier basis function in toroidal direction. Ψj denotespoloidal flux coefficients. The poloidal flux at the interface is given by Ψ =∑j bj Ψj .When a resistive wall is considered, wall currents cannot be eliminated anymore. In this case,the tangential magnetic field is given byBtan =∑ibi∑jM̂ eei,j Ψj +∑kM̂ eyi,k Yk , (2)where wall currents evolve in time according toẎk = −ηwdwM̂yyk,k Yk −∑jM̂yek,j Ψ̇j . (3)Here, M̂ ee, M̂ ey, M̂ye, and M̂yy denote resistive response matrices determined by STARWALL.Indices i and j run over all boundary degrees of freedom of the respective variable and k overall wall current potentials. For cross-checking, a relation between ideal and resistive responsematrices can be derived by letting ηw →∞ (Appendix A).Information about the STARWALL code coupled with CASTOR (sometimes calledSTARWALL C) can be found in Refs. [5, 7]. It is similar to the STARWALL code coupled withJOREK (sometimes called STARWALL J). An article describing more details is in preparationby the STARWALL author Peter Merkel.3. Implementation in JOREKFor the coupling of JOREK and STARWALL, Eqs. (2) and (3) are discretized in time as describedin Section 3.1. The vacuum response then enters into a natural boundary condition in JOREKas discussed in Section 3.2.2 Wall resistivity is normalized the same way as plasma resistivity (see Ref. [6] for JOREK normalizations).3 For simplicity, we will (inaccurately) speak about wall currents when wall current potentials are referred to.4 Btan is the component of (B−Bφ) tangential to the interface, where Bφ denotes the toroidal field component.5 For practical reasons, the algebraic expression is written in terms of the poloidal flux Ψ instead of the normalfield component (which can be calculated from Ψ, of course).3.1. Time-DiscretizationEquation (2) is evaluated at the new timestep n+ 1 and discretizations Ψn+1 = Ψn + δΨn andY n+1 = Y n + δY n are used, where superscripts indicate evaluation at the given timestep. Thetangential magnetic field is thus given byBn+1tan =∑ibi∑jM̂ eei,j(Ψnj + δΨnj)+∑kM̂ eyi,k (Ynk + δYnk ) . (4)The general time-evolution scheme described in Ref. [8],[(1 + ξ)(∂X∂u)n−∆tθ(∂Z∂u)n]δun = ∆t Zn + ξ(∂X∂u)n−1δun−1, (5)is used for the JOREK equations written in the form ∂X(u)/∂t = Z(u), where u denotes thevector of JOREK physical variables like temperature or density, X and Z the left- respectivelyright-hand side terms of the time-evolution equations, and θ and ξ are numerical parameters6.The same general time-evolution method needs to be applied to the wall-current evolution(Eq. (3)) giving(1 + ξ)δY nk +∑jM̂yek,j δΨnj+ ∆t θ ηwdwM̂yyk,k δYnk=−∆t ηwdwM̂yyk,k Ynk + ξδY n−1k +∑jM̂yek,j δΨn−1j .(6)Terms with δY nk are brought to the left hand side. The equation then reads(1 + ξ + ∆tθηwdwM̂yyk,k)δY nk=− (1 + ξ)∑jM̂yek,j δΨnj −∆tηwdwM̂yyk,kYnk + ξδYn−1k + ξ∑jM̂yek,j δΨn−1j .(7)After solving for δY nk , one getsδY nk =∑jÂk,j δΨnj + B̂k,k Ynk + Ĉk,k δYn−1k +∑jD̂k,j δΨn−1j , (8)where some of the “derived response matrices”Ŝk,k = 1 + ξ + ∆tθηwdwM̂yyk,k D̂k,j = ξM̂yek,j/Ŝk,k Ĥi,j = M̂eei,jÂk,j = −(1 + ξ) M̂yek,j/Ŝk,k Êi,j = M̂eei,j +∑k M̂eyi,k Âk,j Ĵi,j =∑k M̂eyi,k D̂k,jB̂k,k = −∆tηwdw M̂yyk,k/Ŝk,k F̂i,k = M̂eyi,k(1 + B̂k,k) K̂k,l = −∆tM̂yyk,kŜl,kĈk,k = ξ/Ŝk,k Ĝi,k = M̂eyi,k Ĉk,k L̂i,l =∑k M̂eyi,k K̂k,l(9)6 Parameter values satisfying θ − ξ = 0.5 are required to guarantee second-order accuracy of the time-evolutionscheme. For instance, (θ = 0.5, ξ = 0) corresponds to a linearized Crank-Nicholson scheme and (θ = 1, ξ = 0.5)to the linearized Gears scheme.have been used. These matrices are computed at the beginning of a JOREK simulation from theSTARWALL response matrices. They need to be updated if one or more parameters enteringthe definitions have changed, e.g., ∆t, ξ, or ηw. Plugging Eq. (8) into Eq. (4) givesBn+1tan =∑ibi[∑j(M̂ eei,j +∑kM̂ eyi,k Âk,j)δΨnj +∑kM̂ eyi,k(1 + B̂k,k)Y nk+∑kM̂ eyi,kĈk,k δYn−1k +∑jM̂ eei,jΨnj +∑j∑kM̂ eyi,kD̂k,jδΨn−1j] (10)and, making use of Eq. (9), we getBn+1tan =∑ibi∑jÊi,j δΨnj +∑kF̂i,k Ynk +∑kĜi,k δYn−1k +∑jĤi,j Ψnj +∑jĴi,j δΨn−1j .(11)3.2. Boundary IntegralThe fixed boundary conditions for poloidal flux and plasma current corresponding to an ideallyconducting wall are removed from JOREK. As a consequence, a boundary integral in the currentdefinition equation resulting from partial integration (which vanishes in fixed-boundary JOREK)needs to be considered now. It can be written in terms of the tangential magnetic field suchthat we can insert Eq. (11) as a natural boundary condition. Details are given in the following.The current definition equation, j = ∆∗Ψ, is written in weak form as∫dVj∗lR2(j −∆∗Ψ) = 0. (12)where j∗l denotes the test-functions taken to be identical with the basis-functions. Using∆∗Ψ ≡ R2 ∇ · (R−2 ∇Ψ), we get∫dVj∗lR2j −∫dV j∗l ∇ ·(1R2∇Ψ)= 0, (13)where the second term of the resulting expression can be integrated by parts (∫dV a ∇ ·b =−∫dV ∇a ·b +∮dA a b · n̂) yielding∫dVj∗lR2j +∫dV1R2∇j∗l · ∇Ψ−∮dAj∗lR(∇Ψ · n̂/R)︸ ︷︷ ︸≡Btan= 0. (14)Here, n̂ denotes the unit vector normal to the interface and the tangential field is identified inthe boundary integral (refer to Appendix B for details). Eq. (11) is inserted into Eq. (14) andafter separating implicit and explicit terms, we have derived the form of the current equationimplemented in JOREK,∑ielem∫dVR2(j∗l δjn +∇j∗l · ∇δΨn)−∑ibnd∮dAj∗lR∑ibi∑jÊi,j δΨnj=−∑ielem∫dVR2(j∗l jn +∇j∗l · ∇Ψn)+∑ibnd∮dAj∗lR∑ibi∑k(F̂i,k Ynk + Ĝi,k δYn−1k)+∑j(Ĥi,j Ψnj + Ĵi,j δΨn−1j) ,(15) 0 0.5 1 1.5 2 2.5 3 3.5 0  0.2  0.4  0.6  0.8  1q-profilePsi normalizedFigure 1Left: JOREK flux-surfaces of a limiter equilibrium (limiter point shown as bluecross) are compared to CEDRES++ last closed flux-surface and magnetic axis.Right: Comparison of the respective q-profiles.with ∫dV =∑ielem∫du dv dφ J2 R, and∮dA =∑ibnd∫du dφ R√(R,u)2 + (Z,u)2, (16)where the integral over u and v is carried out by Gauss-quadrature:∫du X(u) =∑mGwmG X(umG) with weights wmG and positions umG . Coordinate u corresponds to theelement-local coordinate of the 2D Bezier elements oriented poloidally along the interface. Theother element-local coordinate is denoted v. J2 denotes the 2D-Jacobian in the poloidal plane.Wall-currents are updated after each time-step according to Eq. (8) to guarantee consistencywith the implicit time-stepping of JOREK.For the n = 0 component, poloidal field coils need to be taken into account properly. This isdone by replacing the expression for the tangential magnetic field, Eq. (2), byBtan =∑ibi∑jM̂ eei,j (Ψj −Ψcoilj ) +∑kM̂ eyi,k Yk +Bcoiltan,i . (17)Here, Ψcoilj and Bcoiltan,i denote time-independent coil contributions to poloidal flux and tangentialfield at the interface, respectively. This leads to two additional terms in Eq. (15). A similarboundary integral also occurs in the Grad-Shafranov equation giving the plasma equilibrium.Here, only the n = 0 component needs to be considered (axisymmetric equilibria in JOREK)and no time-discretization is required.4. First Benchmarks and ResultsIn this Section, we present first results of the ongoing benchmark effort for code validation.0.01 0.1 1rwall [m] - 1102103growth rate [s-1]CASTORCASTOR + STARWALL_CJOREK + STARWALL_JηJorek=1.e-5ηJorek=1.e-6ηJorek=1.e-7Figure 2Left: Linear growth-rates of a tearing-mode in a circular plasma with concentricideal wall are compared between CASTOR using its own vacuum field module,CASTOR coupled with STARWALL, and JOREK coupled with STARWALL. Verygood agreement is observed for a variety of different plasma resistivities and plasma-wall distances.Right: The current potential distribution (arbitrary units) reflecting the 2/1 tearingmode structure is plotted during the linear phase (ηJorek = 1 · 10−6, rw = 1.2).4.1. Free-boundary equilibriumComputing a free-boundary equilibrium with JOREK and comparing it to the results ofanother equilibrium solver allows to check some parts specific for n = 0 like poloidal field-coil contributions (see Sec. 3.2). We are considering an ITER-like limiter plasma. As seenfrom Figure 1, flux surfaces and q-profile of the JOREK equilibrium agree very well with theequilibrium computed by the CEDRES++ code [9]. Slightly different discretizations of thepoloidal field coils might be responsible for small remaining differences.4.2. Tearing ModeA 2/1 tearing mode in a circular plasma with major radius R = 10, minor radius a = 1 anduniform plasma resistivity surrounded by an ideally conducting wall is considered. Figure 2shows excellent agreement between JOREK and the linear CASTOR code [10]. Repeatingthe simulation for a resistive wall with zero resistivity is physically identical, of course, butnumerically treated in a different way. The exact agreement we get here proves consistency.4.3. Resistive Wall ModeIn this Section, results for 2/1 resistive wall modes (RWMs) in a circular plasma with minorradius a = 1 and major radius R = 10 surrounded by a resistive wall are presented. The leftpart of Figure 3 shows linear growth rates obtained for different wall radii rw as a function of thewall resistivity. Clearly, the modes cannot be stabilized beyond a certain wall distance from theplasma even at very low wall resistivities. The right part of the Figure shows energy time-tracesduring the linear and saturation phase of a typical simulation. 1e-05 0.0001 0.001 1e-07  1e-06  1e-05 0.0001  0.001  0.01  0.1  1linear growth ratewall resistivityrwall=1.1rwall=1.2rwall=1.3rwall=1.4rwall=1.5rwall=1.6rwall=1.8rwall=2.0rwall=4.0rwall=8.0 1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 100000  120000  140000  160000  180000normalized energynormalized timeEmag,01Ekin,00Ekin,01Figure 3Left: Linear RWM growth rates extracted from the simulations are shown for avariety of wall radii and wall resistivities.Right: Energy time-traces are shown for linear and saturation phase of an RWM(rw = 1.2, ηwall = 1 · 10−4). In saturation, wall currents on the high-field side have anamplitude comparable to the low-field side while they are smaller in the linear phase.5. Conclusions and OutlookThe ongoing implementation and verification of a resistive wall model in the non-linear MHD-code JOREK was summarized. Benchmarks for a freeboundary equilibrium and tearing modecases show good agreement. Simulations of the linear and non-linear phase of RWMs werepresented (comparison to analytical theory and linear codes ongoing).Benchmarking will be continued and extended to non-linear comparisons. Realistic X-pointgeometries will be considered requiring a special treatment of grid corners. After completion,the code can be applied to a variety of MHD instabilities interacting with conducting structureslike resistive wall modes, edge localized modes, vertical displacement events or disruptions.AcknowledgementsSimulations were mostly carried out on the HPC-FF computing cluster in Jülich, Germany, andon the IFERC-CSC Helios computing cluster in Rokkasho-Mura, Japan.References[1] Chu M S and Okabayashi M 2010 Plasma Physics and Controlled Fusion 52 123001[2] Igochine V 2012 Nuclear Fusion 52 074010[3] Czarny O and Huysmans G 2008 J. Comput. Phys. 227 7423 – 7445[4] Huysmans G, Pamela S, Beurskens M, Becoulet M and van der Plas E 2010 Proceedings of the 23rd IAEAFusion Energy Conference (Daejon, South Korea) THS/7-1[5] Merkel P and Sempf M 2006 Proceedings of the 21st IAEA Fusion Energy Conference (Chengdu, China)TH/P3-8 URL www-naweb.iaea.org/napc/physics/FEC/FEC2006/papers/th_p3-8.pdf[6] Hölzl M, Günter S, Wenninger R, Müller W C, Huysmans G, Lackner K, Krebs I and the ASDEX UpgradeTeam 2012 Physics of Plasmas (accepted)[7] Strumberger E, Merkel P, Tichmann C, and Günter S 2011 Proceedings of the 38th EPS Conference onPlasma Physics (Strasbourg, France) p5.082 URL http://ocs.ciemat.es/EPS2011PAP/pdf/P5.082.pdf[8] Hirsch C 1989 Numerical Computation of Internal and External Flows, Volume 1, Fundamentals of NumericalDiscretization (Wiley) ISBN 978-0-471-92385-5[9] Hertout P, Boulbe C, Nardon E, Blum J, Brémond S, Bucalossi J, Faugeras B, Grandgirard V and MoreauP 2011” Fusion Engineering and Design 86 1045–1048 Proceedings of the 26th Symposium of FusionTechnology (SOFT-26)[10] Kerner W, Goedbloed J, Huysmans G, Poedts S and Schwarz E 1998 J. Comput. Phys. 142 271 – 303Appendix A. Relation between ideal and resistive response matricesBetween the ideal and resistive response-matrices, the relationM̂ idi,j ≡ M̂ eei,j −∑kM̂ eyi,k M̂yek,j (A.1)must hold which becomes obvious when letting ηw → 0 and inserting Eq. (3) into Eq. (2). Also,the no-wall response can easily be identified asM̂nwi,j ≡ M̂ eei,j (A.2)when letting Yk → 0. Eqs. (A.1) and (A.2) can be used to cross-check some parts of theimplementation as both forms of the wall-treatment with eliminated wall-currents (Eq. (1)) andexplicitly treated wall-currents (Eqs. (2) and (3)) are implemented in JOREK.Appendix B. Magnetic field tangential to the interfaceThe magnetic field in the JOREK reduced-MHD model is given byB =F0Rêφ +1R∇Ψ× êφ. (B.1)Its component tangential to the boundary of the JOREK computational domain (in the poloidalplane) can be determined fromBtan = (B× n̂) · êφ, (B.2)where n̂ denotes the unit normal vector of the interface and êφ the normalized toroidal basisvector. Inserting Eq. (B.1) and using vector identities, this can be written in the formBtan = −1Rêφ · [n̂× (∇Ψ× êφ)]= − 1Rêφ ·(n̂ · êφ)︸ ︷︷ ︸≡0∇Ψ− (n̂ · ∇Ψ)êφ=1Rn̂ · ∇Ψ,(B.3)which can be identified in the boundary integral of Eq. (14) and replaced by the STARWALLresponse.