Plasma Profiles

Overview

In PROCESS the density, temperature and current profiles of the plasma for electrons and ions can take two forms depending on the switch value for ipedestal. Either without a pedestal, ipedestal == 0 or with a pedestal ipedestal == 1. ipedestal == 0 is better suited for modelling L-mode plasmas, while ipedestal == 1 is better suited for modelling H-mode plasmas.

The files responsible for calculating and storing the profiles are plasma_profiles.py and profiles.py. A central plasma profile object is created from the PlasmaProfile class that contains attributes for the plasma density and temperature. The density and temperature profiles are in themselves objects of the Profile abstract base class. Profile, NProfile and TProfile are all defined in profiles.py. PlasmaProfile is exclusively in plasma_profiles.py

Figure 1: UML class breakdown of the plasma profiles

Parabolic Profile | L-mode

If ipedestal == 0 then no pedestal is present and the function describing the profiles is given by:

$$\begin{aligned} \mbox{Density : } n(\rho) & = n_0 \left( 1 - \rho^2 \right)^{\alpha_n} \\ \mbox{Temperature : } T(\rho) & = T_0 \left( 1 - \rho^2 \right)^{\alpha_T} \\ \mbox{Current : } J(r) & = J_0 \left( 1 - \rho^2 \right)^{\alpha_J} \end{aligned}$$

where $\rho = r/a$, and $a$ is the plasma minor radius. This gives volume-averaged values $\langle n \rangle = n_0 / (1+\alpha_n)$, and approximate line-averaged values of $\bar{n} \approx n_0 / \sqrt{(1+\alpha_n)}$, the full equation used can be seen here. These volume- and line-averages are used throughout the code along with the profile indices $\alpha$, in the various physics models, many of which are fits to theory-based or empirical scalings. Thus, the plasma model in PROCESS may be described as ½-D. The relevant profile index variables are alphan, alphat and alphaj, respectively. For full derivation and description of the core and line averaged values please see the parabolic_paramterisation() section which calculates these values for both temperature and density. A table of the the associated variables can be seen below

Profile parameter	Density	Temperature	Current
Plasma centre value	ne0, $n_0$	te0, $T_0$	N/A, $J_0$
Profile index/ peaking parameter	alphan, $\alpha_n$	alphat, $\alpha_T$	alphaj, $\alpha_J$

Plasma current profile

While PROCESS assumes a standard parabolic profile to be the shape of the current profile as per the 1989 ITER physics guidelines² , it does not calculate its values. The profile peaking factor alphaj is calculated in the plasma current calculation relating to iprofile found here. Only the temeprature and density profiles are calculated fully withing the PlasmaProfiles class.

The graph below is for a standard parabolic profile. You can vary the core value (n0) and the profile index (alphan) to see how the function behaves

---

Pedestal Profile | H-mode

If ipedestal == 1 there is now a pedestal present in the profile and they follow the form shown below:

