This notebook estimates several parameters of the plasma in the context of tokamak fusion physics. These parameters include but are not limited to the safety factor, the electron temperature, electron pressure, plasma volume and electron thermal energy and electron energy confinement time. Other more general plasma parameters are calculated as well.
The formulas and explanations are mostly based on the book [1] WESSON, John. Tokamaks. 3. ed. Oxford: Clarendon press, 2004. ISBN 9780198509226. and the reader is encouraged to consult it for details.
The accuracy of these parameters stronly depends on the availability of the plasma position and size reconstruction.
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import integrate, signal, interpolate, constants
import holoviews as hv
hv.extension('bokeh')
import hvplot.pandas
import requests
# 180221 Dirigent
destination='Results/'
os.makedirs(destination, exist_ok=True)
def print_and_save(phys_quant, value, format_str='%.3f'):
print(phys_quant+" = %.5f" % value)
with open(destination+phys_quant, 'w') as f:
f.write(format_str % value)
#update_db_current_shot(phys_quant,value)
def update_db_current_shot(field_name, value):
#os.system('psql -c "UPDATE shots SET '+field_name+'='+str(value)+' WHERE shot_no IN(SELECT max(shot_no) FROM shots)" -q -U golem golem_database')
subprocess.call(["psql -q -U golem golem_database --command='UPDATE shots SET \""+field_name+"\"="+str(value)+" WHERE shot_no IN(SELECT max(shot_no) FROM shots)'"],shell=True)
shot_no = 44267 # 33516 is a good test case
The following analysis makes sense only if a plasma was present in the discharge
def plasma_scalar(name):
r = requests.get(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/PlasmaDetection/Results/{name}')
return float(r.content)
t_plasma_start = plasma_scalar('t_plasma_start')
t_plasma_end = plasma_scalar('t_plasma_end')
if t_plasma_start == t_plasma_end:
raise RuntimeError('no plasma in this discharge, analysis cannot continue')
basic_diagn_signals = ['U_loop', 'Bt', 'Ip']
df = pd.concat([pd.read_csv(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/BasicDiagnostics/Results/{sig}.csv',
names=['time', sig], index_col=0) for sig in basic_diagn_signals], axis='columns')
df = df.loc[t_plasma_start:t_plasma_end] # time slice with plasma
#try:
# df_position = pd.read_csv(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/LimiterMirnovCoils/plasma_position.csv',
# index_col=0)
#except HTTPError:
df_position = None
using_plasma_position = df_position is not None
# TODO load from SQL
R0 = 0.4 # chamber center Major plasma radius [m]
a0 = 0.085 # maximum Minor radius - limiter [m]
If plasma position and size reconstruction is not available, the parameters of the chamber geometry are used for the minor and major plasma radii aa and RR, respectivelly.
def interp_position2basic(name):
interp = interpolate.interp1d(df_position.index, df_position[name], bounds_error=False)
return interp(df.index) * 1e-3 # mm to m
if using_plasma_position:
df['R'] = R0 + interp_position2basic('r')
df['a'] = interp_position2basic('a')
else:
df['R'] = R0
df['a'] = a0
On any given closed flux surface in the plasma in the tokamak the magnetic field line performs qq transits in the toroidal angle ϕϕ per 1 one transit in the poloidal angle θθ. The stronger the toroidal magnetic field is, the more stable the plasma becomes against various instabilities, especially against the kink instability which can occur for q<1q<1. For this reason qq is referred to as the safety factor.
In a simple tokamak with a circular cross-section (such as GOLEM) the poloidal magnetic field can be estimated at least at the very edge of the plasma from the total plasma current IpIp enclosed by the plasma column of minor radius aa and major radius RR as Bθa=μ0Ip2πaBθa=μ0Ip2πa
Typically, in a tokamak the toroidal magnetic field BϕBϕ is several times stronger than the poloidal magnetic field BθaBθa at the egde.
df['B_theta'] = constants.mu_0 * df['Ip'] * 1e3 / (2*np.pi*df['a']) # Ip is in [kA], B_theta will be in [T]
df[['Bt','B_theta']].hvplot(ylabel='B [T]', xlabel='time [ms]', grid=True, logy=True, ylim=(1e-4,None))
For a large aspect ratio tokamak (i.e. the inverse aspect ratio is small ϵ=aR<<1ϵ=aR<<1) such as GOLEM the safety factor at the edge on the last closed flux surface (LCFS) delimited by the limiter ring can be estimated as qa=aBϕRBθqa=aBϕRBθ
df['q_a'] = df['a'] * df['Bt'] / (df['R'] * df['B_theta'])
df['q_a'].hvplot(logy=True, grid=True)
WARNING:param.CurvePlot02008: Logarithmic axis range encountered value less than or equal to zero, please supply explicit lower-bound to override default of 930.290.
To obtain information on qq and BθBθ deeper inside the plasma torus one must have knowledge of or assume a specific profile for the toroidal current density jϕjϕ. A common approximation for a tokamak such as GOLEM is a poloidally symmetric radial profile jϕ(r)=j0(1−(ra))νjϕ(r)=j0(1−(ra))ν where rr is the radius with respect to the plasma center and νν a so called "peaking factor". A common choice is ν=1ν=1 for a "parabolic" profile or ν=2ν=2 for a more peaked profile (likely more realistic). With the average current density defined as ⟨j⟩a=Ipπa2⟨j⟩a=Ipπa2 the maximum current density j0j0 can be estimated from the relation j0⟨j⟩a=ν+1j0⟨j⟩a=ν+1
nu = 2 # probably more realistic than parabolic
df['j_avg_a'] = df['Ip'] *1e3 / (np.pi*df['a']**2)
df['j_0'] = df['j_avg_a'] * (nu+1)
Under this assumption the safety factor in the plasma core (r=0r=0) is reduced according to the relation qaq0=ν+1qaq0=ν+1. which could result in the following profiles for the time when qaqa is the lowest (i.e. closest to an instability).
t_q_a_min = df['q_a'].idxmin()
df_q_a_min = df.loc[t_q_a_min]
print(f'min(q_a)={df_q_a_min["q_a"]:.2f} at t={t_q_a_min:.3f} ms')
min(q_a)=-66054.75 at t=14.252 ms
r = np.linspace(0, df_q_a_min['a'])
df_r = pd.DataFrame({
'j': df_q_a_min['j_0'] * (1-(r/df_q_a_min['a'])**2)**nu,
'B_theta': constants.mu_0 * df_q_a_min['j_0'] * df_q_a_min['a']**2 / (2*(nu+1)) * (1-(1-(r/df_q_a_min['a'])**2)**(nu+1))/r,
'q': 2*(nu+1)/(constants.mu_0 * df_q_a_min['j_0']) * (df_q_a_min['Bt']/df_q_a_min['R']) * (r/df_q_a_min['a'])**2 / (1-(1-(r/df_q_a_min['a'])**2)**(nu+1))
}, index=pd.Index(r, name='r')
)
<ipython-input-1-ec0bca19fa99>:4: RuntimeWarning: invalid value encountered in divide 'B_theta': constants.mu_0 * df_q_a_min['j_0'] * df_q_a_min['a']**2 / (2*(nu+1)) * (1-(1-(r/df_q_a_min['a'])**2)**(nu+1))/r, <ipython-input-1-ec0bca19fa99>:5: RuntimeWarning: invalid value encountered in divide 'q': 2*(nu+1)/(constants.mu_0 * df_q_a_min['j_0']) * (df_q_a_min['Bt']/df_q_a_min['R']) * (r/df_q_a_min['a'])**2 / (1-(1-(r/df_q_a_min['a'])**2)**(nu+1))
df_r.hvplot.line(subplots=True, shared_axes=False, width=300, grid=True)
The plasma is typically as conductive as copper, i.e. is a good conductor with a relatively low resitivity. However, whereas the resitivity of metals increases with temperature, the resitivity of a plasma decreases, because at higher timperatures collisions between particles become less frequent, leading to less resistance to their movement. While with higher particle density the number of collisions increases, the number of charge cariers also increases, so in the end the resistivity does not depend on density.
The simple, unmagnetized plasma resistivity derived by Spitzer ηs=0.51√mee2lnΛ3ϵ20(2πkBTe)32ηs=0.51√mee2lnΛ3ϵ20(2πkBTe)32 with the constants electron mass meme, elementary charge ee, vacuum permitivity ϵ0ϵ0 and kBkB the Boltzmann constant. lnΛlnΛ is the so called Coulomb logarithm which has a weak dependence on density and temperature and for typical GOLEM plasmas can be held at lnΛ∼14lnΛ∼14. The factor 0.51 comes from more precise calculations which show that the parallel resitivity η‖=ηsη∥=ηs (along the magnetic field-line the resistivity is not affected by the field) is halved compared to the classical (analytical) perpendicular resitivity η⊥=1.96η‖η⊥=1.96η∥ though in reality the perpendicular resitivity can be higher due to anomalous transport (turbelence, etc.). If one is interested in the electron temperature TeTe in the units of electron-volts (typically used in the field), the relation is Te[eV]=kBeTe[K]Te[eV]=kBeTe[K].
Additional corrections:
This results in ηmeasured=ηsZeff(1−√ϵ)−2ηmeasured=ηsZeff(1−√ϵ)−2.
These considerations lead to the relation Te[eV]=1e2π(1.96Zeff(1−√ϵ)2ηmeasured3ϵ20√mee2lnΛ)−23Te[eV]=1e2π(1.96Zeff(1−√ϵ)2ηmeasured3ϵ20√mee2lnΛ)−23
def electron_temperature_Spitzer_eV(eta_measured, Z_eff=3, eps=0, coulomb_logarithm=14):
eta_s = eta_measured / Z_eff * (1-np.sqrt(eps))**2
term = 1.96 * eta_s * (3 * constants.epsilon_0**2 /
(np.sqrt(constants.m_e) * constants.elementary_charge**2 * coulomb_logarithm))
return term**(-2./3) / (constants.elementary_charge * 2*np.pi)
To estimate ηmeasuredηmeasured one can use Ohm's law in the form jϕ=σEϕjϕ=σEϕ with the plasma conductivity σ=1ηmeasuredσ=1ηmeasured. The toroidal electric field can be estimated from the loop voltage, but one must take into account inductive effects as well. Neglecting mutual inductances between e.g. the plasma and the chamber, the loop voltage induced in the plasma by the primary winding is "consumed" by the electric field and current inductance as Uloop=2πREϕ+(Li+Le)dIpdtUloop=2πREϕ+(Li+Le)dIpdt where LiLi and LeLe are the internal and external plasma inductances, respectively. The external inductance of a closed toroidal current (assuming a uniform current density) is Le=μ0Rln(8Ra−74)Le=μ0Rln(8Ra−74). The internal plasma inductance is usually parametrized as Li=μ0Rli2Li=μ0Rli2 where lili is the so called normalized internal inductance which depends on the BθBθ (or rather current) profile. For the assumed current profile an accurate estimate is li≈ln(1.65+0.89ν)li≈ln(1.65+0.89ν).
l_i = np.log(1.65+0.89*nu)
df['L_p'] = constants.mu_0 * df['R'] * (np.log(8*df['R']/df['a']) - 7/4. + l_i/2)
dt = np.diff(df.index.values[:2]).item()
n_win = int(0.5 / dt) # use a window of 0.5 ms
if n_win % 2 == 0:
n_win += 1 # window must be odd
# needs SI units: convert current in kA -> A, time step in ms->s
df['dIp_dt'] = signal.savgol_filter(df['Ip']*1e3, n_win, 3, 1, delta=dt*1e-3) # 1. derivative of an order 3 polynomial lsq SG-filter
df['E_phi_naive'] = df['U_loop'] / (2*np.pi*df['R']) # V/m
df['E_phi'] = (df['U_loop'] - df['L_p']*df['dIp_dt']) / (2*np.pi*df['R'])
df[['E_phi_naive', 'E_phi']].hvplot(ylabel='E_phi [V/m]', grid=True)
In the beginning of the discharge the creationg of the poloidal magnetic field by the plasma current diminishes EϕEϕ, and at the end the plasma current and its field dissipates, enhancing EϕEϕ. With the estimated EϕEϕ, one can obtain an average temperature estimate with ⟨j⟩a⟨j⟩a and a (higher) core plasma temperature estimate with j0j0, respectively.
for s in ('0', 'avg_a'):
for k in ('', '_naive'):
df[f'eta_{s}'+k] = df['E_phi'+k] / df[f'j_{s}']
df['eps'] = df['a'] / df['R']
for s in ('0', 'avg_a'):
for k in ('', '_naive'):
df[f'Te_{s}'+k] = electron_temperature_Spitzer_eV(df[f'eta_{s}'+k], eps=df['eps'])
df[['Te_0', 'Te_avg_a', 'Te_0_naive']].hvplot(ylabel='Te [eV]', grid=True)
plt.figure(figsize=(3, 2))
ax = df['Te_0'].plot.line()
ax.set(xlabel='time [ms]', ylabel='Te(r=0) [eV]')
plt.margins(0)
plt.tight_layout(pad=0.1)
plt.savefig('icon-fig.png')
A good estimate of the (line-averaged) electron density (concentration) is typically obtained from the microwave interferoemter. In the absence of this diagnostic an order-of-magntude estimate can be obtained using the ideal gas law applied to the initial inert state of the working gas. Since the whole chamber has a volume of V0≈60lV0≈60l, the working gas with the pre-discharge stationary equilibrium pressure p0p0 at the room temperature T0≈300KT0≈300K will is expected to be composed of NN molecules according to the relation p0V0=NkBT0p0V0=NkBT0. One can assume that for a gven working gas the molecule dissasociates into kaka atoms which can the fully ionaize giving keke electrons. Therefore, one can estimate the order-of-magnitude number of electrons (an upper estimate due to only partial ionaization of the working gas) as Ne≈kakep0V0kBT0Ne≈kakep0V0kBT0
def chamber_parameter(name):
r = requests.get(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Operation/Discharge/{name}')
v = r.content
try:
return float(v)
except ValueError:
return v.strip().decode()
p_0 = chamber_parameter('p_chamber_pressure_predischarge') *1e-3 # from mPa to Pa
working_gas = chamber_parameter('X_working_gas_discharge_request')
working_gas, p_0
('H', 0.019899999999999998)
if working_gas == 'H':
k_a = 2 # binary molecule
k_e = 1
elif working_gas == 'He':
k_a = 2
k_e =2
else:
raise RuntimeError(f'Unknown working gas {working_gas}')
V_0 = 60e-3 # m^3
T_0 = 300
N_e = k_a*k_e* (p_0 * V_0) / (T_0 * constants.k)
N_e
5.765404530767775e+17
To estimate the actual electron density nene , i.e. number of electrons in m−3m−3 one must estimate also the plasma volume VpVp. Assuming a perfect plasma torus, its volume is tha cartesian product of its poloidal cross section (circular - πa2πa2) along the toroidal axis of the torus (length 2πR2πR), together Vp=2π2Ra2Vp=2π2Ra2. The plasma density is then ne≈Ne/Vpne≈Ne/Vp.
df['V_p'] = 2*np.pi**2*df['R']*df['a']**2
df['n_e'] = N_e / df['V_p'] # in m^-3
df['n_e'].mean()
1.0106533066339652e+19
The thermal energy of electrons in the plasma Wth,eWth,e evolves according to the applied heating power PHPH and the (turbulent and radiative) losses summarized by the loss power PLPL as dWth,edt=PH−PLdWth,edt=PH−PL The electron thermal energy can be approximated suing the plasma electron pressure pe=nekBTepe=nekBTe as Wth,e≈TekBneVpWth,e≈TekBneVp.
df['p_e'] = df['n_e'] * df['Te_avg_a'] * constants.elementary_charge # Te is in eV, p_e will be in Pa
df['W_th_e'] = df['p_e'] * df['V_p'] # in Jouls
df[['p_e', 'W_th_e']].mean()
p_e 56.430378 W_th_e 3.219145 dtype: float64
In the absence of auxiliary heating systems such as NBI an ECRH, the only component of the heating power is the resistive (ohmic) heating power density due to the toroidal electric field and current EϕjϕEϕjϕ . Assuming a uniform distribution of this heating density, the total ohmic heating power can be estimated as PH=PΩ=Eϕ⟨jϕ⟩aVpPH=PΩ=Eϕ⟨jϕ⟩aVp. Due to the geometric assumptions used above, this is equivalent to the total induced power with the change of the poloidal magnetic energy subtracted PH=UloopIp−ddt(12(Le+Li)I2p)PH=UloopIp−ddt(12(Le+Li)I2p)
df['P_mag'] = df['L_p'] * df['Ip'] * df['dIp_dt'] # [kW] equivalent after chain rule
df['P_H'] = df['U_loop'] * df['Ip'] - df['P_mag'] # [kW]
A figure of merit critical for thermonuclear fusion is the characteristic time scale at which the thermal energy would be exponentially depleted under the assumption that the loss power is proportional to the stored thermal energy PL∝WthPL∝Wth. This time scale is called the energy confinement time τEτE and for the electron energy it can be estimated from the modified electron thermal energy balance with PL≈Wth,e/τE,ePL≈Wth,e/τE,e dWth,edt=PH−Wth,eτE,edWth,edt=PH−Wth,eτE,e
df['dW_th_e_dt'] = signal.savgol_filter(df['W_th_e'].fillna(0), n_win, 3, 1, delta=dt) # [kW] 1. derivative of an order 3 polynomial lsq SG-filter
df['tau_E_e'] = df['W_th_e'] / (df['P_H'] - df['dW_th_e_dt']) # [ms] -< J/kW
df[['P_H', 'P_mag', 'dW_th_e_dt']].hvplot(grid=True, ylim=(-3, None), ylabel='Power [kW]')
kwd = dict(grid=True, title='', height=250, width=600)
kw = dict(xlabel='', ylim=(0, None), **kwd)
l = df['U_loop'].hvplot(ylabel='loop voltage [V]', **kw) +\
df['Ip'].hvplot(ylabel='plasma current [kA]', **kw) +\
df['a'].hvplot(ylabel='minor radius [m]', **kw) +\
df['q_a'].hvplot(ylabel='edge safety factor [1]', **kw) +\
df['Te_avg_a'].hvplot(ylabel='average elecron temperature [eV]', **kw) +\
df['tau_E_e'].hvplot(ylabel='el. energy confinement time [ms]', xlabel='time [ms]', ylim=(0, 0.5), **kwd)
l.cols(2)
df_overview=df[['q_a', 'E_phi', 'eps', 'V_p', 'Te_0', 'Te_avg_a', 'n_e', 'p_e', 'W_th_e', 'P_H', 'P_mag', 'tau_E_e']].describe().iloc[1:] # skip count
df_overview
q_a | E_phi | eps | V_p | Te_0 | Te_avg_a | n_e | p_e | W_th_e | P_H | P_mag | tau_E_e | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
mean | 17.928906 | 4.692825 | 2.125000e-01 | 5.704631e-02 | 72.490517 | 34.849806 | 1.010653e+19 | 56.430378 | 3.219145 | 30.832803 | -0.000097 | 0.147996 |
std | 1199.535539 | 1.644142 | 4.310641e-14 | 3.101830e-15 | 26.280143 | 12.634175 | 2.490484e+06 | 20.457826 | 1.167044 | 15.965430 | 2.145134 | 1.339362 |
min | -66054.748305 | 2.067672 | 2.125000e-01 | 5.704631e-02 | 0.109054 | 0.052428 | 1.010653e+19 | 0.084893 | 0.004843 | -1.672476 | -5.365057 | -44.281005 |
25% | 7.411773 | 3.360381 | 2.125000e-01 | 5.704631e-02 | 54.497331 | 26.199584 | 1.010653e+19 | 42.423549 | 2.420107 | 18.803222 | -2.311349 | 0.055552 |
50% | 8.548098 | 3.978364 | 2.125000e-01 | 5.704631e-02 | 78.365148 | 37.674034 | 1.010653e+19 | 61.003495 | 3.480024 | 31.966239 | 0.885045 | 0.109408 |
75% | 13.513102 | 6.353391 | 2.125000e-01 | 5.704631e-02 | 95.298185 | 45.814589 | 1.010653e+19 | 74.185049 | 4.231984 | 43.973938 | 1.949422 | 0.179475 |
max | 78566.808542 | 8.374794 | 2.125000e-01 | 5.704631e-02 | 107.942017 | 51.893109 | 1.010653e+19 | 84.027664 | 4.793468 | 60.729293 | 3.315267 | 122.944256 |
for agg in ('mean', 'std','min', 'max'):
for quantity, value in df_overview.loc[agg].iteritems():
print_and_save(quantity+'_'+agg, value)
q_a_mean = 17.92891 E_phi_mean = 4.69283 eps_mean = 0.21250 V_p_mean = 0.05705 Te_0_mean = 72.49052 Te_avg_a_mean = 34.84981 n_e_mean = 10106533066339651584.00000 p_e_mean = 56.43038 W_th_e_mean = 3.21915 P_H_mean = 30.83280 P_mag_mean = -0.00010 tau_E_e_mean = 0.14800 q_a_std = 1199.53554 E_phi_std = 1.64414 eps_std = 0.00000 V_p_std = 0.00000 Te_0_std = 26.28014 Te_avg_a_std = 12.63417 n_e_std = 2490483.57036 p_e_std = 20.45783 W_th_e_std = 1.16704 P_H_std = 15.96543 P_mag_std = 2.14513 tau_E_e_std = 1.33936 q_a_min = -66054.74831 E_phi_min = 2.06767 eps_min = 0.21250 V_p_min = 0.05705 Te_0_min = 0.10905 Te_avg_a_min = 0.05243 n_e_min = 10106533066337161216.00000 p_e_min = 0.08489 W_th_e_min = 0.00484 P_H_min = -1.67248 P_mag_min = -5.36506 tau_E_e_min = -44.28100 q_a_max = 78566.80854 E_phi_max = 8.37479 eps_max = 0.21250 V_p_max = 0.05705 Te_0_max = 107.94202 Te_avg_a_max = 51.89311 n_e_max = 10106533066337161216.00000 p_e_max = 84.02766 W_th_e_max = 4.79347 P_H_max = 60.72929 P_mag_max = 3.31527 tau_E_e_max = 122.94426