diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/NBI/culham_nb.md b/documentation/proc-pages/eng-models/heating_and_current_drive/NBI/culham_nb.md new file mode 100644 index 0000000000..40ecbd1148 --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/NBI/culham_nb.md @@ -0,0 +1,223 @@ +# Culham Neutral Beam Model | `culnbi()` + +- `iefrf/iefrffix` = 8 + + + +This routine calculates Neutral Beam current drive parameters +using the corrections outlined in AEA FUS 172 to the ITER method. +The result cannot be guaranteed for devices with aspect ratios far +from that of ITER (approx. 2.8). + +| Output | Description | +|----------|-------------| +| $\mathtt{effnbss}$ | Neutral beam current drive efficiency in Amperes per Watt | +| $\mathtt{fpion}$ | Fraction of NB power given to ions | +| $\mathtt{fshine}$ | Shine-through fraction of the beam | + +$$ +\mathtt{frbeam} = \frac{R_{\text{tan}}}{R_0} +$$ + +Where $R_{\text{tan}}$ is major radius at which the centre-line of the beam is tangential to the toroidal direction. This can be user defined + +$$ +\left(1+ \frac{1}{A}\right) < \mathtt{frbeam} +$$ + +A quick sanity check is done to make sure no negative roots are formed when calculating $\mathtt{dpath}$ this prevents setups where the NBI beam would miss the plasma + + +$$ +\mathtt{dpath} = R_0 \sqrt{\left(1+\frac{1}{A}\right)^2-\mathtt{frbeam}^2} +$$ + +Beams topping cross section is calculated via $\mathtt{sigbeam}$ found [here](../NBI/nbi_overview.md/#beam-stopping-cross-section-sigbeam). This produces $\mathtt{sigstop}$ + +Calculate number of decay lengths to centre + +$$ +\mathtt{taubeam} = \mathtt{dpath} \times n_{\text{e,0}} \times \mathtt{sigstop} +$$ + +Calculate the shine through fraction of the beam + +$$ +\mathtt{fshine} = e^{\left(-2 \times \mathtt{dpath} \times n_{\text{e,0}} \times \mathtt{sigstop}\right)} +$$ + +Deuterium and tritium beam densities + +$$ +\mathtt{dend} = n_{\text{ion}} \times (1-\mathtt{ftritbm}) +$$ + +$$ +\mathtt{dent} = n_{\text{ion}} \times \mathtt{ftritbm} +$$ + +Power split to the ions and electrons is clauclated with the $\mathtt{cfnbi()}$ method found [here](../NBI/nbi_overview.md/#ion-coupled-power-cfnbi) and outputs $\mathtt{fpion}$ + +## Current drive efficiency | `etanb2()` + +This routine calculates the current drive efficiency in A/W of +a neutral beam system, based on the 1990 ITER model, +plus correction terms outlined in Culham Report AEA FUS 172. + +| Input | Description | +| :---------- | :----------------------------------- | +| $\mathtt{abeam}$ | beam ion mass (amu) | +| $\mathtt{alphan}$, $\alpha_n$ | density profile factor | +| $\mathtt{alphat}$, $\alpha_T$ | temperature profile factor | +| $\mathtt{aspect}$, $A$ | aspect ratio | +| $\mathtt{dene}$, $n_{\text{e}}$ | volume averaged electron density $(\text{m}^{-3})$ | +| $\mathtt{dnla}$, $n_{\text{e,0}}$ | line averaged electron density $(\text{m}^{-3})$ | +| $\mathtt{enbeam}$ | neutral beam energy $(\text{keV})$ | +| $\mathtt{frbeam}$ | R_tangent / R_major for neutral beam injection | +| $\mathtt{fshine}$ | shine-through fraction of beam | +| $\mathtt{rmajor}$, $R$ | plasma major radius $(\text{m})$ | +| $\mathtt{rminor}$, $a$ | plasma minor radius $(\text{m})$ | +| $\mathtt{ten}$ | density weighted average electron temperature $(\text{keV})$ | +| $\mathtt{zeff}$, $Z_{\text{eff}}$ | plasma effective charge | + + +Charge of beam ions +$$ +\mathtt{zbeam} = 1.0 +$$ + +Fitting factor (IPDG89) + +$$ +\mathtt{bbd} = 1.0 +$$ + +Volume averaged electron density ($10^{20} \text{m}^{-3}$) + +$$ +\mathtt{dene20} = n_{\text{e,20}} +$$ + +Line averaged electron density ($10^{20} \text{m}^{-3}$) + +$$ +\mathtt{dnla20} = n_{\text{(e,0) 20}} +$$ + +Critical energy ($\text{MeV}$) (power to electrons = power to ions) (IPDG89) +N.B. ten is in keV + +$$ +\mathtt{ecrit} = 0.01 \times \mathtt{abeam} \times \mathtt{ten} +$$ + +Beam energy in MeV + +$$ +\mathtt{ebmev} = \frac{\mathtt{enbeam}}{10^3} +$$ + +x and y coefficients of function J0(x,y) (IPDG89) + +$$ +\mathtt{xjs} = \frac{\mathtt{ebmev}}{\mathtt{bbd}\times \mathtt{ecrit}} +$$ + +$$ +\mathtt{xj} = \sqrt{\mathtt{xjs}} +$$ + +$$ +\mathtt{yj} = \frac{0.8 \times Z_{\text{eff}}}{\mathtt{abeam}} +$$ + +Fitting function J0(x,y) + +$$ +\mathtt{j0} = \frac{xjs}{(4.0 + 3.0 \times \mathtt{yj} + \mathtt{xjs} \times (\mathtt{xj} + 1.39 + 0.61 \times yj^{0.7}))} +$$ + +Effective inverse aspect ratio, with a limit on its maximum value + +$$ + \mathtt{epseff} = \text{min}(0.2, (0.5 / A)) +$$ + +Reduction in the reverse electron current +due to neoclassical effects + +$$ +\mathtt{gfac} = (1.55 + 0.85 / Z_{\text{eff}}) \times \sqrt{\mathtt{epseff}} - (0.2 + 1.55 / Z_{\text{eff}}) \times \mathtt{epseff} +$$ + +Reduction in the net beam driven current +due to the reverse electron current + +$$ +\mathtt{ffac} = 1.0 - \frac{\mathtt{zbeam}}{Z_{\text{eff}}} \times (1.0 - \mathtt{gfac}) +$$ + +Normalisation to allow results to be valid for +non-ITER plasma size and density: + +Line averaged electron density ($10^{20} \text{m}^{-3}$) normalised to ITER + +$$ +\mathtt{nnorm} = 1.0 +$$ + +Distance along beam to plasma centre + +$$ +\mathtt{r} = \text{max}(R, R \times \mathtt{frbeam}) +$$ + +$$ +\mathtt{eps1} = a / \mathtt{r} +$$ + + +$$ +\mathtt{d} = R \times \sqrt{((1.0 + \mathtt{eps1})^2 - \mathtt{frbeam}^2)} +$$ + +Distance along beam to plasma centre for ITER +assuming a tangency radius equal to the major radius + +$$ +\mathtt{epsitr} = 2.15 / 6.0 +$$ + +$$ +\mathtt{dnorm} = 6.0 \times \sqrt{(2.0 \times \mathtt{epsitr} + \mathtt{epsitr}^2)} +$$ + +Normalisation to beam energy (assumes a simplified formula for +the beam stopping cross-section) + +$$ + \mathtt{ebnorm} = \mathtt{ebmev} \times ((\mathtt{nnorm} \times \mathtt{dnorm}) / (n_{\text{e,0}} \times \mathtt{d})) ^{1.0 / 0.78)} +$$ + +A_bd fitting coefficient, after normalisation with ebnorm + +$$ +\mathtt{abd} = ( + 0.107 + \times (1.0 - 0.35 \times \alpha_n + 0.14 \times \alpha_n^2) + \times (1.0 - 0.21 \times \alpha_T) + \times (1.0 - 0.2 \times \mathtt{ebnorm} + 0.09 \times \mathtt{ebnorm}^2) + ) +$$ + +Normalised current drive efficiency ($\text{A/W} \text{m}^{2}$) (IPDG89) + +$$ +\mathtt{gamnb} = 5.0 \times \mathtt{abd} \times 0.1 \times \mathtt{ten} \times (1.0 - \mathtt{fshine}) \times \mathtt{frbeam} \times \frac{\mathtt{j0}}{0.2} \times \mathtt{ffac} +$$ + +Current drive efficiency (A/W) + +$$ +\text{Current drive efficiency [A/W]} = \frac{\mathtt{gamnb}}{\mathtt{dene20}\times R} +$$ diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/NBI/iter_nb.md b/documentation/proc-pages/eng-models/heating_and_current_drive/NBI/iter_nb.md new file mode 100644 index 0000000000..ee91af1abf --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/NBI/iter_nb.md @@ -0,0 +1,128 @@ +# ITER Neutral Beam Model | `iternb()` + +- `iefrf/iefrffix` = 5 + +| Output | Description | +|----------|-------------| +| $\mathtt{effnbss}$ | Neutral beam current drive efficiency in $\text{A/W}$ | +| $\mathtt{fpion}$ | Fraction of NB power given to ions | +| $\mathtt{fshine}$ | Shine-through fraction of the beam | + +This model calculates the current drive parameters for a neutral beam system, based on the 1990 ITER model.[^1] + + +Firstly the beam access is checked for such that +$$ +\bigg(1+ \frac{1}{A}\bigg) > (R_{\text{tangential}}/R_0) +$$ + +The beam path length to centre is calculated: + +$$ +\underbrace{\mathtt{dpath}}_{\text{Path length to centre}} = R_0 \sqrt{\left(\left(1+\frac{1}{A}\right)^2-\mathtt{frbeam}^2\right)} +$$ + + +Beam stopping cross-section ($\sigma_{\text{beam}}$) is calculated using the `sigbeam` method described [here](nbi_overview.md) : + + +Calculate number of electron decay lengths to centre + +$$ +\tau_{\text{beam}} = \mathtt{dpath}\times n_e \sigma_{\text{beam}} +$$ + +Shine-through fraction of beam: +$$ +f_{\text{shine}} = e^{(-2.0 \times \mathtt{dpath} \times n_e \sigma_{\text{beam}})} \\ +$$ + +Deuterium and tritium beam densities: +$$ +n_D = n_i * (1.0 - \mathtt{ftritbm}) +$$ + +$$ +n_T = n_i * \mathtt{ftritbm} +$$ + +Power split to ions / electrons is calculated via the the `cfnbi` method described [here](nbi_overview.md) + + +## Current drive efficiency | `etanb()` + +This routine calculates the current drive efficiency of +a neutral beam system, based on the 1990 ITER model. +AEA FUS 251: A User's Guide to the PROCESS Systems Code +ITER Physics Design Guidelines: 1989 IPDG89, N. A. Uckan et al, +ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 + +| Input | Description | +|---------|-----------------------------------------------------------| +| $\mathtt{abeam}$, $m_{\text{u,ion}}$ | Beam ion mass ($\text{amu}$) | +| $\mathtt{alphan}$ | Density profile factor | +| $\mathtt{alphat}$ | Temperature profile factor | +| $\mathtt{aspect}$, $A$ | Aspect ratio | +| $\mathtt{dene20}$, $n_{\text{e,20}}$ | Volume averaged electron density ($10^{20} \text{m}^{-3}$) | +| $\mathtt{ebeam}$ | Neutral beam energy ($\text{keV}$) | +| $\mathtt{rmajor}$, R | Plasma major radius ($\text{m}$) | +| $\mathtt{ten}$ | Density weighted average electron temperature ($\text{keV}$) | +| $\mathtt{zeff}$, $Z_{\text{eff}}$ | Plasma effective charge | + +| Output | Description | +|---------|-----------------------------------------------------------| +| $\mathtt{etanb}$ | Neutral beam current drive efficiency in $\text{A/W}$ | + + + + +$$ +\mathtt{zbeam} = 1.0 +$$ + +$$ +\mathtt{bbd} = 1.0 +$$ + + +Ratio of E_beam/E_crit + +$$ +\mathtt{xjs} = \frac{\mathtt{ebeam}}{10 \ m_{\text{u,ion}} \ T_e} +$$ + +$$ +\mathtt{xj} = \sqrt{\mathtt{xjs}} +$$ + +$$ +\mathtt{yj} = 0.8 \frac{Z_{\text{eff}}}{m_{\text{u,ion}}} +$$ + +$$ +\mathtt{rjfunc} = \frac{\mathtt{xjs}}{((4.0 + 3.0\mathtt{yj} + \mathtt{xjs}(\mathtt{xj} + 1.39 + 0.61\mathtt{yj}^{0.7})))} +$$ + +$$ +\mathtt{epseff} = \frac{0.5}{A} +$$ + +$$ +\mathtt{gfac} = \left(1.55 + \frac{0.85}{Z_{\text{eff}}}\right)\left(\sqrt{\mathtt{epseff}}-\left(0.2+\frac{1.55}{Z_{\text{eff}}}\right)\mathtt{epseff}\right) +$$ + +$$ +\mathtt{ffac} = \frac{1}{\mathtt{zbeam}} - \frac{(1-\mathtt{gfac})}{Z_{\text{eff}}} +$$ + +$$ +\mathtt{abd} = 0.107 (1-0.35 \ \mathtt{alphan}+0.14 \ \mathtt{alphan}^2)(1-0.21 \ \mathtt{alphat})(1-0.2\times 10^{-3}\mathtt{ebeam}+0.09\times 10^{-6} \ \mathtt{ebeam}^2) +$$ + +$$ +\text{Current drive efficiency [A/W]} = \mathtt{abd} \times\frac{5}{R_0} \times0.1\frac{\mathtt{ten}}{n_{\text{e},20}} \times \frac{\mathtt{rjfunc}}{0.2}\mathtt{ffac} +$$ + + + +[^1]: N. A. Uckan and ITER Physics Group, *"ITER Physics Design Guidelines: 1989"*, ITER Documentation Series, No. 10, IAEA/ITER/DS/10 (1990) \ No newline at end of file diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/NBI/nbi_overview.md b/documentation/proc-pages/eng-models/heating_and_current_drive/NBI/nbi_overview.md new file mode 100644 index 0000000000..9ececc967e --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/NBI/nbi_overview.md @@ -0,0 +1,165 @@ +# Neutral Beam Injection Heating + +!!! Warning "NBI Models" + At present, the neutral beam models do not include the effect of an edge transport barrier (pedestal) in the plasma profile. + +## Neutral beam access + +If present, a neutral beam injection system needs sufficient space between the TF coils to be able to intercept the plasma tangentially. The major radius `rtanbeam` at which the centre-line of the beam is tangential to the toroidal direction is user-defined using input parameter `frbeam`, which is the ratio of `rtanbeam` to the plasma major radius `rmajor`. The maximum possible tangency radius `rtanmax` is determined by the geometry of the TF coils - see Figure 1, and this can be enforced using constraint equation no. 20 with iteration variable no. 33 (`fportsz`). The thickness of the beam duct walls may be set using input parameter `nbshield`. + +
+![NBI Port Size](../images/portsize.png){ width = "300"} +
Figure 1: Top-down schematic of the neutral beam access geometry. The beam with the maximum possible tangency radius is shown here.
+
+ +## Neutral beam losses + +Input parameter `forbitloss` can be used to specify the fraction of the net injected neutral beam power that is lost between the beam particles' ionisation and thermalisation (known as the first orbit loss). This quantity cannot easily be calculated as it depends on the field ripple and other three-dimensional effects. The power lost is assumed to be absorbed by the first wall. + +The power in the beam atoms that are not ionised as they pass through the plasma (shine-through) is calculated by the code. There are two constraint equations that can be used to control the beam penetration and deposition, as follows: + +- It is necessary to use a beam energy that simultaneously gives adequate penetration of the beam to the centre of the plasma and tolerable shine-through of the beam on the wall after the beam has traversed the plasma. The number of exponential decay lengths, $\tau$, for the beam power to fall before it reaches the plasma centre should be in the region of ~ 4-6[^2],. Constraint equation no. 14 may be used to force $\tau$ to be equal to the value given by input parameter `tbeamin`, and is therefore in effect a beam energy consistency equation. +- Alternatively, constraint equation no. 59 with iteration variable no. 105 (`fnbshineef`) may be used to ensure that the beam power fraction emerging from the plasma is no more than the value given by input parameter `nbshinefmax`. + +It is recommended that only one of these two constraint equations is used during a run. + + +## Beam stopping cross-section | `sigbeam()` + +| Input | Description | +| :---------- | :----------------------------------- | +| $\mathtt{eb}$ | Beam energy $\left(\text{keV}/\text{amu}\right)$ | +| $\mathtt{te}$, $T_{\text{e}}$ | Electron temperature $\left(\text{keV}\right)$ | +| $\mathtt{ne}$, $n_{\text{e}}$ | Electron density $\left(10^{20}\text{m}^{-3}\right)$ | +| $\mathtt{rnhe}$ | Alpha density / $n_{\text{e}}$ | +| $\mathtt{rnc}$, | Carbon density /$n_{\text{e}}$ | +| $\mathtt{rno}$, | Oxygen density /$n_{\text{e}}$ | +| $\mathtt{rnfe}$ | Iron density /$n_{\text{e}}$ | + + +Both the [ITER](./iter_nb.md) and [Culham](culham_nb.md) NBI models both use the `sigbeam` method to calculate the stopping cross section[^1]. It finds a suitable analytic epressing for $\sigma_s^{(Z)}(E,n_{\text{e}},T_{\text{e}},Z_{\text{eff}})$ for fitting $\sigma_s$ data for a single impurity $(\text{Z)}$ plasma: + + + +$$ +\sigma_s^{(Z)}(E,n_{\text{e}},T_{\text{e}},Z_\text{eff}) = \frac{e^{[S_{1} (E,n_{\text{e}}, T_{\text{e}})]}}{E} \times\left [1 +(Z_\text{eff}-1) S_z(E, n_{\text{e}}, T_{\text{e}})\right] \ (\times 10^{-16} \text{cm}^2) +$$ + +where + +$$ +S_{1} = \sum_{i=1}^2 \sum_{j=1}^3 \sum_{k=1}^2 \ \{A_{ijk} \times (\ln E)^{i-1} \ [\ln(n/n_{0})]^{j-1} \ (\ln T_{\text{e}})^{k-1} \} +$$ + +$$ +S_{Z} = \sum_{i=1}^3 \sum_{j=1}^2 \sum_{k=1}^2 \ \{B_{ijk}^{(z)} \times (\ln E)^{i-1} \ [\ln(n/n_{0})]^{j-1} \ (\ln T_{\text{e}})^{k-1} \} +$$ + +with $E, n_e, T$ expressed in units of keV/u, $\text{cm}^3$ and keV, respectively, and $n_0 = 10^{13} \text{cm}^3$. The function $S_1 (E, n_{\text{e}}, T_{\text{e}})$ together with the $E^{-1}$ factor describes the beam stopping in a pure hydrogenic plasma, while the function $(Z_{\text{eff}}- 1)\ S_z (E, n_e, T_e)$ describes the effect of the impurity $Z$ on the beam stopping. + +!!! info "Info" + For the full table of values for $A_{ijk}$ & $B_{ijk}^{(z)}$\ please see the accompanying paper[^1] or `current_drive.py` + +For a plasma having an arbitrary mix of $N$ different types of impurities with densities $n$, and charges $Z_q (q = 1, ..., N)$, the beam stopping cross-section can be represented as the weighted sum of the stopping cross- sections for $N$ reference single-impurity plasmas. In each of these reference plasmas, the electron density and the proton density (including that of deuterium and tritium ions) are the same as in a true plasma. The impurity density, however, is increased in order to satisfy quasi-neutrality. The weighting function is the electron density $n_qZ_q$ associated with the amu impurity (in the true plasma), divided by the sum of these densities. The result is: + +$$ +\sigma_s^{(N)}=\frac{ e^{S_{1}(E, n_{\text{e}}, T_{\text{e}})}}{E} \times\left[1+\frac{1}{n_{\text{e}}} \sum_q n_q Z_q(Z_q-1) S_{Z_q}(E, n_{\text{e}}, T_{\text{e}})\right] +(\times 10^{-16} \mathrm{~cm}^2) +$$ + + +## Ion coupled power | `cfnbi()` +Both the [ITER](./iter_nb.md) and [Culham](culham_nb.md) NBI models both use the `cfnbi` method to calculate the fraction of the fast particle energy coupled to the ions. + +| Input | Description | +| :---------- | :----------------------------------- | +| $\mathtt{afast}$, $m_{\text{u,fast}}$ | Mass of fast particle (units of proton mass) | +| $\mathtt{efast}$, $E_{\text{fast}}$ | Energy of fast particle ($\text{keV}$) | +| $\mathtt{te}$, $T_{\text{e}}$ | Density weighted average electron temperature ($\text{keV}$) | +| $\mathtt{ne}$, $n_{\text{e}}$ | Volume averaged electron density ($\text{m}^{-3}$) | +| $\mathtt{nd}$ | Deuterium beam density ($\text{m}^{-3}$) | +| $\mathtt{nt}$ | Tritium beam density ($\text{m}^{-3}$) | +| $\mathtt{zeffai}$ | Mass weighted plasma effective charge | +| $\mathtt{xlmbda}$ | Ion-electron Coulomb logarithm | + +### Coloumb logarithm | `xlmbdabi()` +Firstly, the Coulomb logarithm for the ion-ion collisions ($\ln \Lambda$) is calculated using the `xlmbdabi` method [^2]. The calculation follows these steps: + +1. Calculate $x_1$ using the formula: + $$ + x_1 = \frac{T_{\text{e}}}{10} E_{\text{fast}} \frac{m_{\text{u,fast}}}{n_{\text{e}}} + $$ + where $T_{\text{e}}$ is the density-weighted average electron temperature (in $\text{keV}$), $E_{\text{fast}}$ is the energy of the fast particle (in $\text{MeV}$), $m_{\text{u,fast}}$ is the mass of the fast particle (in units of proton mass), and $n_{\text{e}}$ is the volume-averaged electron density ($1\times 10^{20}/ \text{m}^3$). + +2. Calculate $x_2$ using the formula: + $$ + x_2 = \frac{m_{\text{u,ion}}}{m_{\text{u,ion}} + m_{\text{u,fast}}} + $$ + where $m_{\text{u,ion}}$ is the mass of the background ions (in units of proton mass). + +3. Calculate $\ln \Lambda$ using the formula: + $$ + \ln \Lambda = 23.7 + \log(x_2 \sqrt{x_1}) + $$ + +--- + + + +--- + +Next, the following calculations are performed: + +1. Calculate `sumln` using the formula: + $$ + \mathtt{sumln} = \mathtt{zeffai} \times \frac{\mathtt{xlmbdai}}{\mathtt{xlmbda}} + $$ + where $\mathtt{zeffai}$ is the mass-weighted plasma effective charge and $\ln\Lambda_{\text{i-e}}$ is the ion-electron Coulomb logarithm. + +2. Calculate `xlnrat` using the formula: + $$ + \mathtt{xlnrat} = 3 \left(\frac{\sqrt \pi}{4}\frac{m_{\text{e}}}{m_{\text{p}} \ \mathtt{sumln}}\right)^{\frac{2}{3}} + $$ + where $m_{\text{e}}$ is the electron mass and $m_{\text{p}}$ is the proton mass. + +3. Calculate $v_{\text{e}}$ using the formula: + $$ + v_{\text{e}} = c \sqrt{\left(\frac{2 \ T_{\text{e}}}{511}\right)} + $$ + where $c$ is the speed of light. + +4. Calculate `ecritfi` using the formula: + $$ + \mathtt{ecritfi} = \frac{m_{\text{u,fast}} m_{\text{p}} v_{\text{e}}^2 \ \mathtt{xlnrat}} {(2000 \times e)} + $$ + where $e$ is the elementary charge. + +5. Calculate `x` using the formula: + $$ + \mathtt{x} = \sqrt{\frac{\mathtt{efast}}{\mathtt{ecritfi}}} + $$ + +6. Calculate `ti` using the formula: + + $$ + \mathtt{ti} = \log{\frac{x^2-x+1}{(x+1)^2}} + $$ + +7. Calculate `thx` using the formula: + $$ + \mathtt{thx} = \frac{(2x-1)}{\sqrt{3}} + $$ + +8. Calculate `t2` using the formula: + $$ + \mathtt{t2} = 2 \sqrt{3} \arctan({\mathtt{thx}})+\frac{\pi}{6} + $$ + +Finally, the fraction of fast particle energy coupled to ions is calculated using the formula: +$$ +\text{Fraction of fast particle energy coupled to ions} = \frac{(\mathtt{t1}+\mathtt{t2})}{3\mathtt{x}^2} +$$ + +[^1]:Janev, R. K., Boley, C. D., & Post, D. E. (1989). *"Penetration of energetic neutral beams into fusion plasmas."* Nuclear Fusion, 29(12), 006. https://doi.org/10.1088/0029-5515/29/12/006 + +[^2]: David R. Mikkelsen & Clifford E. Singer (1983) *"Optimization of Steady-State Beam-Driven Tokamak Reactors"*, Nuclear Technology - Fusion, 4:2P1, 237-252, DOI: 10.13182/FST83-A22816 \ No newline at end of file diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/RF/culham_electron_cyclotron.md b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/culham_electron_cyclotron.md new file mode 100644 index 0000000000..204f6fd053 --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/culham_electron_cyclotron.md @@ -0,0 +1,64 @@ +# Culham Electron Cyclotron Model | `culecd()` + + +- `iefrf/iefrffix` = 7 + +This routine calculates the current drive parameters for a electron cyclotron system, based on the AEA FUS 172 model[^1] + + +1. Local electron temperature $(\mathtt{tlocal})$ is calculated using the `tprofile` method +2. Local electron density $(\mathtt{dlocal})$ is calculated using the `nprofile` method + +3. Calculate the inverse aspect ratio `epsloc`. + +$$ +\mathtt{epsloc} = \frac{1}{3A} +$$ + +4. Calculate the Coulomb logarithm for ion-electron collisions `coulog`.[^2] + +$$ +\mathtt{coulog} = 15.2 - 0.5\log({\mathtt{dlocal}}) + \log({\mathtt{tlocal}}) +$$ + +Calculate normalised current drive efficiency at four different poloidal angles, and average. +cosang = cosine of the poloidal angle at which ECCD takes place = +1 outside, -1 inside. + +## Normalised current drive efficiency + +Uses the `eccdef` model found [here](ec_overview.md) + +5. Calculate the normalised current drive efficiency at four different poloidal angles, and average. + - Set `cosang` to 1.0 and calculate `ecgam1`. + - Set `cosang` to 0.5 and calculate `ecgam2`. + - Set `cosang` to -0.5 and calculate `ecgam3`. + - Set `cosang` to -1.0 and calculate `ecgam4`. + + cosang = 1.0e0 + ecgam1 = self.eccdef(tlocal, epsloc, zlocal, cosang, coulog) + cosang = 0.5e0 + ecgam2 = self.eccdef(tlocal, epsloc, zlocal, cosang, coulog) + cosang = -0.5e0 + ecgam3 = self.eccdef(tlocal, epsloc, zlocal, cosang, coulog) + cosang = -1.0e0 + ecgam4 = self.eccdef(tlocal, epsloc, zlocal, cosang, coulog) + +6. Calculate the normalised current drive efficiency `ecgam` as the average of `ecgam1`, `ecgam2`, `ecgam3`, and `ecgam4`. + +$$ +\mathtt{ecgam} = 0.25(\mathtt{ecgam1} + \mathtt{ecgam2} +\mathtt{ecgam3} + \mathtt{ecgam4}) + $$ + +7. Calculate the current drive efficiency by dividing `ecgam` by `(dlocal * physics_variables.rmajor)`. + +$$ +\text{Current drive efficiency [A/W]} = \frac{\mathtt{ecgam}}{\mathtt{dlocal} \times R_0} +$$ + +Note: The `eccdef` method is called to calculate the current drive efficiency at each poloidal angle. + + + +[^1]: Hender, T.C., Bevir, M.K., Cox, M., Hastie, R.J., Knight, P.J., Lashmore-Davies, C.N., Lloyd, B., Maddison, G.P., Morris, A.W., O’Brien, M.R. and Turner, M.F., 1992. *"Physics assessment for the European reactor study."* AEA FUS, 172. + +[^2]: Wesson, J. and Campbell, D.J., 2011. *"Tokamaks"* (Vol. 149). Oxford university press. \ No newline at end of file diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/RF/culham_lower_hybrid.md b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/culham_lower_hybrid.md new file mode 100644 index 0000000000..60325886d5 --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/culham_lower_hybrid.md @@ -0,0 +1,50 @@ +# Culham Lower Hybrid | `cullhy()` + +- `iefrf/iefrffix` = 6 + + + +This routine calculates the current drive parameters for a +lower hybrid system, based on the AEA FUS 172 model. +AEA FUS 251: A User's Guide to the PROCESS Systems Code +AEA FUS 172: Physics Assessment for the European Reactor Study[^1] + + +1. Call the [`lhrad()`](#lower-hybrid-wave-absorption-radius--lhrad) method to calculate the lower hybrid wave absorption radius, `rratio`. +2. Calculate the penetration radius, `rpenet`, by multiplying `rratio` with the minor radius of the plasma. +3. Calculate the local density, `dlocal`, using the `nprofile()` function from the `profiles_module` module. This function takes into account various plasma parameters such as the density profile, electron density at the edge, pedestal density, separatrix density, and the value of the parameter `alphan`. +4. Similarly, calculate the local temperature, `tlocal`, using the `tprofile()` function from the `profiles_module` module. This function considers parameters such as the temperature profile, electron temperature at the edge, pedestal temperature, separatrix temperature, `alphat`, and `tbeta`. +5. Calculate the local toroidal magnetic field, `blocal`, using the formula `bt * rmajor / (rmajor - rpenet)`. Here, `bt` is the toroidal magnetic field at the magnetic axis, and `rmajor` is the major radius of the plasma. +6. Calculate the parallel refractive index, `nplacc`, which is needed for plasma access. It uses the local density `dlocal` and the local magnetic field `blocal` to calculate a fraction `frac`. `nplacc` is then obtained by adding `frac` to the square root of `1.0 + frac * frac`. +7. Calculate the local inverse aspect ratio, `epslh`, by dividing `rpenet` by `rmajor`. +8. Calculate the LH normalized efficiency, `x`, using the formula `24.0 / (nplacc * sqrt(tlocal))`. +9. Calculate several intermediate terms, `term01`, `term02`, `term03`, and `term04`, using different formulas involving `nplacc`, `physics_variables.zeff`, `tlocal`, `epslh`, and `x`. +10. Calculate the current drive efficiency, `gamlh`, using the formula `term01 * term02 * (1.0e0 - term03 / term04)`. +11. Return the current drive efficiency normalized by the product of `0.1e0 * dlocal` and `physics_variables.rmajor`. + +[^1]: T. C. Hender, M. K. Bevir, M. Cox, R. J. Hastie, P. J. Knight, C. N. Lashmore-Davies, B. Lloyd, G. P. Maddison, A. W. Morris, M. R. O'Brien, M.F. Turner abd H. R. Wilson, *"Physics Assessment for the European Reactor Study"*, AEA Fusion Report AEA FUS 172 (1992) + +## Lower Hybrid wave absorption radius | `lhrad`() + +This routine determines numerically the minor radius at which the damping of Lower Hybrid waves occurs, using a Newton-Raphson method to establish the correct minor radius ratio. The required minor radius ratio has been found when the difference between the results of the two formulae for the energy E given in AEA FUS 172, p.58 is sufficiently close to zero. + +Correction to refractive index (kept within valid bounds) + $\mathtt{drfind} = \min\left(0.7, \max\left(0.1, \frac{12.5}{\text{te0}}\right)\right)$ + +Use Newton-Raphson method to establish the correct minor radius ratio. The required minor radius ratio has been found when the difference between the results of the two formulae for the energy E given in AEA FUS 172, p.58 is sufficiently close to zero. + +Iterate over the following steps to find the minor radius ratio: + +1. Set an initial guess for the minor radius ratio, $\mathtt{rat0}$, to 0.8. +2. Repeat the following steps for a maximum of 100 iterations: + - Calculate the minor radius ratios, $r1$ and $r2$, by subtracting and adding 0.1% of $\mathtt{rat0}$, respectively. + - Evaluate the function $g$ at $\mathtt{rat0}$, $r1$, and $r2$ using the method `lheval(drfind, rat)`. + - Calculate the gradient of $g$ with respect to the minor radius ratio, $\frac{{dg}}{{dr}}$, using the formula $\frac{{g2 - g1}}{{r2 - r1}}$. + - Calculate a new approximation for the minor radius ratio, $\mathtt{rat1}$, using the formula $\mathtt{rat0} - \frac{{g0}}{{\frac{{dg}}{{dr}}}}$. + - Ensure that $\mathtt{rat1}$ is within the bounds of 0.0001 and 0.9999. + - If the absolute value of $g0$ is less than or equal to 0.01, exit the loop. + - Update $\mathtt{rat0}$ with the new approximation, $\mathtt{rat1}$. +3. If the loop completes all 100 iterations without finding a satisfactory solution, report an error and set $\mathtt{rat0}$ to 0.8. +4. Return the final value of $\mathtt{rat0}$. + +[^1]: T. C. Hender, M. K. Bevir, M. Cox, R. J. Hastie, P. J. Knight, C. N. Lashmore-Davies, B. Lloyd, G. P. Maddison, A. W. Morris, M. R. O'Brien, M.F. Turner abd H. R. Wilson, *"Physics Assessment for the European Reactor Study"*, AEA Fusion Report AEA FUS 172 (1992) \ No newline at end of file diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/RF/cutoff_ecrh.md b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/cutoff_ecrh.md new file mode 100644 index 0000000000..1989210b4e --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/cutoff_ecrh.md @@ -0,0 +1,63 @@ +# ECRH with Cutoff + +- `iefrf/iefrffix` = 13 + +| Input | Description | +|-------|-------------| +| `dene`, $n_{\text{e}}$ | Avergae electron temperature $\left[10^{19}\text{m}^{-3}\right]$ | +| `te`, $T_{\text{e}}$ | Avergae electron temperature $\left[\text{keV}\right]$ | +| `rmajor`, $R_0$ | Major radius $\left[\text{m}\right]$ | +| `bt`, $B_{\text{T}}$ | Toroidal magnetic field $\left[\text{T}\right]$ | +| `zeff`, $Z_{\text{eff}}$ | Effective charge | +| `harnum` | Harmonic number | +| `mode` | RF mode | + +---- + + + +$$ +\mathtt{fc} = \frac{\frac{1}{2\pi}eB_{\text{T}}}{m_{\text{e}}} +$$ + +$$ +\mathtt{fp} = \frac{1}{2\pi}\sqrt{\frac{n_{\text{e,19}}e^2}{m_{\text{e}}\epsilon_0}} +$$ + +Apply effective charge correction from GRAY study + +$$ +\mathtt{xi_{CD}} = 0.18\left(\frac{4.8}{2+Z_{\text{eff}}}\right) +$$ + +$$ +\mathtt{effrfss} = \frac{\mathtt{xi_{CD}}T_{\text{e}}}{3.27R_0n_{\text{e,19}}} +$$ + +For the O-mode case: + +$$ +\mathtt{f_{cutoff}} = \mathtt{fp} +$$ + +For the X-mode case: + +$$ +\mathtt{f_{cutoff}} = 0.5\left(\mathtt{fc}+\sqrt{\mathtt{harnum}\times\mathtt{fc}^2+4\mathtt{fp}^2}\right) +$$ + +Plasma coupling only occurs if the plasma cut-off is below the cyclotron harmonic +(a = 0.1). This controls how sharply the transition is reached + +$$ +\mathtt{cutoff_{factor}} = 0.5\left(1+\tanh\left({\left(\frac{2}{a}\right)((\mathtt{harnum}\times \mathtt{fc} -\mathtt{f_cutoff})/\mathtt{fp -a })}\right)\right) +$$ + +$$ +\text{Current drive efficiency [A/W]} = \mathtt{effrfss} \times \mathtt{cutoff_{factor}} +$$ + +
+![ECRH Cutoff](../images/ecrh_cutoff.png){ width = "100"} +
Figure 1: The variation in current drive efficiency as a function of toroidal magnetic field at different harmonics and modes
+
\ No newline at end of file diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ebw_freethy.md b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ebw_freethy.md new file mode 100644 index 0000000000..58c898ee41 --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ebw_freethy.md @@ -0,0 +1,77 @@ +# Electron Bernstein Wave Model + + +The aim of the function below is to produce a 'simple as possible' function for estimating the current drive for an EBW system based solely on the 0D parameters available. It is noted that the efficiency of an actual EBW current drive system depends strongly on how it is used, in particular on the location and spatial distribution of the current and the plasma shape and launcher solution. None of these are included in the below. The model has significant room for improvement, but much more work is required in order to understand how best to do this. + +Electron Bernstein waves (EBW) have the potential to be a source of relatively efficient current drive[^1] and have been predicted to be more efficient at driving current off-axis in spherical tokamaks[^2]. + +Microwave current drive efficiency is often expressed as a normalised quantity, $\xi_{CD}$, which captures its natural variation with temperature and density and is given by[^3]: + +$$ +\xi_{CD} = \frac{3.27 \ I_{\text{P}} [\text{A}] \ R [\text{m}] \ n_{\text{e}}[\times 10^{19} \text{m}^{-3}]}{T_{\text{e}}[\text{keV}] \ P[\text{W}]} +$$ + +We can rearrange this expression to give the current in terms of the other parameters: + +$$ +I_P[\text{A}] = \frac{\xi_{\text{CD}} \ T_{\text{e}}[\text{keV}] \ P[\text{W}]}{3.27 \ R[\text{m}] \ n_{\text{e}} [\times 10^{19} \text{m}^{-3}]} +$$ + + +Experiments of EBW current drive on W7-AS found $\xi_{CD}=0.43\pm0.1$ [^1]. It should be noted that this may be different depending on the plasma and the location within the plasma that the current is driven. A more conservative assumptions might be $\xi_{\text{CD}} = 0.3$, optimistic $\xi_{\text{CD}}$ = 0.5. + +Note that we neglect the coupling of the microwave beam in vacuum to the plasma. + + +Normalized current drive efficiency $\gamma$ + +$$ +\gamma_{\text{CD}}= \frac{\xi_{CD}}{32.7}T_{\text{e}} +$$ + +Absolute current drive efficiency $\gamma$ +$$ +\gamma_{\text{CD}}= \frac{\gamma_{\text{CD}}}{n_{\text{e,20}}R_0} +$$ + + +The EBWs can only couple to the plasma if the cyclotron harmonic is above the plasma density cut-off. In order to capture this behaviour, we introduce an ad-hoc model which reduces the plasma current up to this condition, by way of a tanh function: + +Where $\mathtt{fp}$ is the plasma frequency, $\mathtt{fc}$ is the cyclotron frequency, $\mathtt{harnum}$ is the harmonic number and $a$ is a free parameter which defines the sharpness of the transition. + +The effect of this factor can be seen below: + +$$ +a = 0.1 +$$ + +$$ +\mathtt{fc} = \frac{\frac{1}{2\pi}\times \mathtt{harnum} \times e B_{\text{T}}}{m_{\text{e}}} +$$ + +$$ +\mathtt{fp} = \frac{1}{2\pi}\times \sqrt{\frac{n_{\text{e}}\text{e}^2}{m_{\text{e}} \epsilon_0}} +$$ + +$$ +\mathtt{density factor} = 0.5 \times \left(1 + \tanh{\left(\frac{2}{\mathtt{a}}\times \frac{\mathtt{fp}-\mathtt{fc}}{\mathtt{fp}- \mathtt{a}}\right)}\right) +$$ + +$$ +\text{Current drive efficiency [A/W]} = \gamma_{\text{CD}} \times \mathtt{densityfactor} +$$ + +
+![EBW Coupling](../images/ebw_coupling.png){ width = "100"} +
Figure 1: EBW coupling function as a function of on-axis toroidal field strength
+
+ + + + +[^1]: Laqua HP, Maassberg H, Marushchenko NB, Volpe F, Weller A, Kasparek W; ECRH-Group. Electron-Bernstein-wave current drive in an overdense plasma at the Wendelstein 7-AS stellarator. Phys Rev Lett. 2003 Feb 21;90(7):075003. doi: 10.1103/PhysRevLett.90.075003. Epub 2003 Feb 21. PMID: 12633236. + +[^2]: G. Taylor, P. C. Efthimion, C. E. Kessel, R. W. Harvey, A. P. Smirnov, N. M. Ershov, M. D. Carter, C. B. Forest; Efficient generation of noninductive, off-axis, Ohkawa current, driven by electron Bernstein waves in high +$\beta$⁠, spherical torus plasmas. Phys. Plasmas 1 October 2004; 11 (10): 4733–4739. https://doi.org/10.1063/1.1792635 + +[^3]: Luce, T., Lin-Liu, Y., Harvey, R., Giruzzi, G., Politzer, P., Rice, B., Lohr, J., Petty, C., & Prater, R. (1999). Generation of Localized Noninductive Current by Electron Cyclotron Waves on the DIII-D Tokamak. Phys. Rev. Lett., 83, 4550–4553. \ No newline at end of file diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ebw_overview.md b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ebw_overview.md new file mode 100644 index 0000000000..f6e9fb57ce --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ebw_overview.md @@ -0,0 +1,9 @@ +# Electron Bernstein Wave Heating + +Crawford et al. initially demonstrated the experimental verification of Electron Bernstein Waves (EBWs) in 1964. Subsequent experiments focused on measuring the transmission and emission of EBWs in linear low-temperature plasma devices, with Crawford providing an overview of these foundational experiments. Among these, two notable examples stand out. + +In Landauer's experiments, the strong interaction between EBWs and Electron Cyclotron (EC) waves was showcased. Landauer observed the excitation of EBWs by an unstable electron distribution function, with emission detected even at the 20th harmonic. Furthermore, Landauer confirmed the theoretically anticipated polarization dependence. + +In the experiments led by Leuterer, EBW excitation and detection were achieved using small antennas placed at the plasma edge. Leuterer successfully demonstrated both forward and backward wave propagation. By measuring the phase difference between emitter and detector, the dispersion relation of EBWs could be estimated. + +These fundamental understandings and experimental validations of EBW physics prompted considerations for utilizing EBWs in fusion plasma heating and temperature diagnostics. The appeal of EBWs lies in their ability to access over-dense plasmas, which is limited in conventional plasma heating methods using electromagnetic EC waves. Additionally, EBWs offer high cyclotron absorption, enabling plasma heating at higher resonant frequencies (>2 harmonics), compared to the limited applicability of electromagnetic EC waves primarily at the first harmonic ordinary and extraordinary mode, and the second harmonic extraordinary mode. Higher harmonic heating with electromagnetic EC waves is only feasible under specific high-temperature conditions. \ No newline at end of file diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ec_overview.md b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ec_overview.md new file mode 100644 index 0000000000..b94ead6a88 --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ec_overview.md @@ -0,0 +1,134 @@ +# Electron Cyclotron Heating + +Electron cyclotron resonance heating is the simplest of the radio frequency heating methods. In contrast to ion cyclotron and lower hybrid heating, there is no evanescent region between the antenna and the plasma,although cut-offs can exist within the plasma. +As a result, the antenna can be retracted into a less hostile environment than for the other schemes. Electron cyclotron heating has been made possible by the invention of the gyrotron millimetre wave source. This and related devices have only emerged since the mid-1970s. Gyrotron tubes capable of 0.5-1 MW out- put and with frequencies in the range 100-200 GHz are now under intense development. The frequency required for a reactor would be in the range 100-200 GHz which corresponds to vacuum wavelengths of 1-2 mm. +Since $\omega_{ce}\le \omega_{pe}$ only the electrons respond to electron cyclotron waves and only the electrons are heated directly. However, under reactor- like conditions the ions will be heated collisionally by the electrons. As the density is increased, a limit is encountered above which electron cyclotron waves cannot penetrate to the central regions of a tokamak. Propagation is again described below + +$$ +n_{\perp}^2 = 1 - \frac{\omega_{pe}^2}{\omega^2} \ \ (\text{O-mode}) +$$ + +$$ +n_{\perp}^2=\frac{(1-\frac{\omega_{pe}^2}{\omega^2}-\frac{\omega_{ce}}{\omega})(1-\frac{\omega_{pe}^2}{\omega^2}+\frac{\omega_{ce}}{\omega})}{(1-\frac{\omega_{pe}^2}{\omega^2}-\frac{\omega_{ce}^2}{\omega^2})} \ \ \text{(X-mode)} +$$ + +## Normalised current drive efficiency | `eccdef()` + +One of the methods for calculating the normalised current drive efficiency is the `eccdef` method found below. + +| Input | Description | +| :---------- | :----------------------------------- | +| $\mathtt{tlocal}$ | Local electron temperature (keV) | +| $\mathtt{epsloc}$ | Local inverse aspect ratio | +| $\mathtt{zlocal}$ | Local plasma effective charge | +| $\mathtt{cosang}$ | Cosine of the poloidal angle at which ECCD takes place (+1 outside, -1 inside) +| $\mathtt{coulog}$ | Local coulomb logarithm for ion-electron collisions | + +This routine calculates the current drive parameters for a electron cyclotron system, based on the AEA FUS 172 model. +It works out the ECCD efficiency using the formula due to Cohen quoted in the ITER Physics Design Guidelines : 1989 (but including division by the Coulomb Logarithm omitted from IPDG89). + +We have assumed $\gamma^2$-1 << 1, where gamma is the relativistic factor. +The notation follows that in IPDG89. +The answer ECGAM is the normalised efficiency $n_{\text{e}}IR/P$ with $n_{\text{e}}$ the local density in [$10^{20} / \text{m}^3$], I the driven current in [$\text{MA}$], $R$ the major radius in [$\text{m}$], and $P$ the absorbed power in [$\text{MW}$]. + + +$$ +\mathtt{mcsq} = 9.1095\times10^{-31} \frac{c^2}{1 \text{keV}} +$$ + +$$ +\mathtt{f} = 16\left(\frac{\mathtt{tlocal}}{\mathtt{mcsq}}\right)^2 +$$ + +$\mathtt{fp}$ is the derivative of $\mathtt{f}$ with respect to gamma, the relativistic factor, taken equal to $1 + \frac{2T_{\text{e}}}{(m_{\text{e}}c^2)}$ + +$$ +\mathtt{fp} = 16 \left(\frac{\mathtt{tlocal}}{\mathtt{mcsq}}\right) +$$ + +$\mathtt{lam}$ is IPDG89's lambda. `legend` calculates the Legendre function of order $\alpha$ and argument `lam`, `palpha`, and its derivative, `palphap`. +Here `alpha` satisfies $\alpha(\alpha+1) = \frac{-8}{(1+z_{\text{local}})}$. $\alpha$ is of the form $(-1/2 + ix)$, with x a real number and $i = \sqrt{-1}$. + +$$ +\mathtt{lam} = 1.0 +$$ + +$$ +\mathtt{palpha, palphap} = \text{legend}(\mathtt{zlocal, lam}) +$$ + +$$ +\mathtt{lams} = \sqrt{\frac{2 \times \mathtt{epsloc}}{1 + \mathtt{epsloc}}} +$$ + +$$ +\mathtt{palphas} = \text{legend}(\mathtt{zlocal, lams}) +$$ + +$\mathtt{hp}$ is the derivative of IPDG89's \mathtt{h} function with respect to $\mathtt{lam}$ + +$$ +\mathtt{h} = \frac{-4 \times \mathtt{lam}}{\mathtt{zlocal + 5}} \times \frac{1- (\mathtt{lams \times \mathtt{palpha}})}{\mathtt{lam \times \mathtt{palphas}}} +$$ + +$$ +\mathtt{hp} = \frac{-4}{\mathtt{zlocal}+5} \times \left(1- \mathtt{lams} \times \frac{\mathtt{palphap}}{\mathtt{palphas}} +\right) +$$ + +$\mathtt{facm}$ is IPDG89's momentum conserving factor + +$$ +\mathtt{facm} = 1.5 +$$ + +$$ +\mathtt{y} = \frac{\mathtt{mcsq}}{2 \times \mathtt{tlocal}} \times (1+ \mathtt{epsloc} \times \mathtt{cosang}) +$$ + +We take the negative of the IPDG89 expression to get a positive number + +$$ +\mathtt{ecgam} = \left(\frac{-7.8 \times \mathtt{facm \times \sqrt{\frac{1 + \mathtt{epsloc}}{1- \mathtt{epsloc}}}}}{\mathtt{coulog}}\right) \times (\mathtt{h} \times \mathtt{fp} -0.5 \times \mathtt{y} \times \mathtt{f} \times \mathtt{hp}) +$$ + +---------------------------------------------------------------------------------- +### Legendre function and its derivative | `legend()` + + +| Input | Description | +| :---------- | :----------------------------------- | +| $\mathtt{zlocal}$ | Local plasma effective charge | +| $\mathtt{arg}$ | Argument of Legendre function | + +The `legend()` function is a routine that calculates the Legendre function and its derivative. It takes two input parameters: `zlocal` (local plasma effective charge) and `arg` (argument of the Legendre function). The function returns two output values: `palpha` (value of the Legendre function) and `palphap` (derivative of the Legendre function). + +Here is the explanation of the `legend()` function: + +1. Check if the absolute value of `arg` is greater than `1.0 + 1.0e-10`. If it is, set `eh.fdiags[0]` to `arg` and report an error (error code 18). + +2. Set `arg2` to the minimum value between `arg` and `1.0 - 1.0e-10`. + +3. Calculate `sinsq` as `0.5 * (1.0 - arg2)`. + +4. Calculate `xisq` as `0.25 * (32.0 * zlocal / (zlocal + 1.0) - 1.0)`. + +5. Initialize `palpha` to `1.0`, `pold` to `1.0`, `pterm` to `1.0`, `palphap` to `0.0`, and `poldp` to `0.0`. + +6. Start a loop that iterates up to 10000 times. + +7. Check for convergence every 20 iterations: + - If `n > 1` and `(n % 20) == 1`, calculate `term1` as `1.0e-10 * max(abs(pold), abs(palpha))` and `term2` as `1.0e-10 * max(abs(poldp), abs(palphap))`. + - If the absolute difference between `pold` and `palpha` is less than `term1` and the absolute difference between `poldp` and `palphap` is less than `term2`, return `palpha` and `palphap`. + +8. Update `pold` to `palpha` and `poldp` to `palphap`. + +9. Calculate `pterm` as `pterm * (4.0 * xisq + (2.0 * n - 1.0) ** 2) / (2.0 * n) ** 2 * sinsq`. + +10. Update `palpha` as `palpha + pterm`. + +11. Update `palphap` as `palphap - n * pterm / (1.0 - arg2)`. + +12. If the loop completes without returning, report an error (error code 19). + +[^1]: Abramowitz, Milton. *"Abramowitz and stegun: Handbook of mathematical functions."* US Department of Commerce 10 (1972). \ No newline at end of file diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ecrh_gamma.md b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ecrh_gamma.md new file mode 100644 index 0000000000..399eb1f932 --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ecrh_gamma.md @@ -0,0 +1,15 @@ +# ECRH User Input Gamma Model + +- `iefrf/iefrffix` = 10 + +This model allows the user to input a scaling factor to the current drive efficiency with the variable `gamma_ecrh`. The value of this variable should follow the value and form of the expression below: + +$$ +\gamma_{CD} = \frac{\langle n_{e,20} \rangle I_{CD}R_0}{P_{CD}} +$$ + +The current drive efficiency is then calculated with $\gamma_{CD}$ in the form below: + +$$ +\text{Current drive efficiency [A/W]} = \frac{\gamma_{CD}}{R_{0} n_{e,20}} +$$ \ No newline at end of file diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ehst_lower_hybrid.md b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ehst_lower_hybrid.md new file mode 100644 index 0000000000..15bff14676 --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ehst_lower_hybrid.md @@ -0,0 +1,6 @@ +# Ehst Lower Hybrid +- `iefrf/iefrffix` = 4 +$$ +\text{Current drive efficiency [A/W]} = \frac{T_{\text{e}}^{0.77} (0.034 + 0.196\beta)}{R_0 n_{\text{e},20}}\frac{\frac{32}{5+Z_{\text{eff}}}+2+\frac{\frac{12\left(6+Z_{\text{eff}}\right)}{5+Z_{\text{eff}}}}{3+Z_{\text{eff}}}+\frac{3.76}{Z_{\text{eff}}}}{12.507} +$$ + diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/RF/fenstermacher_electron_cyclotron_resonance.md b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/fenstermacher_electron_cyclotron_resonance.md new file mode 100644 index 0000000000..47165b65b7 --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/fenstermacher_electron_cyclotron_resonance.md @@ -0,0 +1,6 @@ +# Fenstermacher Electron Cyclotron Resonance + +- `iefrf/iefrffix` = 3 +$$ +\text{Current drive efficiency [A/W]} = \frac{0.21 T_{\langle e,n_{\text{e}} \rangle}}{R_0n_{\text{e,20}}(31.0-(\log{n_{\text{e}}}/2)+(\log{(T_\text{e}}\times1000))} +$$ diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/RF/fenstermacher_lower_hybrid.md b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/fenstermacher_lower_hybrid.md new file mode 100644 index 0000000000..a085f17fb3 --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/fenstermacher_lower_hybrid.md @@ -0,0 +1,7 @@ +# Fenstermacher Lower Hybrid + +- `iefrf/iefrffix` = 1: + +$$ +\text{Current drive efficiency [A/W]} = 0.36 \frac{(1+(T_{\text{e}}/25)^{1.16})}{R_{0} n_{\text{e},20}} +$$ \ No newline at end of file diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ic_model.md b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ic_model.md new file mode 100644 index 0000000000..fa8855f893 --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ic_model.md @@ -0,0 +1,7 @@ +# Ion cyclotron model + +- `iefrf/iefrffix` = 2 + +$$ +\text{Current drive efficiency [A/W]} = \frac{\frac{0.063 T_{\langle \text{e}, n_{\text{e}}\rangle}}{2+Z_{\text{eff}}}}{R_0n_{\text{e,20}}} +$$ diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ic_overview.md b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ic_overview.md new file mode 100644 index 0000000000..2d31283e66 --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/ic_overview.md @@ -0,0 +1,16 @@ +# Ion Cyclotron Heating + + +Ion cyclotron heating (ICH) is a technique used in fusion tokamaks to heat ions in a plasma by applying radio frequency (RF) waves. The RF waves used in ICH are typically in the range of tens to hundreds of megahertz. + +To implement ICH in a fusion tokamak, an antenna system is used to couple the RF waves to the plasma. The waves create an oscillating electric field that interacts with the ions in the plasma. As the ions move in the presence of this electric field, they gain energy through a process called resonant absorption. + +Resonant absorption occurs when the frequency of the RF waves matches the cyclotron frequency of the ions. The cyclotron frequency is the frequency at which ions in a magnetic field undergo circular motion. When the RF waves and the ions' cyclotron frequency are in resonance, the ions absorb energy from the waves more efficiently. + +The resonant absorption process leads to the heating of the ions in the plasma. The absorbed energy increases the kinetic energy of the ions, causing them to move faster and collide with other particles in the plasma. These collisions transfer energy to the surrounding particles, resulting in an overall increase in plasma temperature. + +ICH is an important technique in fusion tokamaks and is used to control and manipulate the properties of plasmas by selectively heating specific ions or regions within the plasma. + +In summary, ion cyclotron heating is a method of heating ions in a fusion tokamak by applying RF waves that are in resonance with the ions' cyclotron frequency. This resonant absorption process efficiently transfers energy from the waves to the ions, leading to their heating and increasing the overall plasma temperature. + + diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/RF/lhcd_overview.md b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/lhcd_overview.md new file mode 100644 index 0000000000..84d8cd1a57 --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/lhcd_overview.md @@ -0,0 +1,4 @@ +# Lower Hybrid Current Drive +Lower hybrid heating is generated through the propagation of a lower hybrid wave within a plasma, characterized by a frequency lying between the ion and electron cyclotron frequencies. This wave possesses an electric field component parallel to the magnetic field, enabling the acceleration of electrons traveling along the field lines. + +Lower hybrid resonance heating employs electromagnetic power within the 1-8 GHz frequency range. At 3 GHz, the vacuum wavelength measures 0.1 meters, facilitating the transmission and launching of energy via waveguides operating in their fundamental mode. This feature stands out as a primary advantage of this method. diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/RF/rf_overview.md b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/rf_overview.md new file mode 100644 index 0000000000..cb43a30780 --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/RF/rf_overview.md @@ -0,0 +1,25 @@ +# Radio Frequency Heating + +In the domain of RF heating, electromagnetic energy is introduced into plasma using radio or microwave frequencies, enabling resonant energy transfer from the waves to the plasma. Various RF heating methods are currently being explored, including Electron Cyclotron Resonance Heating (ECRH) operating in the microwave range, Ion Cyclotron Resonance Heating (ICRH) in the radio wave range, and Lower Hybrid Heating (LHH) in the intermediate range between ICRF and ECRH. As a result, RF heating encompasses a diverse range of independent approaches. + +Ion Cyclotron Resonance Heating has demonstrated effectiveness in ion heating within plasma over an extended period. Both slow and fast wave systems have been developed for ion heating in various plasma confinement devices like stellarators, tokamaks, and mirrors. This method has shown efficiency at the first and second harmonic, as well as higher ion cyclotron frequencies. Achieving power levels of several megawatts has been successful, resulting in ion heating efficiencies comparable to those achieved with a neutral beam. + +The transmitter serves as a wave generator, producing the required RF signal in terms of power and frequency, maintaining this signal from reactor startup to the burn phase. The RF wave is then transmitted through a coaxial cable or other suitable transmission line to the tuning network, which optimizes the system impedance for efficient power delivery to the plasma. The actual delivery is carried out by an antenna located in the plasma chamber. This antenna includes a radiating element, often a loop, covered by a Faraday shield. Designing the Faraday shield presents a significant challenge, requiring consideration of polarization requirements, heat loads from the plasma, erosion caused by plasma particles, and induced RF currents. + + +
+![JET 1984 RF Antenna](../images/JET_rf_panel_1984.jpg){ width = "100"} +
Figure 1: The early interior of the Joint European Torus (JET) in 1984. This image shows the vacuum vessel and the first radio frequency antenna being installed.
+
+ + +Electron Cyclotron Resonance Heating (ECRH) utilizes microwave energy within the gigahertz frequency range, generated by a device known as a "gyrotron." This process involves a high-power electron beam moving through a strong magnetic field produced by superconducting magnets. The electron beam follows a spiral path within the magnetic field, emitting microwave radiation. These microwaves are then directed to the plasma chamber via waveguides, passing through a dielectric window typically made from materials like ceramic or diamond, before being directed into the plasma for absorption. Successful ECRH systems, incorporating multiple 1-gigawatt gyrotrons, waveguides, and diamond windows, have been developed and deployed on DIII-D. + +Lower Hybrid Heating (LHH) systems share similarities with ECRH and ICRH systems, employing similar equipment such as transmitters, transmission lines, and wave launchers or antennas. However, the key difference lies in the positioning and design of the launcher. Due to the higher frequency of waves used in ECRH and LHH, a waveguide can replace a loop antenna and Faraday shield. This waveguide allows for strategic placement of critical components, including ceramic windows, within the reactor's shielding, thereby reducing radiation damage to these essential elements. + +In scenarios where the plasma frequency exceeds the microwave frequency, preventing microwaves from penetrating into the core plasma where electron cyclotron resonance (ECR) occurs, an alternative approach is required. In such cases, the electron Bernstein wave (EBW) is utilized due to its electrostatic nature, enabling it to propagate without encountering density limitations. The EBW originates through mode conversion from the slow X-mode at the upper hybrid resonance (UHR) and exhibits resonances at harmonics of the EC frequency. + +In tokamak plasmas, where the magnetic field varies inversely with the major radius, spatial regions for EBW propagation are arranged in a band structure between two adjacent harmonic ECRs. By injecting microwaves from the outboard side and triggering EBW at the UHR layer through the OXB mode-conversion scheme, the EBW is confined to propagate solely between the UHR layer and the harmonic ECR layer on the inboard side. + + + diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/heating-and-current-drive.md b/documentation/proc-pages/eng-models/heating_and_current_drive/heating-and-current-drive.md new file mode 100644 index 0000000000..f6bf579f8d --- /dev/null +++ b/documentation/proc-pages/eng-models/heating_and_current_drive/heating-and-current-drive.md @@ -0,0 +1,55 @@ +# Auxiliary Heating & Current Drive Systems + +## Current Drive + +The use of inductive current drive leads to pulsed plant operation because of the limited flux swing that can be achieved using the central solenoid. This poses problems due to the fact that fatigue failures may result, and there would be a need for thermal storage to maintain output of electricity between pulses, and supply power for starting a new pulse.However, the plasma current can also be produced and maintained (partially or wholly) using non-inductive means which, in principle, removes this restriction. `PROCESS` contains a number of auxiliary current drive schemes, including various RF methods (Lower Hybrid, Electron Cyclotron, Electron Bernstein Wave, and Ion Cyclotron (Fast Wave) current drives) and also Neutral Beam current drive systems. The code calculates the efficiency and the resulting power requirements of the chosen system. + +The fraction of the required plasma current to be produced by non-inductive means, `fvsbrnni`, should be set, and flag `irfcd` should be set to 0 for purely inductive scenarios, or 1 otherwise. The current drive efficiency model to be used in this latter case is defined by the value of switch `iefrf`: + +- `iefrf` = 1: [Fenstermacher Lower Hybrid model](RF/fenstermacher_lower_hybrid.md) +- `iefrf` = 2: [Ion cyclotron model](RF/ic_model.md)[^1], +- `iefrf` = 3: [Fenstermacher electron cyclotron resonance model](RF/fenstermacher_electron_cyclotron_resonance.md) +- `iefrf` = 4: [Ehst Lower Hybrid model](RF/ehst_lower_hybrid.md) +- `iefrf` = 5: [ITER neutral beam model](NBI/iter_nb.md)[^1],[^2], +- `iefrf` = 6: [Culham Lower Hybrid model](RF/culham_lower_hybrid.md)[^2], +- `iefrf` = 7: [Culham electron cyclotron model](RF/culham_electron_cyclotron.md)[^2], +- `iefrf` = 8: [Culham neutral beam model](NBI/culham_nb.md)[^2], +- `iefrf` = 9: Oscillating Field current drive :warning: (OBSOLETE-REMOVED), +- `iefrf` = 10: [ECRH user input gamma](RF/ecrh_gamma.md), +- `iefrf` = 11: ECRH "HARE" model [^3] :warning: (OBSOLETE-REMOVED), +- `iefrf` = 12: [EBW user scaling input. Scaling](RF/ebw_freethy.md) (S. Freethy) +- `iefrf` = 13: [ECRH O-mode cutoff with $Z_{\text{eff}}$ and $T_{\text{e}}$](RF/cutoff_ecrh.md) (S. Freethy) [^4], + +!!! Warning "Warning" + At present, the neutral beam models do not include the effect of an edge transport barrier (pedestal) in the plasma profile. + +It is sometimes useful to adjust artificially the current drive efficiency values produced by these routines. This can be achieved by setting the scaling coefficients `feffcd`. The wall plug to plasma efficiencies can also be adjusted, by changing the relevant variable (`etaech`, `etalh`, `etanbi` or `etaof`). + +### Power limits +The maximum amount of desired heating and current drive power can be set with `pinjalw`. This limit can be enforced by activating constraint equation 30 (`icc=30`). +Similarly the lower bound on required heating and current drive power can be set with `auxmin`. This limit can be enforced by activating constraint equation 40 (`icc=40`). + +### Secondary current drive + +It is possible to have more than one type of heating and current drive system in `PROCESS`. This can be enabled by setting the `iefrffix` switch to the desired current drive scheme, following the same numbered selection for `iefrf`. +The power injected by the secondary current drive scheme has to be set to a fixed value. This value can be set with the `pinjfixmw` variable. + +## Plasma heating only + +In addition to current drive, some auxiliary power can be used to only heat the plasma. The value of input parameters `pheat` determines the amount of auxiliary heating power (in MW) to be applied to the plasma. This variable may be used as an iteration variable (`ixc = 11`). + +### Secondary heating + +Like for a current drive and heating system a fixed amount of heating power that does not drive current can be set with the `pheatfix` variable. + +## Ignited plasma + +Switch `ignite` can be used to denote whether the plasma is ignited, i.e. fully self-sustaining without the need for any injected auxiliary power during the burn. If `ignite` = 1, the calculated injected power does not contribute to the plasma power balance, although the cost of the auxiliary power system is taken into account (the system is then assumed to be required to provide heating etc during the plasma start-up phase only - use `pheat` to indicate the power requirement). If `ignite` = 0, the plasma is not ignited, and the auxiliary power is taken into account in the plasma power balance during the burn phase. Also, constraint equation 28 (`icc = 28`) can be turned on to enforce the fusion gain *Q* to be at least `bigqmin`. + +[^1]: N. A. Uckan and ITER Physics Group, *"ITER Physics Design Guidelines: 1989"*, ITER Documentation Series, No. 10, IAEA/ITER/DS/10 (1990) + +[^2]: T. C. Hender, M. K. Bevir, M. Cox, R. J. Hastie, P. J. Knight, C. N. Lashmore-Davies, B. Lloyd, G. P. Maddison, A. W. Morris, M. R. O'Brien, M.F. Turner abd H. R. Wilson, *"Physics Assessment for the European Reactor Study"*, AEA Fusion Report AEA FUS 172 (1992) + +[^3]: E. Poli, M. Müller, H. Zohm, M. Kovari, *"Fast evaluation of the current driven by electron cyclotron waves for reactor studies"*, Physics of Plasmas 1 December 2018; 25 (12): 122501 + +[^4]: Laqua, H & Maassberg, H & Marushchenko, Nikolai & Volpe, Francesco & Weller, A & Kasparek, W. (2003). *"Electron-Bernstein-Wave Current Drive in an Overdense Plasma at the Wendelstein 7-AS Stellarator"*, Physical review letters. 90. 075003. 10.1103/PhysRevLett.90.075003. \ No newline at end of file diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/images/JET_rf_panel_1984.jpg b/documentation/proc-pages/eng-models/heating_and_current_drive/images/JET_rf_panel_1984.jpg new file mode 100644 index 0000000000..ce4cb98ded Binary files /dev/null and b/documentation/proc-pages/eng-models/heating_and_current_drive/images/JET_rf_panel_1984.jpg differ diff --git a/documentation/proc-pages/eng-models/heating_and_current_drive/images/ebw_coupling.png b/documentation/proc-pages/eng-models/heating_and_current_drive/images/ebw_coupling.png new file mode 100644 index 0000000000..4ead682b3f Binary files /dev/null and b/documentation/proc-pages/eng-models/heating_and_current_drive/images/ebw_coupling.png differ diff --git a/documentation/proc-pages/images/portsize.png b/documentation/proc-pages/eng-models/heating_and_current_drive/images/portsize.png similarity index 100% rename from documentation/proc-pages/images/portsize.png rename to documentation/proc-pages/eng-models/heating_and_current_drive/images/portsize.png diff --git a/mkdocs.yml b/mkdocs.yml index 5c2544ec22..8bbe10cb48 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -55,7 +55,31 @@ nav: - Shield: eng-models/shield.md - Divertor: eng-models/divertor.md - Heat transport: eng-models/power-conversion-and-heat-dissipation-systems.md - - Auxiliary Power Systems: eng-models/heating-and-current-drive.md + - Auxiliary Heating & Current Drive Systems: + - Overview: eng-models/heating_and_current_drive/heating-and-current-drive.md + - Radio Frequency: + - Overview: eng-models/heating_and_current_drive/RF/rf_overview.md + - Lower Hybrid: + - Overview: eng-models/heating_and_current_drive/RF/lhcd_overview.md + - Fenstermacher Model: eng-models/heating_and_current_drive/RF/fenstermacher_lower_hybrid.md + - Ehst Model: eng-models/heating_and_current_drive/RF/ehst_lower_hybrid.md + - Culham Model: eng-models/heating_and_current_drive/RF/culham_lower_hybrid.md + - Electron Cyclotron: + - Overview: eng-models/heating_and_current_drive/RF/ec_overview.md + - Fenstermacher Resonnance Model: eng-models/heating_and_current_drive/RF/fenstermacher_electron_cyclotron_resonance.md + - Culham Model: eng-models/heating_and_current_drive/RF/culham_electron_cyclotron.md + - User Input Gamma Model: eng-models/heating_and_current_drive/RF/ecrh_gamma.md + - Cutoff mode: eng-models/heating_and_current_drive/RF/cutoff_ecrh.md + - Ion Cyclotron: + - Overview: eng-models/heating_and_current_drive/RF/ic_overview.md + - Ion cyclotron model: eng-models/heating_and_current_drive/RF/ic_model.md + - Electron Bernstein Wave: + - Overview: eng-models/heating_and_current_drive/RF/ebw_overview.md + - EBW Model: eng-models/heating_and_current_drive/RF/ebw_freethy.md + - Neutral Beam Injection: + - Overview: eng-models/heating_and_current_drive/NBI/nbi_overview.md + - ITER Model: eng-models/heating_and_current_drive/NBI/iter_nb.md + - Culham Model: eng-models/heating_and_current_drive/NBI/culham_nb.md - Cryostat and vacuum system: eng-models/cryostat-and-vacuum-system.md - Plant Availability: eng-models/plant-availability.md - Power Requirements: eng-models/power-requirements.md diff --git a/process/current_drive.py b/process/current_drive.py index 55280deb5e..ba36e86c2f 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -1876,7 +1876,7 @@ def etanb(self, abeam, alphan, alphat, aspect, dene, ebeam, rmajor, ten, zeff): alphat : input real : temperature profile factor aspect : input real : aspect ratio dene : input real : volume averaged electron density (m**-3) - enbeam : input real : neutral beam energy (keV) + ebeam : input real : neutral beam energy (keV) rmajor : input real : plasma major radius (m) ten : input real : density weighted average electron temp. (keV) zeff : input real : plasma effective charge