A semi-analytical solver for the Grad-Shafranov equation
AA semi-analytical solver for the Grad-Shafranov equation
D. Ciro ∗ and I. L. Caldas † Departamento de F´ısica Aplicada, Universidade de S˜ao Paulo, 05508-090, S˜ao Paulo, Brazil.
In toroidally confined plasmas, the Grad-Shafranov equation, in general a non-linear PDE, de-scribes the hydromagnetic equilibrium of the system. This equation becomes linear when the kineticpressure is proportional to the poloidal magnetic flux and the squared poloidal current is a quadraticfunction of it. In this work, the eigenvalue of the associated homogeneous equation is related withthe safety factor on the magnetic axis, the plasma beta and the Shafranov shift, then, the adjustableparameters of the particular solution are bounded through physical constrains. The poloidal mag-netic flux becomes a linear superposition of independent solutions and its parameters are adjustedwith a non-linear fitting algorithm. This method is used to find hydromagnetic equilibria with nor-mal and reversed magnetic shear and defined values of the elongation, triangularity, aspect-ratio,and X-point(s). The resultant toroidal and poloidal beta, the safety factor at the 95% flux surfaceand the plasma current are in agreement with usual experimental values for high beta dischargesand the model can be used locally to describe reversed magnetic shear equilibria.
I. INTRODUCTION
The solution of the Grad-Shafranov equation [1–3] pro-vides the magnetic field, the current density, and the ki-netic pressure inside an axisymmetric plasma in hydro-magnetic equilibrium. Having analytical solutions to thisequation is convenient to configure physical equilibria asa basis for theoretical studies of transport, waves, andstability. It also allows to estimate the external mag-netic field configuration necessary to confine a toroidalplasma with specified parameters [4].The Grad-Shafranov equation is an elliptic PDE forthe poloidal magnetic flux ψ that labels the magneticsurfaces in an axisymmetric plasma equilibrium. Theequation contains two arbitrary functions p ( ψ ) and F ( ψ )that specify the dependence of the kinetic pressure andthe poloidal plasma current on the magnetic flux ψ .Accordingly, the Grad-Shafranov equation is, in gen-eral, a nonlinear PDE and its solution rely on numeri-cal methods. However, for some choices of the arbitraryfunctions the equation becomes linear and separable, andthe boundary value problem can be solved by superposi-tion of independent solutions. Various classes of analyt-ical solutions have been introduced along the years [5–11], usually involving linear and quadratic dependencesof p and F on ψ . These choices impose some inheritrestrictions on the possible current density profiles, andrelations between the physical parameters, but, in gen-eral, provide good magnetic topology and safety factorprofiles.In this work we study a class of exact solutions result-ing when the pressure is a linear function of ψ and thesquared poloidal current is a quadratic function of ψ [9].These solutions are characterized by an eigenvalue thatis related with the equilibrium parameters of an axisym-metric plasma, specifically, the Shafranov shift ∆, the ∗ [email protected] † [email protected] safety factor at the magnetic axis q and the fraction ofdiamagnetic reduction of the toroidal field. This rela-tion allow us to establish the limits of the model and itsparameters, and, using physical arguments, we define aregion of consistent solutions in the parameters space.Then, we employ the analytical solutions of the Grad-Shafranov equation to build a predictive solver that al-lows to specify the geometrical properties of the toroidalplasma, namely, the aspect ratio (cid:15) , triangularity δ , elon-gation κ and X-point, for a given set of physical param-eters { ∆ , q , f, β } .The treatment presented in this work potentiate theuse of this class of solution for predictive equilibrium cal-culations and to our knowledge the relations introducedhere have not been presented elsewhere and provide avaluable tool for the construction of analytical equilibria.The whole treatment was done in dimensionless variables,so that the results can be properly scaled to any case ofinterest using a couple of machine parameters, the majorradius R and the toroidal vacuum field B .The manuscript is organized as follows. In the Sec-tion II we present a short survey on hydromagnetic equi-librium with the model considered in this work andthe analytical solutions to the Grad-Shafranov equation,then, in the Section III, we study the relations betweenthe physical parameters and the solution parameters. Inthe Section IV, we introduce a numerical method to solvethe boundary value problem, and, in the Section V, wegive some examples of equilibrium calculations performedwith this method and present our conclusions in the Sec-tion VI. II. ANALYTICAL SOLUTIONS
For ideal plasmas, the equilibrium between the kineticand magnetic forces requires that (cid:126)j × (cid:126)B = ∇ p, (1)everywhere inside the plasma. Here, p , (cid:126)j and (cid:126)B are thekinetic plasma pressure, current density and magnetic a r X i v : . [ phy s i c s . p l a s m - ph ] S e p field respectively. Assuming that the system is axisym-metric and using the Amp`ere’s law the equilibrium prob-lem is reduced to the Grad-Shafranov equation [1–3] R ∇ · (cid:18) ∇ ψR (cid:19) = − µ R dpdψ − F dFdψ = − µ Rj φ . (2)Here, R is the distance to the symmetry axis and ψ ( R, z )is the poloidal magnetic flux, calculated though a diskof radius R at the height z . The arbitrary function p ( ψ ) represent the kinetic pressure at the level surface ψ ( R, z ) = const. , and F ( ψ ) is the poloidal plasma cur-rent enclosed by that surface. The magnetic field lines lieon the magnetic surfaces ψ ( R, z ) = const. , that are alsoisobarics. Finally, the second equality in (2) relates thetoroidal (azimuthal) current density with the arbitraryfunctions p ( ψ ) and F ( ψ ).From this point all the calculations are performed indimensionless variables, so that the results can be scaledto any machine size and physical parameters. In dimen-sionless form, the Grad-Shafranov equation becomes ∂ ¯ ψ∂x − x ∂ ¯ ψ∂x + ∂ ¯ ψ∂y = − ϕ (cid:32) ˜ β x d ¯ pd ¯ ψ + ¯ F d ¯ Fd ¯ ψ (cid:33) , (3)where ¯ ψ ( x, y ) = ψ ( R, z ) /ψ is the normalized poloidalmagnetic flux, and ψ is the flux at the magnetic axis ofthe plasma. The variables ( x, y ) = ( R/R , z/R ) are thenormalized cylindrical coordinates, with R the major ra-dius of the plasma measured to the center of the poloidalcross section. The normalized arbitrary functions are¯ p ( ¯ ψ ) = p ( ψ ) /p and ¯ F ( ¯ ψ ) = F ( ψ ) /R B , where p isthe kinetic pressure at the magnetic axis and B is thevacuum toroidal field at R . A characteristic beta was de-fined as ˜ β = 2 µ p /B and the parameter ϕ = R B /ψ characterizes the ratio of toroidal and poloidal magneticfluxes.Setting to zero the poloidal flux at the plasma edge, ¯ ψ grows monotonically towards one at the magnetic axis.If the kinetic pressure is required to have a first orderdependence of ¯ ψ , it has to be simply¯ p ( ¯ ψ ) = ¯ ψ. (4)This guarantees the vanishing of the pressure at theplasma edge, and a maximum value at the magnetic axis.To set the form of the poloidal current F ( ψ ), notice thatthe toroidal magnetic field has the form¯ B φ = ¯ F ( ¯ ψ ) x . (5)This field must tend to its vacuum form B vφ = 1 /x at theplasma edge where the plasma density vanishes. Then,the poloidal current must satisfy ¯ F (0) = 1. Requiringa second order dependence of ¯ F on ¯ ψ leads the generalform ¯ F ( ¯ ψ ) = a ¯ ψ + 2 b ¯ ψ + 1 (6) where the parameters a and b must be related to theequilibrium parameters. Using these arbitrary functionsthe Grad-Shafranov (3) equation takes the linear form ∂ ¯ ψ∂x − x ∂ ¯ ψ∂x + ∂ ¯ ψ∂y = − ϕ (cid:32) ˜ β x + a ¯ ψ + b (cid:33) . (7)The poloidal flux is a superposition of an homogeneoussolution ψ h and a particular solution ψ p . Assuming that ψ p depends only on even powers x , the particular solutionmay take the form ψ p = − ˜ β a x − ba . (8)Other choices lead to infinite series expansions that un-necessarily complicate the analysis or provide solutionsthat appear in the homogeneous solution.To solve the homogeneous equation, define s = aϕ for a > s = − aϕ for a <
0, so that s is always areal number. Then we write the homogeneous equationas an eigenvalue problem ∂ ψ h ∂x − x ∂ψ h ∂x + ∂ ψ h ∂y = ∓ s ψ h . (9)This can be solved by separation of variables with sep-aration constants α, γ related to the eigenvalue through ± α ± γ = ± s .For a > ψ h ( x, y ) = x [ I , K ( αx )][sin , cos( γy )] , γ = α + s [1 , x ][sin , cos( sy )] , α = 0 x [ J , Y ( αx )][sin , cos( γy )] , γ = s − α x [ J , Y ( sx )][1 , y ] , γ = 0 x [ J , Y ( αx )][sinh , cosh( γy )] , γ = α − s (10)Where J , Y , K , I are the Bessel and Bessel modi-fied functions of first order and we have used the ab-breviated notation [ f , f ( x )][ g , g ( y )] = c f ( x ) g ( y ) + c f ( x ) g ( y ) + c f ( x ) g ( y ) + c f ( x ) g ( y ), with c i arbi-trary constants. For the cases with α = 0 or γ = 0 weuse the dominant terms of the Bessel and harmonic func-tions for small arguments. This gives the same solutionsthat solving again the PDE (9) with a single separationconstant.Another possible solution can be obtained withoutseparating variables and assuming spherical symmetry ψ h ( x, y ) = ψ h ( r ) with r = (cid:112) x + y . In this case theeigenvalue problem (9) becomes d dr ψ h ( r ) = − s ψ h ( r ) , (11)with solutions ψ h = sin( sr ) , cos( sr ) . (12)This solution is relevant for the modern small aspect-ratio tokamaks and spheromaks, where the conductingchamber is D-shaped and the magnetic surfaces near theplasma edge are deformed accordingly.In analogy, for a < ψ h ( x, y ) = x [ J , Y ( αx )][sinh , cosh( γy )] , γ = α + s [1 , x ][sinh , cosh( sy )] , α = 0 x [ I , K ( αx )][sinh , cosh( γy )] , γ = s − α x [ I , K ( sx )][1 , y ] , γ = 0 x [ I , K ( αx )][sin , cos( γy )] , γ = α − s (13)and the spherical solutions ψ h = sinh( sr ) , cosh( sr ) . (14)In Fig. 1 we arrange the solutions of (9) in the parameterspace α − γ , this illustrates the relation between the formsof the functions and the possible values of α, γ . FIG. 1. Different forms of the solutions to the eigenvalueproblem respect to the parameter space α − γ . The upperfunctions on each box correspond to the solutions for a > a < In general, to solve a boundary value problem we willexpress the poloidal flux as¯ ψ ( x, y ) = − ba − ˜ β a x + c G ( sr ) + (cid:88) α,γ c α,γ B ( αx ) H ( γy ) , (15)where the form of the spherical solution G and the func-tions B and H depend on the sign of a (Fig. 1). Thevalues of the parameters α , γ and s must be adjustedto satisfy the boundary conditions and the number ofelements in the sum is, in principle, arbitrary. The su-perposition can also be expressed as an integral, but fromthe numerical point of view we only work with discretevalues of α and γ . III. EQUILIBRIUM PARAMETERS
Before dealing with the numerical method to solve theboundary value problem we need to establish relationsbetween physical parameters and the parameters of theanalytical solution (15).From the toroidal magnetic field (5), the requirementfor a diamagnetic plasma is ¯ F ( ¯ ψ ) (cid:46) ≤ ¯ ψ ≤
1. To keep track of this condition wedefine the constant f = ¯ F (1) and use it instead b inthe poloidal flux expansion. Using (6) we obtain 2 b = f − − a , and the squared poloidal current becomes¯ F ( ¯ ψ ) = a ¯ ψ ( ¯ ψ − − (1 − f ) ¯ ψ + 1 . (16)This form is more convenient to define explicitly the dia-magnetic reduction of the toroidal magnetic field insidethe plasma or its increase in paramagnetic cases. Using(16), the toroidal current density in units of B /µ R becomes ¯ j φ = ϕ x [ ˜ βx + a (2 ¯ ψ −
1) + f − . (17)Close to the magnetic axis, the safety factor can be ap-proximated by q ( ρ ) = ρ ¯ B φ (1 + ∆) ¯ B p ( ρ ) . (18)Here, ρ, ∆ are the minor radius of the toroidal magneticsurface and the Shafranov shift in units of R , and ¯ B p , ¯ B φ are the poloidal and toroidal components of the magneticfield in units of B . As ρ →
0, (18) gives the exact valueof the safety factor at the magnetic axis q . For ρ smallthe poloidal magnetic field can be approximated by¯ B p ( ρ ) = ¯ j ρ, (19)where ¯ j is the dimensionless toroidal current density atthe magnetic axis. Setting ¯ ψ = 1 and x = 1 + ∆ in (17)to replace ¯ j in (19), and replacing (16) in (5), the safetyfactor (18) becomes constant and equal to q q = 4 fϕ (1 + ∆)[ ˜ β (1 + ∆) + a + f − . (20)Now, this relation is used to write ϕ in terms of the otherparameters, and is replaced in the eigenvalue equation s = (cid:112) | a | ϕ. (21)This gives s in terms of q , ∆, f and the adjustable pa-rameter as = 4 f (cid:112) | a | q (1 + ∆)[ ˜ β (1 + ∆) + a + f − . (22)For a predictive calculation we can set the values of theparameters f, q , ∆ and ˜ β . For instance, in a usual dia-magnetic configuration f (cid:46) q (cid:38) ≈ − . Toestimate the characteristic beta ˜ β , we use (4), and thedefinition of the toroidal beta β t = 2 µ (cid:104) p (cid:105) /B , leadingto β t = ˜ β (cid:104) ¯ ψ (cid:105) , (23)where (cid:104)(cid:105) denotes a volume average in the plasma domain.In a high-beta plasma the toroidal beta dominates thevalue of the total beta, β ≈ β t [12]. Also, since ¯ ψ = 0 atthe plasma edge and ¯ ψ = 1 on the magnetic axis, we canexpect (cid:104) ¯ ψ (cid:105) ≈ .
5. Then, for a given value of beta, we canapproximate the characteristic beta by˜ β ≈ β (24)To set the eigenvalue of the problem we need to know a in (22). We can define a as an adjustable parameter,but its range of allowed values must be established in aphysical basis. To do this, we require the toroidal currentdensity not to change its sign because it is mainly createdby an inductive electric field. The signs of j φ and ψ mustbe the same, consequently the signs of ¯ j φ and ϕ in (17)are the same, leading to the condition˜ βx + a (2 ¯ ψ −
1) + f − ≥ , (25)in the whole plasma domain. This leads to1 − f − ˜ β (1 + ∆) ≤ a ≤ ˜ β (1 − (cid:15) ) + f − , (26)where (cid:15) is the aspect ratio of the plasma. Now, the con-dition for no poloidal current density inversions, comesfrom requiring the poloidal current ¯ F ( ¯ ψ ) to be a mono-tonic decreasing function of ¯ ψ , i.e. ¯ F (cid:48) <
0. Using (16)we obtain the condition f − < a < − f . (27)This conditions is useful to identify the parameters forpoloidal current density inversions that may be requiredto describe the reversed magnetic shear equilibria emerg-ing in situations with large bootstrap fractions. Follow-ing the conditions (27) and (26) we can identify the re-gions of interest in the parameter space a − ˜ β (Fig. 2). FIG. 2. The light region correspond to pairs ( a, ˜ β ) for whichthere are no toroidal or poloidal current density inversions. Inthe dark region the poloidal current density inverts but thetoroidal do not. In Fig. 2 we depict the regions for poloidal currentinversion ¯ j p , and no-inversions using the definitions β = 2 1 − f (1 + ∆) , β = 2 1 − f (1 − (cid:15) ) (28)and 1 /β = 1 /β + 1 /β . The restrictions over the al-lowed values of a and ˜ β defines through (22) the set ofallowed eigenvalues s that gives physical solutions to theboundary value problem. Given the form of the solutions(15) it is more convenient to adjust the eigenvalue s thanthe parameter a . For this, we invert (22) to write a interms of s and the physical parameters q , ∆ , f, ˜ β . a ± ( s ) = 1 − f − ˜ β (1 + ∆) ± f q s (1 + ∆) × − (cid:115) ± q s (1 + ∆) f [1 − f − ˜ β (1 + ∆) ] , (29)where a + ( s ) is valid for a > a − ( s ) for a < FIG. 3. Like in Fig. 2, the light region correspond to pairs( s, ˜ β ) for which there are no current density inversions andthe darker region lead to poloidal current density inversions.As negative values of s are not allowed the regions for a < a > In Fig. 3, we depict the regions of interest in the pa-rameter space s − ˜ β using the definition A = 2 q (1 + ∆) f (cid:112) − f . (30)As in Fig. 2, the region with a < a >
0. Also, large values of s lead topoloidal current inversions, except for a narrow region of˜ β between β and β . IV. NUMERICAL METHOD
Now that we have characterized the equilibrium solu-tions respect to their position in the parameters space,we can develop a systematic method to build analyticalsolutions with some desired equilibrium properties.Using 2 b = f − − a and (29) the poloidal magneticflux (15) can be casted like¯ ψ ( x, y ) = 12 + 1 − f − ˜ βx a ± ( s ) + c G ( sr )+ (cid:88) i c i B ( α i x ) H ( γ i y ) , (31)with γ i = γ ( α i , s ) (see Fig. 1). The sign of a deter-mines its form in (31) and the functions G, B and H asexplained in the Section II. The sign of a is then kept un-changed during any optimization procedure that modifiesthe eigenvalue s , the coefficients c i and the parameters α i .The Levenberg–Marquardt algorithm [13] is used toadjust the linear and nonlinear parameters involved inthis problem. This method gives good convergence forreasonable choices of the starting parameters. In general,the iterative process consists in the minimization of theerror functional ε ( (cid:126)k ) = N (cid:88) i =1 [ ψ i − ¯ ψ ( p i , (cid:126)k )] , (32)where p i = ( x i , y i ) are points where we know the nu-merical values of the poloidal flux ψ i , and ¯ ψ ( p i , (cid:126)k ) is ourapproximation to that value through (31) for a given setof M parameters (cid:126)k = { s, { c i } , { α i }} . The minimizationof (32) is done by successive variations of (cid:126)k , (cid:126)k = (cid:126)k + (cid:126)δ, (33)where (cid:126)δ must satisfy ε ( (cid:126)k + (cid:126)δ ) < ε ( (cid:126)k ) and is obtained bysolving the linear problem( J T J − λI ) (cid:126)δ ( λ ) = J T [ (cid:126)y − (cid:126)f ( (cid:126)k )] . (34)Here, I is the M × M identity matrix, λ is an adjustableparameter and the vectors are defined by (cid:126)y = ( ψ , ψ , ..., ψ N ) T , (35) (cid:126)f ( (cid:126)k ) = ( ¯ ψ ( p , (cid:126)k ) , ¯ ψ ( p , (cid:126)k ) , ..., ¯ ψ ( p N , (cid:126)k )) T . (36) J is an M × N matrix with entries J i,j = ∂ ¯ ψ ( p j , (cid:126)k ) ∂k i , (37)that in this case can be calculated analytically. To up-date (cid:126)k and λ we calculate the errors for (cid:126)δ ( λ ) and (cid:126)δ ( λ (cid:48) )where λ (cid:48) = rλ and 0 < r <
1. Then (cid:126)k and λ are updatedwith the variation that gives the largest error reduction.If neither reduces the error we do λ → λ/r and repeatthe previous step. Following this procedure we guaran-tee a rapid convergence far from the minimum and morerefined steps close to it.The points where the poloidal flux is known are onthe plasma edge, where ψ i = 0 and the magnetic axiswhere ψ N = 1. To describe the plasma edge we can usea parametric equation containing the relevant geometry x b ( θ ) = 1 + (cid:15) cos( θ + α sin θ ) , (38) y b ( θ ) = (cid:15)κ sin θ, (39) δ = sin α. (40)This describes a D-shape with triangularity δ , elongation κ and minor radius (cid:15) . In the case of a single or double null configuration we can trace straight lines that meetat the X-point at a distance η(cid:15) from the center with adesired angle ξ (see Fig. 4). For given values of ξ and η ,the positions of the X-point p , and the tangency points p , p are uniquely determined and can be found by solv-ing numerically an implicit equation. FIG. 4. The plasma edge is modeled by merging a D-shapewith two tangent straight lines starting at p , p and crossingat the X-point p . The center of the column is at x = 1 andthe magnetic axis is displaced by ∆. V. RESULTS AND DISCUSSION
In the following, the numerical optimization describedin the section IV will be used to find possible equilibriumconfigurations with realistic features in cases with normaland reversed magnetic shear.
A. Normal shear equilibrium
As a first example, the optimization algorithm is usedto describe a shaped plasma with the parameters inTable I. As we can not preset the value of β in our TABLE I. Desired parameters (cid:15) κ δ ξ η ∆ q β . . . . π . . . method, we will set the value of ˜ β following (24), as-suming (cid:104) ¯ ψ (cid:105) ≈ . β ≈ . f ≈ .
95. There is aclose relation between the toroidal field fraction f andthe value of ˜ β , so we perform several runs for differentcombinations of ˜ β and f until we find the combination( ˜ β, f ) = (0 . , . ψ .To choose the basis we first set a > ψ ( x, y ) = 1 / − f − ˜ βx ) / a ( s ) + ( c + c y ) xY ( sx )+ c sin( sy ) + c cos( sy ) + c sin( sr ) + c cos( sr )+ c xJ ( α x ) sin( γ y ) + c xJ ( α x ) cos( γ y )+ c xY ( α x ) sin( γ y ) + c xY ( α x ) cos( γ y )+ c xK ( α x ) sin( γ y ) + c xK ( α x ) cos( γ y )+ c xJ ( α x ) cos( γ y ) + c xY ( α x ) cos( γ y ) . (41)with r = (cid:112) x + y . Using the restrictions (26),(27) over a , we were able to estimate the starting eigenvalue on s ≈ .
0, and the initial values of the starting parame-ters were chosen to be α − = 0 . s , α , = 3 . s and α , = 0 . s . These values evolve independently of s during the optimization process, then, they will spreadin the parameter space α − γ . The initial values of theexpansion coefficients { c i } are calculated by solving thelinear problem of minimizing the error ε ( { c i } ) for s, { α i } fixed on the starting values.In Fig. 5 we can see the evolution of the solution pa-rameters as the error is reduced from 26 . . × − in200 iterations of the method, when the parameters andthe error do not change significantly the run ends.After the minimum is reached and we are satisfied withthe plasma shape we can calculate the relevant plasmaprofiles. In Fig. 6-left we depict the resulting topology ofthe magnetic surfaces and the 50 control points used inthe method, ψ − = 0 for the plasma edge and ψ = 1at the magnetic axis. The plasma edge is in good agree-ment with the desired shape and the magnetic surfacesbehave as expected with the magnetic axis slightly dis-placed from the desired position.We use (5,16,29) to calculate ¯ B φ and compare with thevacuum toroidal field ¯ B vφ = 1 /x . In Fig. 6-right we cansee the reduction of the toroidal field due to the diamag-netic effect controlled by f . The pressure profile is bydefinition the same of ¯ ψ and the toroidal current density¯ j φ is calculated with (17) using the relations (21,22,29).For ˜ β = 0 .
26 the obtained eigenvalue is slightly outsidethe shaded region in Fig. 3, accordingly, there is a moder-ated inversion in j φ coming from the inherit restrictions − − − − − − − − FIG. 5. Traces of the expansion coefficients { c i } and thenonlinear parameters s, { α i } for 200 iterations of the mini-mization method. The eigenvalue s stabilizes at 2 .
69 and theerror at 2 . × − . of the model, and, in addition, ¯ j φ does not vanish at theedge in the low field side.To understand the j φ profile, notice in (17), that thetoroidal current density is not constant at the plasmaedge, where ¯ ψ = 0, but changes with the radial distance x . Consequently, the current density can not vanish inthe low and high field side simultaneously. This is a di-rect consequence of the choice of the profiles of p ( ψ ) and F ( ψ ) that makes linear the Grad-Shafranov equation (3),and such unphysical behavior is acceptable when workingwith analytical models.The safety factor q ( ¯ ψ ) is the constant ratio dφ/dϑ ofthe toroidal and poloidal angles subtended by the mag-netic line as it wanders over its invariant surface. Thepoloidal angle ϑ , is not uniform in the Cartesian space { x, y } , but the toroidal one is, so, we can calculate q ( ¯ ψ )by following the magnetic line until it completes a fullpoloidal cycle, then we use q ( ¯ ψ ) = ∆ φ π . (42)Doing this for a set of initial conditions in the line y = 0 we get the q -profile in Fig. 6-right. The value q min = 1 .
096 is in close agreement with the desired q = 1 .
1, and the safety factor at the 95% flux surface(in our case ¯ ψ = 0 . q = 4 . (cid:104) ψ (cid:105) = 0 .
45, then,from (23) the toroidal beta is β t = 11 . β p = 2 µ (cid:104) p (cid:105) / ¯ B p where . . . . . . . − . − . . . . . . . . . . . . . FIG. 6. In the left, the poloidal magnetic flux contours showgood agreement with the imposed control points ( × ), at theedge and magnetic axis. In the right, the resulting profiles of q ( x ) , B φ ( x ) , ¯ ψ ( x ) , ¯ p ( x ) for a transversal cut y = 0, behave asexpected for a divertor discharge, and the vacuum magneticfield ¯ B vφ = 1 /x stands above B φ , illustrating the diamagneticeffect of the poloidal plasma current. ¯ B p = µ I / πaκ , a is the minor radius, κ the elongationand I the plasma current. In dimensionless variables wecan write this like β p = ˜ β (2 π(cid:15)κ ) (cid:104) ¯ ψ (cid:105) ¯ I , (43)where ¯ I = (cid:82) (cid:82) ¯ j φ ( x, y ) dxdy is the plasma current inunits of R B /µ , namely, the current used to create thevacuum toroidal field. Replacing our values we obtain β p = 2 .
86, then we can calculate the total beta using β − = β − t + β − p , and we obtain β = 11 .
2% that is just1 .
2% above the desired value.
TABLE II. Obtained parameters∆ ˜ β f q min q β t β p β ¯ I .
14 0 .
26 0 .
97 1 .
096 4 .
515 11 .
7% 2 .
86 11 .
2% 1 . Table II summarizes the results presented for this equi-librium. For the chosen poloidal flux expansion (41), theoptimization method led to good agreement with the ex-pected values of Table I, and reasonable prediction forthe values of the plasma current ¯ I p , the poloidal beta β p ,total beta β , and the 95% flux surface safety factor. B. Reversed shear equilibrium
In plasmas with high bootstrap fraction the currentdensity profile is fundamentally changed, presenting acentral minimum, and maximum off-axis. This behav-ior, leads to a non-monotonic safety factor profile withmaximum at the magnetic axis and minimum off-axis. Indivertor discharges the minimum in q is reinforced by the growth of q to the plasma edge, where it diverges. Forthis case we choose a double null equilibrium with theparameters in the Table III. In analogy to the previous TABLE III. Desired parameters (cid:15) κ δ ξ η ∆ q β .
37 1 . . . π .
77 0 .
06 4 . case, we start with the guess ˜ β = 0 .
14 and f = 0 .
95, andperform several runs changing these values for a givenchoice of the poloidal flux expansion. The central re-versed magnetic shear was obtained for f (cid:38)
1, makingthe plasma slightly paramagnetic. Values below one ledto non-monotonic safety factor profiles with several criti-cal points. The parameters that best minimized the errorfunctional were ( ˜ β, f ) = (0 . , . a <
0, correspond-ing to the solution (13) of the Grad-Shafranov equa-tion. Proceeding analogously to the previous case andusing only even functions of y , we consider the followingpoloidal flux expansion¯ ψ ( x, y ) = 1 / − f − ˜ βx ) / a ( s ) + c cosh( sr )+( c + c x ) cosh( sy ) + c xI ( sx ) + c xK ( sx )+ c I x ( α x ) cos( γ y ) + c xK ( α x ) cos( γ y )+ c xJ ( α x ) cosh( γ y ) + c xY ( α x ) cosh( γ y )+ c xI ( α x ) cosh( γ y ) + c xK ( α x ) cosh( γ y )+ c xJ ( α x ) cosh( γ y ) + c xY ( α x ) cosh( γ y ) , (44)the staring eigenvalue was established about s ≈ .
0, andthe parameters were α , = 1 . s , α , = 2 . s , α , =0 . s and α , = 0 .
8. In this case, the optimizationmethod performed 300 cycles and the error stabilized at3 . × − . The eigenvalue stabilized close to the startingvalue at s = 8 . . The resulting profiles are presentedin Fig. 7.For this equilibrium (Fig. 7), the poloidal flux and ki-netic pressure present a stronger drop at the plasma edgeand the toroidal magnetic field is slightly increased (3%)respect to its vacuum value, indicating a paramagneticbehavior. The obtained current density presents a largehole at the plasma center and by inherit model restric-tions it can not develop the off-axis maxima nor decreaseto the plasma edge.The safety factor profile presents the expected maxi-mum at the magnetic axis with q max = 3 .
9, and developsan off-axis minimum near the plasma edge q min = 2 . q close to the separatrix. Thisis a consequence of the inability of the current densityto decrease to the edge. Accordingly, this model is onlyable to represent a global reversed magnetic shear, andshould only be used locally to describe the central regionof the plasma. FIG. 7. In the left, the poloidal magnetic flux contours for adouble null configuration and the control points ( × ), at theedge and magnetic axis. In the right, the resulting profilesof q ( x ) , B φ ( x ) , ¯ ψ ( x ) , ¯ p ( x ) for a transversal cut y = 0, in thiscase we have a hollow current profile and reversed magneticshear. The vacuum magnetic field ¯ B vφ = 1 /x stands below B φ , revealing a paramagnetic behavior of the plasma. The magnetic axis is displaced by 1 .
3% from the de-sired position and the toroidal and poloidal beta are β t = 8 . , β p = 0 .
64, leading to β = 7 .
22% that is0 .
22% above the desired value. The plasma current inunits of R B /µ is ¯ I = 1 . R = 1 . m and B = 2 T , I p = 3 . M A , and is causedby the unavoidable growth of the current density to theplasma edge.
TABLE IV. Obtained parameters∆ ˜ β f q max q min β t β p β ¯ I .
074 0 .
14 1 .
03 3 . .
16 8 .
15% 0 .
64 7 .
22% 1 . For this equilibrium, the values of q max , q min , β t , β p and β are consistent with realistic situations, but themodel is only able to reproduce the internal plasma be- havior, and more flexible profiles for q and j φ requireshigher powers of ¯ ψ on the source functions ¯ p ( ¯ ψ ) and¯ F ( ¯ ψ ), and consequently, the solution of the full non-linearGrad-Shafranov equation. VI. CONCLUSIONS
In this work we have identified the relevant parametersof a linear model of the Grad-Shafranov equation and re-lated them with the plasma parameters of the hydromag-netic equilibrium, revealing the consistency regions in theparameter space. For predictive calculations, the modelparameters can be fixed and the poloidal magnetic fluxbecomes a linear superposition of the solutions to the lin-ear Grad-Shafranov equation. This introduces a numberof free parameters that can be adjusted numerically.We have applied an optimization method to adjust thefree parameters in a situation where the plasma edge andmagnetic axis were established. This method is able toproduce single or double-null equilibrium configurationswith realistic geometry and parameters. It achieves goodmagnetic topology, safety factor profiles, realistic valuesof β , and allows the explicit control the amount of dia-magnetism or even paramagnetism in the plasma. How-ever, it has some inherit limitations in the current densityprofile due to the choice of the arbitrary functions.The presented description of the plasma can be usedglobally to study the equilibrium in usual and high- β dis-charges, and locally to describe the internal plasma in areversed magnetic shear configuration. In both cases theconvergence to the solution is good for reasonable choicesof the basis functions in the poloidal flux expansion.The authors wish to thank Professors R.M.O. Galv˜aoand Z.O. Guimar˜aes-Filho for their useful discussions.This research was depeloped with the financial supportof the Brazilian scientific agencies: CAPES, CNPq andthe S˜ao Paulo Research Foundation (FAPESP), grants2012/18073-1 and 2011/19269-11. [1] V. Hain, R. L¨ust, and A. Schl¨uter, Z. Naturforschg ,833 (1957)[2] Grad and Rubin, Proc. 2nd UN Conf. on the PeacefulUses of Atomic Energy , 190 (1958)[3] V. Shafranov, Sov. Phys. JETP , 545 (1958)[4] V. Shafranov and L. Zakharov, Nucl. Fusion , 599(1972)[5] L. Solov’ev, Sov. Phys. JETP , 400 (1968)[6] E. K. Maschke, Plasma Phys. , 535 (1973)[7] S. Zheng, A. Wootton, and E. Solano, Phys. Plasmas ,1176 (1996)[8] G. Ludwig, Plasma Phys. Control. Fusion , 2021 (1997) [9] P. McCarthy, Phys. Plasmas , 3554 (1999)[10] L. Guazzotto and J. Freidberg, Phys. Plasmas , 112508(2007)[11] A. J. Cerfon and J. P. Freidberg, Phys. Plasmas ,032502 (2010)[12] J. Freidberg, Ideal Magnetohydrodynamics (PleniumPress, NY, 1987) p. 71[13] D. Marquardt, J. Soc. Indust. Appl. Math. , 431 (1963)[14] E. Strait, L. Lao, M. Mauel, B. Rice, T. Taylor, K. Bur-rell, M. Chu, E. Lazarus, T. Osborne, S. Thompson, andA. Turnbull, Phys. Rev. Lett.75