$$\begin{aligned} \mbox{Density:} \qquad n(\rho) = \left\{ \begin{aligned} & \text{n}_{\text{ped}} + (n_0 - \text{n}_{\text{ped}}) \left( 1 - \frac{\rho^2}{\rho_{\text{ped},n}^2}\right)^{\alpha_n} & 0 \leq \rho \leq \rho_{\text{ped},n} \\ & \text{n}_{\text{sep}} + (\text{n}_{\text{ped}} - \text{n}_{\text{sep}})\left( \frac{1- \rho}{1-\rho_{\text{ped},n}}\right) & \rho_{\text{ped},n} < \rho \leq 1 \end{aligned} \right. \end{aligned}$$

$$\begin{aligned} \mbox{Temperature:} \ T(\rho) = \left\{ \begin{aligned} & \text{T}_{\text{ped}} + (T_0 - \text{T}_{\text{ped}}) \left( 1 - \frac{\rho^{\beta_T}} {\rho_{\text{ped},T}^{\beta_T}}\right)^{\alpha_T} & 0 \leq \rho \leq \rho_{\text{ped},T} \\ & \text{T}_{\text{sep}} + (\text{T}_{\text{ped}} - \text{T}_{\text{sep}})\left( \frac{1- \rho}{1-\rho_{\text{ped},T}}\right) & \rho_{\text{ped},T} < \rho \leq 1 \end{aligned} \right. \end{aligned}$$

Subscripts $0$, $\text{ped}$ and $\text{sep}$, denote values at the centre ($\rho = 0$), the pedestal ($\rho = \rho_{\text{ped}}$) and the separatrix ($\rho=1$), respectively. The density and temperature peaking parameters $\alpha_n$ and a $\alpha_T$ as well as the second exponent $\beta_T$ (input parameter tbeta, not to be confused with the plasma beta) in the temperature profile can be chosen by the user, as can the pedestal heights and the values at the separatrix (neped, nesep for the electron density, and teped, tesep for the electron temperature); the ion equivalents are scaled from the electron values by the ratio of the volume-averaged values).

Pedestal setting

If ipedestal == 1 then the pedestal density neped is set as a fraction fgwped of the Greenwald density (providing fgwped >= 0). The default value of fgwped is 0.8².

A table of the the associated variables can be seen below

Profile parameter	Density	Temperature
Pedestal radius (r/a)	rhopedn, $\rho_{\text{ped},n}$	rhopedt, $\rho_{\text{ped},T}$
Plasma centre value	ne0, $n_0$	te0, $T_0$
Pedestal value	neped, $n_{\text{ped}}$	teped, $T_{\text{ped}}$
Separatrix value	nesep, $n_{\text{sep}}$	tesep, $T_{\text{sep}}$
Profile index/ peaking parameter	alphan, $\alpha_n$	alphat, $\alpha_T$
Profile index $\beta$		tbeta, $\beta_T$

The graph below is for a standard pedestal profile. You can vary its attributes given in the table above to see how the function behaves.

Un-realistic profiles

If setting ipedestal == 1 it is highly recommended to make sure constraint equation 81 (icc=81) is active. This enforces solutions in which $n_0$ has to be greater than $n_{\text{ped}}$. Negative $n_0$ values can also arise during iteration, so it is important to be weary on how low the lower bound for $n_{\text{e}} (\mathtt{dene})$ is set. More info can be found here

---

Plasma Profile Class | PlasmaProfile

Initialization | __init__()

The parent plasma profile class is PlasmaProfile. Initialization sets the profile class size and neprofile and teprofile to NProfile & TProfile objects from profiles.py

Note

Profile sizes are set to 501 point by default. this can be varied in the __init__ of PlasmaProfile. Changing this will affect the values when doing Simpsons rule integration on the profiles.

---

Runner function | run()

This function is used to run the appropriate function relating to the calculation of the PlasmaProfile class. At the moment it only runs the pramterise_plasma() function

---

Plasma Parameterization | parameterise_plasma()

The input variable tratio can be used to define the ratio between the volume averaged ion and electron temperatures. By default it is set so that the volume averaged ion and electron temperatures are equal. The value of $T_{\text{i}}$ is set such that:

$$T_{\text{i}} = \mathtt{tratio}\times T_{\text{e}}$$

Depending on the value of ipedestal different functions will be ran, the different run cases are given below:

Parabolic Profile | ipedestal == 0

parabolic_paramterisation()

If pedestal profile values are set they are reset to have values that agree with the original form of the parabolic profiles. Such that $\rho_{\text{ped}} = 1$ and that pedestal and separatrix densities and temepratures are zero. This will then warn the user in the terminal.

The density and temperature profile runner function TProfile/NProfile.run() is then called to re-calculate the profile and core values.

Ratio of density-weighted to volume-averaged temperature factor is calculated:

$$\mathtt{pcoef} = \frac{(1+\alpha_n)(1+\alpha_T)}{(1+\alpha_T + \alpha_n)}$$

---

The line averaged density is then calculated for the profile paramaters

Line averaged density derivation

Line averaged electron density is calculated by integrating the profile across the normalised width of the profile and then dividing by the width of the integration bounds

$$\mathtt{dnla}\ | \ \bar{n_{\text{e}}} = \frac{\int^1_0 n_0(1-\rho^2)^{\alpha_n} \ d\rho}{\rho}$$

This can be more easily represented in radial coordinates and returning $\rho = \frac{r}{a}$

$$\mathtt{dnla} \ | \ \bar{n_{\text{e}}} = \frac{\int^a_0 n_0(1-r^2/a^2)^{\alpha_n} \ dr}{a}$$

$$\mathtt{dnla} \ | \ \bar{n_{\text{e}}} = \frac{\left[\frac{an_0}{2}\frac{\Gamma(1/2)\Gamma(\alpha_n+1)}{\Gamma(\alpha_n+3/2)}\right]}{a}$$

$\Gamma$ is the gamma function and is calculated in the code with scipy.special.gamma()

$$\mathtt{dnla} / \bar{n_{\text{e}}} = \frac{n_0}{2}\frac{\Gamma(1/2)\Gamma(\alpha_n+1)}{\Gamma(\alpha_n+3/2)}$$

This is in agreement with the derivation from the ITER Physics Design 1989 ²

$\blacksquare$

---

The density weighted temperatures are set:

$$\begin{aligned} \mathtt{ten} = \mathtt{pcoef} \times T_\text{e} \\ \mathtt{tin} = \mathtt{pcoef}\times T_\text{i} \end{aligned}$$

---

The central values for temperature and density is then calculated for the profile paramaters

Core temperature & density derivation

We calculate the volume integrated profile and then divide by the volume of integration to get the volume average density $\langle n \rangle$. If we assume the plasma to be a torus of circular cross-section then we can use spherical coordinates. We can simplify the problem by representing the torus as a cylinder of height equal to the circumference of the torus equal to $2\pi R$ where $R$ is the major radius of the torus.

The cylindrical volume element is given by:

$$V = \int \int \int dV = \int^{2\pi R}_0 \int^{2\pi}_0 \int^a_0 r \ dr \ d\theta \ dz$$

Inserting our density function in the form where $\rho$ is expanded as $\rho = r/a$ we get:

$$\int^{2\pi R}_0 \int^{2\pi}_0 \int^a_0 r \left(n_0(1-r^2/a^2)^{\alpha_n}\right) \ dr \ d\theta \ dz$$

Since our density function is only a function of $r$, and the torus is symmetric around its center, the integration simplifies to integrating over $r$ and the $d\theta ,\ dz$ integrals are solved to give values for the full poloidal angle and cylindrical height and torus length, leading to:

$$4\pi^2R \int^a_0 r \left(n_0(1-r^2/a^2)^{\alpha_n}\right) \ dr$$

In the form of volume average density where the volume integrated density function has to be divided by the volume of the cylinder / torus we get:

$$\langle n \rangle = \frac{4\pi^2R \int^a_0 r \left(n_0(1-r^2/a^2)^{\alpha_n}\right) \ dr}{2\pi^2Ra^2}$$

Simplify:

$$\langle n \rangle = \frac{2 \int^a_0 r \left(n_0(1-r^2/a^2)^{\alpha_n}\right) \ dr}{a^2}$$

Integrate:

$$\langle n \rangle = \frac{2}{a^2} \frac{a^2n_0}{2(\alpha_n+1)}$$

$$\therefore \langle n \rangle = \frac{n_0}{\alpha_n+1}$$

This is in agreement with the derivation from the ITER Physics Design 1989. ²

Since all parabolic profiles are of the same form this proof holds for the density and temperature.

A similar derivation is found in normalised coordinates, $\rho$ if the line averaged integral is swept around the poloidal circumference then divided by the poloidal area.

$$\langle n \rangle = \frac{2\pi \int^1_0 \rho \left(n_0(1-\rho^2)^{\alpha_n}\right) \ d\rho}{\pi\rho^2}$$

$\blacksquare$

---

The core values for the temperature and density are then calculated for the electrons and ions

$$\begin{aligned} \mathtt{te0} = T_\text{e} \times (\alpha_T+1) \\ \mathtt{ti0} = T_\text{i} \times (\alpha_T+1) \end{aligned}$$

$$\begin{aligned} \mathtt{ne0} = n_{\text{e}} \times (\alpha_n+1) \\ \mathtt{ni0} = \mathtt{dnitot} \times (\alpha_n+1) \end{aligned}$$

---

calculate_profile_factors()

The central plasma pressure is calculated from the ideal gas law.

$$p_0 = (\mathtt{ne0} \times \mathtt{te0}+\mathtt{ni0}\times \mathtt{ti0})\times (1000 \times \text{e})$$

With the coefficients used to turn the temperature from $\text{keV}$ back to Joules.

Since we multiply the density and temperature profiles together to get the pressure we can find the pressure profile factor by adding the profile factors for temperature and density. $$ \alpha_p = \alpha_n + \alpha_T $$

Pressure profile factor

The calculation of $\alpha_p$ is only valid assuming a parabolic profile case. The calculatio of $p_0$ is still true as the core values are calculated independantly for each profile type.

$p_0$ is NOT equal to $\langle p \rangle \times (1 + \alpha_p)$, but $p(\rho) = n(\rho)T(\rho)$ and $\langle p \rangle = \langle n \rangle$ $T_n$, where $T_n$ is the density-weighted temperature.

---

calculate_parabolic_profile_factors()

This is used by the stellerator module to find temperature gradients. We first need to find the normalised position in the plasma where the change in the pressure gradient is at a minima or maxima. This is just when the second derivative of the profile is zero. Using the density profile as an example, though this is true for temperature also.

$$T(\rho) = n_{0} \left(1 - {\rho}^{2}\right)^{\alpha_T}$$

$$\frac{{\mathrm{d}T}}{\mathrm{d}\rho} = -2T_{0} {\alpha_T}{\rho} \left(1 - {\rho}^{2}\right)^{{\alpha_T} - 1}$$

$$\frac{\mathrm{d}^2T}{\mathrm{d}\rho^2} = \frac{2T_{0} {\alpha_T} \left(1 - {\rho}^{2}\right)^{\alpha_T} \left(\left(2{\alpha_T} - 1\right) {\rho}^{2} - 1\right)}{\left({\rho}^{2} - 1\right)^{2}}$$

The positive root for the 2^nd derivative above:

$$\mathtt{rho\_te\_max} = \frac{1}{\sqrt{2{\alpha_T} - 1}}$$

Substituting the root above into the first derivative as $\rho$ to get the minima or maxima gradient value:

$$-2T_0\alpha_T\left(\frac{1}{\sqrt{2{\alpha_T} - 1}}\right)\left(1-\left(\frac{1}{\sqrt{2{\alpha_T} - 1}}\right)^2\right)^{\alpha_T-1}$$

$$\frac{2T_{0} {\alpha_T} \left(1 - \frac{1}{2{\alpha_T} - 1}\right)^{{\alpha_T} - 1}}{\sqrt{2{\alpha_T} - 1}}$$

Simplify fully:

$$\mathtt{dtdrho\_max}= -T_0\alpha_T(2^{\alpha_T})(-1+\alpha_T)^{(-1+\alpha_T)} (-1+2\alpha_T)^{(-0.5-\alpha_T)}$$

Find the temperature on the profile where the gradient change is largest by substituting into the original profile function:

$$\mathtt{te\_max} = T_0\left(1-\mathtt{rho\_te\_max}^2\right)^{\alpha_T}$$

This solution for the profile gradient only holds true if $\alpha_T \ge 1$.

In the region $0 \le \alpha_T \le 1$ when substituting the roots of the second derivative into the first derivative the function diverges into the complex number solution space.

To overcome this we can assume the second derivative root to be a value. In this case we assume a default value of $\mathtt{rho\_te\_max}$ = 0.9. Then we just substitute this value into the first derivative to

$$-2T_0\alpha_T(0.9)(1-(0.9)^2)^{-1+\alpha_T}$$

The value of $\mathtt{rho\_te\_max}$ can be manually changed in the code depending on the derivative profile required. You can use the slider in the graph below to experiment with the value of $\mathtt{rho\_te\_max}$ and how this changes the profile solutions.

We can now set the normalised gradient length:

$$\mathtt{gradient\_length\_te} = -(\mathtt{dtdrho\_max})a \times \frac{\mathtt{rho\_te\_max}}{\mathtt{te\_max}}$$

Temperature & Density Case

Calculations are identical for the temperature and density case, just with the substitution of ne and te in the variable names

---

H-mode profile | ipedestal == 1

pedestal_parameterisation()

The density and temperature profile runner function TProfile/NProfile.run() is firstly called to re-calculate the profile and core values.

Perform integrations to calculate ratio of density-weighted to volume-averaged temperature, etc. Density-weighted temperature = $\frac{\int{nT \ dV}}{\int{n \ dV}}$, which is approximately equal to the ratio $\frac{\int{\rho \ n(\rho) T(\rho) \ d\rho}}{\int{\rho \ n(\rho) \ d\rho}}$

Density weighted temperatures are thus set as such

$$\mathtt{ten} = \frac{\int_0^1{\rho \ n(\rho) T(\rho) \ d\rho}}{\int_0^1{\rho \ n(\rho) \ d\rho}} \\ \mathtt{tin} = \mathtt{ten}\times \frac{T_\text{i}}{T_\text{e}}$$

The above is done numerically with scipy.integrate.simpson integration.

Set profile factor which is the ratio of density-weighted to volume-averaged temperature

$$\mathtt{pcoef} = \frac{\mathtt{ten}}{T_\text{e}}$$

Calculate the line averaged electron density by integrating the normalised profile using the class integrate_profile_y() function

$$\mathtt{dnla} = \int_0^1{n(\rho) \ d\rho}$$

A divertor variable prn1 is set to be equal to the separatrix density over the mean density:

$$\mathtt{prn1} = \frac{n_{\text{e,sep}}}{n_{\text{e}}}$$

---

calculate_profile_factors()

The same function is run from the ipedestal == 0 profile case, found here

---

Key Constraints

---

Setting pedestal values as fractions of the Greenwald limit

By default, the values of $n_{\text{ped}}$ and $n_{\text{sep}}$ are set as fractions of the Greenwald limit such as:

$$n_{\text{ped}} = \mathtt{fgwped} \times \frac{I_p [\text{A}]}{\pi a^2 [\text{m}^2]} \times 10^{14}$$

$$n_{\text{sep}} = \mathtt{fgwsep} \times \frac{I_p [\text{A}]}{\pi a^2 [\text{m}^2]} \times 10^{14}$$

To set the values of $n_{\text{ped}}$ and $n_{\text{sep}}$ directly, the user can input the value of $\mathtt{fgwped}$ or $\mathtt{fgwsep}$ to be less than 0.0 (i.e negative) to prevent the Greenwald fraction value being set.

$\mathtt{fgwped}$ and $\mathtt{fgwsep}$ can be set as iteration variables respectively by using ixc = 45 and ixc = 152 respectively

---

Pedestal Density Upper limit

This constraint can be activated by stating icc = 81 in the input file

To prevent unrealistic profiles when iterating the values of ne0, nesep, neped etc, this constraint ensures that the value of ne0 is always higher that neped to get a converged solution. This can be scaled with fne0

---

Eich Critical Separatrix Density Model

This constraint can be activated by stating icc = 76 in the input file

Eich et.al has given a usable formula for the critical separatrix density in terms of fundamental plasma parameters⁴. A function for the ballooning parameter at the separatrix $(\alpha_{\text{sep}}^{\text{crit}})$ is seen to increase about linearly with the separatrix density normalised to Greenwald density. This is seen in JET and ASDEX based on Thomson-scattering measurements, a clear correlation of the density limit of the tokamak H-mode high-confinement regime with the approach to the ideal ballooning instability threshold at the periphery of the plasma.

$$\alpha_{\text{sep}}^{\text{crit}} \approx \kappa^{1.2}(1+1.5\delta) \approx 2.5$$

$$\alpha_{\text{sep}}^{\text{crit}} = \kappa^{1.2}(1+1.5\delta)$$

Where $\kappa$ is the plasma elongation and $\delta$ is the plasma triangularity.

$$n_{\text{sep}}^{\text{crit}} = 5.9 \times \alpha_{\text{crit}} A^{-\frac{2}{7}}\left(\frac{1+\kappa^2}{2}\right)^{-\frac{6}{7}} P_{\text{sep}}^{-\frac{11}{70}} \times n_{\text{GW}}$$

Where $A$ is the plasma aspect ratio, $P_{\text{sep}}$ is the power crossing the separatrix and is the heating power minus the radiated power in this case, $n_{\text{GW}}$ is the Greenwald density limit.

The value of nesep is check to make sure that it does not go higher than $n_{\text{sep}}^{\text{crit}}$. This can be scaled with fnesep

---

1. M. Bernert et al. Plasma Phys. Control. Fus. 57 (2015) 014038 ↩

2. N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989', ↩↩↩↩

3. Johner Jean (2011) HELIOS: A Zero-Dimensional Tool for Next Step and Reactor Studies, Fusion Science and Technology, 59:2, 308-349, DOI: 10.13182/FST11-A11650 ↩

4. T. Eich et al 2018 Nucl. Fusion, 'Correlation of the tokamak H-mode density limit with ballooning stability at the separatrix' 58 034001 ↩