# Tokamak GOLEM - electron energy confinement time calculation

Strategic flowchart

## Initialization

First we clear matlab workspace and define fundamental constants and parameters

clc; clear; close all

VacuumShotNo = 22475; % Vacuum shot number

PlasmaShotNo=22471; % Plasma shot number

baseURL='http://golem.fjfi.cvut.cz/utils/data/';

RogowskiCalibration=5.3e6; % Calibration for Rogowski coil [A/V(10us)] 03/2010, CASTOR calibration (Mr. Zacek), 10 us sampling, Mr. Zacek states 5.3e6 A/Vs

UloopCalibration=5.5; % Calibration for U_loop measurement

Mu0=4*pi*1e-7; % [H/m] Permeability

kB=1.3806503e-23; % [J/K] Botlzmann constant

eV=1.60217646e-19; % [J/eV] 1 eV equivalent to J

LimiterPosition=0.085; % [m] max plasma radius 0.05 for ring with Mach probes

## Chamber resistence calculation from vacuum shot

Plot raw data from NIstandard DAS (Data Acquisition System)

First let's get the data from the DAS in the columns specified by the given identifiers

Let's show raw data

figure(1);

h1 = subplot(2,1,1);

plot(Ul(:,1)*1e3,Ul(:,2))

ylabel('channel #1 [V]'); set(gca,'XTickLabels',[]); ylim([0 2]); xlim([0 40])

legend('ID:uloop')

title(['Vacuum shot: #' num2str(VacuumShotNo) ' - raw data from NIstandard DAS'])

h2 = subplot(2,1,2);

p = get(h1, 'pos'); p(2) = p(2) - p(4)*1.1; set(h2,'pos',p) % proper positon of subplot

plot(Ir(:,1)*1e3,Ir(:,2))

ylabel('channel #3 [V]'); xlabel('Time [ms]'); ylim([-.1 0.3]); xlim([0 40])

legend('ID:irog')

Get real physical quantities from raw data

Physical quantities can be obtained by multiplying the raw data by appropriate calibration constants. In the case of the current, the derivative is measured, so the signal has to be numerically integrated.

from ID:uloop needs only multiplication

U_l = Ul(:,2)*UloopCalibration;

Chamber current needs offset correction, integration and calibration

offset = mean(Ir(Ir(:,1)<=.004,2)); % Up to 4 ms no signal, good zero specification

I_ch = cumtrapz(Ir(:,1),Ir(:,2)-offset);

I_ch = I_ch*RogowskiCalibration;

And let´s plot it

figure(2);

h1 = subplot(2,1,1);

plot(Ul(:,1)*1e3,U_l)

ylabel('U_l [V]'); set(gca,'XTickLabels',[]); ylim([-4 12]); xlim([0 40])

legend('Loop voltage')

title(['Vacuum shot: #' num2str(VacuumShotNo) ' - real data'])

h2 = subplot(2,1,2);

p = get(h1, 'pos'); p(2) = p(2) - p(4)*1.1; set(h2,'pos',p) % proper positon of subplot

plot(Ir(:,1)*1e3,I_ch)

ylabel('I_{ch} [A]'); xlabel('Time [ms]'); ylim([-1e3 1.5e3]); xlim([0 40])

legend('Chamber current I_{ch}')

Final chamber resistance specification

The chanber resistnace should be calculated in the time span where the physical quantities are reasonably stationary.

t_start = find(Ul(:,1)>=0.014,1,'first');

t_end = find(Ul(:,1)>=0.015,1,'first');

Calculation via Ohm's law

R_ch = U_l(t_start:t_end)./I_ch(t_start:t_end);

R_mean = mean(R_ch);

R_err = std(R_ch);

disp(['Chamber resistance = (' num2str(R_mean*1e3,'%.2f') ' +/- ' num2str(R_err*1e3,'%.2f') ') mOhm'])

Chamber resistance = (9.38 +/- 0.02) mOhm

## Energy confinement time calculation from plasma shot

Plot raw data from NIstandard DAS and interferometer

Get data

Interpolate the density on the same time axis as the rest of the quantities

n_e = interp1(ne(:,1),ne(:,2),Ul(:,1),'pchip');

Let´s show raw data from NIstandard

figure(3);

h1 = subplot(2,1,1);

plot(Ul(:,1)*1e3,Ul(:,2))

ylabel('channel #1 [V]'); set(gca,'XTickLabels',[]); ylim([0 2]); xlim([0 40])

legend('ID:uloop')

title(['Plasma shot: #' num2str(PlasmaShotNo) ' - raw data from NIstandard DAS'])

h2 = subplot(2,1,2);

p = get(h1, 'pos'); p(2) = p(2) - p(4)*1.1; set(h2,'pos',p) % proper positon of subplot

plot(Ir(:,1)*1e3,Ir(:,2))

ylabel('channel #3 [V]'); xlabel('Time [ms]'); ylim([-.1 0.3]); xlim([0 40])

legend('ID:irog')

Let´s show raw data from interferometer diagnostics

figure(4)

plot(Ul(:,1)*1e3,n_e)

ylabel('n_e [m^{-3}]'); ylim([0 6e18]); xlim([0 40])

legend('ID:electron\_density')

title(['Plasma shot: #' num2str(PlasmaShotNo) ' - raw data from TektronixDPO'])

