Lower Hybrid Simulation Code Manual

D. W. Ignat,
S. C. Jardin, D. C. McCune and E. J. Valeo
Plasma Physics Laboratory
Princeton University
 
November 2000

Abstract

To read the mathematics in this version, which was produced by TtH, under Unix and Netscape, put the following line in .Xdefaults or .Xresources:
 
Netscape*documentFonts.charset*adobe-fontspecific: iso-8859-1

Nov 19, 2000

Contents

1  Introduction

Lower hybrid heating has been of interest for some time. Computational analysis perhaps started with an approach to calculating the spectrum of waves launched from a coupler [1]. The eikonal method [2], or geometrical optics [3], also referred to as ray tracing, was known to be available to problems in plasma physics, and several workers started investigating the consequenses of the ray method [4,5,6,7,8]. Success at driving the toroidal current in a tokamak [9,10] following theoretical predictions [11] caused increased interest in theory [12,13], modelling [14,15,16,17] and related considerations [18,19,20,21,22]. A natural extension was to add to the comprehensive analysis tools such as TRANSP [23,24] and TSC [25] the the best models of lower hybrid heating (LHH) and lower hybrid current drive (LHCD). Such an extension was particularly desirable becasue of a possible investment in a major new project to operate a tokamak in a steady state mode [26,27].

A static but magnetically accurate model [28,29] plus studies of wave chaos and accessibility [18,22,30,31] revealed many of the important issues.

The Lower hybrid Simulation Code (LSC) is a computational model of lower hybrid current drive in the presence of an electric field [32] derived from direct forerunners in references [8] and [15]. Heating of electrons and ions is also computed. Details of geometry and the plasma profiles are treated, so that the accessibility problems are treated as fully as possible in the ray model [32,21]. Two-dimensional velocity space effects are approximated in a one-dimensional Fokker-Planck treatment [12]. LSC contains a model of an xray camera to help compare results of the model with the physical camera [33,34,35,36].

Applications include work in references [32], [37] and [38].

The LSC was originally written to be a module for lower hybrid current drive and heating called by the Tokamak Simulation Code (TSC) , which is a numerical model of an axisymmetric tokamak plasma and the associated control systems [25,39,40]. The TSC simulates the time evolution of a free boundary plasma by solving the MHD equations on a rectangular computational grid. The MHD equations are coupled to the external circuits (representing poloidal field coils) through the boundary conditions. The code includes provisions for modeling the control system, external heating, and fusion heating.

The LSC module can also be called by the TRANSP [23,24] code. TRANSP represents the plasma with a axisymmetric, fixed-boundary model and focuses on calculation of plasma transport to determine transport coefficients from data on power inputs and parameters reached.

This manual covers the basic material needed to use the LSC. If run in conjunction with TSC, the ``TSC Users Manual'' should be consulted [39,40]. If run in conjunction with TRANSP, on-line documentation will be helpful.

The few general mathematical subroutines (root finder, matrix invertor, random number generator) were taken from the the standard ``Numerical Recepies'' textbook [41]. X ray cross sections were entered from the literature [42,43]. Therefore no external libraries are needed for LSC.

A theoretical background of the governing equations and numerical methods is given in section 2. Information on obtaining, compiling and running the code is given in section 3.

This manual can be found at w3.pppl.gov through one or all of ~ignat, NTCC or topdac.

At such locations there will be links to the source code .

The manual is available in html written by Ian Hutchinson's program TtH and in html written by latex2html. For viewing TtH output with netscape under Unix you will need to put the following line in either .Xdefaults or .Xresources:

Neither html version is meant for printing.

The manual is also available in PostScript as a .ps.gz file and in the Portable Document Format pdf with activated internal links as a .pdf file.

2  Theoretical Background

2.1  Geometrical Configuration

For LSC and TSC the poloidal flux is minimum in the center of the plasma and rises toward the plasma boundary. For LSC the B_Z at the outer midplane is positive for a positive current.

2.2  Wave Propagation

2.2.1  Dispersion relation and energy absorption

E = Sj Ej(r ) expi Fj  ,

(1)

where the rapid space and time variation occurs through the exponential

Fj º ó õ kj(r ) ·dr  -  wt  ,

(2)

which depends on the local wave vector k_j and the frequency w. Additional time behavior is ignored because the transit time of a wave is short compared to the time scale of plasma evolution. The mode amplitudes E_j vary much less rapidly. Suppressing the index j, each mode locally satisfies the matrix equation, [44]

[ kk - I k2 + k02 K(r,k, w) ]·E = 0

(3)

where k₀ = w/c represents the free-space wavenumber, and K is the dielectric tensor, which has a nontrivial solution for a vanishing determinant

e(w, k, r) º |[ kk - I k2 + k02 K(r,k, w) ]| = 0     .

(4)

The dispersion relation  (4) connects the vector k and frequency w at the spatial location r (explicit time variation of the dielectric properties are assumed negligible from here on). To proceed it is convenient to choose a local Cartesian coordinate system [44], such that [^(z)] ×B = 0, and k is contained in the x-z plane:

k = k||   ^ z   + k^  ^ x    .

(5)

For lower hybrid waves, the ion cyclotron frequency is much less than the wave frequency which is much less than the electron cyclotron frequency: w_ci << w << w_ce, an inequality assumed throughout. In the Hermitian part of the dielectric tensor we keep only cold plasma terms, except that the dominant warm-plasma term is carried to guard against singularity near the lower hybrid resonance.

The anti-Hermitian part of the tensor is retained as a perturbation. For the case at hand, the principal such term enters as an imaginary correction to K_zz, and describes the interaction between the component of the wave electric field parallel to B and electrons whose speed along B matches that of the wave (Landau damping ). Thus, the plasma dielectric behavior is described by the following tensor elements, all other elements being zero:

Kxx = Kyy = S - ak^2     , (6) -Kyx = Kxy = i D = i wpe2 wwce    , (7) Kzz = P     + i Kzz,i    , (8)

where

S = 1 + wpe2 wce2 - å j  wpi,j2 w2     , (9) a = 3 4 wpe2 wce4 vTe2 + 3 å j  wpi,j2 vTi,j2 w4     , (10) P = 1 - wpe2 w2     , (11) Kzz,i = - p wpe2 w ó õ dv||v|| ¶fe ¶v|| d(w- k||v||)     , (12)

and the new symbols are as follows: w_pe is the electron plasma frequency n_e e²/(e₀ m_e), w_pi,j is the ion plasma frequency of the j^th species, v_Te and v_Ti,j are electron and ion thermal velocities ( º kT/m), and f_e(v_||) is a one-dimensional electron velocity distribution function normalized such that

ó õ dv||fe(v||) = 1  .

(13)

The permittivity of free space is e₀. The thermal term designated by a would be important near lower hybrid resonance , S = 0. However, in cases treated in this paper S does not approach zero, and a is not important. This would be true whenever w² >> w_cew_ci. The wave-particle interaction responsible for electron heating and current drive is in K_zz,i. In the event the lower hybrid resonance should become important, a thermal term representing the ion wave-particle interaction would have to be added to K_xx and K_yy, but here we assume such terms vanish.

In the coordinate system of equation (5), we decompose the dispersion relation  (4) into its real and imaginary parts using equations (6) - (12),

e = er + i ei = 0  ,

(14)

where

er = -ak^6 + k^4 S + k^2   é ë (P + S)(k||2-k02S) +k02D2 ù û +P   é ë (k||2-k02S)2-k04D2 ù û     , (15) ei = ¶er ¶P Kzz,i     . (16)

If n º k c / w >> 1 a simplified dispersion relation can be found from the matrix equation leading up to equation (4) by asymptotically expanding

^ n   ~ ^ n   0   + 1 n2 ^ n   1   + O( 1 n4 )  , E ~ E0 + 1 n2 E1 + O( 1 n4 )  . (17)

To lowest order, I_^·E⁰ = 0, where I_^ = I - [^(n)] [^(n)], so that E⁰ = [^(n)]⁰ E⁰. In first order, there arises the solubility condition

^ n   0   ·K · ^ n   0   = 0  ,

(18)

which gives the electrostatic limit of the dispersion relation :

er = - ak^4 + k^2 S + k||2 P     .

(19)

The consistency condition |E¹| = |K ·E⁰| < n^-2 |E⁰| is satisfied for n_||² > 1 + (w_pe² / w_ce²).

The extension of the local solutions to a spatially inhomogeneous plasma is accomplished by the eikonal method with the useful result that an initial wave field at r with an initial propagation vector k evolves according to Hamiltonian equations which preserve the local dispersion relation e_r = 0 along the ray trajectories

dr dt = - ¶er ¶k / ¶er ¶w (20) dk dt = + ¶er ¶r / ¶er ¶w (21)

As in Hamiltonian mechanics, the spatial coordinates denoted by r are canonically conjugate to the wave number coordinates denoted by k. For example, in a Cartesian frame one would have (x,y,z) and (k_x,k_y,k_z) and fields varying as expi ò(k_x dx +k_y dy + k_z dz); and in a cylindrical frame one would have (R,Z,f) and (k_R,k_Z,l) where R, Z have dimensions of length, k_R, k_Z have dimensions of inverse length (wave number), and l is the dimensionless toroidal mode number. In the study of the axisymmetric tokamak, it is the cylindrical coordinate system that is the most natural since the toroidal wave number l is constant along the ray path ( d l / dt µ ¶e_r / ¶f = 0).

In formulating e_r it is necessary to use the canonical forms by writing k_|| = k ·B/B and k_^² = k² - k_||² where B and B are magnitude and vector quantities of the local magnetic field.

A spectral component of power W experiences a change in power DW over time interval Dt:

DW = -2  ei / ( ¶er ¶w ) W Dt = -2  ¶er ¶P Kzz,i /( ¶er ¶w )  W Dt = 2p wpe2 w   ó õ dv||v||  ¶fe ¶v||  d(w- k||v||)     ¶er ¶P /( ¶er ¶w )   W Dt     . (22)

Note that equations (20), (21), and (22) also carry an index numbering the individual spectral component, which has been suppressed to improve readability. Given the velocity distribution and the profiles of macroscopic plasma parameters, the absorption of a lower hybrid spectrum can be computed. An actual incident wave spectrum is a continuous function of the parallel wave number. Computationally, this continuous spectrum is approximated by assigning the input power to a number of discrete rays as in equation (1), each ray having a definite initial k_|| and launched power. The number of rays and their distribution over k_|| can be chosen to adequately resolve the spectrum in cases of reasonably strong single-pass damping. On the other hand, if many passes are required to absorb the power, a very dense packing of rays into the spectrum might be required to obtain well behaved results.

The plasma is divided up into discrete shells between flux surfaces, typically 20 - 100. For each ray, equations (20) through (22) allow calculation of the needed quantities. The width of the shells and the detailed trajectory determine Dt in equation (22). As a practical matter, one cannot ensure that the shells are fine enough to keep |DW| / W < 1, so this quantity must be monitored and Dt adjusted downward the amount necessary to keep |DW| < W . The new Dt is recorded for use in equation (26) of section 2.3.1.

2.2.2  Changes to k_||

In applications of the model to many devices the following features are significant:

1. toroidal effects move the value of k_|| up and down during propagation - a rise helping absorption, and a fall retarding absorption
2. the rise in k_|| is a stronger effect as density, poloidal field inside the plasma increase, and as aspect ratio (R/a) decreases

These effects are well established in the literature [14,18,19,20,21,32].

An important approximate relationship [18,19,21], exact under the electrostatic form of the cold plasma dispersion relation, between k_||0, the the initial k_|| of a ray when launched, and k_||max, the maximum possible value at a particular point in the plasma is the following.

k||max £ k||0 ×  R0/R 1-(Bp/Bf)( wpe/w) S-1/2     .

(23)

Here R₀ is the initial major radius of the LHCD ray of frequency w, and R, B_p, B_f, w_pe are the major radius, poloidal and toroidal magnetic field, and plasma frequency at some point of interest. The perpendicular dielectric constant is represented by S, a slowly varying quantity of order unity. It has been shown, without benefit of approximations [21], that this maximum tends to be reached in elongated plasmas after many reflections from the wall, whereas the maximum may not be reached in similar circular plasmas. It is worth noting that Kupfer, Moreau, and Litaudon [22] have used physics closely related to the basis of Eq. (23) to formulate a statistical theory of current drive with lower hybrid waves under conditions of multi-pass absorption.

2.3  Fokker-Planck Equation

2.3.1  Quasilinear Diffusion Coefficient

Dql(v||) = p 2 æ ç è e me ö ÷ ø 2    E||2 d(w- k|| v||),

(24)

where E_|| represents the amplitude of the wave field parallel to the local static magnetic field, and the electron charge and mass are e and m_e [44]. One instructive way to find the relationship between field E_|| and wave power W is to equate P_ql, the energy per unit time per unit volume going into electrons and out of the wave from the quasilinear point of view,

Pql = - ne me ó õ dv||v||Dql (v||) ¶fe ¶v||

(25)

with the similar quantity from the ray point of view as assembled from equations (12), (16), (22) to obtain

Dql(v||) = 2 æ ç è p e0 ö ÷ ø æ ç è e me ö ÷ ø 2    W   ¶er /¶P w¶er/¶w   æ ç è Dt DV ö ÷ ø  d(w- k|| v||)

(26)

for the incremental D_ql from a wave of power W traversing a flux shell of volume DV in time Dt.

2.3.2  Kinetic Equation

¶fe ¶t    =   æ ç è ¶fe ¶t ö ÷ ø c   +   æ ç è ¶fe ¶t ö ÷ ø w

(27)

Here the ()_c term is the Coulomb collision operator, and the ()_w term is the wave diffusion (quasilinear) operator. There is no term involving a steady electric field. The wave diffusion operator is the one-dimensional divergence of the rf-induced flux:

æ ç è ¶fe ¶t ö ÷ ø w  = ¶ ¶v|| Dql(v||) ¶fe ¶v||  ,

(28)

where now the term D_ql signifies a sum over all waves in existence on a flux surface, with the appropriate powers and velocities. Note that the use of a simple sum means that we assume there are no interference effects.

We employ a one-dimensional collision operator of the form,

æ ç è ¶fe ¶t ö ÷ ø c  = ¶ ¶v|| é ê ë æ ç è  Dc(v||) ¶ ¶v|| + nc(v||)v||  ö ÷ ø fe(v||) ù ú û     .

(29)

The collisional diffusion and drag coefficients are given by

Dc(v||) = bZ  G/vTe  [1 + (v||/ vTe)2]-3/2     , (30) nc(v||) = bZ  G/vTe3  [1 + (v||/ vTe)2]-3/2     , (31)

where G = lnL  n_e e⁴ / (4pe₀² m_e² ) and the normalization coefficient b_Z is chosen to yield the correct value for the electrical conductivity in the absence of rf. Approximately, b_Z = (1 + Z) / 5, with Z the effective ion charge. The collision operator has several correct properties, such as conservation of particles, producing a Maxwellian (with thermal velocity v_Te) in the absence of wave power, and correct asymptotic velocity dependence of the coefficients for high speeds v >> v_Te.

In solving for f_e we set ¶/¶t = 0 because the time for equilibration between rf power and the electron distribution is short compared to the time for plasma to evolve. Then the solution for f_e is an integral in velocity space.

fe(v||) = 1 Ö 2pvTe2  exp æ ç è - ó õ v|| 0  nc(v¢)  v¢ dv¢ Dc(v¢) + Dql(v¢) ö ÷ ø

(32)

2.3.3  Some Numerical Details

The rays do not move monotonically in ``radius'' (which in our calculation is measured by poloidal flux y ) and in fact typically make many in-out changes of direction. To compute needed quantities at each flux surface, the motion of ray number i_r must be recorded by a zone number i_z, which is sequential with the intersection of that ray with any flux surface. For each zone index, the index i_p of the y surface intersected is kept, where i_p is 1 at the center of the plasma, increasing toward the edge. For each trajectory the i_z starts at 1 and increases to a maximum value; but i_p starts near its maximum, decreases initially, and exhibits individual behavior thereafter.

Changes in direction can be of several types. An inward-going ray can change to an outward-going ray, or vice-versa, because:

1. the roots of the dispersion relation converge (mode conversion);

2. speaking by analogy to a mechanical Hamiltonian system, the radial penetration limit is reached owing to the inpenetrability of a centrifugal barrier;

3. the root of the dispersion relation approaches zero (cutoff);

4. the ray is captured and specularly reflected near the plasma edge. This is accomplished by logic in the calculation which senses that a cutoff is approaching.

The quantities D_ql;iv,ip f_e;iv,ip, W_iz,ir are all made consistent by iteration, beginning with very small W_{iz = 1, ir} and ramping up toward the full power launched. We use an under-relaxation algorithm in the iteration of f_e, in that the past f_e is averaged with the newly computed f_e to get the new function propagated to the next cycle.

