This script generates the data file of an magnetic equilibrium created by the HELENA magnetic equlibrium solver (helmex77 version) By Y. Peysson (CEA/DSM/IRFM, yves.peysson@cea.fr) and J. Decker (CEA/DSM/IRFM, joan.decker@cea.fr)
0001 % 0002 % This script generates the data file of an magnetic equilibrium created 0003 % by the HELENA magnetic equlibrium solver (helmex77 version) 0004 % 0005 % By Y. Peysson (CEA/DSM/IRFM, yves.peysson@cea.fr) and J. Decker (CEA/DSM/IRFM, joan.decker@cea.fr) 0006 % 0007 clear all 0008 close all 0009 % 0010 id = 'JETh77';% scenario identification 0011 % 0012 nit = 0;%number of iteration for equilibrium consistency (jmoy <-> psi) 0013 % 0014 % Equilibrium parameters 0015 % 0016 % Here find the definitions for the shape of the separatrix 0017 % 0018 % - R_sep = (R_sep_max + R_sep_min)/2, for the guess. 0019 % - a_sep = (R_sep_max - R_sep_min)/2, for the guess. 0020 % - Z_sep = Z(R_sep_max), for the guess. 0021 % - elongation_sep = (Z_sep_max - Z_sep_min)/2 (always positive) (xpoint(1) + xpoint(5) is approximately elongation_sep*a_sep) 0022 % - upper_triangularity_sep = (R_sep - R(Z_sep_max))/2 (always positive) xpoint(2) = upper_triangularity_sep*a_sep 0023 % - lower_triangularity_sep = (R_sep - R(Z_sep_min))/2 (always positive) xpoint(6) = lower_triangularity_sep*a_sep 0024 % 0025 % 0026 % if no upper X-point exists, xpoint(2) = xpoint(3) = 0. In this case,xpoint(1) is the upper altitude of the separatrix in minor radius unit. 0027 % if no lower X-point exists, xpoint(7) = xpoint(8) = 0. In this case,xpoint(5) is the lower altitude of the separatrix in minor radius unit (always positive). 0028 % if an upper X-point exists, xpoint(2) + xpoint(3) close to 90 degrees 0029 % if an lower X-point exists, xpoint(7) + xpoint(8) close to 90 degrees 0030 % 0031 R_sep = 3.0;% Geometrical major radius of the separatrix (m) 0032 a_sep = 0.9;% Geometrical minor radius of the separatrix (m) 0033 Z_sep = 0.0;% Geometrical vertical shift of the separatrix (m) 0034 % 0035 xpoint(1) = 1.687;%Upper altitude X point in minor radius unit 0036 xpoint(2) = 0.466;%Upper triangularity in minor radius unit 0037 xpoint(3) = 0;%Upper separatrix angle (R,X) (LFS, degrees) 0038 xpoint(4) = 0;%Upper separatrix angle (-R,X) (HFS, degrees) 0039 xpoint(5) = 2.001;%Lower altitude X point in minor radius unit (always positive) 0040 xpoint(6) = 0.568;%Lower triangularity in minor radius unit 0041 xpoint(7) = 22.46;%Lower separatrix angle (R,X) (LFS, degrees) 0042 xpoint(8) = 67.92;%Lower separatrix angle (-R,X) (HFS, degrees) 0043 xpoint(9) = a_sep;%Lower separatrix angle (-R,X) (HFS, degrees) 0044 xpoint(10) = R_sep;%Lower separatrix angle (-R,X) (HFS, degrees) 0045 xpoint(11) = Z_sep;%Lower separatrix angle (-R,X) (HFS, degrees) 0046 % 0047 npsi = 101;%ONLY 101 radial profile grid points for HELMEX77 0048 ntheta = 65;%ONLY 65 poloidal grid points for HELMEX77 0049 % 0050 Ip = 2;% Plasma current (MA) 0051 % 0052 qmin = 1;% Safety factor q0 at plasma center 0053 eq = 3;% Exponent for q radial profile (q = (qmax - qmin)*(r/a).^eq + qmin, qmax calaculated by the Ampere's theorem at the plasma edge with a circular plasma cross-section) 0054 qopt = 1;%Option for q profile. (0): constant (default, uniform current, psi=rho^2), (1): from qmin and eq 0055 % 0056 Bt = 3.5;% Toroidal magnetic field on the magnetic axis (T) 0057 % 0058 % Kinetic parameters (profiles) 0059 % 0060 Zi = [1,1,1,6];% Ion types: (1) H/D/T, (2) He, ..., (6) C [1,p] (WARNING: Zi must be [1,1,1,imp1,imp2] for hydrogen plasmas) 0061 mi = [1,2,3,12];% Ion mass (uma) [1,p] (WARNING: Zi must be [1,2,3,mimp1,mimp2] for hydrogen plasmas) 0062 fi = [0,1,0];% Hydrogen isotopic fraction (H/D/T) [1,3] (WARNING: only used when hydrogen plasmas are considered) 0063 % 0064 Te0 = 5;% Core electron temperature (keV) 0065 Tea = 0.5;% Edge electron temperature (keV) 0066 eTe1 = 4;% Exponent #1 for Te profile 0067 eTe2 = 2;% Exponent #2 for Te profile 0068 % 0069 ne0 = 2.0e19;% Core electron density (m-3) 0070 nea = 2.0e17;% Edge electron density (m-3) 0071 ene1 = 4;% Exponent #1 for ne profile 0072 ene2 = 0.5;% Exponent #2 for ne profile 0073 % 0074 Ti0 = 2;% Core ion temperature (keV) 0075 Tia = 0.2;% Edge ion temperature (keV) 0076 eTi1 = 4;% Exponent #1 for Ti profile 0077 eTi2 = 2;% Exponent #2 for Ti profile 0078 % 0079 Zeff0 = 1.5;% Core effective charge (a.u.) 0080 Zeffa = 1.5;% Edge effective charge (a.u.) 0081 eZeff = 0;% Define the slope at the radial point where the derivative is maximum (must be >> 1) 0082 xZeff = 0.5;%Radial normalized position at which the derivative is maximum 0083 % 0084 % Initial calculations to set up the guess for HELMEX77 (cylindrical magnetic equilibrium with circular concentric flux surfaces) 0085 % 0086 % Magnetic data 0087 % 0088 [ppsin,ppsi_apRp,theta,x,y,Bx,By,BPHI,pBpp,pq_Rpap,pj,pmag,Ip_test,prho_mag,pspsin,pspsinT,ppsiT] = idealequilcyl_yp(a_sep,R_sep,Z_sep,Bt,Ip,qmin*R_sep/a_sep,eq,qopt,npsi,ntheta);%Cylindrical magnetic equilibrium 0089 % 0090 [prho,pTe,pne,pzTi,pzni,zZi,zmi,fi,pkin] = idealprof_yp(Zi,mi,fi,Te0,Tea,eTe1,eTe2,ne0,nea,ene1,ene2,Ti0,Tia,eTi1,eTi2,Zeff0,Zeffa,eZeff,xZeff,npsi);%Profiles 0091 % 0092 figure('Name','Guess magnetic equilibrium'), 0093 subplot(2,2,1),plot(prho,pj/1e6),xlabel('\rho');ylabel('j (MA/m^2)');title(['I_p = ',num2str(round(Ip_test)/1e6),' MA']) 0094 subplot(2,2,2),plot(prho,pBpp),xlabel('\rho');ylabel('B_p(T)');title(['B_t ',num2str(Bt),' T']) 0095 subplot(2,2,3),plot(prho,pq_Rpap*a_sep/R_sep),xlabel('\rho');ylabel('q') 0096 subplot(2,2,4),plot(prho,ppsin),xlabel('\rho');ylabel('\psi_n') 0097 % 0098 pj_min = min(pj); 0099 if pj_min < 0,pj = pj - pj_min;end,%pj must be non-negative. Only the shape is considered and Ip for the total plasma current 0100 [equil_magnetic,prho_out,pspsin_out,pspsinT_out,pp_out,li_out] = equil_magnetic_helmex77_yp(xpoint,Bt,Ip,ppsi_apRp*R_sep/a_sep,pj,pkin,ntheta,nit,2); 0101 % 0102 [qe] = pc_dke_yp; 0103 figure('Name','Pressure profile'), 0104 subplot(121),plot(prho,pmag*qe*1000/100000),xlabel('\rho');ylabel('Magnetic pressure (Bar)');title(['Guess magnetic']),axis('square') 0105 subplot(122),plot(prho,pkin*qe*1000/100000,'b',prho_out,pp_out*qe*1000/100000,'r--'),xlabel('\rho');ylabel('Kinetic pressure (Bar)');legend('Guess kinetic','HELENA'),axis('square') 0106 % 0107 equil_prof.pTe = pTe; 0108 equil_prof.pne = pne; 0109 equil_prof.pzTi = pzTi; 0110 equil_prof.pzni = pzni; 0111 equil_prof.zZi = zZi; 0112 equil_prof.zmi = zmi; 0113 % 0114 equil = conc_struct_jd(equil_magnetic,equil_prof); 0115 equil.id = id; 0116 % 0117 eval(['save EQUIL_',id,'.mat equil']); 0118 % 0119 equilibrium_jd(equil,'',1) 0120 0121 0122 0123 0124 0125