Get real physical data

from ID:uloop needs only multiplication

U_l = Ul(:,2)*UloopCalibration;

Plasma current needs offset correction, integration calibration and chamber current subtraction (getting from Ohm´s law)

offset = mean(Ir(Ir(:,1)<=.004,2));

I_chANDp = cumtrapz(Ir(:,1),Ir(:,2)-offset);

I_chANDp = I_chANDp*RogowskiCalibration;

I_p = I_chANDp-U_l/R_mean;

And let´s plot it

figure(5);

h1 = subplot(3,1,1);

plot(Ul(:,1)*1e3,U_l)

ylabel('U_l [V]'); set(gca,'XTickLabels',[]); ylim([-4 12]); xlim([0 40])

legend('Loop voltage')

title(['Plasma shot: #' num2str(PlasmaShotNo) ' - real data'])

h2 = subplot(3,1,2);

p = get(h1, 'pos'); p(2) = p(2) - p(4)*1.1; set(h2,'pos',p) % proper positon of subplot

plot(Ir(:,1)*1e3,I_chANDp/1e3)

ylabel('I_{ch} [kA]'); xlabel('Time [ms]'); ylim([-1 5]); xlim([0 40])

legend('Chamber + plasma current I_{ch+p}')

h3 = subplot(3,1,3);

p = get(h2, 'pos'); p(2) = p(2) - p(4)*1.1; set(h3,'pos',p) % proper positon of subplot

plot(Ir(:,1)*1e3,I_p/1e3)

ylabel('I_p [kA]'); xlabel('Time [ms]'); ylim([-1 5]); xlim([0 40])

legend('Plasma current I_p')

## Data correlation ...

Select only a signal with plasma

i_start = find(Ul(:,1)>=tstart,1,'first');

i_end = find(Ul(:,1)>=tend,1,'first');

n_e = n_e(i_start:i_end);

I_p = I_p(i_start:i_end);

U_l = U_l(i_start:i_end);

time = Ul(i_start:i_end,1);

And let´s plot it

figure(6);

h1 = subplot(3,1,1);

plot(time*1e3,U_l)

ylabel('U_l [V]'); set(gca,'XTickLabels',[]); ylim([-4 12]); xlim([tstart tend]*1e3)

legend('Loop voltage U_l')

title(['Plasma shot: #' num2str(PlasmaShotNo) ' - real data'])

h2 = subplot(3,1,2);

p = get(h1, 'pos'); p(2) = p(2) - p(4)*1.1; set(h2,'pos',p) % proper positon of subplot

plot(time*1e3,I_p/1e3)

ylabel('I_p [kA]'); set(gca,'XTickLabels',[]); ylim([-1 5]); xlim([tstart tend]*1e3)

legend('Plasma current I_p')

h3 = subplot(3,1,3);

p = get(h2, 'pos'); p(2) = p(2) - p(4)*1.1; set(h3,'pos',p) % proper positon of subplot

plot(time*1e3,n_e)

ylabel('n_e [m^{-3}]'); ylim([0 7e18]); xlabel('Time [ms]'); xlim([tstart tend]*1e3)

legend('Electron density n_e')

## Final calculations

Plasma resistivity adopted to avoid problems @ T_e calculation

R_p = U_l./I_p;

R_p(R_p<0 & isnan(R_p)) = 0;

Plasma heating power

P_OH = U_l.*I_p;

Central electron temperature (the GOLEM specific case)

T_e = 0.9.*R_p.^(-2/3);

Energy content in the plasma

W_p = PlasmaVolume.*n_e.*eV.*T_e./3;

Let´s make final calculations in the stationary phase to avoid complex derivation calculations of plasma energy balance

imin = find(time>=0.015,1,'first');

imax = find(time>=0.020,1,'first');

P_OH = P_OH(imin:imax);

W_p = W_p(imin:imax);

time = time(imin:imax);

tau = W_p./P_OH;

Plot

figure(7);

h1 = subplot(3,1,1);

plot(time*1e3,P_OH/1e3)

ylabel('P_{OH} [kW]'); set(gca,'XTickLabels',[]); ylim([20 27]);

legend('Ohmic heating power P_{OH}')

title(['Plasma shot: #' num2str(PlasmaShotNo) ' - energy balance'])

h2 = subplot(3,1,2);

p = get(h1, 'pos'); p(2) = p(2) - p(4)*1.1; set(h2,'pos',p) % proper positon of subplot

plot(time*1e3,W_p)

ylabel('W_p [J]'); set(gca,'XTickLabels',[]); ylim([0.6 1.3]);

legend('Plasma energy W_p')

h3 = subplot(3,1,3);

p = get(h2, 'pos'); p(2) = p(2) - p(4)*1.1; set(h3,'pos',p) % proper positon of subplot

plot(time*1e3,tau*1e6)

ylabel('n_e [m^{-3}]'); ylim([25 50]); xlabel('Time [ms]')

legend('Energy confinement time \tau_E')

grid on

Final calculation

tau_mean = mean(tau);

tau_err = std(tau);

disp(['Energy confinement time = (' num2str(tau_mean*1e6,'%.0f') ' +/- ' num2str(R_err*1e6,'%.0f') ') us'])

Energy confinement time = (41 +/- 17) us