2.4  RF-Driven Current

Jrf = -e ne nr ó õ dv||Dql(v||) ¶fe(v||) ¶v|| · ¶Ws(u) ¶u     .

(33)

In the above, n_r = G/|v_r|³, u = v_||/v_r, and v_r = - sign(eE_dc) Ö{m_eG/|eE_dc|}, as explained in the reference [12]. These definitions give as J_rf the current density of stopped electrons, and use the function W_s(u)/u which is given as a tabulation of coefficients fitting simple algebraic terms to solution of the complete Fokker-Planck equation in two velocity dimensions (v_^ and v_||).

The key quantity is W_s(u), the energy (normalized to m_e v_r²/2) imparted to the electric field E_dc by an electron as it slows down. The local electric field E_dc is either prescribed for a static simulation (usually set to zero) or supplied by TSC as part of an iteration to be described.

For very small electric field, W_s(u) can be represented in a power series. From the first two terms, one finds

Jrf = -e ne G ó õ dv||Dql(v||) ¶fe(v||) ¶v||   v||3  4 5 + Z   é ê ë m- 1 + Z/2 + 3m2/2 3+Z   v||2 vr2 ù ú û     ,

(34)

where m = -1 for cooperative current drive, and m = +1 for current drive into an opposing field.

The function W_s(u) behaves in such a way that for u approaching -1, corresponding to rf current drive cooperating with the ohmic current driven by the ambient electric field E_dc at a velocity near the runaway velocity, the calculated current grows extremely large. Because of the runaway condition, this is to be expected. Whenever u » -1, the rf-driven current is not properly calculable from our approach.

It is often necessary to iterate between TSC and LSC for the field and current at the exact time that LSC is called for power and current information, as follows. Suppose at time ^(n-1) the local total current density is J^(n-1), the rf-driven current density is J_rf^(n-1), and the local electric field is E_dc^(n-1). We want to advance quantities to the new time, ⁽ⁿ⁾, and be consistent with Ohm's Law: J = sE_dc + J_rf, where s is the parallel neoclassical conductivity. An initial guess at what the new electric field E_dc^(n-) must be is formed from the new total current, as advanced by TSC, and the old rf-driven current and old electric field.

J(n) = sEdc(n-)+Jrf(n-1)(Edc(n-1))     ,

(35)

From the E_dc^(n-) obtained from equation (35) one can find the estimated new rf current density J_rf^(n-) from equation (33). Its derivative with respect to E_dc,  ( ¶J_rf/ ¶E_dc)^(n-), is formed by numerical differentiation. Algebra yields estimates for the new electric field and for the new rf-driven current density,

Edc(n) = Edc(n-) + Jrf(n-1) - Jrf(n-) s+ ( ¶Jrf/ ¶Edc)(n-)     ,

(36)

Jrf(n) = Jrf(n-) + ( ¶Jrf/ ¶Edc )(n-)( Jrf(n-1) - Jrf(n-) ) s+ ( ¶Jrf/ ¶Edc)(n-)     .

(37)

The new values E_dc⁽ⁿ⁾ and J_rf⁽ⁿ⁾ may not be consistent. This is checked by assigning them to the previous step labeled ^(n-1) and repeating the process until the result is stable to a chosen accuracy. At that point, it is desired that the electric field is small enough that the runaway situation does not exist for electrons interacting with the waves. An inequality expressing this condition, as suggested by equation (20) of reference 12, is

tr-loss . n   er  ne = tr-lossDql(vr) ê ê ê ¶fe(vr) ¶v|| ê ê ê   <<  1

(38)

where t_r-loss is a loss time for electrons pushed by rf waves into a runaway region and [n\dot]_er is rate of increase of the density of runaway electrons owing to the rf diffusion from a lower velocity. A simple estimate for t_r-loss would be the confinement time ( ~ 10 msec). If equation (38) is not satisfied, then the calculation is inappropriate.

Iteration is not necessary if when the ambient electric field is low and the density is high, such that w << |v_r k_|||. In that case the evolution of the electric field and the current is handled without special attention as the transport calculation evolves.

2.5  Broadening of rf-Driven Current

Diffusing a velocity distribution is beyond the scope of the LSC, but we have introduced a heuristic approach similar to that used by V. Fuchs and others [17].

The idea is that the actual, diffused, RF-driven current J_d is modified from the computed RF-driven current J_rf by slowing-down, characterized by an inverse-time n_slow, and a cross-field diffusion, characterized by D_Jrf.

Assuming for the moment slab geometry with x the spatial variable the diffusion-like equation we use is:

¶Jd/ ¶t = nslow ( Jrf- Jd)  +  ¶/¶x  ( DJrf ¶Jd/ ¶x )

(39)

This can be transformed to a normalized flux coordinate [^x] º (y- y_min)/(y_max - y_min) with standard transformations for curvilinear coordinates. Under the reasonably good assumption that y ~ r² the result is:

¶Jd/ ¶t = nslow ( Jrf- Jd)  +  (4/a2)  ¶/¶ ^ x    ( ^ x    DJrf ¶Jd/ ¶ ^ x   )    ,

(40)

where a represents the nominal minor radius of the plasma.

Equation 40 clearly has desirable properties: if D_Jrf is small then the J_d moves to J_rf on the n_slow time scale; and if D_Jrf is large then J_d and J_rf can be quite different, in particular J_d will be smooth compared to J_rf over a length of order Ö{ D_Jrf /  n_slow}.

