ELSEVIER
Journal of Nuclear Materials 220-222 (1995) 1024-1027
Techniques and results of tokamak-edge simulation G.R. Smith a, P.N. Brown a, R.B. Campbell b, D.A. Knoll c, P.R. McHugh c, M.E. Rensink a, T.D. Rognlien a ;t Lawrence Livermore National Laboratory, University of California, Livermore, CA 94551, USA b Sandia National Laboratories - New Mexico, Albuquerque, N M 87185, USA c Idaho National Engineering Laboratory, Idaho Falls, ID 83415, USA
Abstract This paper describes recent development of the U E D G E code in three important areas. (1) Non-orthogonal grids allow accurate treatment of experimental geometries in which divertor plates intersect flux surfaces at oblique angles. (2) Radiating impurities are included by means of one or more continuity equations that describe transport and sources and sinks due to ionization and recombination processes. (3) Advanced iterative methods that reduce storage and execution time allow us to find fully converged solutions of larger problems (i.e., finer grids). Sample calculations are presented to illustrate these developments.
1. Introduction Fluid models of the tokamak edge plasma are developing rapidly through the efforts of several groups and are reviewed in Ref. [1]. The U E D G E simulation code is used to find equilibrium and time-dependent solutions of a two-dimensional fluid model of the tokamak edge as described in Ref. [2]. All of the 2-D models solve in the poloidal plane and assume classical parallel (and thus poloidal) transport and anomalous radial transport. The solution algorithms do differ, and the U E D G E code utilizes a fully implicit method with a numerical Jacobian. The need to interpret experiments and design new divertors has led to ever more realistic models. This paper describes recent improvements made to U E D G E and presents some new results.
other set locally parallel to the plates or material surfaces. The numerical task of forming the finitevolume equations to be solved requires finite-difference representation of the fluxes passing through the faces of the grid-cell volumes. A 9-point difference stencil is used. Further details of the differencing procedure are given in Ref. [3]. In Fig. 1 we show an example of a non-orthogonal grid used for calculations described in Section 3. The grid is distorted to accommodate protruding shelves that are intended to direct recycling neutrals back toward the separatrix. In the design, a reentrant or baffled region also exists under the outer shelf to facilitate the pumping of gas. Our grid does not include the reentrant region where neutrals escape to the pump. 2.2. Impurities
2. Improvements in the UEDGE simulation code 2.1. Non-orthogonal grid
Divertor plates and other material surfaces that interact with the edge plasma are generally not orthogonal to the magnetic flux surfaces. To treat such situations effectively, we use a non-orthogonal grid that has one set of grid lines parallel to flux surfaces and the
The average-ion and multi-charge impurity models developed previously and described in Refs. [4,5] have been incorporated into UEDGE. Continuity equations are added to describe the spatial variation of the impurities; one equation is added for the average-ion model, and an equation for each charge state is added for the multi charge impurity model. Ionization and recombination are included. The velocities normal to
0022-3115/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0022-3115(94)00466-8
G.R. Smith et aL /Journal of Nuclear Materials 220-222 (1995) 1024-1027 the flux surfaces are described by the same anomalous diffusion process as for the hydrogenic species. In the poloidal direction, the impurity velocity U~z is obtained from parallel momentum balance neglecting inertia and viscosity,
( Uxz
-
-
10
°
,
i
,
,
I
~
1025
i
,
i
I
~
I
,
i
,
~
I
i
10-2 -~
10-4
'~ 10.s
Uxh )miniV hz
n-
10-a
= I-if]
[ - - - Ox + aznz-'~x + ~ z n z ~ x
-Zenz
10-1°
,
(1)
where the various terms are hydrogenic velocity Uxh, interspecies friction rate Vh~, ratio of poloidal to total magnetic field B i B , impurity pressure P~, electron and ion temperatures Te,i, and electrostatic potential q~. The a and/3 are numerical coefficients [5]. The average-ion and multi-charge impurity models are available as alternatives to the simple fixed-fraction model used previously [6]. For the new models, the impurity density distribution is calculated along with hydrogenic-ion density and poloidal flow velocity, electron and ion temperature, and gas density, during iterations to an equilibrium solution.
2.3. Advanced iterative solution methods Our fully implicit formulation of the fluid equations gives us access to the strong convergence properties of
1.0
~
i
i
1'.4
1'.6
0.8
N 0.6
~ > 0.4
0.2
110
500
1.'2
'
MAJOR RADIUS (m) Fig. 1. The region near the X-point of a proposed configuration for the DIII-D experiment at General Atomics. Solid
lines show the grid, which is distorted to conform at the poloidal grid boundaries to protruding material surfaces. Dashed lines show the edges of the shelves and the baffled entrances to the pumping ducts; the baffled region is not simulated here. The X's indicate the location of the boundary of the orthogonal grid that is used for comparisons with the non-orthogonal grid in Section 3.
,
~
,
N
,
300 .~
200
::= 100 O
0 0
2 4 6 8 Index of nonlinear iterations
10
Fig. 2. Performance during Newton-Krylov iterations of three solution methods for the edge-plasma equations. See the text in Section 2.3 for a description of the methods. (a) Reduction of the norm of the residuals as a function of iteration number. (b) Number of evaluations of the residuals, summed from iteration 1 to iteration i, as a function of iteration number i. Newton's method for the entire finite-differenced equation set. For fine grids the associated Jacobian matrix can require significant computational resources (large storage and long inversion times) [2]. NewtonKrylov methods allow reduction in the storage and time requirements, because an approximate factorization of the Jacobian matrix is sufficient and the sparsity of the Jacobian and its factors can be exploited [7]. Furthermore, the approximation involved in the factorization process can be improved considerably by reordering the rows and columns of the Jacobian matrix [8]. The Jacobian is recalculated only every fifth iteration for added efficiency. In Fig. 2 we illustrate the strong convergence and the benefits of reordering during Newton-Krylov iteration of a problem described by a total of 5518 equations (5 equations in each cell of a 62 x 18 grid). The problem being modeled here is gas puffing in the DIII-D tokamak as described in Ref. [9]. Three different approximate factorization methods are compared. The first method (curves A in Fig. 2) stores only the matrix elements of the factors L and U that occur in the diagonals that had non-zero elements in the Jacobian. The second method (curves B) does no reordering but allows as many as 100 additional non-zero elements ("fill-in") in each row of L and U compared to the number of non-zero elements in the Jacobian. The third method (curves C) does reordering but im-
1026
G.R. Smith et al. /Journal of Nuclear Materials 220-222 (1995) 1024-1027
poses a tighter limit (35) on the fill-in. Part (a) of Fig. 2 shows the performance of Newton-Krylov iteration in reducing the norm of the residuals after an equilibrium solution was perturbed by changing a boundary condition, T¢ at the core-edge boundary, from 200 to 250 eV. Curve A does not show the strong convergence of the other two methods. Curves B and C show similarly good performance in reducing the residuals. However, Fig. 2b shows another aspect of performance, in which reordering (curve C) is superior. That figure shows, for each method, how the total number of evaluations of the residual (right-hand-side) function accumulates during the outer, nonlinear iteration loop of the Newton-Krylov algorithm. The number of residual evaluations is nearly equal to the number of linear iterations in the inner loop and is approximately proportional to the required execution time. Fig. 2b shows that reordering reduces total execution time by a factor of about 2. Reordering reduces residual evaluations, because the exact L and U factors have fewer non-zero elements and L and U can therefore be represented more accurately with a given limit on the fill-in. More accurate L and U lead to more rapid convergence of the linearized problem that is solved for each of the nonlinear iterations. We also take advantage of advanced factorization methods during time-dependent relaxation towards equilibrium, which is of great utility in practice because of its greater radius of convergence than Newton methods.
3. Results with a non-orthogonal grid
The effect of a non-orthogonal grid is illustrated for a case corresponding to the proposed radiative divertor design for DIII-D [10]. The grid for the lower half of this double null configuration is shown in Fig. 1. Grid lines that measure poloidal position are not orthogonal to magnetic flux surfaces but are distorted to conform to the outer protruding shelf and the shaped inner plate surface. The X's show the outline of the orthogonal grid used for a comparison calculation that omits the protruding structures. We consider a case where the electron and ion temperatures (Te and Ti) on the core boundary are 175 eV and the boundary density is 3 × 1019 m 3. The anomalous radial diffusion coefficients are D = 0.25 m2/s, ~(e = 1.5 me/s, and Xi = 1.0 m2/s. The resulting power crossing the separatrix is 8.8 MW, with half going to the lower portion being simulated here. The particle recycling coefficient on the divertor plates is set to 1.0. We use these parameters for two calculations, one with an orthogonal grid and normal divertor plates, and one with the non-orthogonal grid and protruding structures. The non-orthogonal grid does not
lal ~
1.5
~
.0
_~ 0.5
-T
~"
60
~
2O
. . . . .
,
',b-,,
tlJ 1.00
1.04
1.08
Normalized poloidal flux
Fig. 3. Comparison of results computed with orthogonal and non-orthogonal grids. (a) Ion density in m 3 plotted versus normalized poloidal flux, t0., on the outer divertor plate for the case in Fig. 1; solid line, orthogonal grid, and dashed line, non-orthogonal grid with protruding shelf. The separatrix has ~,~ = 1.0. (b) The corresponding electron temperature profiles.
degrade the convergence properties of the numerical algorithm except that about 50% more linear iterations are required for each nonlinear step. One of the goals of the divertor design is to prevent gas from escaping the plate region and fueling the core. While the outer protruding shelf is planned to redirect the recycling gas back into the hot separatrix region where it is ionized, there is the concern that some of the gas may more easily penetrate to the core because the recycling surface is now closer to the core. Our non-orthogonal grid calculations show that the net effect of the protruding shelf is to reduce the plasma and gas density on the outer flux surfaces while increasing these densities near the separatrix. The results are shown in Fig. 3a for the ion density; the gas density profiles are very nearly the same. The profiles are plotted versus the normalized flux qJn, with unity corresponding to the separatrix. The non-orthogonal grid begins to rise from the horizontal for ~, = 1.01, just where the density peaks. The increase in density near the separatrix is caused by the neutrals being directed away from the shelf and toward the separatrix on the angled portion of the shelf. The ionization of these neutrals enhances the plasma density near the separa-
G.R. Smith et aL /Journal of Nuclear Materials 220-222 (1995) 1024-1027
trix. The second effect of the protruding shelf is to increase Te along it, as shown in Fig. 3b. Because the neutrals are directed away from the shelf, the plasma density decreases there, and approximate pressure balance causes the electron temperature to increase.
1027
ated by the carbon; the hydrogenic radiation is only 0.14 M W for this low-recycling case. Note that the largest radiation is on the inner divertor leg, resulting in the lowest Te there. This is qualitatively similar to experimental observations, but this initial illustrative calculation needs be extended to higher recycling and a more detailed carbon source from sputtering.
4. Results with the average-impurity-ion model The average-ion model is illustrated with a calculation corresponding again to the D I I I - D tokamak. The M H D equilibrium is the same as that used for Fig. 1 and corresponds to the gas-puffing case in Ref. [9]. However, to highlight the impurities, we do not gas puff here and reduce the recycling coefficient to the low value of 0.8. The core-boundary values are Te,i = 100 eV and n i = 2 × 1019 m -3. The 200 A source of the carbon impurities is a uniform flux distributed around the outer-wall boundary. The resulting contours of carbon radiation in the divertor region are shown in Fig. 4. The total power crossing the separatrix is 2 MW, and 0.4 M W is radi-
Acknowledgements W e gratefully acknowledge discussions with R.H. Cohen, M.E. Fenstermacher, A.C. Hindmarsh, and G.D. Porter. This work was performed under the auspices of the US D e p a r t m e n t of Energy by L L N L under contract no. W-7405-Eng-48, by Sandia-NM under contract no. D E AC04-76DP00789, and by I N E L under contract no. DE-AC07-76ID01570.
References
0.6
Z ~9
0.4
,,=, >
0.2
__,
1~2
i
1'.4
'
116
MAJOR RADIUS (m)
Fig. 4. Carbon radiation contours for the average-ion model obtained for a DIII-D simulation with carbon injected around the outer wall. Numbers on the contours are multiples of 0.1 MW m -3.
[1] D. Reiter, J. Nuci. Mater. 196-198 (1992) 80. [2] T.D. Rognlien, J.L. Milovich, M.E. Rensink and G.D. Porter, J. Nucl. Mater. 196-198 (1992) 347. [3] G.R. Smith, P.N. Brown, R.B. Campbell et al., Technical Report UCRL-JC-115909, Lawrence Livermore National Laboratory, 1994. [4] D.A. Knoll, R.B. Campbell and P.R. McHugh, Contrib. Plasma Phys. 34 (1994) 386. [5] M. Keilhacker, R. Simonini, A. Taroni and M.L. Watkins, Nucl. Fusion 31 (1991) 535. [6] S.L. Allen, M.E. Rensink, D.N. Hill et al., J. Nucl. Mater. 196-198 (1992) 804. [7] P.N. Brown and Y. Saad, SIAM J. Sci. Stat. Comput. 11 (1990) 450. [8] Y. Saad, Technical Report UMSI 92/38, University of Minnesota Supercomputer Institute, 1992. [9] T.D. Rognlien, P.N. Brown, R.B. Campbell et al., Contrib. Plasma Phys. 34 (1994) 362. [10] S.L. Allen, these Proceedings (PSI-11), J. Nucl. Mater. 220-222 (1995) 336.