Ralph Dux STRAHL User Manual IPP 10/30 September, 2006
Contents
1. Ar 1 37 Transformation of the transport equation to p yields for the radial derivatives 0 Or p 9 0p and 0 0r p 202 Op p 9 0p Using the symmetric three point formulas for the p derivatives leads to an increase of the truncation error as compared to uniform grids Thompson et al 1985 Fletcher 1997 however the increase is small as long as the change of Ar for two adjacent grid intervals i e d2r dp p p is low It is appropriate to keep the change of Ar below 10 1 4 Detailed Numerical Scheme The discretisation of the coupled impurity transport equations 1 12 and the numerical solution scheme as given in equations 1 31 1 32 is described in detail We rewrite Eq 1 12 for the ion stage z On 7 ld dn De an Og Mele me PRO 1 38 l The parallel loss in the SOL is approximated by a volume loss term with characteristic decay time 7jj The volume sources due to recombination and ionisation between ions with z z 1 are contained in the reaction term Rz R ne Sz 1Nz 1 Nz11 N Ne S 0 1 1 39 Q describes an external source e g fuelling of singly ionised ions by ionisation of neutrals or generation of Het ions by fusion of D and T The mesh points shall be equidistant in the coordinate p r Transformation of the transport equation to p yields with dp dr p dD On Sen o Uz Bee dp The change of the density at the radial mesh point i during a time2. a atanh y atanh y 1 ef ops 4 l a An example of a geometry input file prepared with grid circular is given on the following pages It includes small headers for all quantities and shows all the units which have to be used 20 cv cv cv cv cv cv cv rho volume LCFS cm number of grid points sqrt 0 00000 0 44207 0 78404 1 00448 rho 00000 35000 70000 00714 Or major 00000 11667 23333 33571 bere major 00000 88333 76667 66429 O OC O l safetv 1 00000 1 24500 1 98000 50 0 28 Psi Psi ax 0 06659 0 49783 0 82429 1 00893 volume rho 0 05000 0 40000 0 75000 1 01429 radius low field side Psi sep Psi ax 19786 60231 89935 20 174 0 13269 0 55125 0 86269 101335 R_axis cm 150 0 HOOG U_loop V 0 5000 21 volume LCFS 10000 45000 80000 02143 eooo 1 01667 LLIS 1 25000 1 33810 103333 1 15000 1 26667 1 34048 radius high 0 98333 0 86667 0 75000 0 66190 factor 1 00500 1 32000 2 12500 field side R_axis 0 96667 0 85000 0 73333 0 65952 bh 02000 40500 28000 DH 0 0 0 a PPpPeeEN 0 0 0 0 15000 50000 85000 02857 R_axis 05000 16667 28333 34286 95000 83333 71667 65714 04500 50000 44500 re AN eS a Ho 00 ROOO oooo points up to separtrix 26169 65104 93439 02212
3. 20000 55000 90000 03571 06667 18333 30000 34524 93333 81667 70000 65476 08000 60500 62000 CHAPTER 2 INPUT FILES time s 2 0000 fourier coefficients 3 32384 69752 96790 02646 E OGOOGO 25000 60000 95000 04286 HM O0 0 08333 20000 31667 34762 Rae 91667 80000 68333 65238 oooo 1 12500 1 72000 2 80500 Oi Orie leet lil ra RP 00 HH OO Hr 38403 74182 00000 03078 30000 65000 00000 05000 10000 21667 33333 35000 90000 78333 66667 65000 18000 84500 00000 2 2 CV CV CV CV CV GEOMETRY N 0 0 0 4 6 N w 67868 43526 29223 91828 96971 04964 00229 01266 01455 01420 10825 23182 25036 25715 27607 00000 00000 00000 17E 03 94E 02 63E 02 45E 01 51E 01 fraction of circulating particles 1 00000 0 81250 0 73624 0 51965 0 48928 0 46127 0 34660 0 32750 0 30940 Integral dl p B p m T 4 71238 4 73529 4 80396 5 82685 6 16481 6 54599 9 07297 9 69584 10 35518 lt B total gt IT 2 00000 2 00028 2 00107 2 00888 2 01036 2 01163 2 01463 2 01471 2 01467 lt Btotal 2 gt T1 2 4 00000 4 00166 4 00650 4 06336 4 07796 4 09295 4 17396 4 19217 4 21141 lt 1 B_total 2 gt 1 T 2 0 25000 0 25004 0 25015 0 25285 0 25403 0 25546 0 26650 0 26946 0 27265 lt RB T gt MAI 3 00000 3 00000 3 00000 3 00000 3 00000
4. e usr home strahl atomdat adf15 pec96 c xyz c2 dat The first name follows strictly the ADAS conventions while the second name omits the directory name in the file name The c2 after the last is constructed from the character after in pec96 c and the charge of the ion 2 Please note the added adf15 in the directory path Having found the right file the program looks for all lines which fall into the specified wavelength window A AA and adds up all the selected contributions excitation recombination CX of all lines Thus it is possible to add several lines of a pec file which are not resolved by the spectrometer Note that the wavelengths in the pec files are sometimes not the same as the measured wavelengths On the next page C_ atomdat is shown as an example 2 4 MAIN PARAMETER FILE 25 cv main ion brems SXR spectral brems thermal CX 1 1 1 0 cv diagnostic lines 1 cc begin atomicData acd acd96 c dat recombination scd scd96 c dat ionisation prb prb96 c dat continuum radiation plt plt96 c dat line radiation ccd ccd89 c dat thermal charge exchange prce prc89_c dat thermal charge exchange contiuum radiation pls plsxi c dat line radiation in SXR range prs prsxl_c dat continuum radiation in SXR range fis sxrfill dat sensitivity of SXR brs brs10400 dat spectral bremsstrahlung cc end atomic Data REEAKARRRERKARE REREAD LAGNOSTLC Lines kK KKK KK KK KK KKK KK KK cd excitation re
5. nn e the flux surface averages B B 1 B RBr RB B 1 RP e the fourier coefficients cos m B sin m B cos m0 B In B sin mO B In B The maximum order of m can be set in the header part of the grid file 2 2 GEOMETRV 19 The fraction of circulating particles f can be calculated with Be 3 B ae AA az VI AB For circular concentric flux surfaces the above quantities can be expressed analytically as a function of the major radius of the magnetic axis Ro the toroidal field at the magnetic axis B o the safety factor q and the inverse aspect ratio of the flux surface e r Ro The toroidal poloidal and total field vary as eda ON an ee P cos6 t l Hecosf 1 ecos9 with Bi o e p 0 l a an 0 to 1 FEED One finds for the flux surface averages di 2rr 2 Ba 2 1 3 5 B B B 1 B l B Bio 0 Vie BE 35 RB RoB 322 32 Po 3 ya t 02t 0 gt p Bl 2 R l 2 sin m8 B 0 cos m0 B 2 A vl e Bo sin m6 B ln B 0 cos m6 B In B To The program grid circular prepares a geometry file for circular concentric flux surfaces assuming a parabolic profile for the safety factor q where the values of q on axis and at the LCFS are specified by the user For a q profile which is given by q r qo 1 ae the polidal flux label ppor is given by atanh y atanh y V1 e th pa gt mm WI
6. Lackner et al 1982 Here the densities and radial derivatives of the transport part are again treated like in the time centred Crank Nicolson scheme however for the ionisation and recombination part an alternation between an implicit and explicit 1 3 NUMERICAL SOLUTION 7 method is chosen Thus there is a time step which calculates the ionisation rate per volume with the density 7 implicit and the recombinations with the densities nil explicit ltt ql Ai D dD Onit1 2 dos Sy pett A jii tat PET er ee ANA i PUN 1 1 2 At Or HG A dr Or Gua nit Rr d 1 31 followed by a time step where jil Fl is used for the recombinations and fil for the ionisations 2141 _ pl 271 1 2 D D l41 2 A 5 a On D dD an db y AA E y a Pat a STA mr 2 p A jal 1 22 At Or 7 i dr r Gr l Rn d 1 32 The method is again unconditionally stable Lackner et al 1982 and has the great advantage to lead to Z tridi agonal matrix equations on the number of mesh points which can very effectively be solved using the Thomas algorithm For the first time step Eq 1 31 the equations are solved in increasing order of the charge state so that the densities of the lower charge state at the new time point 1 are known when needed for the calculation of the ionisation source into the present charge state For the second time step Eq 1 32 the equations for the charge states are solved in decreasing order since
7. f 2r Raris Ve With this choice the averaged continuity equation is casted in the familiar cylindrical form Oni z 1 dt a rly z Q z 1 9 For the flux density perpendicular to the magnetic surface an ansatz with a diffusive and a convective part is used D 8 is the radial diffusion coefficient and v 8 is the radial drift velocity where the dependence on the poloidal angle 6 is explicitly denoted The radial covariant component T 7 z of the flux density is Prz D 0 Vnr z v n z D0 ZE Vr v On z 1 10 Here the density gradient has been seperated into a poloidally dependent and independent part using Vn7 z Onr z Or Wr The radial contravariant flux density is obtained by multiplying with Vr an I z T 2 Vr D Vr Be IZ y 9 Vriniz 1 11 Thus we obtain the radial impurity transport equation Onr Z i 0 anr Z lt Du 1 12 Ot T Far r Or ki nrz Qrz Here D and v are the flux surface averaged transport coefficients 27 DX Vr jp VridS Vr ao DIVE an DIVrdS a DON 1 1 a odo a Z 1 1 v v Vr TE f vas We Rast fi 8 78 do 1 13 For the collisional transport parameters we know the dependence of D and v on the poloidal angle and the correct surface averages are calculated in the neoclassical transport part of the code Fig 1 1 depicts the behaviour for the classical part of the collisional impurity flux The
8. velocitv or parameter cv f of interpolation points 30 CHAPTER 2 INPUT FILES n rho cv rho poloidal grid for interpolation rho 1 rho n_rho cv D m 2 s v m s alphall VO t 1 y 1 t 1 y n_rho t 1 VO t nt y l t nt y n_rho t nt funct invokes a gauss function with different decav lengths inside and outside some radial position Ti PAXTLCFS For r lt ra the function is r r y r yo y yo exp 5 Po For r gt r the function is 2 Dones a This function is quite flexible to produce various shapes The data block looks like cv function only for Drift funct velocity or parameter cv D m 2 s v m s alpha 1 pill yO t 1 yl t 1 y2 t 1 p0 t 1 pl t 1 p2 t 1 yO t n yl t n y2 t n po t n pl t n p2 t n funet2 This option uses the sum of two functions of the type funct Thus for every time point two lines have to be specified Each line has the same order and meaning as for the previous option funct const c This option just applies to the drift The code looks up the electron density profile in the middle of the time intervall and adjusts the drift velocity in such a way that the impurity concentration nimp ne in the source free region for ppor lt 0 9 is constant Thus the drift velocity for ppo lt 0 9 is 1 dne D 3 ne dr For Ppor gt 0 9 v is set to zero In the c
9. 01 in Vol gt 1 10e 00 in gt 1 65e 00 in gt 1 83e 00 in gt 2 09e 00 in gt 2 20e 00 in gt 4 55e 00 in gt 2 10e 00 in gt 1 93e 00 in gt 1 77e 00 in gt 1 85e 00 in gt 2 57e 00 in gt 3 40e 00 in gt 8 11e 00 in gt 4 72e 01 in gt 1 23e 01 in gt 4 76e 00 in 19 all stages gt 100 in Vol tot tot Vol Vol Vol Vol Vol Vol Vol Vol Vol Vol Vol Vol Vol Vol tot tot tot tot tot tot tot tot tot STRAHL 1999 Mon Jul 12 15 53 46 t 2 000008 a 61 6cm Z A plasm 1 2 imp 18 40 lt ne gt 6 90e 19m Te 0 3 53keV ne 0 8 58e 19m Zeff 0 1 12 for rho 0 1 0 4 0 9 D 0 10 0 16 2 20 m s v 0 0 0 0 0 0 m s neocl 0 CEX 0 influx s valve 1 90e 20 wall 0 00e 00 div 0 00e 00 div main 64 7 tau ms sol 2 5 div Inf pump Inf separatrix Te 1 01e 02eV Ne 3 49e 19m falloff Te 1 5cm ne 1 9cm w SOl 2 5cm lon Length 0 02cm IAI NOA SANILNOY ONILLO Id Iv Bibliographv Behringer 1987 Behringer K 1987 Description of the impurity transport code STRAHL Description of the impurity transport code STRAHL JET R 87 08 JET Joint Undertaking Culham Fletcher 1997 Fletcher C A J 1997 Computational Techniques for Fluid Dynamics Springer Berlin Hei delberg New York 2nd edition Fu mann 1992 FuBmann G 1992 Teilchentransport in magnetisch eingeschlossenen plasmen Technical Report 1 273 IPP Garching German
10. 1 For the magnetic flux surface whose poloidal projection is drawn on the left with a thick line the poloidal variation of the toroidal magnetic field red the gradient blue and the contributions to the classical part of the collisional radial impurity fluxes black are shown in the polar plot on the right All values are normalized to their maximum value the volume V inside a flux surface which leads to O 3 3 5 rza phi 2d5 f Qrzav 1 2 The flux surface shall be labeled with p where p rises from centre to edge Thus V p V pl is the unit vector perpendicular to the flux surface pointing towards the outside of the enclosed volume and the flux density integral becomes Vods as p I zd Iy 1 3 gr LZ pr N f ZT Tp 1 3 where Pr z is the radial contravariant component of the particle flux density PS z TizVp Using the flux surface average of an arbitrary scalar quantity F AVNTI dS 2 37 f IVpl and r av f Fae I doy 1 5 0 the continuitv equation 1 2 becomes p p OV OV OV Lobi ee se en il a f doz 3 Ciz Wp Q1z 1 6 0 0 1 1 IMPURITY TRANSPORT EQUATION 3 which can be differentiated with respect to p For a flux surface geometry which is independent of time we arrive at nrz NTa NV si ha ap ap ha 0 i We choose the flux surface label r which is calculated from the enclosed volume V The volume V depends on r as the volume of a cylinder with length 27 Razis V 1 8
11. 3 Te leV Til ev VO t 1 y 1 t 1 y n rho t 1 y0 t nt y l t nt y n_rho t nt decay lambda t 1 lambda t nt e interpa also causes quadratic interpolation between values at given radial positions however the interpolation is used for the whole radaial grid i e also outside the LCFS The data block looks the same as for the previous option The only difference is that the decay length block is left out cv cv cv cv cv cv times EXI a EICHE function interpa radial coordinate poloidal rho or volume rho of interpolation points n_rho radial grid for ne interpolation rho 1 rho n rho ne cm 3 Te eV Ti eV yO t 1 y 1 t 1 y n rho t 1 yO t nt y 1 t nt y n_rho t nt e exppolo exponential of an even polynom y p yo exp pop pip p2p pap 15 16 CHAPTER 2 INPUT FILES cv times nt tl su Ent cv function exppol0 cv radial coordinate poloidal rho or volume rho cv ne cm 3 Te eV Ti eV y0 t 1 pO t 1 pl t 1 p2 t 1 p3 t 1 yO t nt pO t nt pl t nt p2 t nt p3 t nt cv decay lambda t 1 lambda t nt exppoll exponential of a polynom with zero slope on the axis y p yo exp pop pip pap psp cv times nt t 1 t nt cv function exppoll cv radial coordinate poloidal rho or volume rho cv ne c
12. 3 00000 3 00000 3 00000 3 00000 lt R 2 B_p 2 B 2 gt m 2 0 00E 00 6 19E 04 2 41E 03 2 03E 02 2 37E 02 2 68E 02 3 52E 02 3 58E 02 3 61E 02 lt 1 R 2 gt 1 m 2 4 44E 01 4 45E 01 4 45E 01 4 A8E 01 4 48E 01 4 50E 01 4 57E 01 4 59E 01 4 61E 01 63E 01 WW A A A o N w 63095 41097 21991 07805 43518 17176 00382 01346 01436 02426 12390 329355 25071 25911 27972 00000 00000 00000 64E 03 15E 02 63E 02 45E 01 52E 01 66E 01 Ww o Ww 58960 38821 26038 28299 94154 93199 00552 01403 01412 03614 13998 27676 25122 26132 28359 00000 00000 00000 25E 02 31E 02 62E 02 46E 01 54E 01 69E 01 Ww o Ww 55285 36680 24560 53274 48782 32864 00724 01442 01384 04932 15662 30197 25192 26379 28767 00000 00000 00000 64E 02 44E 02 60E 02 47E 01 95E 01 71E 01 21 22 CHAPTER 2 INPUT FILES cv lt cos m theta B_total 2 gt T 2 0 00E 00 3 33E 02 6 68E 02 1 00E 01 1 34E 01 1 68E 01 2 03E 01 2 38E 01 2 73E 01 3 09E 01 3 45E 01 3 81E 01 4 18E 01 4 56E 01 4 94E 01 5 32E 01 5 72E 01 6 12E 01 6 53E 01 6 95E 01 7 38E 01 0 00E 00 2 78E 04 1 11E 03 2 51E 03 4 48E 03 7 03E 03 1 02E 02 1 39E 02 1 83E 02 2 33E 02 2 89E 02 3 5
13. 5 30E 04 6 92E 04 8 85E 04 1 11E 03 1 38E 03 1 68E 03 2 03E 03 2 43E 03 2 88E 03 3 39E 03 cv sin m theta B total 1In B_total gt TIn T 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 2 3 Atomic Data The atomic data are read from files which contain the according rates or emissivity coefficients in tabulated form Most of these data are stored according to ADAS conventions Since the ADAS tables are written in fixed format editing of these files has to be done with great care The interface between the data tables and STRAHL is a file named Xx atomdat It is located in the current directory and Xx is the element symbol e g Ar Ne For elements with a one character symbol like C and O C and O are used 2 3 ATOMIC DATA 23 The first two lines of Xx atomdat contain some switches e Bremsstrahlung due to electron scattering at the main ions is not calculated e Emission in the spectral range of
14. RECYCLING cv decay length of impurity outside last grid point cm 2 0 28 CHAPTER 2 INPUT FILES cv Rec ON 1 0FF 0 wall rec Tau div gt SOL ms Tau pump ms 0 thes 99 240 cv SOL width cm 5 0 N div Tdiv gt SOL divertor puff Ndiv Tpump Pout LE pump Figure 2 1 Recycling fluxes in STRAHL Impurity behaviour outside the LCFS in the regime of open field lines is not very well modelled in a one dimen sional code local effects parallel transport The following edge and recycling parameters are needed The flux of ions from the last mesh point towards the wall is specified by the diffusion constant on the last mesh point and a gradient length outside the grid which has to be given in the first line The radial flux density towards the wall is then Ty Dedgenz edge A The parallel loss is specified by a time constant TsOL div Which has to be given in the transport section see below A sketch of the recycling of impurities is given in figure 2 1 Wall recycling is specified by the factor R which is in the range from O no wall recycling to 1 complete wall recycling The particles that are lost parallel to the field lines go into the divertor reservoir From there the particles recycle with a time constant Tdiv sozr Or are pumped by the pumps with a time constant Tyump The particles Nq in the divertor reservoir cause a recycling flux rate Naiy Taiv sor and a rate of pumped partic
15. Soft X ray diagnostic is not calculated e Spectral emissivity of bremsstrahlung for one wavelegth is not calculated e Charge Exchange with thermal neutral hydrogen is not calculated e Line emission of special lines diagnostic lines is not calculated Then the names of the needed atomic data files are specified The names are given in a block which starts and ends with cc Only the first fifteen characters of the names are read by the code Each name is given after a three character label followed by a colon e g acd The label tells the program which data set is given in this line Data sets which are not needed for the current setting of the switches do not have to be specified The meaning of the labels are e acd rate coefficients for recombination e scd rate coefficients for ionisation e prb power loss coefficients for continuum radiation e plt power loss coefficients for line radiation e ccd rate coefficients for thermal charge exchange e prc power loss coefficients for thermal charge exchange contiuum radiation e pls power loss coefficients for line radiation in SXR range e prs power loss coefficients for continuum radiation in SXR range e fis sensitivity of Soft X ray diagnostic as function of photon energy e brs spectral power loss coefficients for bremsstrahlung All these data files are in the sub directory atomdat newdat In the diagnostic line section of Xx atomdat the data line identifier changes
16. and temporal step sizes 3 Ary s D dD Bow Tz l 3 S R 1 30 We arrive at a linear system of N x Z algebraic equations which yield the densities at the new time point The equation for the density of an ion stage at mesh point k contains the densities of the neighbouring radial mesh points as well as the neighbouring ion stages Thus the difference formulation is rather similar to that of a two dimensional transport equation The coefficient matrix of the linear system can be ordered such that it has tridiagonal form plus additional terms that are displaced from the main diagonal by the number of ion stages Z which result from ionisation and recombination An iterative solution of the matrix equation has been used by Behringer Behringer 1987 and is one of the two implemented solution methods in STRAHL To this end the full matrix equation is split into Z matrix equations on the number of radial mesh points which are solved with the fast Thomas algorithm These Z matrix equations are solved in ascending order of the ion charge and the densities of the lower charge state at the time point 1 are known for the calculation of the time centred ionisation rate into the present charge state The time centred recombination rate into the present charge state uses the density of the next higher charge state at l 1 as calculated by the preceeding iteration step An alternative method can be used in STRAHL which has been proposed by Lackner et al
17. calc 1 take distr from shot at time 0 0 0 0 On Ur TT PUNT cv save all cycles 1 save final and start distribution 0 al TIMES EPS cv number of changes start tim Fstop time 2 cv time dt at start increase of dt after cycle steps per cycle 1 0 l e 3 1 20 10 2 0 1 e 3 1 20 10 S UR CE cv position cm constant rate 1 s time dependent rate from file 1000 0 2 5e21 0 cv divertor puff source width in cm source width out cm 0 0 EDGE R ti CYCLING Tau div gt SOL ms 0 decay length of impurity outside last grid point cm Tau pump ns 39 240 cv 2 50 cv Rec ON 1 0FF 0 wall rec 0 1 cv SOL width cm 5 0 DENS ITY TEMPERATURE AND NEUTRAL Hr R GE N FOR CX 2 4 MAIN PARAMETER FILE 33 CV take from file with shot index 1 5 N amp LA 505 LACA G TRANSPORT method 0 off gt 0 of Drift l approx cv lt 0 figure out but dont use 2 3 NEOART neoclassics for rho_pol lt 0 2 95 ANOMALOUS TRANS PORT cv f of changes for transport 1 cv time vector 0 0000 cv parallel loss times ms DD cv Diffusion m 2 s funct cv D D SOL Dmax beta gamma r_max r LCFS 0 1 1 0 245 15 6 0 0 8 cv Drift function Drift Parameter Velocity const c velocitv cv of sawteeth inversion radius cm 0 1 5 cv times of sawteeth 0 0 Chapter 3 Execution and Output Once vou prepared all the necessarv input files vou finallv can run STRA
18. step At i e from t to t At is calculated by taking the spatial derivatives of the density at time t At 2 This approach leads to the Crank Nicolson difference scheme Denoting the densities at the old time with superscript O and the densities at the new time step with 1 the derivatives at the ith mesh point are given by z 2 2 1 BS HE ef AD U n P 12 Pip JD p an dn Neg Nz i dt At dn Ne ay EN MALA SRA dp 4Ap den N qa tiga MIL 2n NL tne AAA A ADS 1 41 dp 2A p The following notation for the derivatives of diffusion coefficient and drift velocity is used dD _ Desira Daina AD dp 2Ap 2Ap dv 0 Uzyitl 7 Uzi Av 1 42 dp 2Ap 2Ap 1 4 DETAILED NUMERICAL SCHEME 9 This yields the discretisation of the transport equation At D p nii tie A Gi ay Head Qn i nia niia At p p T Ff p Te 2 Nps Ap Mein FMg Ma ni 1 At v p 1 A AtR At 1 43 A on aes nz Q 1 43 With p p 2Ap and q p 2Ap the equation reads 1 0 At na a Tap D n Ne ig t ni 1 mn In ne Mia At 22 3 q D PAD puz n neta Fra nia 1 Ma 1 At ea se a pAv nes n i AUR AtQ 1 44 The reaction term changes it form each time step according to equations 1 31 1 32 A time step I with implicit treatment of ionisation and explicit treatment of recombination is fo
19. term Q7 z couples the transport equation of each ionisation stage with the neighbouring stages MER C Qiz neSiz ne01 2 HNHAL z NI Z H meSIZ 1N11 Z 1 near z ENHAT Z41 NI Z 1 1 14 Sr z is the rate coefficient for ionisation of impurity species J in ionisation stage Z ar z is the recombination coefficient for radiative and di electronic recombination from ionisation stage Z and af is the respective recom bination coefficient due to charge exchange The sum of the transport equations for all ions not the neutral of one species is a very useful simplification since all source terms cancel except for the ionisation and recombination between the neutral atom and the single ionised ion Ont 10 On r D 1 15 Ot T e Or on r D and v are average values of the transport coefficients for each ion stage weighted with the fractional abundance of each ion stage at a certain flux surface 1 2 Neoclassical Transport Two neoclassical code packages are now coupled with STRAHL In both toroidal rotation of the plasma is not taken into account and equal temperatures of main ion D and impurity I Tp Tr T are assumed In the first package the neoclassical transport coefficients are calculated with the numerical code NEOART by Arthur Peeters It solves the set of linear coupled equations for the parallel velocities in arbitrary toroidally sym metric geometry for all collision regimes The classical fluxes are giv
20. the recombination source with the densities at the new time point are now needed The truncation error of the two steps E a On E By ee on 5 a gt gt At A A On Es E S 1 2 CN Ty S R Er 1 33 are similar to the time centred scheme however additional terms linear in At appear which change sign from step to step The numerical solution is found to oscillate from time step to time step and has to be taken after an even number of time steps A detailed treatment of the discretisation scheme 1 31 and 1 32 can be found in section 1 4 At the inner boundary the plasma axis a vanishing density gradient and zero drift velocity are requested For r 0 the transport equations reads Oni r 0 ao d A A gt TY MA BIE 9 Rite 0 4 1 34 Ot Or r 0 drir 0 vile Here the second derivative of the densities is replaced by 2 r o Ar while the gradient of the drift velocities at r 0 has to be approximated with the asymmetric two point formula Dirichlet boundary conditions at the outermost grid point are easily incorporated by setting the density of the last grid point to the given value More commonly a decay length is used as boundary condition In order to maintain a second order truncation error in the radial step size this boundary condition is implemented with the symmetric three point formula using the fictitious grid point N 1 i e On tiny TiN 1 fi
21. time steps For the first half time step the respective equations are solved in ascending order of z and for the second half time step the deschending order is used When using a value less or equal 10 and greater than 20 for the threshold parameter the Lackner method is used without splitting the time step into two half time steps This method is more for educational purposes to see the oscillation of the solution between even and uneven number of time steps START COUN DAL EA ONS cv start new 0 from old calc 1 take distr from shot at time 0 12345 1 50 During a STRAHL run the impurity distribution can be saved and can be used in a later run as start conditions You can either start from zero impuritv densitv choose 0 or from the distribution of a previous run choose 1 Shot number time and the element svmbol is used to find the file with the old impuritv distribution In this example STRAHL would search for the file result C 12345t1 00 2 4 MAIN PARAMETER FILE 27 OUT PUT cv save all cycles 1 save final and start distribution 0 1 The time steps in a STRAHL run have to be defined before run time A number of equal time steps is called a cycle After each cycle STRAHL stores the current impurity distribution and related quantities choose 1 Sometimes it 1s not necessary to store the complete time evolution and one just wishes to have the result at the end of the run In this case choose O and only the start distribution a
22. to ionisation the density of neutral atoms decays with decreasing radius pedge r nS no r no r 19 exp o gr 1 58 T T and the radial distribution of the ionisation source from neutrals to singly ionised ions is Qo 1 r no r ne r So r 1 59 The source distribution Qo 1 r is needed for the solution of the transport equations 1 12 and is calculated from equations 1 57 and 1 58 Volume integration of the source distribution yields the total number of ionised neutrals per unit time pedge Qo41 anh no r ne r So r rdr 1 60 It is assumed that the plasma is sufficiently hot and dense to cause complete ionisation of the neutrals Thus the total number of ionised neutrals per unit time equals the total number of neutrals entering the plasma per unit time However the mesh points at the edge are usually not sufficently dense and the numerically integrated total source Qo_41 might be substantially different from Thus we multiply the source distribution Qo 1 7 with the factor Qo 1 which is 1 for correct numerical integration to insure the correct particle balance even for very strongly radially decaying profiles of the neutral impurity Chapter 2 Input Files Everv STRAHL run needs several input files to specifv the parameters of the plasma background e g electron den sitv electron temperature the geometrv of the plasma the transport parameters and the atomic data for ionisation recombina
23. which are just different by a few parameters The pp file has the following structure e data block for the electron density profiles e data block for the electron temperature profiles e data block for the ion temperature profiles optional e data block for the neutral hydrogen density optional An example for a data block describing the electron density for the time points t 2 0s and t 4 0s looks like cv time vector 2 2 0 4 0 CV ne function 14 CHAPTER 2 INPUT FILES interp cv radial coordinate poloidal rho cv f of interpolation points 19 cv radial grid for ne interpolation 0 0000 0 1000 0 2000 0 3000 0 4000 0 5000 0 6000 0 7000 0 8000 0 8200 0 8400 0 8600 0 8800 0 9000 0 9200 0 9400 0 9600 0 9800 1 0000 cv ne cm 3 5 000E 13 1 000E 00 9 752E 01 9 076E 01 8 088E 01 6 977E 01 6 019E 01 5 464E 01 5 240E 01 5 043E 01 4 959E 01 4 826E 01 4 673E 01 4 477E 01 4 195E 01 3 827E 0 3 430E 01 3 080E 01 2 74IE 01 2 355E 01 6 000E 13 1 000E 00 9 752E 01 9 076E 01 8 088E 01 6 977E 01 6 019E 01 5 464E 01 5 240E 01 5 043E 01 4 959E 01 4 826E 01 4 673E 01 4 477E 01 4 195E 01 3 827E 01 3 430E 01 3 080E 0 2 747E 01 2 355E 01 cv ne decay length cm in rho_volume 4 2 4 2 Each block contains the following data e number of time points followed by the time points e type of parametrization e radial coordinate of parametrization e one profile for every time point The code
24. 2E 02 4 22E 02 5 00E 02 5 84E 02 6 76E 02 7 77E 02 8 85E 02 1 00E 01 1 13E 01 1 27E 01 0 00E 00 2 32E 06 1 86E 05 6 28E 05 1 50E 04 2 93E 04 5 10E 04 8 15E 04 1 22E 03 1 76E 03 2 43E 03 3 26E 03 4 27E 03 5 48E 03 6 91E 03 8 59E 03 1 05E 02 1 28E 02 1 54E 02 1 84E 02 2 17E 02 cv lt sin m theta B_total 2 gt T 2 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 0 00E 00 cv lt cos m theta B total ln B_total gt Tln T 0 00E 00 1 67E 02 3 34E 02 5 01E 02 6 69E 02 8 37E 02 1 01E 01 1 18E 01 1 35E 01 1 52E 01 1 69E 01 1 86E 01 2 03E 01 2 21E 01 2 38E 01 2 56E 01 2 74E 01 2 91E 01 3 09E 01 3 27E 01 3 46E 01 0 00E 00 6 95E 05 2 78E 04 6 27E 04 1 12E 03 1 75E 03 2 52E 03 3 44E 03 4 51E 03 5 72E 03 7 09E 03 8 60E 03 1 03E 02 1 21E 02 1 41E 02 1 63E 02 1 86E 02 2 11E 02 2 37E 02 2 66E 02 2 96E 02 0 00E 00 3 86E 07 3 09E 06 1 04E 05 2 48E 05 4 86E 05 8 43E 05 1 34E 04 2 01E 04 2 88E 04 3 97E 04
25. 3 NUMERICAL SOLUTION 5 dinnp HP dinT dr KPS dr where the term in the first row of eq 1 19 reduces to 247 B in the case of large aspect ratio and circular geometry The factors KFS and H are functions of the impurity strength parameter a n Z np and the main ion collisionalitv v vpp vpr Rog vine viz DE ZJ 1 20 JARA 0 52a ii 0 59 a 1 34 3 2v 2 1 2 grs y 0 29 0 680 1 21 2 059 0 134 3 2p5 2 The banana plateau coefficients are RBr ke Tuy pe ee E 1 22 LZ Gare Zn as dinn dinT vfz Dpp z HBP 1 23 where the term with the flux surface averages in eq 1 22 reduces to q e BZ in the case of large aspect ratio and circular geometry e is the inverse aspect ratio u y is the weighted viscosity coefficient uoo of impurity and en MDNDHDoo MINILIO lin 1 1 24 MPDNDHDOO MINIHIOO The factor HFF in the drift term is connected to the ratio of the viscosity coefficients 401 Hoo of impurity and main ion Brei PU ais ae 1 25 uno Z H100 All viscosity coefficients are calculated using the expressions given by Kim Kim etal 1991 For the poloidal plasma cross section of ASDEX Upgrade with elongation x 1 6 and average triangularity 6 0 2 the ratios of the true geometry factor and the large aspect ratio circular geometry approximation are 1 1 1 3 at pro 0 3 0 9 for the CL term 0 75 0 42 for the PS term and 0 75 0 44 for the BP
26. Chapter 1 Introduction STRAHL is a code to calculate the radial transport and the emission of impurities in the plasma bulk The code solves the radial continuity equation for each ionisation stage of the impurity in 1 D geometry For the transport an ansatz of anomalous diffusivities and radial drift velocities is used a full neoclassical treatment of impurity transport can be switched on if desired The model focuses on the impurity transport and radiation while the parameters of the background plasma are taken from the experiment In this chapter the main aspects of STRAHL are introduced The first version of STRAHL was written around 1980 Over the years the code has gone through several mod ifications concerning the treatment of atomic data neoclassical transport and the general scheme of input output processing Even though the description of the numerical algorithm Behringer 1987 is still valid a manual has become necessary to introduce the new user to the modifications and to describe the technical details of data input output in STRAHL 1 1 Impurity Transport Equation A short introduction to the radial transport equation in the plasma bulk and the definition of the transport parame ters will be given first The law of particle conservation for the particle density nr z of impurity J in ion stage Z may be written as onr z _vf 1 1 ET VIr z Qr z 1 1 where 1 z is the flux density of the impurity and Qr z represents
27. HL and look at the results with several IDL plotting routines When using STRAHL and the plotting routines vou should know two things about the input on the screen There are alwavs default values for the control variables which are shown in brackets at the end of the input line The default values can be accepted by entering a slash Very often there is a list of key words and one of the characters is displayed in inverse mode white on black You activate the option by entering the highlighted character lower or upper case 3 1 Executing STRAHL An interactive STRAHL run starts by just typing st rah1 in the current directory The code comes back and asks for the name of the main input files Give the file name without the path specification i e without param files Now you are asked whether you want to halt the calculation at some time in the time intervall you were defining in the input file By default STRAHL runs to the end of the simulation time After the transport calculations have finished the radiation is calculated and then STRAHL asks you whether you want to save the current distribution option A or whether you want to save the time development option S If you are not at the end of the simulation time you can check what you got so far by using one of the plotting routines and come back to STRAHL to continue the run STRAHL stores the result in two files e Every run is saved in result strahl_result dat Thus the next run o
28. IPP Report Ralph Dux STRAHL User Manual IPP 10 30 September 2006 Dieser IPP Bericht ist als Manuskript des Autors gedruckt Die Arbeit entstand im Rahmen der Zusammenarbeit zwischen dem IPP und EURATOM auf dem Gebiet der Plasmaphysik Alle Rechte vorbehalten This IPP Report has been printed as author s manuscript elaborated under the col laboration between the IPP and EURATOM on the field of plasma physics All rights reserved STRAHL User Manual Ralph Dux Max Planck Institut fiir Plasmaphysik D 85748 Garching b M nchen Assoziation EURATOM IPP September 2006 Contents 1 Introduction 1 1 Impurity Transport Equation 1 2 Neoclassical Transport 5 4 3 28 A KOR A ER Ban WH WA SS 1 3 Numerical Solution crre Sack ee m fan Bo os tae Ra es HB ak e wl a 1 4 Detailed Numerical Scheme Como nun 1 5 Neutral Impurities a2 22 22 he ee Sea ER ER ae Dre oe eee a 2 Input Files 2 1 Plasma Backeround ars ts ee Oe ee E Be REE U 2 27 SGEOMeIhy en ok a at tebe te as Seat Doe et ew eel Gee gd ee 2 3 Atomic Data sies ea scat ete er a ace Bock ee ee ores le ee SD 2 4 Main Parameter File 4 4 3 4 2 4 4 eee PR aS Eh a A eee Gg A SR ee ET 3 Execution and Output 341 Exe utine STRAHL y 2 0 20 205 a 8 8 2 a bb Ble et ee DA ee IL The Result Files is 2 au ma a td ee du fe SNe as od hrs ii 3 3 Plotting Routines for IDL 2 2 4 0 eA se we 2 PE PE A AG OR A BA A 0o fe m
29. ase of the drift velocity the line which contains the identifier of the functional form has to be completed with velocity or parameter to notify whether v or a shall be used For the option const_c the code wants to read this parameter even though it is not needed CV CV of sawteeth inversion radius cm 0 15 times of sawteeth 0 0 2 4 MAIN PARAMETER FILE 31 The last two lines give the possibilty to invoke a very simple sawtooth model You just give the inversion radius and the times for sawtooth collapses At the crash times the code redistributes the impurity density inside r V2 x Tinversion to completely flat profiles keeping the total number of impurity particles constant The time intervall A is reduced to the settings which are in the last line of the time step block Also the number of time steps and the parameter how to increase At are taken from there Now you might have a look at the complete input file and ask yourself whether you still remember everything 32 CHAPTER 2 INPUT FILES ELEMENT cv element atomic weight amu energy of neutrals eV ASES 12 1 0 cv background ion atomic weight charge 25 1 G Res De FP 15 LE cv shot index 99999 0 GRID POINTS AND ITERATION cv K number of grid points dr center cm dr edge cm 6 0 101 2 0 0 1 max iterations stop iteration cv at fixed time if change below 100 0 02 START CONDITIONS cv start new 0 from old
30. bes the source rate by ionisation of neutral impurity atoms An unconditionally stable and effective discretisation method for diffusion convection equations is the Crank Nicolson scheme Here the change of the density during the time step At from time point to 1 is calculated by taking the average values of densities and spatial derivatives at the time points and 1 for the terms which appear on the right hand side of Eq 1 26 With the notation al gl 1 pi Y n 1 27 3 1 27 the Crank Nicolson scheme has the following form gr jil poe 2 5 dD 5 on 1 2 ii e ta din jap At Or r dr Or r dr SHH RR 4d 1 28 An equidistant mesh of N radial grid points with radial step size Ar shall be considered The spatial derivatives at mesh point k in Eq 1 28 are replaced by the symmetric three point formulas which have a truncation error quadratic in Ar Of _ fe fr Ar rof A le one a tar Of Fur 2 fo fr Ar r tf i El SEE TA EN 0 Ar 1 29 The truncation error of the finite difference formulation can be calculated by expanding the densities and transport coefficients at the different nodes in a Taylor series and forming the difference with the original transport equation When expanding around t At 2 the truncation error reads a me ale Or r dr y Or dr Or ars At on d PR D db R Ott EAS E EEE A case e A a 24 5 r T dr Ot r T dr 0 Ordot a The truncation error is quadratic in the radial
31. cal varift D radial distribution or time development for one radius e collisionality for main ion and impurity plus Pfirsch Schliiter and banana plateau limit radial distribution or time development for one radius e various volume integrals of total power loss density time development e impurity particle rates for flux through valve recycling fluxes perpendicular and parallel losses time devel opment e number of impurity particles in the confined plasma at the wall in the divertor and in the pump time development e energy content of the plasma time development e emissivity due to line radiation of diagnostic lines radial distribution or time development for one radius e emissivity due to line radiation of diagnostic lines following corona ionisation equilibrium radial distribu tion or time development for one radius e Intensity of diagnostic lines when integrating the emissivity along the major radius at the mid plane l e R dR time development e the same for the total line radiation of each stage time development e the same for the total soft X ray radiation time development e mean poloidal flux label of emission shells of diagnostic lines L a E P Pro R AR T e R dR time development 40 CHAPTER 3 EXECUTION AND OUTPUT e mean electron temperature of emission shells of diagnostic lines ie e R dR time development e mean impurity density in emission shells of diagnostic lines ji
32. classical diffusion coefficient Doy has a magnetic field dependence Dcr x 1 B and the polar plot of Fig 1 1 depicts the integrand which appears in Eq 1 13 All terms i e the magnetic field dependence the gradient and the larger surface on the outboard favour the radial flux contribution from the outboard side being about a factor of 5 larger than the inboard contribution The fluxes connected to the classical drift velocity vcz have just the same dependence vcz is proportional to Dcr times a combination of gradients of temperature and density which leads to vc Vr B For the turbulent transport it is also expected that the main contribution is on the outboard side however the exact poloidal dependence is not known It is important to note that D and v depend on the choice of the flux surface label Multiplying r with any constant number k would blow up or shrink our cylinder and yield transport cofficients D k D and k kv The best choice is a flux surface label with surface averages of V p close to 1 which is fulfilled by r The superscript x will be omitted from now on Outside the LCFS in the regime of open field lines the one dimensional model is not very well suited Parallel transport to the divertor limiter can be considered here in the form of a parallel loss time 7 Thus in the scrape off layer SOL a term n7 z T is added on the right hand side of equation 1 12 4 CHAPTER 1 INTRODUCTION The source sink
33. combination charge exchange 0 0 cd of lines cd charge of ion wavelength A half width of window A file extension 1 903 7 3 pju 2 977 0 3 pju 2 1175 7 3 pju 2 4650 1 3 pju 3 1549 1 3 pju 3 312 4 3 pju 3 384 1 3 pju 3 419 7 8 pju 4 40 3 idi pju 5 33 8 1 pju 2 4 Main Parameter File The main parameter file is located in the sub directorv param files and has arbitrarv name At run time STRAHL asks for the name of the main parameter file It contains information about the anomalous transport the impuritv source and recvcling and some code specific data like the time step evolution We now go through an example file line bv line The complete example is given at the end of this section ELEMENT cv element atomic weight amu energy of neutrals eV DEL 12 1 0 cv main ion atomic weight amu charge 26 CHAPTER 2 INPUT FILES The first line contains the element symbol the atomic weight amu and the energy of the neutrals eV when starting at the wall The element symbol is just used to construct the file name of Xx atomdat and the same name conventions apply here Then follows the atomic weight amu of the main plasma ion and the charge GR E D cE DE cv shot index 99999 0 The next line gives the number and the index of the geometry file Here the geometry is taken from Mmete grid 99999 0 GRID POCEN TOS AND OLE ETR A TEJO N cv K number of grid points
34. dr_center cm dr_edge cm 6 0 101 2 0 O L max iterations stop iteration Gv at fixed time if change below 100 0 02 Then follows the grid parameters STRAHL uses a grid which is not equidistant in r You may choose from the following two options see chapter 1 e Setdr center to a negative number specify K with K gt 1 and the number of grid points Here the grid is equidistant in rX e Specify dr_center dr_edge and K The distance of the grid points is then given by equation 1 37 and the number of grid points is calculated by STRAHL In the next line the numerical treatment of the differential equation solver can be influenced When the equations are solved in ascending order of the charge state Z the density n 41 Of the next higher stage at the new time point is not known and similar in descending order n _ is unknown STRAHL iteratively solves the equations in ascending order by calculating the recombination from n 1 in the first iteration In the following iterations ni 1 1s taken from the previous iteration step The iteration stops when the relative change of the densities at the new time point is below a certain threshold which can be set here When using a value less than O and greater than 10 for the threshold parameter second parameter of the second line the iterative method is switched off and the Lackner method according to eqautions 1 31 and 1 32 is used To this end the time step is split into two half
35. e e R dR time development e spectrum of the diagnostic lines for a given line width You should be able to find out how the program works There is also an option to print the data on a PostScript printer One option should be noted You might call the programs with the keyword save_file and input file You start by using pstrahl save_file my_plots Now you continue with pstrahl and produce some plots on the printer After changing something in the input parameter sets you repeat the STRAHL run and want to look at the same type of plots for the new result Now you just enter pstrahl in put_file my_plots The previos session of pstrahl was saved in the file my_plots and is now repeated by taking the commands from the file When having a number of plots this might be a little fast on the screen and it is only useful to combine this option with plots that are printed An example plot of the impurity density distribution of argon as produced by pst rah1 is shown on the following page Impurity Densities Ar00001t1 00_2 00_0 5 0 1016 4 0 1018 3 0 1016 16 _ 1 2 0 10 1 0 1016 0 16 AA oe 49___ 19_ 0 00 0 10 0 20 0 40 rhO os m 0 Z 0 1 Z 1 2 Z 2 3 Z 3 4 Z 4 5 Z 5 6 Z 6 7 Z 7 8 Z 8 9 Z 9 10 Z 10 11 Z 11 12 Z 12 13 Z 13 14 Z 14 jem p u SN UI IN IN IN oo NOOO 18 Z 18 gt 4 51e 02 in Vol gt 1 57e 01 in Vol gt 4 49e
36. e usage of read strahl pro is very easy when you have some basic knowledge of IDL The header of the routine hopefully gives enough explanations pro read_strahl file_name var_name data time time rho rho stage stage gatr gatr all all verbose verbos ERKRANKEN ER BRN BRE IRB DEBRA KR IR RBS IR ARE RK RENTEN Reads variables and attributes from STRAHL result files T Input gt file name a Give full name including path b name without path takes file from SHome strahl result 7 where Home is the home directory as stated in the UNIX environment variable Home gt var name Name of the variable attribute for an attribute set the keyword gatr Keywords gt time 0 gt all times are retrieved H t gt 0 gt only result at time t is retrieved gt rho 0 gt al radii are retrieved H r gt 0 gt only result at normalized poloidal flux label H rho r is retrieved gt stage 0 gt all stages lines are retrieved H n gt 0 gt only stage line n is retrieved 38 CHAPTER 3 EXECUTION AND OUTPUT gt all 1 gt get all data of the variable gt gatr 0 gt read a variable i 1 gt read an attribute gt verbose write information to the screen Output lt data data H kkikkkkkkkkkkkkkkkkkkAkkkAkkAkkAkkkAkkAkkAkkAkkAkAAkkAkkkkAAAkkAkAkAkAAAAX Generallv the meaning of the different dimensions of the variabl
37. en by eq 5 9 and 5 10 in section 5 of Hirshman and Sigmar Hirshman and Sigmar 1981 The equations for the banana plateau contribution are equal to that used by Houlberg Houlberg et al 1997 The Pfirsch Schl ter contribution is calculated from the coupled equations 6 1 2 and 6 14 15 of Hirshman and Sigmar Hirshman and Sigmar 1981 as described in Ref Peeters 2000 For all contributions a reduced charge state formalism is applied In the second package approximative analytical expressions according to Hirshman and Sigmar Hirshman and Sigmar 1981 are used We follow the description given by FuBmann Fu mann 1992 The classical diffusion coefficient and drift velocity is given by 1 HEBE mikgplvip Dee le ee 1 16 Da ETEA B ez oe dinnp 1dlnT CL _ pOL MAL of DFG2 Ir 5 dr 1 17 Here W is the poloidal magnetic flux B is the poloidal magnetic field and vrp is the binary collision frequency of impurity and main ion The binary collision frequency of species a with species b is for equal temperatures Ay2re 1 MaMbp 472 InA Vab n a 3 Ma m m kg M3 R 1 18 The brackets denote a flux surface average as defined in equation 1 4 The product of the first two terms in equation 1 16 reduces to 1 B for the large aspect ratio circular cross section case B is the vacuum toroidal field on axis The Pfirsch Schl ter coefficients are DIE e B KPS mikgTvip MW 1 19 1
38. es is obvious However the radiation stages need some extra explanation e impurity_radiationand sxr radiation index 1 total line radiation of neutral impurity index 2 total line radiation of singly ionised impurity index n total line radiation of hydrogen like ion index n 1 bremsstrahlung due to electron scattering at main ion index n 2 total continuum radiation of impurity bremsstrahlung and recombination continua index n 3 bremsstrahlung due to electron scattering at impurity index n 4 total radiation of impurity and main ion if set in Xx atomdat e spectral bremsstrahlung index 1 0 index 2 bremsstrahlung due to electron scattering at singly ionised impurity index n 1 bremsstrahlung due to electron scattering at fully ionized impurity index n 2 bremsstrahlung due to electron scattering at main ion index n 3 total bremsstrahlung of impurity and main ion if set in Xx atomdat 3 3 Plotting Routines for IDL There are two plotting routines available which run under IDL They are called pstrahl and pdiag All the routines which are needed to run these two programs are in the directory dl The following selection of plots is implemented so far pstrahl e impurity densities of each ion stage and total impurity density radial distribution or time development for one radius e impurity densities neglecting the transport effect i e in corona ionisation equilibrium radial distribution or time development for one
39. flux surface at the height of the magnetic axis i e for Z Zmag e major radius R for z 0 at low field side unit cm e major radius R for z 0 at high field side unit cm The first two flux surface labels p and r are used to map from the r to the ppo grid which is usualy taken by the other tokamak diagnostics The pro grid is not essential for the code if you do not use ppo in your densitv temperature profile definition The major radius R at z 0 is only used in the plotting routines to calculate line integrals of the radiation at mid plane All the following quantities are just needed for the calculation of neo classical transport and may be set to zero if this feature is not needed Many of these quantities are flux surface averages The definition of the flux surface average was given in eq 1 4 An alternative definition of the identical quantity is the following expression Adlp Bp y dl B where di is the length along the flux surface in poloidal direction B is the poloidal magnetic field and the integration is over one poloidal turn It is only defined inside the LCFS A NEOART uses a fourier expansion of some quantities in the generalized poloidal angle This angle is defined as _ ondo Bal Bo O 27 Bal Bp We need to specify the following quantities e the loop voltage e the safety factor q e the fraction of circulating particles fe BL This quantity is also related to dW dr 211 Raxis
40. from cc to cd_ In the first line there are switches which indicate whether emission of diagnostic lines following excitation recombinaion and or charge exchange with thermal neutral hydrogen is considered The number of diagnostic lines is given in the next line Then follows the block of diagnostic lines which contains the information to choose the emission lines from specific files e the charge state of the ion e the central wavelength A of the wavelength window in Angstrom e the half width of the wavelength window A in Angstrom e a three character file name identifier Some of the information in this data block is needed to construct the name of the file which contains the photon emissivity coefficients pec files The construction of the file name is somewhat complicated and is explained using an example of a CIII line The file name is build from the nuclear charge of the element 6 the ion charge 2 the file name identifier e g xyz and additional information about the directory path in a file called pec files in the current directory In pec files the code goes to line number 6 nuclear charge of C and reads the main directory e g usr home strahl atomdat and the path name for the subdirectory of the element e g pec96 c From this information the code looks for a file with one of the following two names 24 CHAPTER 2 INPUT FILES e usr home strahl atomdat adf15 pec961tc pec96tc xvztc2 dat
41. l rho or volume rho cv f of interpolation points n_rho cv radial grid for interpolation rho 1 rho n_rho cv neutral hydrogen density cm 3 y0 t 1 y 1 t 1 y n_rho t 1 yO t nt y 1 t nt y n_rho t nt 2 2 Geometry The plasma geometry is assumed to be independent of time The parameters for the plasma geometry are in the directory nete The filename is nete grid nnnnn i where nnnnn is a 5 digit number shot number and i is a l digit number index The program reads all parameters with a free format read statement in Fortran77 read unit Each item in the data block is preceded by a line which starts with cv see 2 1 The grid file gives some flux functions or flux surface averages on a common mesh of flux surfaces Most quantites are needed for the calculation of neoclassical transport The quantities which have to be specified on the mesh inside and outside the LCFS are the following 18 CHAPTER 2 INPUT FILES e polidal flux label ppo defined inside and outside LCFS y Waris Ppol ow a re Viors Waris e r is the radius of a circular torus with same volume V as the flux surface contour unit cm Inside the LCES r is well defined V r 4 gt gt 5 Qn Raxis Outside the LCFS r is not defined For diverted ASDEX Upgrade discharges I blow up the last close flux surface to match for a given ppo the the mean value and the difference of the major radii of the true
42. les Naiw Tpump All recycling particles enter again as neutrals The recycling can be switched on and off with the first parameter in that line The last parameter in that block specifies the radial width in units of r of the open field line region SOL DENS ITY TEMPERATURE AND NEUTRAL HYDROGEN FOR CX CV take from file with shot index 1 5 This line specifies which file for the plasma background shall be used In this example it would be the file mete pp00001 5 N E07 Col A 5 5 LC 4 L TRANSPORT method 0 off 20 of Drift 1 approx 2 4 MAIN PARAMETER FILE 29 cv lt 0 figure out but dont use 2 3 NEOART neoclassics for rho_pol lt 0 2 95 The rest of the input file deals with transport The neo classical transport can be switched on by using a positive number for the first parameter 7 In the code the neo classical drift velocity is multiplied with 7 100 Forn lt 0 the code calculates the neo classical transport but does not use it for the solution of the differential equations For n 0 neo classical transport is switched off completely The different methods to figure out the neo classical values were explained in section 1 and can be chosen next 1 analytical expressions 2 NEOART with one effective impurity species in NEOART with density total density of all ion stages and z mean impurity charge 3 NEOART with each ion stage as a seperate fluid The neo classical transport e
43. ll the input files and denote the parameter that you want to change by the program with a cv in the appropriate input file These parameter are then transferred to STRAHL via the file ext parameter dat and STRAHL is started with option a 3 2 The Result File The result files are written in the Network Common Data Form netCDF The netCDF format is self describing and allows direct access of a small subset of the data The netCDF software is freely available for C C Fortran and other languges It is also included in IDL which I used to write readout and plotting routines Documentation about netCDF can be get from http www unidata ucar edu and the software for various platforms is available via anonymous FTP from the directory ftp ftp unidata ucar edu pub netcdf The fastest way to inspect the output file is the program ncdump which comes with the netCDF software The best is to pipe the result through more or less by entering ncdump filename more Partof the ncdump output is shown here for a run with carbon netcdf strahl_result dimensions grid_points 101 ionisation_stages 7 time UNLIMITED 18 currently variables float time time time units s time _FillValue O time C_format 4g float rho_poloidal_grid grid_points rho_poloidal_grid units 1 rho_poloidal_grid _FillValue 0 f rho_poloidal_grid C_format 4g float impurity_density time ionisation_stages grid_point
44. llowed by a time step II with explicit treatment of ionisation and implicit treatment of recombination 1 no 3 R ei ti 795 202 1 D 1 45 ne S2 nd_1 nl md 5 ni 0 z 1 ID This difference method leads to the solution of a tridiagonal matrix equation for each ion stage on the number of mesh points any b ob nk ni jai d d 1 46 The coefficients b b are on the main diagonal a stands in the lower c in the upper diagonal and d d on the right hand side The coefficients with superscript R are due to ionisation recombination reactions At p a 2p D At a gt a TID p AD pvz At 1 b 1 4pD At 5 pAv tr c 4PD At a an i 2 biz end AtQ 1 47 e acne D neQz 1 ID d At ne S2 nz 1 aend 1 9201 D 1 48 Ne Spin a A n Sz ID For case I the m matrix equations are solved in ascending order of z and the N ni of the next lower z 1 i stage at the new time point is known when needed in the d term of the equation for ni For case ID nl 41 61 needed and the matrix equations are solved in descending order of z At the axis r 0 the density gradient and the drift velocity is zero Expanding density and drift velocity in the vicinity of r 0 yields the transport equation for r 0 On din du TAL 2D 0 Po et 2 poz 0 R Qz 1 49 10 CHAPTER 1 INTRODUCTION The derivative of v at r 0 is replaced by v 1
45. m 3 Tel leV Ti eV yO t 1 pO t 1 pl t 1 p2 t 1 p3 t 1 yO t nt pO t nt pl t nt p2 t nt p3 t nt cv decay lambda t 1 lambda t nt ratfun rational function of the form cv y p yo 1 po 1 pP P po times ELE ern EG function ratfun radial coordinate 2 2 GEOMETRV 17 poloidal rho or volume rho cv ne cm 3 Te eV Tilev y0 t 1 pO t 1 pl t 1 p2 t 1 p3 t 1 vO t nt pO t nt pl t nt p2 t nt p3 t nt cv decay lambda t 1 lambda t nt The parameter pz is not relevant but must be given in the data block The ion temperature T is just needed when calculations of neo classical transport are performed If Tr is not needed or if Ty is equal T the number of time points for the ion temperature must be set to zero and the rest of the profile information can be left out STRAHL assumes equal temperature of electrons and ions in this case Neutral hydrogen density is needed for the calculation of charge exchange CX with thermal neutral hydrogen The inclusion of CX is optional It can be set in Xx atomdat as discussed below When CX is switched on the neutral hydrogen density has to be specified in the pp file The only possible parametrization is quadratic interpolation for the whole radial grid The input block looks like cv times nt EICH tos TELE CV function interpa cv radial coordinate poloida
46. n ls A 1 35 Or 2Ar 15392 The densities y 1 are then replaced by fin 1 2Ar A rin in equations 1 28 1 31 and 1 32 respectively The numerical implementation usually uses a non uniform mesh of radial grid points Close to the last closed flux surface strong gradients of electron temperature and density appear which can have decay lengths as low as a few mm The electron temperature might rise from a few eV to keV within 2 cm and strong gradients in the densities of the different ionisation stages can be present Therefore the radial step size used in the numerical calculations needs to be around 1 mm in the boundary region while in the central part of the plasma a grid resolution around 1 cm is usually sufficient It is advantageous to change the radial grid spacing from centre to edge in order to keep the number of radial mesh points as low as possible The mesh points shall be equidistant in the coordinate 8 CHAPTER 1 INTRODUCTION p p r STRAHL has two options for the grid Both choices give a higher density of grid points at the plasma edge For the first option p r with k gt 1 and the radial distance of two neighbouring grid points is proportional to 1 r 1 The second option is p _ 1 36 Arcentre k 1 Com Ar centre Tedge This choice produces a radial step size Ar which decreases from Arcentre at the axis to Areage at the boundary 1 1 1 T 1 Ar centre Af edge Ar centre Tedge
47. nd the final distribution is saved TIMES TEP S cv number of changes start time stop time 2 cv time dt at start increase of dt after cycle steps per cycle 1 0 1 e 3 1 20 10 2 0 1 e 3 1 20 10 Now we reallv come to the time step definition In the first line vou give the number of time step changes including the start and stop time Then you give as many lines as you specified above The first number is the time at which the time steps shall change in this example it is just the start and the stop time Then comes the time increment At with which the program shall start after the change The last number gives the number of time steps n of one cycle After the completion of a cycle the impurity distribution is stored in a result file Then the time increment is multiplied with a constant number a and in the new cycle we have Atyew Ato a X a The third parameter in the line defines a In this example only the stop time from the last input line is relevant and the rest of the line can be filled with dummies They are only needed when the primitive saw tooth model is switched on see below SO URCE cv position cm constant rate 1 s time dependent rate from file 1000 0 2 5e21 0 cv divertor puff source width in cm source width out cm 0 0 0 The puff rate of neutral impurities Py is given in the next line If is constant in time it is defined by the pa rameter constant rate and the switch under time dependent ha
48. quations apply only to the region inside the LCFS The radial region where neo classical values shall be calculated and used can be set by the last parameter ANOMALOUS TRANS PORT cv f of changes for transport 1 cv time vector 0 0000 cv parallel loss times ms 2 9 cv Diffusion m 2 s hunt cv D D_SOL Dmax beta gamma r_max r LCFS 051 1 0 229 LS 6 0 0 8 cv Drift function Drift Parameter Velocity const c velocitv The anomalous transport can be given as a function of time The code uses the settings of the first time point for times which are before the first time point and the settings of the last time point for times which are after the last time point For the times which are inside the specified time intervall a linear interpolation is used The first line gives the number of time points and the second line lists all times It is followed by the parallel loss times in the region of open field lines Tso L aiy ms The diffusion constant D has unit m s and the drift velocity v has unit m s The drift velocity can be specified by using the drift parameter a In this case the drift velocity is calculated using 2Dr v Q 3 LOFS The diffusion constant and drift velocity can be parametrized in several ways and the upper example is just one possibilty The possible options are e interp causes quadratic interpolation between values at given radial positions The data block looks like cv function only for Drift interp
49. radius e impurity concentrations of each ion stage and and total impurity concentration radial distribution or time development for one radius e AZ of each ion stage total AZ y of impurity and Zepp radial distribution or time development for one radius e emissivity due to total line radiation of each ion stage bremsstrahlung and recombination radiation The total emissivity is also compared with the emissivity which would follow from corona ionisation equilibrium radial distribution or time development for one radius e the same for the spectral range of SXR 3 3 PLOTTING ROUTINES FOR IDL 39 e bremsstrahlung due to impurity and main ion radial distribution or time development for one radius e emissivity due to line radiation of diagnostic lines radial distribution or time development for one radius e total power loss density due to radiation ionisation of impurity radial distribution or time development for one radius e electron density main ion density main ion loss due to impurity radial distribution or time development for one radius e electron temperature main ion temperature radial distribution or time development for one radius e anomalous and neo classical total and parts diffusion coefficent radial distribution or time development for one radius e anomalous and neo classical total and parts drift velocity radial distribution or time development for one radius e anomalous and neo classi
50. s impurity_density units 1 cm 3 impurity_density _FillValue 0 f 3 2 THE RESULT FILE 37 impuritv densitv C format 4g global attributes mass of impuritv 12 f maximum charge 6 data time 1 1 009 1 02 1 033 1 05 1 069 1 093 1 121 1 155 1 196 1 244 15308 L 313 143 15958 15068 e825 2 5 The structure of the different data arrays is very easy to understand from the output of ncdump The dimen sions of the different variables are shown first then comes a block with the variables followed by a block with global attributes and finally there is the long data section The definition of three varibles is shown in this ex ample The first is time with one dimension time and unit s The second is rho poloidal grid with dimension grid points and unit 1 and the third is impurity_density with three dimesions time ionisation_stages grid points and unit cm There are also global attributes defined likemass of impuritvandmaximum charge I think the names and the structure of the variables need no long explanations An alternative way is a public domain IDL program which is called ncbrowse With ncbrowse you can inspect the structure of the variables and make plots When you want to write your own program to read and analyse the data you might use the routine read_strahl pro which runs under IDL With read_strahl pro you can access all data in a conve nient way Th
51. s to be set to zero When setting the switch to 1 a time dependent rate can be specified in the file nete Xxfixnnnnn dat where Xx is the ele ment symbol and nnnnn is the same number as is used for the plasma background the file nete ppnnnmn i In mete Xxfixnnnmn dat the first line contains the number of time points n and in the following n lines you give the times and the according rates The radial position where the neutrals enter the plasma usually is the plasma edge and f source parameter position has to be set to a value which is greater or equal the maximum r It is also possible to choose f souree lt Tmax to simulate pellet ablation Another more convenient option for pellet simula tionare source width inandsource width out Forsource width inorsource width out greater than zero the code assumes a neutral source distribution with a gaussian shape which has the maximum at f source and a full width at half maximum FWHM of source width in to the inside and a FWHM of source width out to the outside For source width in and source width out less than zero the code assumes a neutral source being located at the single grid point which is closest to rsource The op tion divertor puff puffs the neutral impurities in a seperate reservoir the divertor The particles Nq in that reservoir enter the plasma at source at a rate Nav Taiw gt soL Which is specified by the decay time constant Tdiv SOL see the recycling options below EDGE
52. term 1 3 Numerical Solution The impurity transport equations 1 12 for the different ion stages of one impurity species can only be solved numerically Similarly the solution of the transport equation for the sum of all ion stages 1 15 is only known if the transport coefficients have simple radial dependences e g D r const and v r var a or D r Do Da r a and v r var a Therefore computer codes which provide a numerical solution of Eq 1 12 and 1 15 are indispensable A brief description of the preferred numerical schemes shall be given here We use a matrix notation for the Z transport equations 1 12 of all ion stages of an impurity with atomic number Z on on 10 A x gt A gt D gt T a op pe ee an D dD sn 0 dd a a gt i zz an is a vector in the space of of the ionisation stages containing the number densities of the different ionisation stages D and 6 are diagonal matrices with the transport coefficients for each ion stage and R are matrices with the ionisation and recombination rates i e the elements of are the product of n and the rate coefficients for electron impact ionisation while the elements of R are the sum of the products of n and the rate coefficients for 6 CHAPTER 1 INTRODUCTION radiative and di electronic recombination has a non zero diagonal and lower diagonal while R has a non zero diagonal and upper diagonal d has only one non zero element at Z 1 and descri
53. the sources and sinks due to ionisation recombi nation and charge exchange The source and sink term connects neighbouring charge states With an ansatz for the density dependence of the transported flux density ii 1 z the particle conservation equation shall be transformed into an impurity transport equation Here we are interested in the radial transport perpendicular to the magnetic flux surfaces since the large transport coefficients parallel to the magnetic field cause a practically constant density ny z on the magnetic flux surface Even though ny z and also the source sink term 7 z are constant on the flux surface the radial component of r I z shows a variation with the poloidal angle 9 One reason is the poloidally varying distance of the flux surfaces being shorter on the outboard than on the inboard This leads to a poloidal variation of the gradients of density and temperature which are the driving terms for the radial impurity flux Another reason is the 1 R dependence of the toroidal magnetic field Br In Fig 1 1 the poloidal projection of a few flux surfaces are sketched For one flux surface the dependence of the gradient and of Br on are shown in a polar plot Both quantities vary by about a factor of 2 The correct flux surface average of Eq 1 1 is derived by the following steps see Ref Hinton and Hazeltine 1976 Inside the last closed flux surface LCFS Eq 1 1 is integrated over 2 CHAPTER 1 INTRODUCTION Figure 1
54. tion and emission In this chapter everv input file is described All path names directory names are relative to the directory where STRAHL is executed which is called the current directory The needed files are located in the current directory or in sub directories of the current directory sub directorv All input files are in ASCII format 2 1 Plasma Background The parameters for the plasma background are in the directory nete The filename is nete ppnnnnn i where nnnnn is a 5 digit number shot number and iis a 1 digit number index The program reads all parameters with a free format read statement in Fortran77 read unit Thus it is not necessary to obey strict formatting rules Character constants have to be given within single quotes e g text Each item in the data block is preceded by a line which starts with cv and the programm scans through the input file to search for the next appearance of cv When it finds the next cv it starts to read the data values from the next line The advantage of this method is the possibility to include descriptions in the input file which help to understand the structure or memorize changes Old input data can become comments by replacing cv by anything else in the preceding line When replacing cv_ by cv the next data block is read from the input file with fixed name ext_parameter dat in the current directory This alternative input might simplify scans of different STRAHL runs
55. uses the profile of the first time point for times which are before the first time point and the profile of the last time point for times which are after the last time point For the times which are inside the specified time intervall a linear interpolation in time is used The different types of possible parametrizations are identical for the electron density the electron temperature and the ion temperature In the radial range outside the last closed flux surface LCFS the profile is just described by an exponential decay length A unit cm The profile outside the LCFS with radius r Lcrs is then given by T FLORS y r y LCFS exp For the profiles inside the LCFS four different functions are accepted Each parametrization can be used with ppot or r for the radial coordinate The formulas are given here with the dummy p Each profile contains a scale factor Yo which transforms the mathematical function into a physical qunatity The scale factor for the electron density has the unit cm7 while the scale factor for the temperatures has the unit eV e interp causes quadratic interpolation between values at given radial positions The data block looks like 2 1 PLASMA BACKGROUND CV cv cv cv cv cv cv times HL ee EME function interp radial coordinate poloidal rho or volume rho of interpolation points n_rho radial grid for ne interpolation rho 1 rho n rho ne cm
56. verwrites the old results e If you ask for a backup of the result option S STRAHL invokes a system call to copy re sult strahl_result dat to a file with name result XXnnnnntx xx_y yy_i XX is the element symbol nnnnn is the number of the file for the background plasma x xx and y yy give the time intervall and i is an index you enter at run time For times larger the 10s the format xx x is used There are several options which can be supplied with the call of STRAHL e no option strahl All input parameters you enter at run time are stored in the file strahl control e strahl a All input parameters are read from strahl control e strahl n No radiation results are calculated 36 CHAPTER 3 EXECUTION AND OUTPUT trahl d nly the impurity distribution and the diagnostic line radiation is calculated ou strahl s Only the total impurity density and the soft X ray radiation is saved e strahl w Some atomic data are saved in directoy tmp e strahl v STRAHL produces a lot of output on the screen This option might be helpful in the case of runtime errors or to detect errors in the input files The best way to use this option is to combine it with the a option and enter strahl a v gt strahl trace Now you can search the error in strahl trace e strahl h All options are displayed If you want to include STRAHL into other programs e g to perform scans or fitting you might use the following scheme Setup a
57. vz 9 Ar and the discretisation is nsi Ne o 2A nd Mra De no an Yn o AtR AtQ 1 50 The coefficients are a 0 b 1 2A 5 D c 215 d 2 b nd end AtQ 1 51 At the outermost radial mesh point k the density shall decay with decay length A dn Nz k Se 0 1 52 dr en The boundarv condition defines the density at the virtual mesh point k 1 in terms of the densities at k and k 1 dn Nz Nlz k dr P N k 1 nz k 1 TA Nak Mak l EST 1 53 Eq 1 44 is modified for the outermost grid point k by replacing the densities at k 1 with Eq 1 53 Furthermore derivatives of D and v are neglected At 1 hex n 8p D 1 ni nk N fi p1 2 2pr Atl p ya TID poz nz nz At 1 FEA tbs AR AtQ 1 54 II The coefficients are a 4pD At 1 At 1 p At v 1 b 14 4pD At ED ll A 4p D nt mar z Aa ea c 0 and pi 2 b n y AtQ 1 55 For the hard boundary case nz edge 0 the coefficients are a 0 bi 1 c 0 d d 1 1 56 1 5 NEUTRAL IMPURITIES 11 1 5 Neutral Impurities The transport equations 1 12 are only solved for the ions while the neutrals only act as source for the first ionisation stage The neutrals are assumed to enter the plasma with a given uniform radial speed vo and the total number of neutrals entering the plasma per unit time shall be given by B 4r Rr vpng r 1 57 Due
58. y Habilitationsschrift Hinton and Hazeltine 1976 Hinton F L and Hazeltine R D 1976 Theory of plasma transport Reviews of Modern Physics 48 2 239 308 Hirshman and Sigmar 1981 Hirshman S P and Sigmar D J 1981 Tokamak impurity transport Nuclear Fusion 21 9 1079 1201 Houlberg et al 1997 Houlberg W A Shaing K C Hirshman S P and Zarnstorff M C 1997 Bootstrap current and neoclassical transport in tokamaks of arbitrary collisionality and aspect ratio Physics of Plasmas 4 9 3230 3241 Kim et al 1991 Kim Y B Diamond P H and Groebner R J 1991 Neoclassical poloidal and toroidal rotation in tokamaks Physics of Fluids B Plasma Physics 3 8 2050 2060 Lackner et al 1982 Lackner K Behringer K Engelhardt W and Wunderlich R 1982 An algorithm for the solution of impurity diffusion under finite reaction rates Zeitschrift fiir Naturforschung 37a 5 93 1 938 Peeters 2000 Peeters A G 2000 Reduced charge state equaitons that describe pfirsch schl ter impurity transport in tokamak plasma Physics of Plasmas 7 1 268 275 Thompson et al 1985 Thompson J F Warsi Z U A and Mastin C W 1985 Numerical Grid Generation Foundations and Applications North Holland New York 1st edition 43
Download Pdf Manuals
Related Search