In LSC D_Jrf on input is one number, constant across the plasma, in units of meter² per second. Internally, D_Jrf is an array, so generalization would be simple. The n_slow are found from equation 31 for a particular phase speed, corresponding to n_|| = 2, also constant across the plasma. The diffusion process is imagined to be complete at each call to LSC, so that the time derivative is set to zero. Under that condition, the solution for J_d comes from one inversion of a tridiagonal matrix.

The boundary condition on J_d at the outer boundary is a zero value. At the inner boundary, the center, the condition of zero flux (zero derivative with respect to radius) translates to the derivative with respect to [^x] being better behaved than 1/Ö{[^x]}, which cannot be expressed in a manner practical for the finite difference equations. Instead, we force dJ_d/d[^x] to be constant over the first two grid spacings. This is equivalent to specifying a zero value for a certain linear combination of J_d and dJ_d/d[^x] at the first interior grid point.

The rf-driven current found J_rf leads directly to the smoothed (diffused) J_d after the inversion of a tri-diagonal matrix [41].

Note that the integral of J_d over the plasma cross section, or total diffused rf-driven current, is not constrained to have any particular relationship with the same integral over J_rf, the total undiffused current. In fact, it seems reasonable that diffusion can in some circumstances increase current by moving fast electrons to regions of lower collision rate. However, it appears that the diffused current is less than the undiffused current for many actual situations, typically having deposition at mid-radius and mildly peaked density profile.

The heuristic smoothing method of V. Fuchs [17] is similar to the method described here, but apparently treats D_Jrf and n_slow as constants, sets the change in total current from an estimate incorporating the formula for n_slow, and finds the value needed for D_Jrf as an eigenvalue.

Our model is not a simulation of energetic electron diffusion in that it does not transport fast particles or change the rf wave damping because of their motion.

2.6  Heuristic Power Diffusion Estimate

The undiffused current J_rf was computed from the ray tracing, quasi-linear development of an electron distribution function in the parallel velocity giving the local power deposition per unit volume P₀, and finally the Karney-Fisch [12] electric field adjustment to current density.

As explained in the previous section, the diffused current density J_d can be quite different from J_rf, and therefore quite different in radial distribution from the P₀ found from ray tracing. This situation is contrasts with the intuition that power density is proportional to current density times background number density. Therefore, we added an option to spread the deposited rf power according to

P( ^ x   ) = as    |Jd( ^ x   )|  ne( ^ x   ) ó õ |Jd( ^ x   )|  ne( ^ x   )  dV( ^ x   )    ó õ P0 dV  + (1-as) P0( ^ x   )     ,

(41)

where P is the diffused power, P₀ is the power deposited from the ray tracing and quasilinear calculation, a_s ranges from 0 (no spreading of P₀) to 1 (full spreading), dV is a volume element.

2.7  X-ray Camera Image

The PBX-M experiment to which we have applied our calculation the most has a 2-dimensional x-ray camera with an axis roughly tangent to a mid-plane circle just inside the major radius of the the plasma. The images found by this camera reveal clues to the behavior of the fast electrons in the plasma. At the same time, such images can be computed from the electron velocity distribution we establish in our model.

Taking account of obstructions, the actual location of pinhole, focal plane, etc., we evaluate the 2-d image with the following formula.

Ix-p = 1 dAp ó õ dl dA(l) dWp(l)  ni(y(l))  ne(y(l)) × ó õ dv||v||fe(v||,y( l)) × ó õ E(v||) 0  dE¢E¢ ¶2 sei ¶E¢¶W [E(v||), E¢, ^ n   · ^ b   ] T(E¢) (42)

In the above, I_x-p is the x-ray intensity on pixel p of area dA_p at the camera focal plane, dl is an element of path length in the plasma on a line passing through the pinhole and the pixel, dA is an area element for emission into the pixel p, dW_p is the solid angle subtended by the pinhole as seen from the emitting volume element, n_i(e) is the ion (electron) density at the emitting location, E¢ denotes photon energy emitted by an electron of energy E, ¶² s_ei/¶E¢¶W is the cross section for electron-ion bremsstrahlung differential in energy and solid angle, [^(n)] is a unit vector connecting the emitting volume with the camera element, and [^(b)] is a unit vector along the magnetic field, and T(E¢) is the transmission factor through vessel windows and other absorbers. The differential cross section used in our calculation is that given by formula 2BN of Koch and Motz [42]. The significance of ``2'' is the differential in energy and angle; ``B'' is Born approximation; ``N'' is no screening of the ions. Note that the distribution function f_e is that computed from the one-dimensional Fokker-Planck model, and is itself not influenced by the electric field E_dc. The absorption from vessel and windows is computed from fits to experimental absorption data [43].

3  LSC Usage

There are graphics calls embedded in the code for PPPL's ``Scientific Graphics,'' or SGlib , or just SG, and for ``gnuplot.''

Information and source for SGlib is at:
w3.pppl.gov/rib/repositories/NTCC/catalog/Asset/sglib.html.

Information and source for gnuplot is at:
www.cs.dartmouth.edu/gnuplot_info.html.

The combination LSC/gnuplot should work on any Unix, VMS, MS-DOS, or Windows system with a Fortran 77 compiler. Gnuplot for he Macintosh has been reported, but not found on the Internet in September 2000.

The combination LSC/sglib should work on any Unix or VMS system with a Fortran 90 compiler and a C compiler.

The source for LSC is constructed with calls for both sglib and gnuplot imbedded, but there are exceptions, as follows. Contour plots and the 3-D-like isometric sketch of ray paths in the toroid require sglib to be loaded; and, none of the graphs used by the X ray camera have been converted to gnuplot. Converting contour plots to gnuplot does not appear to be possible in a satisfactory way with a modest amount of work, and will probably never be done. Fortunately, these plots never seemed to be particularly useful. Converting the 2-D plots of the X ray camera to gnuplot is probably easy enough, and can be done if there is a demand for these graphs under gnuplot.

