jbootrap_th_cyl_yp

 
 Calculates the bootstrap current from neoclassical models by Hinton & Hazeltine, Hirshman and Sauter 
 in the limit of large aspect ratio ia << 1.  
 Note: The calculation is based on the model of circular concentric flux-surface, with parameters derived from the magnetic geometry.

   INPUTS: 

       - tokname: Name of the tokamak

   OUTPUTS:

       - prho: radial grid [1,np]
       - pcurr_hir_ref_fsav: Electron bootstrap current from Hirshman model, low collisionality limit (MA/m2) [1,np]
       - pcurr_sau_ref_fsav: Electron bootstrap current from Sauter model (MA/m2) [1,np]
       - pcurr_hin_ref_fsav: Electron bootstrap current from Hinton & Hazeltine model (MA/m2) [1,np]
       - pcurr_i_ref_fsav: Ion contribution to total bootstrap current for Hirschman, Sauter and Hinton & Hazeltine models (MA/m2) [1,np]
       - pcurr_tag_e_ref_fsav: Electron bootstrap current from Taguchi model (anistropic electron distribution) (MA/m2) [1,np] 
       - pcurr_tag_i_ref_fsav: Ion bootstrap current from Taguchi model (anistropic electron distribution) (MA/m2) [1,np] 

by Y. Peysson <yves.peysson@cea.fr> and Joan Decker <joan.decker@epfl.ch>

