subroutine ip_chi(RLT,RLN,RLNe,q,kappa, . shat,zth,nbeam,tau,eps,gnu,beta, . g_perp,rmajor,rho_i,v_ti,rho_e,v_te,etg,rlte, . RLTcrit,RLTcritz,chi_0,gfac,chi_i,chi_e,gamma) c c The formulas embodied in this subroutine are documented in the Physics c of Plasmas article entitled "Quantitative Predictions of Tokamak Energy c Confinement from First-Principles Simulations with Kinetic Effects", c by M. Kotschenreuther, W. Dorland, M.A. Beer, and G.W. Hammett, c Vol. 2, p. 2381, (1995). Extensions to non-circular cross-sections are c described below. c c There is a significant typographical error in that paper. R/Ln* is c defined to be max(6,R/Ln); it should be min(6,R/Ln) as defined in c this subroutine. c c Also, note that in deriving these formulas, we assumed that the c density gradient scale lengths for the different species were equal. c This is an approximation that needs to be relaxed. c c As emphasized in the paper, these formulas were derived numerically and c are therefore not trustworthy outside a particular region of parameter c space. For example, we did not parameterize the heat flux in the weak c magnetic shear limit; thus, one should not use the model in this limit. c I have attempted to reduce related strange numerical behaviors by c limiting some inputs to be roughly within their range of validity. c c Questions, problems, errors, etc. should be reported to c bdorland@zonker.ph.utexas.edu or bdorland@pppl.gov. c c I will reply as quickly as possible. c c Stiffness: c ********** c For many cases that we have simulated, the transport equations c tend to be very stiff. That is, the plasma temperature gradient scale c length tends to adjust itself to be close to the critical gradient scale c length over some region of the plasma, because chi becomes very large c very fast for shorter temperature gradient scale lengths. Typically, c we have had to be very careful with the numerical algorithm used in the c transport equation solver with this experience in mind. The details c of our implementation are available to anyone that is interested. c c Geometry: c ******** c c The nonlinear simulations that were done to obtain these formulas were c mostly done in a simplified geometry, using a shifted circle, low beta, c high aspect ratio expansion. Some modifications due to more sophisticated c geometrical models have been calculated and have been included here, c but should be considered preliminary. There are two important issues c that must be noted. First, we derived our formulas using a different c radial coordinate. Second, since we are actually calculating the c transport coefficients in general geometry, we require less assumptions c for the form of the transport equation to be solved. c c Let me describe the <|grad rho|> issue first: c c The database standard modeling assumptions that were agreed upon for c this exercise include the assumption that the anomalous fluxes for c non-circular cross-sections are simply related to the anomalous c fluxes for related circular cross-section plasmas. That is, in order c to get the factor of < |grad rho|**2 > that appears as a coefficient c of chi in the energy transport equations, one assumes that c c chi_anom_general = chi_anom_circular * (grad rho). c c One need not make this assumption; one can just calculate the quantity c chi_anom_general directly. One would then have a transport equation c of the form c c (3/2) d(n T)/dt = (1/V') d/drho [V' <|grad rho|> n chi d/drho(T)] + ... c c in which (grad rho) appears to the first power, rather than the second, c and chi is the thermal diffusivity from a general geometry theory. c c This is arguably the better way to proceed, since c c V' <|grad rho|> = A, c c where A is the surface area. In this form, the quantity c c -n chi dT/drho c c can be identified as the heat flux per unit area, a natural c quantity from a theoretical/simulation point of view. This chi is c the quantity returned by this subroutine. c c If you are solving the transport equations in the form that c the ITER Expert Group agreed upon, e.g., c c (3/2) d(n T)/dt = (1/V') d/drho [V' <|grad rho|**2> n chi d/drho(T)] + ... c c then you need to multiply the chi_i and chi_e reported by this c subroutine by the factor <|grad rho|>/<|grad rho|**2>. This should c result in only small corrections to the predicted profiles. c c The choice of radial coordinate is more difficult to resolve: c c We did not use the sqrt(toroidal flux) radial coordinate in our c non-circular cross-section simulations. Instead, we used "rho_d", c where rho_d is defined to be the average horizontal minor radius c at the elevation of the magnetic axis, normalized to the value of c this quantity at the LCFS. c c In other words, denote by "d" the horizontal diameter of a given flux c surface measured at the elevation of the magnetic axis. Denote by c "D" the horizontal diameter of the last closed flux surface at the c elevation of the magnetic axis. Then rho_d = d/D. I believe this c is variable number 67 (RMINOR) in the ITER Profile Database Standard c List. c c It is not difficult to allow for an arbitrary radial coordinate c in one's transport code. One must obtain all of the radial c quantities as functions of rho_d rather than rho via interpolation. c c However, I do not expect everyone to go to this length to test our c model, since you agreed to use the sqrt(toroidal flux) definition c of the radial coordinate. Thus, I suggest the following alternative: c Simply use the rho_d coordinate to define the scale lengths that c appear in the formulas below. For most quantities (such as R/LT), c this simply amounts to including an additional factor d rho/d rho_d c in the expressions passed to this subroutine. While not completely c correct, this workaround captures the dominant effect, related to c evaluating the flux near the critical gradient. c c ****** Summary of comments: *********** c c (1) The general geometry extensions to the IFS/PPPL model were derived c using rho = d/D = RMINOR as the radial coordinate. To be c most accurate, the transport equation should be solved using d/D c as the radial coordinate. c c If you use rho proportional to sqrt(toroidal flux) instead of rho=d/D c as your radial coordinate, you should at least carefully define the c scale lengths as indicated below (using rho=d/D). c c (2) This routine should be used to return the thermal transport c coefficients (chi_i, chi_e) for energy transport equations of the form c c (3/2) d(n T)/dt = (1/V') d/drho [V' <|grad rho|> n chi d/drho(T)] + ... c c Note that <|grad rho|> only appears to the first power according c to this definition of chi. If your code is hardwired to solve an c equation of the form c c (3/2) d(n T)/dt = (1/V') d/drho [V' <|grad rho|**2> n chi d/drho(T)] + ... c c then multiply the chi_i and chi_e obtained from this routine by c the factor <|grad rho|>/<|grad rho|**2>. c c ***************************************************************** c c RLT R/L_Ti, where R = the major radius and c 1/L_Ti = -1/T_i dT_i/drho_d c RLN R/L_ni, where 1/L_ni = -1/n_i dn_i/drho_d c RLNe R/L_ne, where 1/L_ne = -1/n_e dn_e/drho_d c q The safety factor. c kappa The elongation, defined here to be the c kappa = maximum height/maximum diameter c shat == rho_d/q dq/drho_d c zth Thermal Z_eff. The simulations that were carried out to c generate the formulae in this subroutine assumed the plasma c was composed of a thermal hydrogenic species, thermal carbon, c a hydrogenic beam species, and electrons. We found that low-Z c impurities primarily act to dilute the main ion concentration, c and can accounted for to first order by modifying the c definition of Z_eff. Some of the more important effects of c the fast ions in the plasma are also partially accounted c for by this parameter, which is: c zth == (n_i + 36 n_C)/(n_e - n_beam) c c nbeam == local fast ion (beam) density normalized to the electron c density. c tau == T_i/T_e. Note that this is opposite to a widely c used convention. c eps == rho_d/R, the local minor radius normalized to the major radius. c gnu Dimensionless collisionality parameter. c gnu == 2.5e-7 * n_e / (T_e**1.5 T_i**0.5) * rmajor c where n_e is in units of cm**-3, T_e and T_i are in eV, and c rmajor is in units of m. For an R = 2.4 m, 100 eV, 1e13 plasma, c gnu=600. c beta Local total beta; presently ignored. c g_perp == velocity shear parameter. Use g_perp=0 (actual value c discussed in the Waltz paper cited below). c rmajor == major radius of the plasma c rho_i == local thermal gyroradius of thermal hydrogenic species. c v_ti == sqrt(T_i/m_i) where T_i and m_i are the local thermal c hydrogenic temperature and average thermal hydrogenic mass. c c rho_e,v_te,etg,rlte: dummy parameters at present. Ignore. c c Units: The only dimensional parameters in the inputs are the major c radius, rho_i, and v_t. Their units should be consistent; the chi's c that are returned will be in units of rho_i**2 v_ti / rmajor. c c OUTPUT: c RLTcrit: R/L_Tcrit for ITG mode c RLTcritz: R/L_Tcrit for carbon branch c chi_0: normalized chi (ignore) c gfac: L_Tc/L_T, where L_Tc is the critical temperature gradient c scale length for the deuterium branch of the ITG mode. c chi_i: Anomalous ion thermal diffusivity from toroidal ITG mode. c chi_e: Anomalous electron thermal diffusivity from toroidal ITG mode. c c This parameterization of chi is not complete. There are significant c neglected physical processes that are known to be important in c many operational regimes. c c The most significant problems are: c c (1) Trapped ion/long wavelength ITG modes. These modes are known c to be unstable for typical edge tokamak parameters. However, until c we have nonlinear estimates of the associated thermal diffusivity, c these modes are ignored, leading to overly optimistic predictions of c edge thermal confinement. c (2) Trapped electron modes, which can alter the stability boundary c significantly for low collisionality. At high collisionality these c modes are generally stable and thus largely irrelevant. When present, c they are associated most strongly with particle transport, although c there is also an associated heat transport. c (3) Minority ion density gradients, which can strongly change chi c and LT_crit. c (4) Sheared flows, which are stabilizing. This includes diamagnetic c and ExB shear flows. c (5) Finite beta effects, generally stabilizing. c implicit none real RLT,RLN,RLNe,shat,zth,tau,rmajor,chi_i,q, . nbeam,taub,rho_i,v_ti,eps,chi_e,nu,chi_0,g_perp,gnu real RLTcrit,RLTcritz,RLTcritg,chi0,f_0,f_z,chiz,esqrt real trln,trlne,tshat,kappa,rho_e,v_te,rlte,beta real chie1,chie2,gamma,tgamma,rot,g_fac1,g_facz,gfac real etg_rltecrit,f_e real reduced logical etg real RLT_mult external reduced data RLT_mult/1.0/ ! multiplier on R/L_Tcrit, for sensitivity studies taub=tau/(1.-nbeam) nu=gnu*0.84 tRLN=min(abs(RLN),6.0)*sign(1.0,RLN) tRLNe=min(abs(RLNe),6.0)*sign(1.0,RLNe) c c Formula is not applicable for shat<0.5; does not matter most of the time, c since low shear regions of the plasma tend to be in the region affected c by sawteeth. We will include the effects of weak or negative c magnetic shear in a future release. c tshat=min(max(shat,0.25),2.5) c C critical ion temperature gradient: C (If this number is negative, check the input variable definitions, c especially gnu. Then ask me about it, at c bdorland@hagar.ph.utexas.edu) c if(q.le.0.) write(*,*) 'q= ',q if(taub.le.0.) write(*,*) 'taub= ',taub if(eps.le.0.) write(*,*) 'eps= ',eps if(zth.le.0.) write(*,*) 'zrh= ',zth if(nu.le.0.) write(*,*) 'nu= ',nu RLTcrit=2.46*(1.+2.78/q**2)**0.26*(zth/2.)**0.7*taub**0.52 . *( (0.671+0.570*tshat-0.189*tRLNe)**2 . +0.335*tRLNe+0.392-0.779*tshat+0.210*tshat**2) . *( 1.-0.942*(2.95*eps**1.257/nu**0.235-0.2126) . *zth**0.516/abs(tshat)**0.671) RLTcrit=RLTcrit*RLT_mult c if(RLTcrit.le.0) then c write(*,*) 'rltcrit,q,zth,taub,tshat,trlne,eps,nu' c write(*,*) rltcrit,q,zth,taub,tshat,trlne,eps,nu c write(*,*) 1.-0.942*(2.95*eps**1.257/nu**0.235-0.2126) c . *zth**0.516/abs(tshat)**0.671 c write(*,*) (0.671+0.570*tshat-0.189*tRLNe)**2 c . +0.335*tRLNe+0.392-0.779*tshat+0.210*tshat**2 c endif RLTcritz=0.75*(1.+taub)*(1.+tshat)*max(1.,3.-2.*tRLNe/3.)* . (1.+6.*max(0.,2.9-zth)) RLTcritz=RLTcritz*RLT_mult f_0=11.8 . *min(1.,(3./zth)**1.8) . *q**1.13 . /(1.+tshat**0.84) . /taub**1.07 . *(1.+6.72*eps/nu**0.26/q**0.96) . /(1.+((kappa-1.)*q/3.6)**2) f_z=7.88 . /(1.+tshat) . *max(0.25,zth-3) . /taub**0.8 . /(1.+((kappa-1.)*q/3.6)**2) chi0=f_0 * rho_i**2*v_ti/rmajor chiz=f_z * rho_i**2*v_ti/rmajor g_fac1=max(0.,min(esqrt(RLT-RLTcrit),RLT-RLTcrit)) g_facz=max(0.,min(esqrt(RLT-RLTcritz),RLT-RLTcritz)) chi_i=max(chi0*g_fac1,chiz*g_facz) gfac=RLT/RLTcrit chi_0=chi0*sqrt(abs(RLTcrit)) chie1=chi0*g_fac1*1.44*tau**0.4*(q/tshat)**0.3*nu**0.14 . *max(0.16667,eps) chie1=0.5*chie1*(1.+tRLNe/3.) chie2=0.5*max(2.,(1.+RLNe/3.))*0.526*tau*nu**0.22*chiz*g_facz c c Correction for n_i/n_e and ratio of heat fluxes rather than chi's: c chi_e=max(chie1,chie2)*(7.-zth)/6. c if(chi_i.gt.0.) then c write(*,*) chi_i/chi_e,eps,(q/tshat)**0.3 c endif c c **Preliminary** model of shear flow stabilization. Based on c work described in Waltz, et al., Phys. of Plasmas, Vol. 1, p. 2229 c (1992), and in Hahm and Burrell, Phys. Plasmas 2, 1648 (1995). c The fundamental concept of sheared-flow stabilization of turbulence c in plasmas originated in the seminal paper by c Biglari, Diamond, Terry, Phys. Plasmas B2, 1 (1990), which was c particularly fundamental in understanding H-mode experiments. c [There are some roots to this concept in general fluid turbulence c also, such as the shearing rate or eddy turn-over time in Orszag's c EDQNM, general concepts of Kolmogorov scaling, etc. But there are c also important differences such as the 2-D vs. 3-D nature of the c plasma vs. fluid nonlinearities, and the existence of the c Kelvin-Helmholtz instability in sheared flows in normal fluids. c This instability is different in magnetized plasmas.) c c gperp is essentially the shearing rate d/dr(v_perp), and approximate c stabilization is assumed when the shearing rate is of order the linear c growth rate (note that gperp has units of 1/sec, i.e. the units of a c growth rate). In the absence of full-torus nonlinear simulations, c gperp should include the shear in the perpendicular phase velocity of c the plasma turbulence waves, but we will assume this is primarily the c shear in ExB flow for now, and ignore shear in omega/k_theta. Also, c we will ignore the Kelvin-Helmholtz destabilization which Waltz points c out can occur due to parallel flow shear. In the limit where v_parallel c combines with v_perp to make the net flow purely toroidal (a common c assumption), gperp can be written as (Waltz p. 2235): c c gperp=(r/q) d/dr(v_phi/R) c c Waltz didn't write down a more general form, but the ballooning c coordinates he uses implies that it is the shear in the perpendicular c flow which matters. For general geometry and general flows, the above c expression still applies if r/q is interpreted as the local value of c R*B_theta/B, if d/dr is interpreted as a local gradient, and c and if v_phi/R is replaced by (B/B_theta)*v_ExB/R. [This may not be c obvious to some, but it is a consequence of the field-line-following c ballooning coordinates we use, which involve a parallel coordinate and c a toroidal angle coordinate, i.e., the toroidal angle (at fixed parallel c coordinate) is the effective perpendicular coordinate. This is related c to the way one changes the "slab" definition of omega_* into a general c geometry definition involving k_phi=n/R instead of k_theta, using m=n*q c etc, and is related to the fact that the properly expressed omega_* is a c flux function.] Note that v_phi/R and (B/B_theta)*v_ExB/R are both flux c functions, i.e., constant on a flux surface. c c Hahm and Burrell derive an expression directly in general geometry, c using a more formal 2-pt correlation function theory. It is essentially c equivalent to the above expression, if one makes the correct assumptions c above implied by the ballooning coordinates, and if one reduces the Hahm c expression for practical experimental comparisons by assuming the radial c and poloidal correlation lengths are equal. In flux coordinates c (dpsi = R*B_theta*dr), and in CGS units (where v_ExB = c*E_r/B), c this can be written as c c gperp = (R**2*B_theta**2/B) c d**2(Phi)/dpsi**2 c c where Phi(psi) is the electrostatic potential. Note that the factor out c front gives significant poloidal variation, particularly when the c Shafranov shift is large so that B_theta is compressed. This is good, c since theta=0 is where the bad-curvature drive of the toroidal ITG mode c is worst. Thus gperp should be evaluated at theta=0. c c Rough growth rate parameterization. Crudely assume the form of the c growth rate for the carbon and deuterium branches is the same, with c just a different threshold: RLTcritg=RLTcrit if (chiz*g_facz .gt. chi0*g_fac1) RLTcritg=RLTcritz gamma=0.25/(1+0.5*tshat**2)/tau . *(1+3.*max(0.,eps-0.16667)) . /(1+max(0.,q-3)/15.) . *max(0.0,RLT-RLTcritg)*v_ti/rmajor tgamma=max(gamma,1.e-8*v_ti/rmajor) rot=1.-abs(g_perp)/tgamma rot=min(1.,max(0.0,rot)) chi_i=max(0.,chi_i*rot) chi_e=max(0.,chi_e*rot) if(etg) then etg_rltecrit=2.184*sqrt(.5+1./q)*(zth/tau)**.61 . *(1.-0.85*eps/tshat**.25) . *(abs(.1976-.4550*tshat+.1616*RLNe)**.769 . +.7813+.2762*tshat+.3967*tshat**2) etg_rltecrit=etg_rltecrit*RLT_mult f_e=14.*(1.+eps/3.)*(1+0.1*tau)*tau*q/(2+tshat) if(RLTe.gt.etg_rltecrit) then chi_e=chi_e+f_e*rho_e**2*v_te/rmajor . *min((RLTe-etg_RLTecrit)**0.5,(RLTe-etg_RLTecrit)) endif endif c chi_i=reduced(eps) c chi_e=reduced(eps) return end real function esqrt(x) real x esqrt=0. if(x.gt.0.) esqrt=sqrt(x) return end real function reduced(chi) reduced=chi return end