By default all calls to sglib are compiled into the executable code as no-operations with the cpp command #ifndef USE_SGLIB. In a similar way the X ray camera is not compied via the cpp command #ifndef USE_X_RAY_CAMERA.

As mentioned in the introduction, this manual and links to the source code can be found at w3.pppl.gov, through one of ~ignat, NTCC or topdac . The source code should be accompanied by at least one ``Makefile'' which works in the development environment, Linux Fujitsu Fortran, for the stand-alone code. Linux GNU Fortran 77 (g77) was also used in the development, however, the g77 executable is somewhat slower.  LSC is constructed as a subroutine to be called as a module of a transport code such as TSC or TRANSP, but stand-alone use with a simple driver program is easily done. The distribution provides a driver, lscdrive.F , and a simplifed mock-equilibrium writer, WrCircEq.F . Stand-alone runs can be extremely useful in studying the behavior and assuring reasonable convergence with the parameters chosen.

Naturally, stand-alone computations cannot show the interaction of the current driven with the electric field generated, reacting to the current, in a complete calculation.

3.1  Input Files

1. The plasma (geometry, fields, density, temperature); This specification is done through a large, strictly formatted, equilibrium file. The name given this file is arbitrary, but stand-alone runs have commonly used jardin.d and TSC runs have commonly used lhcdoua . For TRANSP the ``file'' is in memory; in other words there are no disk read/writes for the equilibrium. (Im vague on this; needs checking.)
2. LH parameters such as the frequency, the spectrum, and current/power diffusion. This specification is done through short NAMELIST (list-directed) file input.lhh under the namelist name inpvalue
3. Computational paramers such as the step size and the number of steps for each ray, strategy for smoothing in velocity space , number of iterations in constructing the distribution function. This specification is done through the same short NAMELIST file input.lhh under the namelist name inpexprt . (Putting the two namelists in different files might be a good idea in some cases, and doing so would be simple enough.)
4. Flags to specify graphs and printouts. This specification is also under the namelist name inpexprt .
5. LHCD power, i/o unit numbers, a flag on the strategy for re-tracing rays, and a flag which can turn off all plots . This specification comes in the argument list in the statement calling LSC.

Details follow in the following sub-sections.

3.1.1  Plasma

A header line (informational comments) equhd is read first followed by integers npsitm, nspc, kcycle, a plasma-time in the calling program, times (which is ignored by LSC) plasma parameters anecc, tekev, anicc, tikev, amass, achrg, the electric field (loop voltage) voltlp, and specifications of the equilibrium, rho, vptemp, xsv, pary, ppary, gary, gpary, nx, nz, isym, iplim, psimin, psilim, xmag, zmag, rgzero, bgzero, apl.

The meaning of the variables, and units used, are documented in the source.

3.1.2  LH parameters

The specification is done through the NAMELIST (list-directed) file input.lhh under the namelist name inpvalue

The namelist variables are documented in source file Inpval.inc .

The number of rays to be launched is nrays, which decomposed into the product of ntors and npols, so each toroidal spectral component can have a number of poloidal components, each of which behave somewhat differntly. Defaults are in params.inc .

There are two methods of representing the spectrum, a simple model of Gaussian shapes, and a calculation stemming from Brambilla paper [1].

Three tables list the LH parameters according to the following plan.

Table 1 lists the most important input paramters, and also the ones needed to specify a Gaussian spectrum.

Table 2 lists the input parameters to choose a Brambilla calculation of the waveguide spectrum. Note that some of the parameters of Table 1 are ``recycled'' for the Brambilla case, and therefore have a somewhat different meaning: nGrps, powers.

Table 3 lists the input parameters for the heuristic model of current and power diffusion.