0001 function [prho,pcurr_hir_ref_fsav,pcurr_sau_ref_fsav,pcurr_hin_ref_fsav,pcurr_i_ref_fsav,L31] = jbootrap_th_cyl_yp(equil)
0002 %
0003 % Calculates the bootstrap current from neoclassical models by Hinton & Hazeltine, Hirshman and Sauter
0004 % in the limit of large aspect ratio ia << 1.
0005 % Note: The calculation is based on the model of circular concentric flux-surface, with parameters derived from the magnetic geometry.
0006 %
0007 %   INPUTS:
0008 %
0009 %       - tokname: Name of the tokamak
0010 %
0011 %   OUTPUTS:
0012 %
0013 %       - prho: radial grid [1,np]
0014 %       - pcurr_hir_ref_fsav: Electron bootstrap current from Hirshman model, low collisionality limit (MA/m2) [1,np]
0015 %       - pcurr_sau_ref_fsav: Electron bootstrap current from Sauter model (MA/m2) [1,np]
0016 %       - pcurr_hin_ref_fsav: Electron bootstrap current from Hinton & Hazeltine model (MA/m2) [1,np]
0017 %       - pcurr_i_ref_fsav: Ion contribution to total bootstrap current for Hirschman, Sauter and Hinton & Hazeltine models (MA/m2) [1,np]
0018 %       - pcurr_tag_e_ref_fsav: Electron bootstrap current from Taguchi model (anistropic electron distribution) (MA/m2) [1,np]
0019 %       - pcurr_tag_i_ref_fsav: Ion bootstrap current from Taguchi model (anistropic electron distribution) (MA/m2) [1,np]
0020 %
0021 %by Y. Peysson <yves.peysson@cea.fr> and Joan Decker <joan.decker@epfl.ch>
0022 %
0023 [qe,me,mp,mn,e0,mu0,re,mc2,clum,alpha] = pc_dke_yp;%Universal physics constants
0024 %
0025 zZi = equil.zZi;
0026 zmi = equil.zmi;
0027 theta = equil.theta;
0028 pTe = equil.pTe;
0029 pne = equil.pne;
0030 pzTi = equil.pzTi;
0031 pzni = equil.pzni;
0032 ptx = equil.ptx;
0033 pty = equil.pty;
0034 psi_apRp = equil.psi_apRp;
0035 ptBx = equil.ptBx;
0036 ptBy = equil.ptBy;
0037 BPHI = equil.ptBPHI;
0038 Rp = equil.Rp;
0039 Zp = equil.Zp;
0040 %
0041 Bt = BPHI(1,1);
0042 ap = ptx(end,1);
0043 %
0044 R = Rp + ptx;
0045 Z = Zp + pty;
0046 ppsi = psi_apRp*Rp*ap;
0047 %
0048 BR = ptBx;
0049 BZ = ptBy;
0050 %
0051 [npsi,ntheta] = size(R);
0052 np = length(ppsi);
0053 %
0054 % FFT interpolation on R,Z, BR, BZ, BPHI as a function of theta (very precise, works for uniform theta grids only)
0055 %
0056 mfactor = 10;
0057 mtheta = 0:2*pi/((ntheta-1)*mfactor):2*pi*((ntheta-1)*mfactor-1)/((ntheta-1)*mfactor); %super theta grid
0058 %
0059 mR = interpft(R(:,1:(ntheta-1))',(ntheta-1)*mfactor)'; %FFT interpolation
0060 mZ = interpft(Z(:,1:(ntheta-1))',(ntheta-1)*mfactor)'; %FFT interpolation
0061 mBR = interpft(BR(:,1:(ntheta-1))',(ntheta-1)*mfactor)'; %FFT interpolation
0062 mBZ = interpft(BZ(:,1:(ntheta-1))',(ntheta-1)*mfactor)'; %FFT interpolation
0063 mBPHI = interpft(BPHI(:,1:(ntheta-1))',(ntheta-1)*mfactor)'; %FFT interpolation
0064 %
0065 mtheta = [mtheta,2*pi]; %2*pi point retrieved (for future interpolation)
0066 mR = [mR,mR(:,1)]; %2*pi point retrieved
0067 mZ = [mZ,mZ(:,1)]; %2*pi point retrieved
0068 mBR = [mBR,mBR(:,1)]; %2*pi point retrieved
0069 mBZ = [mBZ,mBZ(:,1)]; %2*pi point retrieved
0070 mBPHI = [mBPHI,mBPHI(:,1)]; %2*pi point retrieved
0071 %
0072 nmtheta = length(mtheta);
0073 %
0074 Rp = mR(1,1);%Major radius of magnetic Axis
0075 Zp = mZ(1,1);%Elevation of magnetic Axis
0076 %
0077 Rmax = mR(npsi,1);%Major radius where the last close flux-surface crosses the LFS midplane
0078 ap =  Rmax - Rp;%Plasma minor radius on the LFS midplane
0079 prho = (R(:,1) - Rp)'/ap;
0080 %
0081 psia = ppsi(npsi);%Poloidal flux at the edge
0082 ppsin = ppsi/psia;%Normalized poloidal flux
0083 %
0084 pR = mR;
0085 pZ = mZ;
0086 pBR = mBR;
0087 pBZ = mBZ;
0088 pBphi = mBPHI;
0089 pBT = abs(pBphi);%Toroidal B-field
0090 pBp = sqrt(pBR.^2 + pBZ.^2);%Poloidal B-field
0091 pB = sqrt(pBphi.^2 + pBR.^2 + pBZ.^2);%Total B-field
0092 %
0093 %index 0 refers to the minimum of the magnetic field
0094 %
0095 [pB0,piB0] = min(pB'); %Minimum of B on flux surface
0096 pBmax = max(pB'); %Maximum of B on flux surface
0097 %
0098 pmhu0T2 = 1 - pB0./pBmax; %Pitch-angle (squared) at the trapped-passing boundary
0099 %
0100 piB0 = (piB0 - 1)*npsi + (1:npsi);%Indices in (psi,theta) grid corresponding of B=Bmin
0101 pBT0 = pBT(piB0);%Toroidal field at B=Bmin
0102 pBp0 = pBp(piB0);%Poroidal field at B=Bmin
0103 pR0 = pR(piB0);%Major radius at B=Bmin
0104 %
0105 pBp0(pBp0==0) = eps;%To avoid singularities
0106 %
0107 % note: model with one ion species only. to be completed
0108 %
0109 pTi = pzTi(1,:);
0110 Zi = zZi(1);
0111 %
0112 drho = diff(prho);
0113 %
0114 dTedr = diff(pTe)./drho/ap;
0115 dTedr = [dTedr(1),(dTedr(1:np-2) + dTedr(2:np-1))/2,dTedr(np-1)]; 
0116 dnedr = diff(pne)./drho/ap;
0117 dnedr = [dnedr(1),(dnedr(1:np-2) + dnedr(2:np-1))/2,dnedr(np-1)]; 
0118 dTidr = diff(pTi)./drho/ap;
0119 dTidr = [dTidr(1),(dTidr(1:np-2) + dTidr(2:np-1))/2,dTidr(np-1)]; 
0120 %
0121 dlnTedr = 1./(pTe+eps).*dTedr;
0122 dlnnedr = 1./(pne+eps).*dnedr;
0123 dlnTidr = 1./(pTi+eps).*dTidr;
0124 dlnnidr = dlnnedr;
0125 %
0126 ppe = pne.*pTe;
0127 dlnpedr = dlnTedr + dlnnedr;
0128 dlnpidr = dlnTidr + dlnnidr;
0129 %
0130 pia = prho*ap/Rp;
0131 %
0132 px = (1.46*pia.^(0.5) + 2.40*pia)./(1 - pia).^(1.5);%Small inverse aspect ratio expansion
0133 %[pft,pftu,pftl] = ft_llm_dke_yp(pia,1000);
0134 %px = pft./(1-pft);
0135 %
0136 Dx = 1.414*Zi + Zi^2 + px*(0.754 + 2.657*Zi + 2*Zi^2) + px.^2*(0.348 + 1.243*Zi + Zi^2);
0137 %
0138 L31 = px.*(0.754 + 2.21*Zi + Zi^2 + px*(0.348 + 1.243*Zi + Zi^2))./Dx;
0139 L32 = -px*(0.884 + 2.074*Zi)./Dx;
0140 L34 = L31;
0141 alpha = -1.172./(1 + 0.462*Zi);
0142 %
0143 %From Sauter paper
0144 %
0145 %pft
0146 %pnhu_e_s
0147 %X = pft./(1 + (1-0.1*pft).*sqrt(pnhu_e_s) + 0.5*(1-pft).*pnhu_e_s/Zi);
0148 %L31_s = (1 + 1.4/(Zi+1))*X - 1.9/(Zi+1)*X.^2 + 0.3/(Zi+1)*X.^3 + 0.2/(Zi+1)*X.^4;
0149 %X = pft./(1 + 0.26*(1-pft).*sqrt(pnhu_e_s) + 0.18*(1-0.37*pft).*pnhu_e_s/sqrt(Zi));
0150 %F32_ee_s = (X-X.^4)*(0.05+0.62*Zi)/(Zi+0.44*Zi^2) + (X.^2-X.^4 - 1.2*(X.^3-X.^4))/(0.22*Zi+1) + 1.2/(0.5*Zi+1)*X.^4;
0151 %Y = pft./(1 + (1+0.6*pft).*sqrt(pnhu_e_s) + 0.85*(1-0.37*pft).*pnhu_e_s/(1+Zi));
0152 %F32_ei_s = -(Y-Y.^4)*(0.56+1.93*Zi)/(Zi+0.44*Zi^2) + 4.95*(Y.^2-Y.^4 - 0.55*(Y.^3-Y.^4))/(2.48*Zi+1) - 1.2/(0.5*Zi+1)*Y.^4;
0153 %L32_s = F32_ee_s + F32_ei_s;
0154 %X = pft./(1 + (1-0.1*pft).*sqrt(pnhu_e_s) + 0.5*(1-0.5*pft).*pnhu_e_s/Zi);
0155 %L34_s = (1 + 1.4/(Zi+1))*X - 1.9/(Zi+1)*X.^2 + 0.3/(Zi+1)*X.^3 + 0.2/(Zi+1)*X.^4;
0156 %
0157 % Sauter Correction for L32 in case of large aspect ratio limit (S.Schultz thesis p 91-92)
0158 %
0159 px_s = 1 - (1 - pia).^2./sqrt(1-pia.^2)./(1 + 1.46*sqrt(pia));
0160 
0161 %L31_s =
0162 
0163 L32_s = -(0.51 + 1.31*Zi)/Zi/(1 + 0.44*Zi)*(px_s - px_s.^4)...
0164         +(5.95 + 3.57*Zi)/(1 + 2.70*Zi + 0.546*Zi^2)*(px_s.^2 - px_s.^4)...
0165         -(3.92 + 3.57*Zi)/(1 + 2.70*Zi + 0.546*Zi^2)*(px_s.^3 - px_s.^4);
0166 %L34_s = L31_s;
0167 %
0168 A1 = dlnpedr + pTi./(Zi*pTe).*dlnpidr;
0169 A2 = dlnTedr;
0170 A4 = alpha.*pTi./(Zi*pTe).*dlnTidr;
0171 %
0172 Jbavnorm = -ppe.*pBT0./pBp0.*(L31.*A1 + L32.*A2 + L34.*A4);%WARNING: pBT0./pBp0 not equal to pq./pia when pia becomes large (i.e. > 0.2)
0173 Jbavnorm_s = -ppe.*pBT0./pBp0.*(L31_s.*A1 + L32_s.*A2 + L34_s.*A4); 
0174 Jbavnorm_hh = -ppe.*pBT0./pBp0.*(0.69*dlnTedr + 2.44*dlnnedr); 
0175 %
0176 pcurr_hir_ref_fsav = Jbavnorm/Bt*qe*1e-6*1e3; 
0177 pcurr_sau_ref_fsav = Jbavnorm_s/Bt*qe*1e-6*1e3;
0178 pcurr_hin_ref_fsav_e = Jbavnorm_hh/Bt*qe*1e-6*1e3;
0179 %
0180 % Ion Flux to be added in the fluid limit: Hinton & Hazeltine, Hirshman and Sauter relations use the same term
0181 %
0182 A1_i = pTi./(Zi*pTe).*dlnpidr;
0183 Jbavnorm_i = -ppe.*pBT0./pBp0.*(L31.*A1_i + L34.*A4);
0184 pcurr_i_ref_fsav = Jbavnorm_i/Bt*qe*1e-6*1e3; 
0185 %
0186 pcurr_hin_ref_fsav = pcurr_hin_ref_fsav_e + pcurr_i_ref_fsav;
0187 %
0188 %Taguchi's model for anisotropic electron distribution function (Ref: J. Physical Soc. Japan, 65,11 (1996) 3530)
0189 %
0190 pcurr_tag_e_ref_fsav = NaN;
0191 pcurr_tag_i_ref_fsav = NaN;