Regarding Table 2, the `USRSPECn' waveguides are set equal to TFTRLHCD in the code as distributed. The intent is for the data statements to be replaced as needed. The format of the specification of the waveguide coupler is the following: az(i) is the toroidal distance to the near edge of i^th guide in cm; bz(i) is the toroidal distance to the far edge of i^th guide in cm; so that az(2)-bz(1) is the the width of the septum between the first two guides.

3.1.3  Computational parameters

Table 4 gives the computational parameters, and lists the defaults for each one.

3.1.4  Flags for graphs and print files

Typical input namelists, to be found in input.lhh, are found in Table 7.

 &inpvalue 
 fghz = 4.6, 
 nGrps = 1,
 powers(1) = 1.00, powers(2) = 1.00,
 phaseDeg(1) = 130., phaseDeg(2) = 130., nslices=301,
 ntors=15, npols=1,
 DoBram = 1, 
 couplers(1)='TFTRLHCD',
 couplers(2)='TFTRLHCD',
 DiffuJrf = 0.00, 
 /
 &inpexprt
 HstpLH = .01,  nstep = 20000,  nfreq = 50,  npsi = 100,
 nzones = 2000,
 nv = 401, nsmoo = 9, nsmw = 3,  weghtitr = .2,
 plflg( 1) = 0, 1, 0, 
 plflg( 4) = 1, 0, 0, 0, 
 prflg( 1) = 0, 0, 0,
 incLabls = 1,
 /

3.1.5  Calling arguments

Table 8 lists the arguments that must be supplied by the calling program, and is largely self-explanatory. The power is the most significant parameter. For use of LSC with a transport code and determination of both current and electric field in a self-consistent way iRayTrsi is crucial, and almost always equal to 0 because of the number of iterations needed.

3.2  Output files

3.3  Output graphs

Figure 1: ``Graphed'' input parameters to LSC for keeping track of the nature of the calculation by putting important information in the graphics output file.

Figure 2: A launched spectrum from a Brambilla calculation, including the negative-going minor components. Upper-right: Bar graph showing relative powers versus v/c . Lower-left: Same, but represented versus n-parallel. Lower-right: The launched spectrum as smoothed in velocity by the same smoothing function used in constructing the quasi-linear distribution function.

Figure 3: A typical result of the calculation for one ray. Left side: Evolution of n-parallel (top) and n-perpendicular (bottom) versus square root of normalized poloidal flux. The dotted line on the top frame indicates the region of strong damping, ie, if the magnitude of n-parallel is larger than the magnitude shown dotted, then the linear damping is strong. Quasi-linear damping may not be strong, owing to quasi-linear burn-through. Right side: Ray trajectory projected onto a poloidal cross section of the plasma. The ``roundness'' comes from an idealized test case; in general the cross section is elongated with the correct Shafranov shift of flux surfaces.

Figure 4: Power and current deposited without diffusing the current or the power. More rays would produce a smoother result. Horizontal axis is the square root of normalized poloidal flux. Left two: Volumetric and integrated power from the ray point of view for 1 MW launched power. Center two: As on left, but from the quasi-linear power deposition point of view. Right two: Volumetric and integrated current. The negative current stems from negative-going components of a realistic launched spectrum. The power and current found represented volumetrically and integrated from the center outwards.

Figure 5: Power and current deposited with diffusing the current and then the power. The current diffusion coefficient was set to 0.05 meter-squared per second, and the weighting of the diffusion effect on the power was 90 percent. Note that the center frames show un-diffused power.

Figure 6: A typical electron distribution and quasilinear diffusion coefficient

4  Appendix: Brief Description of Source Files

 &inpxry
 E_max = 0.20,  E_min = 0.01,  E_ph_min = 1.0e-6, 
 FoilCode = 'AG', 
 PusherMajor = 1.152, PusherMinor = 0.190 
 R_bound_max = 2.0, R_bound_min = 1.1,
 R_tangent = 1.525, z_tangent = 0.0,  
 Rpinhole = 2.661, Zpinhole = 0.0, pinhole_size = 0.00625,
 pinhole_size = 0.00625,
 dE_ph = 0.01, dFoilTCM = 0.0, 
 dlxray = 0.005,  focal_length = 0.495, 
 iAbsXray = 1,
 screen_d = 0.215, 
 /

Acknowledgments

F. W. Perkins devised the approach taken in the consideration of scattering. M. Ono , F. Paoletti , and H. Takahashi contributed discussions regarding the limited up-shift of parallel index of refraction. S. von Goeler , S. E. Jones , and D. P. Enright provided valuable help regarding the x-ray camera. V. Fuchs encouraged the approximate treatment of rf current diffusion.

Operation of the code by M. Shoucri (of the Tokamak de Varennes) in conjunction with TSC and by S. M. Kaye , R. V. Budny (of PPPL) and P. Stubberfield (of JET) in conjunction with TRANSP led to early improvements.

Operation by M. Ju (of KSTAR) led to a recent and significant improvement in the stability of the contruction of the electron distribution function from ray data. S. Bernabei and M. Okabayashi supported, encouraged, and advised this project over the past several years.

An earlier version of this manual depended on the work of A. J. Redd . The present manual re-uses many of his contributions.

References

[1]

M. Brambilla, ``Slow-wave launching at the lower hybrid frequency,'' Nucl. Fusion 16 47 (1976).

[2]

S. Weinberg, ``Eikonal method in magnetohydrodynamics,'' Phys. Rev. 126 1899 (1962).

[3]

I. B. Bernstein, ``Geometric optics in space- and time- varying plasmas,'' Phys. Fluids 18 320 (1975).

[4]

J. M. Wersinger, E. Ott, and J. M. Finn, ``Ergodic behavior of lower hybrid decay wave trajectories,'' Phys. Fluids 21 2263 (1978).

[5]

Yu. F. Baranov and V. I. Fedorov, ``Lower hybrid wave propagation in tokamaks,'' Soviet Physics Technical Physics Letters 4 332 (1978); Nucl. Fusion 20 1111 (1980).

[6]

P. T. Bonoli and E. Ott, ``Accessibility and energy deposition of lower hybrid waves in a tokamak with density fluctuations,'' Phys. Rev. Lett. 46 424 (1981).

[7]

P. T. Bonoli and E. Ott, ``Toroidal and scattering effects on lower hybrid propagation,'' Phys. Fluids 25 359 (1982).

[8]

D. W. Ignat, ``Toroidal effects on propagation, damping, and linear mode conversion of lower hybrid waves,'' Phys. Fluids 24 1110 (1981).

[9]

F. C. Jobes, et al., ``Current ramp-up by lower hybrid waves in the PLT tokamak,'' Phys. Rev. Lett. 55 1295 (1985).

[10]

C. F. F. Karney, N. J. Fisch, and F. C. Jobes, ``Comparison of the theory and practice of lower hybrid current drive,'' Phys. Rev. A 32 2554 (1985).

[11]

C. F. F. Karney and N. J. Fisch, ``Numerical studies of current generation by rf traveling waves,'' Phys. Fluids 22 1817 (1979).

[12]

C. F. F. Karney and N. J. Fisch, ``Current in wave-driven plasmas,'' Phys. Fluids 29 180 (1986).

[13]

N. J. Fisch, ``Theory of current drive in plasmas,'' Rev. Mod. Phys. 59 175 (1987).

[14]

P. T. Bonoli and R. C. Englade, ``Simulation model for lower hybrid current drive,'' Phys. Fluids 29 2937 (1986).

[15]

E. J. Valeo and D. C. Eder, ``Numerical modeling of lower hybrid heating and current drive,'' J. Comput. Phys. 69 341 (1987).

[16]

P. T. Bonoli, M. Porkolab, Y. Takase, and S. F. Knowlton, ``Numerical Modeling of the lower hybrid heating and current drive experiments in Alcator C,'' Nucl. Fusion 28 991 (1988).

[17]

V. Fuchs, I. P. Shkarofsky, R. A. Cairns, and P. T. Bonoli, ``Simulations of lower hybrid current drive and ohmic transformer recharge,'' Nucl. Fusion 29 1479 (1989).

[18]

K. Kupfer and D. Moreau, ``Wave chaos and the dependence of LHCD efficiency on temperature,'' Nucl. Fusion 32 1845 (1992).

[19]

H. Takahashi, et al., ``Accessibility for lower hybrid waves in PBX-M,'' in 1993 Conference on Controlled Fusion and Plasma Physics (Proc. Conf. Lisbon, 1994) Vol 17C, Part III, European Physical Society, Geneva (1993) 901.

[20]

J. Kesner, et al., ``Simulations of lower hybrid current drive in shaped tokamaks,'' Nucl. Fusion 34 619 (1994).

[21]

F. Paoletti, et al., ``LHCD Accessibility Study with Reconstructed Magnetic Equilibria in PBX-M,'' Nucl. Fusion 34 771-776 (1994).

[22]

K. Kupfer, et al., ``Statistical theory of wave propagation and multi-pass absorption for current drive in tokamaks,'' Phys. Fluids B 5 4391 (1993).

[23]

R. Hawryluk, in Physics of Plasmas Close to Thermonuclear Conditions, (Varenna, 1979), Commission of the European Communities, Brussels, I 61 (1979).

[24]

R. J. Goldston, D. C. McCune, H. H. Towner, et al., J. Comput. Phys. 43 61 (1981).

[25]

S. C. Jardin, N. Pomphrey, and J. DeLucia, ``Dynamic modeling of transport and positional control of tokamaks,'' J. Comput. Phys. 66 481 (1986).

[26]

W. M. Nevins, et al., ``Mission and Design of the Tokamak Physics Experiment / Steady State Advanced Tokamak,'' Fourteenth International Conference on Plasma Physics and Controlled Nuclear Fusion Research, Würzburg, Germany, 30 September - 7 October 1992, paper IAEA-CN-56/F1-5.

[27]

J. A. Schmidt, et al., ``The Design of the Tokamak Physics Experiment,'' J. Fusion Energy 12 455 (1993).

[28]

K. Tani, M. Asumi, and R. S. Devoto, ``Numerical analysis of 2D MHD equilibrium with non-inductive plasma current in tokamaks,'' J. Comput. Phys. 98 332 (1992).

[29]

R. S. Devoto, et al., ``Modeling of lower hybrid current drive in self-consistent elongated tokamak equilibria,'' Nucl. Fusion 32 773 (1992).

[30]

J. P. Bizarro and D. Moreau, ``On ray stochasticity during lower hybrid current drive in tokamaks,'' Phys. Fluids B 5 1227 (1993).

[31]

J. P. Bizarro, ``On the dynamics of the launched power spectrum during lower hybrid current drive in tokamaks,'' Nucl. Fusion 33 831 (1993).

[32]

D. W. Ignat, E. J. Valeo, and S. C. Jardin, ``Dynamic Modeling of Lower Hybrid Current Drive,'' Nucl. Fusion 34 837-852 (1994).

[33]

S. von Goeler, et al., ``Initial lower hybrid current drive results and first two-dimensional images of hard x-ray emission in PBX-M,'' in Proceedings of the Nineteenth European Physical Society Conference on Controlled Fusion and Plasma Physics, (Innsbruck, Austria, 1992) II 949 (1992).

[34]

S. E. Jones, et al., ``Calculation of an upper limit for an effective fast electron diffusion constant using the hard x-ray camera on PBX-M,'' Plasma Phys. Control. Fusion 35 1003 (1993).

[35]

J. P. Bizarro, et al., ``On self-consistent ray-tracing and Fokker-Planck/ modeling of the hard x-ray emission during lower-hybrid current drive in tokamaks,'' Phys. Fluids B 5 3276 (1993).

[36]

S. von Goeler, et al., ``A camera for imaging hard x-rays from suprathermal electrons during lower hybrid current drive on PBX-M,'', Rev. Sci. Instrum. 65 1621 (1994).

[37]

D. W. Ignat, R. Kaita, S. C. Jardin, and M. Okabayashi, ``Spreading of Lower Hybrid Wave Driven Currents in PBX-M,'' Nucl. Fusion 36 1733-1742 (1996)

[38]

Myunghee Ju, Jinyong Kim, KSTAR Team, ``A predictive study of non-inductive current driven KSTAR tokamak discharge modes using a new transport-heating simulation package,'' Nucl. Fusion 40 1859-1866 (2000)

[39]

W. Daughton, ``TSC Users Manual,'' internal PPPL working document dated August 1993; updated January 1997 by W. Daughton, S. C. Jardin, C. Kessel, N. Pomphrey

[40]

S. C. Jardin, ``TSC Users Manual,'' 1999.

[41]

William H. Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling, Numerical Recipes; The Art of Scientific Computing - FORTRAN Version, (Cambridge University Press, New York, 1989).

[42]

H. W. Koch and J. W. Motz, ``Bremsstrahlung Cross-Section Formulas and Related Data,'' Rev. Mod. Phys. 31 920 (1959).

[43]

W. H. McMaster, et al., ``Compilation of X-ray Cross Sections,'' University of California, Livermore, California, Report UCRL 50174, Sec. II rev. 1, May 1969.

[44]

T. H. Stix, Waves in Plasmas, (American Institute of Physics, New York, 1992).