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 index 992f32b917..d9524b9c2d 100644 --- 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 @@ -1,4 +1,4 @@ -# Neutral Beam Injection Heating +# Neutral Beam Injection Heating | `NeutralBeam` !!! Warning "NBI Models" At present, the neutral beam models do not include the effect of an edge transport barrier (pedestal) in the plasma profile. 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 index f6e9fb57ce..491f7cedc2 100644 --- 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 @@ -1,4 +1,4 @@ -# Electron Bernstein Wave Heating +# Electron Bernstein Wave Heating | `ElectronBernstein` 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. 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 index 0dee08f3f7..86a22bf4a2 100644 --- 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 @@ -1,4 +1,4 @@ -# Electron Cyclotron Heating +# Electron Cyclotron Heating | `ElectronCyclotron` 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. 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 index d9d8bbce47..c841f0b07a 100644 --- 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 @@ -1,5 +1,8 @@ -# Ehst Lower Hybrid +# Ehst Lower Hybrid | `lower_hybrid_ehst()` + + - `i_hcd_primary/i_hcd_secondary` = 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 index 0cbdc459a1..2e3c16a7ae 100644 --- 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 @@ -1,6 +1,9 @@ -# Fenstermacher Electron Cyclotron Resonance +# Fenstermacher Electron Cyclotron Resonance | `electron_cyclotron_fenstermacher()` + +- `i_hcd_primary/i_hcd_secondary` = 3 [^1] -- `i_hcd_primary/i_hcd_secondary` = 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))} $$ + +[^1]: T.C. Hender et al., 'Physics Assessment of the European Reactor Study', AEA FUS 172, 1992. 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 index 11c72f2714..31b54e6c36 100644 --- 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 @@ -1,7 +1,13 @@ -# Fenstermacher Lower Hybrid +# Fenstermacher Lower Hybrid | `lower_hybrid_fenstermacher()` - `i_hcd_primary/i_hcd_secondary` = 1: +This forumla was originally in the Oak RidgeSystems Code[^1], attributed to Fenstermacher and is used in the AEA FUS 172 report[^2]. + $$ \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 +$$ + +[^1]: R.L.Reid et al, Oak Ridge Report ORNL/FEDC-87-7, 1988 + +[^2]: T.C. Hender et al., 'Physics Assessment of the European Reactor Study', AEA FUS 172, 1992. \ 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 index c1cff1ae44..880a68e797 100644 --- 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 @@ -1,7 +1,11 @@ -# Ion cyclotron model +# IPDG89 ion cyclotron model | `ion_cyclotron_ipdg89()` -- `i_hcd_primary/i_hcd_secondary` = 2 +- `i_hcd_primary/i_hcd_secondary` = 2 [^1] [^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}}} $$ + +[^1]: N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989', https://inis.iaea.org/collection/NCLCollectionStore/_Public/21/068/21068960.pdf + +[^2] T.C. Hender et al., 'Physics Assessment of the European Reactor Study', AEA FUS 172, 1992. 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 index 2d31283e66..581e8a3805 100644 --- 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 @@ -1,4 +1,4 @@ -# Ion Cyclotron Heating +# Ion Cyclotron Heating | `IonCyclotron` 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. 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 index 84d8cd1a57..098de4ba7b 100644 --- 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 @@ -1,4 +1,5 @@ -# Lower Hybrid Current Drive +# Lower Hybrid Current Drive | `LowerHybrid` + 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/heating-and-current-drive.md b/documentation/proc-pages/eng-models/heating_and_current_drive/heating-and-current-drive.md index a3025ec314..08f732084a 100644 --- 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 @@ -1,10 +1,10 @@ -# Auxiliary Heating & Current Drive Systems +# Auxiliary Heating & Current Drive Systems | `CurrentDrive` ## 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 `i_hcd_calculations` 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 `i_hcd_primary`: +The fraction of the required plasma current to be produced by non-inductive means, `f_c_plasma_non_inductive`, should be set, and flag `i_hcd_calculations` 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 `i_hcd_primary`: - `i_hcd_primary` = 1: [Fenstermacher Lower Hybrid model](RF/fenstermacher_lower_hybrid.md) - `i_hcd_primary` = 2: [Ion cyclotron model](RF/ic_model.md)[^1], diff --git a/documentation/proc-pages/physics-models/plasma_current/inductive_plasma_current.md b/documentation/proc-pages/physics-models/plasma_current/inductive_plasma_current.md index 90f00358ee..415cab6fa0 100644 --- a/documentation/proc-pages/physics-models/plasma_current/inductive_plasma_current.md +++ b/documentation/proc-pages/physics-models/plasma_current/inductive_plasma_current.md @@ -1,6 +1,6 @@ -Currently in `PROCESS` the inductive current fraction from the CS is not calculated directly but is just equal to ($1 - \mathtt{fvsbrnni}$). Where $\mathtt{fvsbrnni}$ is the sum of the fractions of current driven by non inductive means. +Currently in `PROCESS` the inductive current fraction from the CS is not calculated directly but is just equal to ($1 - \mathtt{f_c_plasma_non_inductive}$). Where $\mathtt{f_c_plasma_non_inductive}$ is the sum of the fractions of current driven by non inductive means. -This calculated fraction (`inductive_current_fraction`) is then used in the `calculate_volt_second_requirements()` and `burn()` functions to calculate the volt-second requirements and the burn time for a pulsed machine. +This calculated fraction (`f_c_plasma_inductive`) is then used in the `calculate_volt_second_requirements()` and `burn()` functions to calculate the volt-second requirements and the burn time for a pulsed machine. !!! info "Inductive plasma current fraction refactor" diff --git a/examples/data/csv_output_large_tokamak_MFILE.DAT b/examples/data/csv_output_large_tokamak_MFILE.DAT index b9cbd5faaf..82ed4e7d65 100644 --- a/examples/data/csv_output_large_tokamak_MFILE.DAT +++ b/examples/data/csv_output_large_tokamak_MFILE.DAT @@ -127,9 +127,9 @@ bore____________________________________________________________________ (itvar035)____________________ 2.0588E+00 bore_(final_value/initial_value)________________________________________ (xcm035)______________________ 1.0294E+00 bore_(range_normalised)_________________________________________________ (nitvar035)___________________ 1.9786E-01 - fvsbrnni________________________________________________________________ (itvar036)____________________ 4.1495E-01 - fvsbrnni_(final_value/initial_value)____________________________________ (xcm036)______________________ 1.0374E+00 - fvsbrnni_(range_normalised)_____________________________________________ (nitvar036)___________________ 4.1437E-01 + f_c_plasma_non_inductive________________________________________________________________ (itvar036)____________________ 4.1495E-01 + f_c_plasma_non_inductive_(final_value/initial_value)____________________________________ (xcm036)______________________ 1.0374E+00 + f_c_plasma_non_inductive_(range_normalised)_____________________________________________ (nitvar036)___________________ 4.1437E-01 tdmptf__________________________________________________________________ (itvar037)____________________ 1.9607E+01 tdmptf_(final_value/initial_value)______________________________________ (xcm037)______________________ 7.8427E-01 tdmptf_(range_normalised)_______________________________________________ (nitvar037)___________________ 1.9526E-01 @@ -545,14 +545,14 @@ Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 4.1436E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 Pfirsch-Schlueter_fraction______________________________________________ (f_c_plasma_pfirsch_schluter)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 5.9038E-04 - Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 5.8505E-01 - Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 - Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 4.1495E-01 ITV - Electron_cyclotron_injected_power_(MW)__________________________________ (p_ecrh_injected_mw)______________________ 7.5213E+01 OP + Auxiliary_current_drive_fraction________________________________________ (f_c_plasma_auxiliary)_______________________ 5.9038E-04 + Inductive_fraction______________________________________________________ (f_c_plasma_inductive)_______________________ 5.8505E-01 + Total___________________________________________________________________ (plasipf+f_c_plasma_auxiliary+f_c_plasma_inductive)_________ 1.0000E+00 + Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (f_c_plasma_non_inductive)____________________ 4.1495E-01 ITV + Electron_cyclotron_injected_power_(MW)__________________________________ (p_hcd_ecrh_injected_total_mw)______________________ 7.5213E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (p_hcd_injected_max)_____________________ 2.0000E+02 ECH_wall_plug_efficiency________________________________________________ (eta_ecrh_injector_wall_plug)______________________ 5.0000E-01 - ECH_wall_plug_power_(MW)________________________________________________ (echwpow)_____________________ 1.5043E+02 OP + ECH_wall_plug_power_(MW)________________________________________________ (p_hcd_ecrh_electric_mw)_____________________ 1.5043E+02 OP Total_V-s_capability_of_Central_Solenoid/PF_coils_(Wb)__________________ (abs(vs_cs_pf_total_pulse))__________________ 5.8939E+02 Required_volt-seconds_during_start-up_(Wb)______________________________ (vssoft)______________________ 2.8401E+02 Available_volt-seconds_during_burn_(Wb)_________________________________ (vsmax)_______________________ 2.8571E+02 @@ -1158,7 +1158,7 @@ Total_(MW)______________________________________________________________ ______________________________ 2.1348E+03 Net_electric_power_output(MW)___________________________________________ (pnetelmw.)___________________ 4.0000E+02 Required_Net_electric_power_output(MW)__________________________________ (pnetelin)____________________ 4.0000E+02 - Electric_power_for_heating_and_current_drive_(MW)_______________________ (pinjwp)______________________ 1.5043E+02 + Electric_power_for_heating_and_current_drive_(MW)_______________________ (p_hcd_electric_total_mw)______________________ 1.5043E+02 Electric_power_for_primary_coolant_pumps_(MW)___________________________ (htpmw)_______________________ 1.8742E+02 Electric_power_for_vacuum_pumps_(MW)____________________________________ (vachtmw)_____________________ 5.0000E-01 Electric_power_for_tritium_plant_(MW)___________________________________ (trithtmw)____________________ 1.5000E+01 @@ -1432,9 +1432,9 @@ ixc = 29 boundl(29) = 0.1 dr_bore = 2.0 -* fvsbrnni +* f_c_plasma_non_inductive ixc = 44 -fvsbrnni = 0.4 +f_c_plasma_non_inductive = 0.4 * tdmptf [s] ixc = 56 diff --git a/examples/data/large_tokamak_1_MFILE.DAT b/examples/data/large_tokamak_1_MFILE.DAT index e69ca39509..820cd132c3 100644 --- a/examples/data/large_tokamak_1_MFILE.DAT +++ b/examples/data/large_tokamak_1_MFILE.DAT @@ -128,9 +128,9 @@ bore____________________________________________________________________ (itvar035)____________________ 1.9854E+00 bore_(final_value/initial_value)________________________________________ (xcm035)______________________ 9.9271E-01 bore_(range_normalised)_________________________________________________ (nitvar035)___________________ 1.9045E-01 - fvsbrnni________________________________________________________________ (itvar036)____________________ 4.3358E-01 - fvsbrnni_(final_value/initial_value)____________________________________ (xcm036)______________________ 1.0839E+00 - fvsbrnni_(range_normalised)_____________________________________________ (nitvar036)___________________ 4.3301E-01 + f_c_plasma_non_inductive________________________________________________________________ (itvar036)____________________ 4.3358E-01 + f_c_plasma_non_inductive_(final_value/initial_value)____________________________________ (xcm036)______________________ 1.0839E+00 + f_c_plasma_non_inductive_(range_normalised)_____________________________________________ (nitvar036)___________________ 4.3301E-01 tdmptf__________________________________________________________________ (itvar037)____________________ 1.8429E+01 tdmptf_(final_value/initial_value)______________________________________ (xcm037)______________________ 7.3716E-01 tdmptf_(range_normalised)_______________________________________________ (nitvar037)___________________ 1.8347E-01 @@ -543,14 +543,14 @@ Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 4.1920E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 Pfirsch-Schlueter_fraction______________________________________________ (f_c_plasma_pfirsch_schluter)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 1.4376E-02 - Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 5.6642E-01 - Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 - Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 4.3358E-01 ITV - Electron_cyclotron_injected_power_(MW)__________________________________ (p_ecrh_injected_mw)______________________ 8.0143E+01 OP + Auxiliary_current_drive_fraction________________________________________ (f_c_plasma_auxiliary)_______________________ 1.4376E-02 + Inductive_fraction______________________________________________________ (f_c_plasma_inductive)_______________________ 5.6642E-01 + Total___________________________________________________________________ (plasipf+f_c_plasma_auxiliary+f_c_plasma_inductive)_________ 1.0000E+00 + Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (f_c_plasma_non_inductive)____________________ 4.3358E-01 ITV + Electron_cyclotron_injected_power_(MW)__________________________________ (p_hcd_ecrh_injected_total_mw)______________________ 8.0143E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (p_hcd_injected_max)_____________________ 2.0000E+02 ECH_wall_plug_efficiency________________________________________________ (eta_ecrh_injector_wall_plug)______________________ 5.0000E-01 - ECH_wall_plug_power_(MW)________________________________________________ (echwpow)_____________________ 1.6029E+02 OP + ECH_wall_plug_power_(MW)________________________________________________ (p_hcd_ecrh_electric_mw)_____________________ 1.6029E+02 OP Total_V-s_capability_of_Central_Solenoid/PF_coils_(Wb)__________________ (abs(vs_cs_pf_total_pulse))__________________ 5.7728E+02 Required_volt-seconds_during_start-up_(Wb)______________________________ (vssoft)______________________ 2.8601E+02 Available_volt-seconds_during_burn_(Wb)_________________________________ (vsmax)_______________________ 2.7014E+02 @@ -1153,7 +1153,7 @@ Total_(MW)______________________________________________________________ ______________________________ 2.1752E+03 Net_electric_power_output(MW)___________________________________________ (pnetelmw.)___________________ 4.0001E+02 Required_Net_electric_power_output(MW)__________________________________ (pnetelin)____________________ 4.0000E+02 - Electric_power_for_heating_and_current_drive_(MW)_______________________ (pinjwp)______________________ 1.6029E+02 + Electric_power_for_heating_and_current_drive_(MW)_______________________ (p_hcd_electric_total_mw)______________________ 1.6029E+02 Electric_power_for_primary_coolant_pumps_(MW)___________________________ (htpmw)_______________________ 1.9053E+02 Electric_power_for_vacuum_pumps_(MW)____________________________________ (vachtmw)_____________________ 5.0000E-01 Electric_power_for_tritium_plant_(MW)___________________________________ (trithtmw)____________________ 1.5000E+01 @@ -1426,9 +1426,9 @@ ixc = 29 boundl(29) = 0.1 dr_bore = 2.0 -* fvsbrnni +* f_c_plasma_non_inductive ixc = 44 -fvsbrnni = 0.4 +f_c_plasma_non_inductive = 0.4 * tdmptf [s] ixc = 56 diff --git a/examples/data/large_tokamak_2_MFILE.DAT b/examples/data/large_tokamak_2_MFILE.DAT index 76c83a6731..17c6b98a33 100644 --- a/examples/data/large_tokamak_2_MFILE.DAT +++ b/examples/data/large_tokamak_2_MFILE.DAT @@ -128,9 +128,9 @@ bore____________________________________________________________________ (itvar035)____________________ 1.9854E+00 bore_(final_value/initial_value)________________________________________ (xcm035)______________________ 9.9271E-01 bore_(range_normalised)_________________________________________________ (nitvar035)___________________ 1.9045E-01 - fvsbrnni________________________________________________________________ (itvar036)____________________ 4.3358E-01 - fvsbrnni_(final_value/initial_value)____________________________________ (xcm036)______________________ 1.0839E+00 - fvsbrnni_(range_normalised)_____________________________________________ (nitvar036)___________________ 4.3301E-01 + f_c_plasma_non_inductive________________________________________________________________ (itvar036)____________________ 4.3358E-01 + f_c_plasma_non_inductive_(final_value/initial_value)____________________________________ (xcm036)______________________ 1.0839E+00 + f_c_plasma_non_inductive_(range_normalised)_____________________________________________ (nitvar036)___________________ 4.3301E-01 tdmptf__________________________________________________________________ (itvar037)____________________ 1.8429E+01 tdmptf_(final_value/initial_value)______________________________________ (xcm037)______________________ 7.3716E-01 tdmptf_(range_normalised)_______________________________________________ (nitvar037)___________________ 1.8347E-01 @@ -543,14 +543,14 @@ Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 4.1920E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 Pfirsch-Schlueter_fraction______________________________________________ (f_c_plasma_pfirsch_schluter)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 1.4376E-02 - Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 5.6642E-01 - Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 - Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 4.3358E-01 ITV - Electron_cyclotron_injected_power_(MW)__________________________________ (p_ecrh_injected_mw)______________________ 8.0143E+01 OP + Auxiliary_current_drive_fraction________________________________________ (f_c_plasma_auxiliary)_______________________ 1.4376E-02 + Inductive_fraction______________________________________________________ (f_c_plasma_inductive)_______________________ 5.6642E-01 + Total___________________________________________________________________ (plasipf+f_c_plasma_auxiliary+f_c_plasma_inductive)_________ 1.0000E+00 + Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (f_c_plasma_non_inductive)____________________ 4.3358E-01 ITV + Electron_cyclotron_injected_power_(MW)__________________________________ (p_hcd_ecrh_injected_total_mw)______________________ 8.0143E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (p_hcd_injected_max)_____________________ 2.0000E+02 ECH_wall_plug_efficiency________________________________________________ (eta_ecrh_injector_wall_plug)______________________ 5.0000E-01 - ECH_wall_plug_power_(MW)________________________________________________ (echwpow)_____________________ 1.6029E+02 OP + ECH_wall_plug_power_(MW)________________________________________________ (p_hcd_ecrh_electric_mw)_____________________ 1.6029E+02 OP Total_V-s_capability_of_Central_Solenoid/PF_coils_(Wb)__________________ (abs(vs_cs_pf_total_pulse))__________________ 5.7728E+02 Required_volt-seconds_during_start-up_(Wb)______________________________ (vssoft)______________________ 2.8601E+02 Available_volt-seconds_during_burn_(Wb)_________________________________ (vsmax)_______________________ 2.7014E+02 @@ -1153,7 +1153,7 @@ Total_(MW)______________________________________________________________ ______________________________ 2.1752E+03 Net_electric_power_output(MW)___________________________________________ (pnetelmw.)___________________ 4.0001E+02 Required_Net_electric_power_output(MW)__________________________________ (pnetelin)____________________ 4.0000E+02 - Electric_power_for_heating_and_current_drive_(MW)_______________________ (pinjwp)______________________ 1.6029E+02 + Electric_power_for_heating_and_current_drive_(MW)_______________________ (p_hcd_electric_total_mw)______________________ 1.6029E+02 Electric_power_for_primary_coolant_pumps_(MW)___________________________ (htpmw)_______________________ 1.9053E+02 Electric_power_for_vacuum_pumps_(MW)____________________________________ (vachtmw)_____________________ 5.0000E-01 Electric_power_for_tritium_plant_(MW)___________________________________ (trithtmw)____________________ 1.5000E+01 @@ -1426,9 +1426,9 @@ ixc = 29 boundl(29) = 0.1 dr_bore = 2.0 -* fvsbrnni +* f_c_plasma_non_inductive ixc = 44 -fvsbrnni = 0.4 +f_c_plasma_non_inductive = 0.4 * tdmptf [s] ixc = 56 diff --git a/examples/data/large_tokamak_3_MFILE.DAT b/examples/data/large_tokamak_3_MFILE.DAT index 0a09b73653..7f9443783c 100644 --- a/examples/data/large_tokamak_3_MFILE.DAT +++ b/examples/data/large_tokamak_3_MFILE.DAT @@ -128,9 +128,9 @@ bore____________________________________________________________________ (itvar035)____________________ 1.1854E+00 bore_(final_value/initial_value)________________________________________ (xcm035)______________________ 9.9271E-01 bore_(range_normalised)_________________________________________________ (nitvar035)___________________ 1.9045E-01 - fvsbrnni________________________________________________________________ (itvar036)____________________ 4.5358E-01 - fvsbrnni_(final_value/initial_value)____________________________________ (xcm036)______________________ 1.0839E+00 - fvsbrnni_(range_normalised)_____________________________________________ (nitvar036)___________________ 4.3301E-01 + f_c_plasma_non_inductive________________________________________________________________ (itvar036)____________________ 4.5358E-01 + f_c_plasma_non_inductive_(final_value/initial_value)____________________________________ (xcm036)______________________ 1.0839E+00 + f_c_plasma_non_inductive_(range_normalised)_____________________________________________ (nitvar036)___________________ 4.3301E-01 tdmptf__________________________________________________________________ (itvar037)____________________ 1.2429E+01 tdmptf_(final_value/initial_value)______________________________________ (xcm037)______________________ 7.3716E-01 tdmptf_(range_normalised)_______________________________________________ (nitvar037)___________________ 1.8347E-01 @@ -543,14 +543,14 @@ Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 4.1920E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 Pfirsch-Schlueter_fraction______________________________________________ (f_c_plasma_pfirsch_schluter)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 1.4376E-02 - Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 5.6642E-01 - Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 - Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 4.3358E-01 ITV - Electron_cyclotron_injected_power_(MW)__________________________________ (p_ecrh_injected_mw)______________________ 8.0143E+01 OP + Auxiliary_current_drive_fraction________________________________________ (f_c_plasma_auxiliary)_______________________ 1.4376E-02 + Inductive_fraction______________________________________________________ (f_c_plasma_inductive)_______________________ 5.6642E-01 + Total___________________________________________________________________ (plasipf+f_c_plasma_auxiliary+f_c_plasma_inductive)_________ 1.0000E+00 + Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (f_c_plasma_non_inductive)____________________ 4.3358E-01 ITV + Electron_cyclotron_injected_power_(MW)__________________________________ (p_hcd_ecrh_injected_total_mw)______________________ 8.0143E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (p_hcd_injected_max)_____________________ 2.0000E+02 ECH_wall_plug_efficiency________________________________________________ (eta_ecrh_injector_wall_plug)______________________ 5.0000E-01 - ECH_wall_plug_power_(MW)________________________________________________ (echwpow)_____________________ 1.6029E+02 OP + ECH_wall_plug_power_(MW)________________________________________________ (p_hcd_ecrh_electric_mw)_____________________ 1.6029E+02 OP Total_V-s_capability_of_Central_Solenoid/PF_coils_(Wb)__________________ (abs(vs_cs_pf_total_pulse))__________________ 5.7728E+02 Required_volt-seconds_during_start-up_(Wb)______________________________ (vssoft)______________________ 2.8601E+02 Available_volt-seconds_during_burn_(Wb)_________________________________ (vsmax)_______________________ 2.7014E+02 @@ -1153,7 +1153,7 @@ Total_(MW)______________________________________________________________ ______________________________ 2.1752E+03 Net_electric_power_output(MW)___________________________________________ (pnetelmw.)___________________ 4.0001E+02 Required_Net_electric_power_output(MW)__________________________________ (pnetelin)____________________ 4.0000E+02 - Electric_power_for_heating_and_current_drive_(MW)_______________________ (pinjwp)______________________ 1.6029E+02 + Electric_power_for_heating_and_current_drive_(MW)_______________________ (p_hcd_electric_total_mw)______________________ 1.6029E+02 Electric_power_for_primary_coolant_pumps_(MW)___________________________ (htpmw)_______________________ 1.9053E+02 Electric_power_for_vacuum_pumps_(MW)____________________________________ (vachtmw)_____________________ 5.0000E-01 Electric_power_for_tritium_plant_(MW)___________________________________ (trithtmw)____________________ 1.5000E+01 @@ -1426,9 +1426,9 @@ ixc = 29 boundl(29) = 0.1 dr_bore = 2.0 -* fvsbrnni +* f_c_plasma_non_inductive ixc = 44 -fvsbrnni = 0.4 +f_c_plasma_non_inductive = 0.4 * tdmptf [s] ixc = 56 diff --git a/examples/data/large_tokamak_4_MFILE.DAT b/examples/data/large_tokamak_4_MFILE.DAT index 71d77629a8..5f3c9f43df 100644 --- a/examples/data/large_tokamak_4_MFILE.DAT +++ b/examples/data/large_tokamak_4_MFILE.DAT @@ -128,9 +128,9 @@ bore____________________________________________________________________ (itvar035)____________________ 1.9854E+00 bore_(final_value/initial_value)________________________________________ (xcm035)______________________ 9.9271E-01 bore_(range_normalised)_________________________________________________ (nitvar035)___________________ 1.9045E-01 - fvsbrnni________________________________________________________________ (itvar036)____________________ 4.3358E-01 - fvsbrnni_(final_value/initial_value)____________________________________ (xcm036)______________________ 1.0839E+00 - fvsbrnni_(range_normalised)_____________________________________________ (nitvar036)___________________ 4.3301E-01 + f_c_plasma_non_inductive________________________________________________________________ (itvar036)____________________ 4.3358E-01 + f_c_plasma_non_inductive_(final_value/initial_value)____________________________________ (xcm036)______________________ 1.0839E+00 + f_c_plasma_non_inductive_(range_normalised)_____________________________________________ (nitvar036)___________________ 4.3301E-01 tdmptf__________________________________________________________________ (itvar037)____________________ 1.8429E+01 tdmptf_(final_value/initial_value)______________________________________ (xcm037)______________________ 7.3716E-01 tdmptf_(range_normalised)_______________________________________________ (nitvar037)___________________ 1.8347E-01 @@ -543,14 +543,14 @@ Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 4.1920E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 Pfirsch-Schlueter_fraction______________________________________________ (f_c_plasma_pfirsch_schluter)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 1.4376E-02 - Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 5.6642E-01 - Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 - Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 4.3358E-01 ITV - Electron_cyclotron_injected_power_(MW)__________________________________ (p_ecrh_injected_mw)______________________ 8.0143E+01 OP + Auxiliary_current_drive_fraction________________________________________ (f_c_plasma_auxiliary)_______________________ 1.4376E-02 + Inductive_fraction______________________________________________________ (f_c_plasma_inductive)_______________________ 5.6642E-01 + Total___________________________________________________________________ (plasipf+f_c_plasma_auxiliary+f_c_plasma_inductive)_________ 1.0000E+00 + Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (f_c_plasma_non_inductive)____________________ 4.3358E-01 ITV + Electron_cyclotron_injected_power_(MW)__________________________________ (p_hcd_ecrh_injected_total_mw)______________________ 8.0143E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (p_hcd_injected_max)_____________________ 2.0000E+02 ECH_wall_plug_efficiency________________________________________________ (eta_ecrh_injector_wall_plug)______________________ 5.0000E-01 - ECH_wall_plug_power_(MW)________________________________________________ (echwpow)_____________________ 1.6029E+02 OP + ECH_wall_plug_power_(MW)________________________________________________ (p_hcd_ecrh_electric_mw)_____________________ 1.6029E+02 OP Total_V-s_capability_of_Central_Solenoid/PF_coils_(Wb)__________________ (abs(vs_cs_pf_total_pulse))__________________ 5.7728E+02 Required_volt-seconds_during_start-up_(Wb)______________________________ (vssoft)______________________ 2.8601E+02 Available_volt-seconds_during_burn_(Wb)_________________________________ (vsmax)_______________________ 2.7014E+02 @@ -1153,7 +1153,7 @@ Total_(MW)______________________________________________________________ ______________________________ 2.1752E+03 Net_electric_power_output(MW)___________________________________________ (pnetelmw.)___________________ 4.0001E+02 Required_Net_electric_power_output(MW)__________________________________ (pnetelin)____________________ 4.0000E+02 - Electric_power_for_heating_and_current_drive_(MW)_______________________ (pinjwp)______________________ 1.6029E+02 + Electric_power_for_heating_and_current_drive_(MW)_______________________ (p_hcd_electric_total_mw)______________________ 1.6029E+02 Electric_power_for_primary_coolant_pumps_(MW)___________________________ (htpmw)_______________________ 1.9053E+02 Electric_power_for_vacuum_pumps_(MW)____________________________________ (vachtmw)_____________________ 5.0000E-01 Electric_power_for_tritium_plant_(MW)___________________________________ (trithtmw)____________________ 1.5000E+01 @@ -1426,9 +1426,9 @@ ixc = 29 boundl(29) = 0.1 dr_bore = 2.0 -* fvsbrnni +* f_c_plasma_non_inductive ixc = 44 -fvsbrnni = 0.4 +f_c_plasma_non_inductive = 0.4 * tdmptf [s] ixc = 56 diff --git a/examples/data/large_tokamak_IN.DAT b/examples/data/large_tokamak_IN.DAT index b374171e95..dbe7a0b1b1 100644 --- a/examples/data/large_tokamak_IN.DAT +++ b/examples/data/large_tokamak_IN.DAT @@ -236,9 +236,9 @@ ixc = 29 boundl(29) = 0.1 dr_bore = 2.0 -* fvsbrnni +* f_c_plasma_non_inductive ixc = 44 -fvsbrnni = 0.4 +f_c_plasma_non_inductive = 0.4 * tdmptf [s] ixc = 56 diff --git a/examples/data/large_tokamak_once_through_IN.DAT b/examples/data/large_tokamak_once_through_IN.DAT index 48235fad95..2e5b0f5782 100644 --- a/examples/data/large_tokamak_once_through_IN.DAT +++ b/examples/data/large_tokamak_once_through_IN.DAT @@ -58,7 +58,7 @@ ixc = 29 * bore boundl(29) = 0.1 ixc = 37 * coheof ixc = 41 * fcohbop -ixc = 44 * fvsbrnni +ixc = 44 * f_c_plasma_non_inductive ixc = 56 * tdmptf ixc = 57 * dr_tf_nose_case ixc = 58 * thwcndut @@ -321,7 +321,7 @@ dene = 7.796223900029837e+19 * electron density (/m3) (`iteration variable 6 beta_norm_max = 3.0 * Troyon-like coefficient for beta scaling calculated fgwsep = 0.5 * fraction of Greenwald density to set as separatrix density; If `<0`; separatrix fkzohm = 1.02 * Zohm elongation scaling adjustment factor (`ishape=2; 3`) -fvsbrnni = 0.4242184436680697 * fraction of the plasma current produced by non-inductive means (`iteration variable 44`) +f_c_plasma_non_inductive = 0.4242184436680697 * fraction of the plasma current produced by non-inductive means (`iteration variable 44`) ejima_coeff = 0.3 * Ejima coefficient for resistive startup V-s formula hfact = 1.185971818905028 * H factor on energy confinement times; radiation corrected (`iteration variable 10`); i_bootstrap_current = 4 * switch for bootstrap current scaling diff --git a/examples/data/scan_MFILE.DAT b/examples/data/scan_MFILE.DAT index 8c7ce00364..6e648b3317 100644 --- a/examples/data/scan_MFILE.DAT +++ b/examples/data/scan_MFILE.DAT @@ -398,14 +398,14 @@ Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 Pfirsch-Schlueter_fraction______________________________________________ (f_c_plasma_pfirsch_schluter)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 - Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV - Electron_cyclotron_injected_power_(MW)__________________________________ (p_ecrh_injected_mw)______________________ 5.1000E+01 OP + Auxiliary_current_drive_fraction________________________________________ (f_c_plasma_auxiliary)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (f_c_plasma_inductive)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+f_c_plasma_auxiliary+f_c_plasma_inductive)_________ 1.0000E+00 + Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (f_c_plasma_non_inductive)____________________ 3.8724E-01 ITV + Electron_cyclotron_injected_power_(MW)__________________________________ (p_hcd_ecrh_injected_total_mw)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (p_hcd_injected_max)_____________________ 5.1000E+01 ECH_wall_plug_efficiency________________________________________________ (eta_ecrh_injector_wall_plug)______________________ 4.0000E-01 - ECH_wall_plug_power_(MW)________________________________________________ (echwpow)_____________________ 1.2750E+02 OP + ECH_wall_plug_power_(MW)________________________________________________ (p_hcd_ecrh_electric_mw)_____________________ 1.2750E+02 OP Total_V-s_capability_of_Central_Solenoid/PF_coils_(Wb)__________________ (abs(vs_cs_pf_total_pulse))__________________ 6.2115E+02 Required_volt-seconds_during_start-up_(Wb)______________________________ (vssoft)______________________ 3.6189E+02 Available_volt-seconds_during_burn_(Wb)_________________________________ (vsmax)_______________________ 2.3054E+02 @@ -980,7 +980,7 @@ Total_(MW)______________________________________________________________ ______________________________ 2.6502E+03 Net_electric_power_output(MW)___________________________________________ (pnetelmw.)___________________ 5.0000E+02 Required_Net_electric_power_output(MW)__________________________________ (pnetelin)____________________ 5.0000E+02 - Electric_power_for_heating_and_current_drive_(MW)_______________________ (pinjwp)______________________ 1.2750E+02 + Electric_power_for_heating_and_current_drive_(MW)_______________________ (p_hcd_electric_total_mw)______________________ 1.2750E+02 Electric_power_for_primary_coolant_pumps_(MW)___________________________ (htpmw)_______________________ 2.3440E+02 Electric_power_for_vacuum_pumps_(MW)____________________________________ (vachtmw)_____________________ 5.0000E-01 Electric_power_for_tritium_plant_(MW)___________________________________ (trithtmw)____________________ 1.5000E+01 @@ -1393,14 +1393,14 @@ Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 Pfirsch-Schlueter_fraction______________________________________________ (f_c_plasma_pfirsch_schluter)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 - Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV - Electron_cyclotron_injected_power_(MW)__________________________________ (p_ecrh_injected_mw)______________________ 5.1000E+01 OP + Auxiliary_current_drive_fraction________________________________________ (f_c_plasma_auxiliary)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (f_c_plasma_inductive)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+f_c_plasma_auxiliary+f_c_plasma_inductive)_________ 1.0000E+00 + Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (f_c_plasma_non_inductive)____________________ 3.8724E-01 ITV + Electron_cyclotron_injected_power_(MW)__________________________________ (p_hcd_ecrh_injected_total_mw)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (p_hcd_injected_max)_____________________ 5.1000E+01 ECH_wall_plug_efficiency________________________________________________ (eta_ecrh_injector_wall_plug)______________________ 4.0000E-01 - ECH_wall_plug_power_(MW)________________________________________________ (echwpow)_____________________ 1.2750E+02 OP + ECH_wall_plug_power_(MW)________________________________________________ (p_hcd_ecrh_electric_mw)_____________________ 1.2750E+02 OP Total_V-s_capability_of_Central_Solenoid/PF_coils_(Wb)__________________ (abs(vs_cs_pf_total_pulse))__________________ 6.2115E+02 Required_volt-seconds_during_start-up_(Wb)______________________________ (vssoft)______________________ 3.6189E+02 Available_volt-seconds_during_burn_(Wb)_________________________________ (vsmax)_______________________ 2.3054E+02 @@ -1975,7 +1975,7 @@ Total_(MW)______________________________________________________________ ______________________________ 2.6502E+03 Net_electric_power_output(MW)___________________________________________ (pnetelmw.)___________________ 5.0000E+02 Required_Net_electric_power_output(MW)__________________________________ (pnetelin)____________________ 5.0000E+02 - Electric_power_for_heating_and_current_drive_(MW)_______________________ (pinjwp)______________________ 1.2750E+02 + Electric_power_for_heating_and_current_drive_(MW)_______________________ (p_hcd_electric_total_mw)______________________ 1.2750E+02 Electric_power_for_primary_coolant_pumps_(MW)___________________________ (htpmw)_______________________ 2.3440E+02 Electric_power_for_vacuum_pumps_(MW)____________________________________ (vachtmw)_____________________ 5.0000E-01 Electric_power_for_tritium_plant_(MW)___________________________________ (trithtmw)____________________ 1.5000E+01 @@ -2388,14 +2388,14 @@ Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 Pfirsch-Schlueter_fraction______________________________________________ (f_c_plasma_pfirsch_schluter)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 - Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV - Electron_cyclotron_injected_power_(MW)__________________________________ (p_ecrh_injected_mw)______________________ 5.1000E+01 OP + Auxiliary_current_drive_fraction________________________________________ (f_c_plasma_auxiliary)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (f_c_plasma_inductive)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+f_c_plasma_auxiliary+f_c_plasma_inductive)_________ 1.0000E+00 + Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (f_c_plasma_non_inductive)____________________ 3.8724E-01 ITV + Electron_cyclotron_injected_power_(MW)__________________________________ (p_hcd_ecrh_injected_total_mw)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (p_hcd_injected_max)_____________________ 5.1000E+01 ECH_wall_plug_efficiency________________________________________________ (eta_ecrh_injector_wall_plug)______________________ 4.0000E-01 - ECH_wall_plug_power_(MW)________________________________________________ (echwpow)_____________________ 1.2750E+02 OP + ECH_wall_plug_power_(MW)________________________________________________ (p_hcd_ecrh_electric_mw)_____________________ 1.2750E+02 OP Total_V-s_capability_of_Central_Solenoid/PF_coils_(Wb)__________________ (abs(vs_cs_pf_total_pulse))__________________ 6.2115E+02 Required_volt-seconds_during_start-up_(Wb)______________________________ (vssoft)______________________ 3.6189E+02 Available_volt-seconds_during_burn_(Wb)_________________________________ (vsmax)_______________________ 2.3054E+02 @@ -2970,7 +2970,7 @@ Total_(MW)______________________________________________________________ ______________________________ 2.6502E+03 Net_electric_power_output(MW)___________________________________________ (pnetelmw.)___________________ 5.0000E+02 Required_Net_electric_power_output(MW)__________________________________ (pnetelin)____________________ 5.0000E+02 - Electric_power_for_heating_and_current_drive_(MW)_______________________ (pinjwp)______________________ 1.2750E+02 + Electric_power_for_heating_and_current_drive_(MW)_______________________ (p_hcd_electric_total_mw)______________________ 1.2750E+02 Electric_power_for_primary_coolant_pumps_(MW)___________________________ (htpmw)_______________________ 2.3440E+02 Electric_power_for_vacuum_pumps_(MW)____________________________________ (vachtmw)_____________________ 5.0000E-01 Electric_power_for_tritium_plant_(MW)___________________________________ (trithtmw)____________________ 1.5000E+01 @@ -3383,14 +3383,14 @@ Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 Pfirsch-Schlueter_fraction______________________________________________ (f_c_plasma_pfirsch_schluter)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 - Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV - Electron_cyclotron_injected_power_(MW)__________________________________ (p_ecrh_injected_mw)______________________ 5.1000E+01 OP + Auxiliary_current_drive_fraction________________________________________ (f_c_plasma_auxiliary)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (f_c_plasma_inductive)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+f_c_plasma_auxiliary+f_c_plasma_inductive)_________ 1.0000E+00 + Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (f_c_plasma_non_inductive)____________________ 3.8724E-01 ITV + Electron_cyclotron_injected_power_(MW)__________________________________ (p_hcd_ecrh_injected_total_mw)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (p_hcd_injected_max)_____________________ 5.1000E+01 ECH_wall_plug_efficiency________________________________________________ (eta_ecrh_injector_wall_plug)______________________ 4.0000E-01 - ECH_wall_plug_power_(MW)________________________________________________ (echwpow)_____________________ 1.2750E+02 OP + ECH_wall_plug_power_(MW)________________________________________________ (p_hcd_ecrh_electric_mw)_____________________ 1.2750E+02 OP Total_V-s_capability_of_Central_Solenoid/PF_coils_(Wb)__________________ (abs(vs_cs_pf_total_pulse))__________________ 6.2115E+02 Required_volt-seconds_during_start-up_(Wb)______________________________ (vssoft)______________________ 3.6189E+02 Available_volt-seconds_during_burn_(Wb)_________________________________ (vsmax)_______________________ 2.3054E+02 @@ -3965,7 +3965,7 @@ Total_(MW)______________________________________________________________ ______________________________ 2.6502E+03 Net_electric_power_output(MW)___________________________________________ (pnetelmw.)___________________ 5.0000E+02 Required_Net_electric_power_output(MW)__________________________________ (pnetelin)____________________ 5.0000E+02 - Electric_power_for_heating_and_current_drive_(MW)_______________________ (pinjwp)______________________ 1.2750E+02 + Electric_power_for_heating_and_current_drive_(MW)_______________________ (p_hcd_electric_total_mw)______________________ 1.2750E+02 Electric_power_for_primary_coolant_pumps_(MW)___________________________ (htpmw)_______________________ 2.3440E+02 Electric_power_for_vacuum_pumps_(MW)____________________________________ (vachtmw)_____________________ 5.0000E-01 Electric_power_for_tritium_plant_(MW)___________________________________ (trithtmw)____________________ 1.5000E+01 @@ -4378,14 +4378,14 @@ Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 Pfirsch-Schlueter_fraction______________________________________________ (f_c_plasma_pfirsch_schluter)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 - Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV - Electron_cyclotron_injected_power_(MW)__________________________________ (p_ecrh_injected_mw)______________________ 5.1000E+01 OP + Auxiliary_current_drive_fraction________________________________________ (f_c_plasma_auxiliary)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (f_c_plasma_inductive)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+f_c_plasma_auxiliary+f_c_plasma_inductive)_________ 1.0000E+00 + Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (f_c_plasma_non_inductive)____________________ 3.8724E-01 ITV + Electron_cyclotron_injected_power_(MW)__________________________________ (p_hcd_ecrh_injected_total_mw)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (p_hcd_injected_max)_____________________ 5.1000E+01 ECH_wall_plug_efficiency________________________________________________ (eta_ecrh_injector_wall_plug)______________________ 4.0000E-01 - ECH_wall_plug_power_(MW)________________________________________________ (echwpow)_____________________ 1.2750E+02 OP + ECH_wall_plug_power_(MW)________________________________________________ (p_hcd_ecrh_electric_mw)_____________________ 1.2750E+02 OP Total_V-s_capability_of_Central_Solenoid/PF_coils_(Wb)__________________ (abs(vs_cs_pf_total_pulse))__________________ 6.2115E+02 Required_volt-seconds_during_start-up_(Wb)______________________________ (vssoft)______________________ 3.6189E+02 Available_volt-seconds_during_burn_(Wb)_________________________________ (vsmax)_______________________ 2.3054E+02 @@ -4960,7 +4960,7 @@ Total_(MW)______________________________________________________________ ______________________________ 2.6502E+03 Net_electric_power_output(MW)___________________________________________ (pnetelmw.)___________________ 5.0000E+02 Required_Net_electric_power_output(MW)__________________________________ (pnetelin)____________________ 5.0000E+02 - Electric_power_for_heating_and_current_drive_(MW)_______________________ (pinjwp)______________________ 1.2750E+02 + Electric_power_for_heating_and_current_drive_(MW)_______________________ (p_hcd_electric_total_mw)______________________ 1.2750E+02 Electric_power_for_primary_coolant_pumps_(MW)___________________________ (htpmw)_______________________ 2.3440E+02 Electric_power_for_vacuum_pumps_(MW)____________________________________ (vachtmw)_____________________ 5.0000E-01 Electric_power_for_tritium_plant_(MW)___________________________________ (trithtmw)____________________ 1.5000E+01 @@ -5373,14 +5373,14 @@ Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 Pfirsch-Schlueter_fraction______________________________________________ (f_c_plasma_pfirsch_schluter)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 - Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV - Electron_cyclotron_injected_power_(MW)__________________________________ (p_ecrh_injected_mw)______________________ 5.1000E+01 OP + Auxiliary_current_drive_fraction________________________________________ (f_c_plasma_auxiliary)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (f_c_plasma_inductive)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+f_c_plasma_auxiliary+f_c_plasma_inductive)_________ 1.0000E+00 + Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (f_c_plasma_non_inductive)____________________ 3.8724E-01 ITV + Electron_cyclotron_injected_power_(MW)__________________________________ (p_hcd_ecrh_injected_total_mw)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (p_hcd_injected_max)_____________________ 5.1000E+01 ECH_wall_plug_efficiency________________________________________________ (eta_ecrh_injector_wall_plug)______________________ 4.0000E-01 - ECH_wall_plug_power_(MW)________________________________________________ (echwpow)_____________________ 1.2750E+02 OP + ECH_wall_plug_power_(MW)________________________________________________ (p_hcd_ecrh_electric_mw)_____________________ 1.2750E+02 OP Total_V-s_capability_of_Central_Solenoid/PF_coils_(Wb)__________________ (abs(vs_cs_pf_total_pulse))__________________ 6.2115E+02 Required_volt-seconds_during_start-up_(Wb)______________________________ (vssoft)______________________ 3.6189E+02 Available_volt-seconds_during_burn_(Wb)_________________________________ (vsmax)_______________________ 2.3054E+02 @@ -5955,7 +5955,7 @@ Total_(MW)______________________________________________________________ ______________________________ 2.6502E+03 Net_electric_power_output(MW)___________________________________________ (pnetelmw.)___________________ 5.0000E+02 Required_Net_electric_power_output(MW)__________________________________ (pnetelin)____________________ 5.0000E+02 - Electric_power_for_heating_and_current_drive_(MW)_______________________ (pinjwp)______________________ 1.2750E+02 + Electric_power_for_heating_and_current_drive_(MW)_______________________ (p_hcd_electric_total_mw)______________________ 1.2750E+02 Electric_power_for_primary_coolant_pumps_(MW)___________________________ (htpmw)_______________________ 2.3440E+02 Electric_power_for_vacuum_pumps_(MW)____________________________________ (vachtmw)_____________________ 5.0000E-01 Electric_power_for_tritium_plant_(MW)___________________________________ (trithtmw)____________________ 1.5000E+01 @@ -6368,14 +6368,14 @@ Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 Pfirsch-Schlueter_fraction______________________________________________ (f_c_plasma_pfirsch_schluter)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 - Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV - Electron_cyclotron_injected_power_(MW)__________________________________ (p_ecrh_injected_mw)______________________ 5.1000E+01 OP + Auxiliary_current_drive_fraction________________________________________ (f_c_plasma_auxiliary)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (f_c_plasma_inductive)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+f_c_plasma_auxiliary+f_c_plasma_inductive)_________ 1.0000E+00 + Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (f_c_plasma_non_inductive)____________________ 3.8724E-01 ITV + Electron_cyclotron_injected_power_(MW)__________________________________ (p_hcd_ecrh_injected_total_mw)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (p_hcd_injected_max)_____________________ 5.1000E+01 ECH_wall_plug_efficiency________________________________________________ (eta_ecrh_injector_wall_plug)______________________ 4.0000E-01 - ECH_wall_plug_power_(MW)________________________________________________ (echwpow)_____________________ 1.2750E+02 OP + ECH_wall_plug_power_(MW)________________________________________________ (p_hcd_ecrh_electric_mw)_____________________ 1.2750E+02 OP Total_V-s_capability_of_Central_Solenoid/PF_coils_(Wb)__________________ (abs(vs_cs_pf_total_pulse))__________________ 6.2115E+02 Required_volt-seconds_during_start-up_(Wb)______________________________ (vssoft)______________________ 3.6189E+02 Available_volt-seconds_during_burn_(Wb)_________________________________ (vsmax)_______________________ 2.3054E+02 @@ -6950,7 +6950,7 @@ Total_(MW)______________________________________________________________ ______________________________ 2.6502E+03 Net_electric_power_output(MW)___________________________________________ (pnetelmw.)___________________ 5.0000E+02 Required_Net_electric_power_output(MW)__________________________________ (pnetelin)____________________ 5.0000E+02 - Electric_power_for_heating_and_current_drive_(MW)_______________________ (pinjwp)______________________ 1.2750E+02 + Electric_power_for_heating_and_current_drive_(MW)_______________________ (p_hcd_electric_total_mw)______________________ 1.2750E+02 Electric_power_for_primary_coolant_pumps_(MW)___________________________ (htpmw)_______________________ 2.3440E+02 Electric_power_for_vacuum_pumps_(MW)____________________________________ (vachtmw)_____________________ 5.0000E-01 Electric_power_for_tritium_plant_(MW)___________________________________ (trithtmw)____________________ 1.5000E+01 @@ -7363,14 +7363,14 @@ Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 Pfirsch-Schlueter_fraction______________________________________________ (f_c_plasma_pfirsch_schluter)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 - Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV - Electron_cyclotron_injected_power_(MW)__________________________________ (p_ecrh_injected_mw)______________________ 5.1000E+01 OP + Auxiliary_current_drive_fraction________________________________________ (f_c_plasma_auxiliary)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (f_c_plasma_inductive)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+f_c_plasma_auxiliary+f_c_plasma_inductive)_________ 1.0000E+00 + Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (f_c_plasma_non_inductive)____________________ 3.8724E-01 ITV + Electron_cyclotron_injected_power_(MW)__________________________________ (p_hcd_ecrh_injected_total_mw)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (p_hcd_injected_max)_____________________ 5.1000E+01 ECH_wall_plug_efficiency________________________________________________ (eta_ecrh_injector_wall_plug)______________________ 4.0000E-01 - ECH_wall_plug_power_(MW)________________________________________________ (echwpow)_____________________ 1.2750E+02 OP + ECH_wall_plug_power_(MW)________________________________________________ (p_hcd_ecrh_electric_mw)_____________________ 1.2750E+02 OP Total_V-s_capability_of_Central_Solenoid/PF_coils_(Wb)__________________ (abs(vs_cs_pf_total_pulse))__________________ 6.2115E+02 Required_volt-seconds_during_start-up_(Wb)______________________________ (vssoft)______________________ 3.6189E+02 Available_volt-seconds_during_burn_(Wb)_________________________________ (vsmax)_______________________ 2.3054E+02 @@ -7945,7 +7945,7 @@ Total_(MW)______________________________________________________________ ______________________________ 2.6502E+03 Net_electric_power_output(MW)___________________________________________ (pnetelmw.)___________________ 5.0000E+02 Required_Net_electric_power_output(MW)__________________________________ (pnetelin)____________________ 5.0000E+02 - Electric_power_for_heating_and_current_drive_(MW)_______________________ (pinjwp)______________________ 1.2750E+02 + Electric_power_for_heating_and_current_drive_(MW)_______________________ (p_hcd_electric_total_mw)______________________ 1.2750E+02 Electric_power_for_primary_coolant_pumps_(MW)___________________________ (htpmw)_______________________ 2.3440E+02 Electric_power_for_vacuum_pumps_(MW)____________________________________ (vachtmw)_____________________ 5.0000E-01 Electric_power_for_tritium_plant_(MW)___________________________________ (trithtmw)____________________ 1.5000E+01 @@ -8358,14 +8358,14 @@ Bootstrap_fraction______________________________________________________ (bootipf)_____________________ 3.8470E-01 Diamagnetic_fraction____________________________________________________ (diaipf)______________________ 0.0000E+00 Pfirsch-Schlueter_fraction______________________________________________ (f_c_plasma_pfirsch_schluter)_______________________ 0.0000E+00 - Auxiliary_current_drive_fraction________________________________________ (aux_current_fraction)_______________________ 2.5418E-03 - Inductive_fraction______________________________________________________ (inductive_current_fraction)_______________________ 6.1276E-01 - Total___________________________________________________________________ (plasipf+aux_current_fraction+inductive_current_fraction)_________ 1.0000E+00 - Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (fvsbrnni)____________________ 3.8724E-01 ITV - Electron_cyclotron_injected_power_(MW)__________________________________ (p_ecrh_injected_mw)______________________ 5.1000E+01 OP + Auxiliary_current_drive_fraction________________________________________ (f_c_plasma_auxiliary)_______________________ 2.5418E-03 + Inductive_fraction______________________________________________________ (f_c_plasma_inductive)_______________________ 6.1276E-01 + Total___________________________________________________________________ (plasipf+f_c_plasma_auxiliary+f_c_plasma_inductive)_________ 1.0000E+00 + Fraction_of_the_plasma_current_produced_by_non-inductive_means__________ (f_c_plasma_non_inductive)____________________ 3.8724E-01 ITV + Electron_cyclotron_injected_power_(MW)__________________________________ (p_hcd_ecrh_injected_total_mw)______________________ 5.1000E+01 OP Maximum_allowable_ECRH_power_(MW)_______________________________________ (p_hcd_injected_max)_____________________ 5.1000E+01 ECH_wall_plug_efficiency________________________________________________ (eta_ecrh_injector_wall_plug)______________________ 4.0000E-01 - ECH_wall_plug_power_(MW)________________________________________________ (echwpow)_____________________ 1.2750E+02 OP + ECH_wall_plug_power_(MW)________________________________________________ (p_hcd_ecrh_electric_mw)_____________________ 1.2750E+02 OP Total_V-s_capability_of_Central_Solenoid/PF_coils_(Wb)__________________ (abs(vs_cs_pf_total_pulse))__________________ 6.2115E+02 Required_volt-seconds_during_start-up_(Wb)______________________________ (vssoft)______________________ 3.6189E+02 Available_volt-seconds_during_burn_(Wb)_________________________________ (vsmax)_______________________ 2.3054E+02 @@ -8940,7 +8940,7 @@ Total_(MW)______________________________________________________________ ______________________________ 2.6502E+03 Net_electric_power_output(MW)___________________________________________ (pnetelmw.)___________________ 5.0000E+02 Required_Net_electric_power_output(MW)__________________________________ (pnetelin)____________________ 5.0000E+02 - Electric_power_for_heating_and_current_drive_(MW)_______________________ (pinjwp)______________________ 1.2750E+02 + Electric_power_for_heating_and_current_drive_(MW)_______________________ (p_hcd_electric_total_mw)______________________ 1.2750E+02 Electric_power_for_primary_coolant_pumps_(MW)___________________________ (htpmw)_______________________ 2.3440E+02 Electric_power_for_vacuum_pumps_(MW)____________________________________ (vachtmw)_____________________ 5.0000E-01 Electric_power_for_tritium_plant_(MW)___________________________________ (trithtmw)____________________ 1.5000E+01 @@ -9033,7 +9033,7 @@ ixc = 41 * f_j_cs_start_pulse_end_flat_top ixc = 42 * dr_cs_tf_gap boundl(42) = 0.05 boundu(42) = 0.1 -ixc = 44 * fvsbrnni +ixc = 44 * f_c_plasma_non_inductive ixc = 48 * fstrcase ixc = 49 * fstrcond ixc = 50 * fiooic @@ -9235,7 +9235,7 @@ aspect = 3.1 * Aspect ratio (iteration variable 1) dene = 7.983e+19 * Electron density (/m3) (iteration variable 6) beta_norm_max = 3.0 * (troyon-like) coefficient for beta scaling; fkzohm = 1.0245 * Zohm elongation scaling adjustment factor (i_plasma_geometry=2; 3) -fvsbrnni = 0.4434 * Fraction of the plasma current produced by +f_c_plasma_non_inductive = 0.4434 * Fraction of the plasma current produced by ejima_coeff = 0.3 * Ejima coefficient for resistive startup v-s formula hfact = 1.1 * H factor on energy confinement times (iteration variable 10) i_bootstrap_current = 4 * Switch for bootstrap current scaling; @@ -9331,7 +9331,7 @@ t_burn = 1.0d4 * Burn time (s) (calculated if i_pulsed_plant=1) fjohc0 = 5.3923E-01 f_j_cs_start_pulse_end_flat_top = 9.3176E-01 dr_cs_tf_gap = 5.0000E-02 - fvsbrnni = 3.9566E-01 + f_c_plasma_non_inductive = 3.9566E-01 fstrcase = 1.0000E+00 fstrcond = 9.2007E-01 fiooic = 6.3437E-01 diff --git a/examples/data/scan_example_file_IN.DAT b/examples/data/scan_example_file_IN.DAT index 051a50be64..0154ccc43b 100644 --- a/examples/data/scan_example_file_IN.DAT +++ b/examples/data/scan_example_file_IN.DAT @@ -236,9 +236,9 @@ ixc = 29 boundl(29) = 0.1 dr_bore = 2.0 -* fvsbrnni +* f_c_plasma_non_inductive ixc = 44 -fvsbrnni = 0.4 +f_c_plasma_non_inductive = 0.4 * tdmptf [s] ixc = 56 diff --git a/process/costs.py b/process/costs.py index 6ef2b1ffb0..a8dbc66222 100644 --- a/process/costs.py +++ b/process/costs.py @@ -1807,7 +1807,8 @@ def acc223(self): self.c2231 = ( 1.0e-6 * cost_variables.ucech - * (1.0e6 * current_drive_variables.p_ecrh_injected_mw) ** exprf + * (1.0e6 * current_drive_variables.p_hcd_ecrh_injected_total_mw) + ** exprf ) if cost_variables.ifueltyp == 1: @@ -1820,13 +1821,15 @@ def acc223(self): self.c2232 = ( 1.0e-6 * cost_variables.uclh - * (1.0e6 * current_drive_variables.plhybd) ** exprf + * (1.0e6 * current_drive_variables.p_hcd_lowhyb_injected_total_mw) + ** exprf ) else: self.c2232 = ( 1.0e-6 * cost_variables.ucich - * (1.0e6 * current_drive_variables.plhybd) ** exprf + * (1.0e6 * current_drive_variables.p_hcd_lowhyb_injected_total_mw) + ** exprf ) if cost_variables.ifueltyp == 1: @@ -1835,7 +1838,7 @@ def acc223(self): # Account 223.3 : Neutral Beam - # self.c2233 = 1.0e-6 * cost_variables.ucnbi * (1.0e6*pnbeam)**exprf + # self.c2233 = 1.0e-6 * cost_variables.ucnbi * (1.0e6*p_hcd_beam_injected_total_mw)**exprf # #327 self.c2233 = ( @@ -2471,7 +2474,7 @@ def acc26(self): if cost_variables.ireactor == 0: pwrrej = ( physics_variables.fusion_power - + heat_transport_variables.pinjwp + + heat_transport_variables.p_hcd_electric_total_mw + tfcoil_variables.tfcmw ) else: diff --git a/process/current_drive.py b/process/current_drive.py index 9d066259eb..028e9aafd6 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -6,7 +6,6 @@ from process.exceptions import ProcessError, ProcessValueError from process.fortran import ( constants, - cost_variables, current_drive_variables, heat_transport_variables, physics_variables, @@ -17,1400 +16,944 @@ from process.plasma_profiles import PlasmaProfile -class CurrentDrive: +class NeutralBeam: def __init__(self, plasma_profile: PlasmaProfile): self.outfile = constants.nout self.plasma_profile = plasma_profile - def cudriv(self, output: bool): - """Routine to calculate the current drive power requirements + def iternb(self): + """Routine to calculate ITER Neutral Beam current drive parameters author: P J Knight, CCFE, Culham Science Centre - outfile : input integer : output file unit - iprint : input integer : switch for writing to output file (1=yes) - This routine calculates the power requirements of the current - drive system, using a choice of models for the current drive - efficiency. + effnbss : output real : neutral beam current drive efficiency (A/W) + f_p_beam_injected_ions : output real : fraction of NB power given to ions + fshine : output real : shine-through fraction of beam + This routine calculates the current drive parameters for a + neutral beam system, based on the 1990 ITER model. + ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, + ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 """ - - current_drive_variables.p_ecrh_injected_mw = 0.0e0 - current_drive_variables.pnbeam = 0.0e0 - current_drive_variables.plhybd = 0.0e0 - current_drive_variables.c_beam_total = 0.0e0 - beam_current_fix = 0.0e0 - current_drive_variables.p_beam_orbit_loss_mw = 0.0e0 - porbitlossmwfix = 0.0e0 - - pinjmw1 = 0.0 - pinjmwfix = 0.0 - pinjimw1 = 0.0 - pinjemw1 = 0.0 - pinjemwfix = 0.0 - pinjimwfix = 0.0 - auxiliary_cdfix = 0.0 - aux_current_fraction_fix = 0.0 - gamcdfix = 0.0e0 - - # To stop issues with input file we force - # zero secondary heating if no injection method - if current_drive_variables.i_hcd_secondary == 0: - current_drive_variables.p_hcd_secondary_extra_heat_mw = 0.0 - - # check for unphysically large heating in - # secondary injected power source - if ( - current_drive_variables.p_hcd_secondary_extra_heat_mw - > current_drive_variables.p_hcd_secondary_injected_mw - ): - current_drive_variables.p_hcd_secondary_extra_heat_mw = ( - current_drive_variables.p_hcd_secondary_injected_mw + # Check argument sanity + if (1 + physics_variables.eps) < current_drive_variables.frbeam: + raise ProcessValueError( + "Imminent negative square root argument; NBI will miss plasma completely", + eps=physics_variables.eps, + frbeam=current_drive_variables.frbeam, ) - # current_drive_variables.i_hcd_calculations | switch for current drive calculation - # = 0 | turned off - # = 1 | turned on - if current_drive_variables.i_hcd_calculations != 0: - # put electron density in desired units (10^-20 m-3) - dene20 = physics_variables.dene * 1.0e-20 - - # If present we must calculate second current drive - # efficiencies in units of Amps/Watt using the fixed - # values from user input - # current_drive_variables.i_hcd_secondary | switch for fixed current drive efficiency model - - # Fenstermacher Lower Hybrid model - if current_drive_variables.i_hcd_secondary == 1: - effrfssfix = ( - (0.36e0 * (1.0e0 + (physics_variables.te / 25.0e0) ** 1.16e0)) - / (physics_variables.rmajor * dene20) - * current_drive_variables.feffcd - ) - eta_cd_hcd_secondary = effrfssfix - # Ion-Cyclotron current drive - elif current_drive_variables.i_hcd_secondary == 2: - effrfssfix = ( - 0.63e0 - * 0.1e0 - * physics_variables.ten - / (2.0e0 + physics_variables.zeff) - / (physics_variables.rmajor * dene20) - * current_drive_variables.feffcd - ) - eta_cd_hcd_secondary = effrfssfix - # Fenstermacher Electron Cyclotron Resonance model - elif current_drive_variables.i_hcd_secondary == 3: - effrfssfix = ( - 0.21e0 - * physics_variables.ten - / (physics_variables.rmajor * dene20 * physics_variables.dlamee) - * current_drive_variables.feffcd - ) - eta_cd_hcd_secondary = effrfssfix - # Ehst Lower Hybrid / Fast Wave current drive - elif current_drive_variables.i_hcd_secondary == 4: - effrfssfix = ( - physics_variables.te**0.77e0 - * (0.034e0 + 0.196e0 * physics_variables.beta) - / (physics_variables.rmajor * dene20) - * ( - 32.0e0 / (5.0e0 + physics_variables.zeff) - + 2.0e0 - + (12.0e0 * (6.0e0 + physics_variables.zeff)) - / (5.0e0 + physics_variables.zeff) - / (3.0e0 + physics_variables.zeff) - + 3.76e0 / physics_variables.zeff - ) - / 12.507e0 - * current_drive_variables.feffcd - ) - eta_cd_hcd_secondary = effrfssfix - elif current_drive_variables.i_hcd_secondary == 5: - ( - effnbss, - current_drive_variables.f_p_beam_injected_ions, - current_drive_variables.f_p_beam_shine_through, - ) = self.iternb() - effnbssfix = effnbss * current_drive_variables.feffcd - eta_cd_hcd_secondary = effnbssfix - # Culham Lower Hybrid current drive model - elif current_drive_variables.i_hcd_secondary == 6: - effrfss = self.cullhy() - effrfssfix = effrfss * current_drive_variables.feffcd - eta_cd_hcd_secondary = effrfssfix - # Culham ECCD model - elif current_drive_variables.i_hcd_secondary == 7: - effrfss = self.culecd() - effrfssfix = effrfss * current_drive_variables.feffcd - eta_cd_hcd_secondary = effrfssfix - # Culham Neutral Beam model - elif current_drive_variables.i_hcd_secondary == 8: - ( - effnbss, - current_drive_variables.f_p_beam_injected_ions, - current_drive_variables.f_p_beam_shine_through, - ) = self.culnbi() - effnbssfix = effnbss * current_drive_variables.feffcd - eta_cd_hcd_secondary = effnbssfix - # ECRH user input gamma - elif current_drive_variables.i_hcd_secondary == 10: - # Normalised current drive efficiency gamma - current_drive_variables.eta_cd_norm_hcd_secondary = ( - current_drive_variables.eta_cd_norm_ecrh - ) - - # Absolute current drive efficiency - effrfssfix = current_drive_variables.eta_cd_norm_hcd_secondary / ( - dene20 * physics_variables.rmajor - ) - eta_cd_hcd_secondary = effrfssfix - # EBW scaling - elif current_drive_variables.i_hcd_secondary == 12: - # Scaling author Simon Freethy - # Ref : PROCESS issue 1262 - # Normalised current drive efficiency gamma - current_drive_variables.eta_cd_norm_hcd_secondary = ( - current_drive_variables.xi_ebw / 32.7e0 - ) * physics_variables.te - - # Absolute current drive efficiency - effrfssfix = current_drive_variables.eta_cd_norm_hcd_secondary / ( - dene20 * physics_variables.rmajor - ) - eta_cd_hcd_secondary = effrfssfix - - # EBWs can only couple to plasma if cyclotron harmonic is above plasma density cut-off; - # this behaviour is captured in the following function (ref issue #1262): - # current_drive_variables.n_ecrh_harmonic = cyclotron harmonic number (fundamental used as default) - # constant 'a' controls sharpness of transition - a = 0.1e0 - - fc = ( - 1.0e0 - / (2.0e0 * np.pi) - * current_drive_variables.n_ecrh_harmonic - * constants.electron_charge - * physics_variables.bt - / constants.electron_mass - ) - fp = ( - 1.0e0 - / (2.0e0 * np.pi) - * np.sqrt( - physics_variables.dene - * constants.electron_charge**2 - / (constants.electron_mass * constants.epsilon0) - ) - ) - - density_factor = 0.5e0 * ( - 1.0e0 + np.tanh((2.0e0 / a) * ((fp - fc) / fp - a)) - ) - - eta_cd_hcd_secondary = eta_cd_hcd_secondary * density_factor - effrfssfix = effrfssfix * density_factor - elif current_drive_variables.i_hcd_secondary == 13: - # ECCD model for O-mode cut-off with added Te and Zeff dependance - # Scaling author: Simon Freethy - # Ref : PROCESS issue #2994 - - fc = ( - 1 - / (2 * np.pi) - * constants.electron_charge - * physics_variables.bt - / constants.electron_mass - ) - fp = ( - 1 - / (2 * np.pi) - * np.sqrt( - (physics_variables.dene / 1.0e19) - * constants.electron_charge**2 - / (constants.electron_mass * constants.epsilon0) - ) - ) - - xi_CD = 0.18e0 # This is tuned to the results of a GRAY study - xi_CD = xi_CD * ( - 4.8e0 / (2 + physics_variables.zeff) - ) # Zeff correction - effrfssfix = ( - xi_CD - * physics_variables.te - / ( - 3.27e0 - * physics_variables.rmajor - * (physics_variables.dene / 1.0e19) - ) - ) + # Calculate beam path length to centre + dpath = physics_variables.rmajor * np.sqrt( + (1.0 + physics_variables.eps) ** 2 - current_drive_variables.frbeam**2 + ) - # O-mode case - if current_drive_variables.i_ecrh_wave_mode == 0: - f_cutoff = fp + # Calculate beam stopping cross-section + sigstop = self.sigbeam( + current_drive_variables.e_beam_kev / physics_variables.m_beam_amu, + physics_variables.te, + physics_variables.dene, + physics_variables.f_nd_alpha_electron, + physics_variables.rncne, + physics_variables.rnone, + physics_variables.rnfene, + ) - # X-mode case - elif current_drive_variables.i_ecrh_wave_mode == 1: - f_cutoff = 0.5 * ( - fc - + np.sqrt( - current_drive_variables.n_ecrh_harmonic * fc**2 + 4 * fp**2 - ) - ) + # Calculate number of decay lengths to centre + current_drive_variables.n_beam_decay_lengths_core = ( + dpath * physics_variables.dene * sigstop + ) - # 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 - cutoff_factor = 0.5 * ( - 1 - + np.tanh( - (2 / (a)) - * ( - (current_drive_variables.n_ecrh_harmonic * fc - f_cutoff) - / fp - - a - ) - ) - ) - eta_cd_hcd_secondary = effrfssfix * cutoff_factor - elif current_drive_variables.i_hcd_secondary != 0: - raise ProcessValueError( - f"Current drive switch is invalid: {current_drive_variables.i_hcd_secondary = }" - ) + # Shine-through fraction of beam + fshine = np.exp(-2.0 * dpath * physics_variables.dene * sigstop) + fshine = max(fshine, 1e-20) - if current_drive_variables.i_hcd_secondary in [1, 2, 4, 6]: - # Injected power - pinjemwfix = current_drive_variables.p_hcd_secondary_injected_mw + # Deuterium and tritium beam densities + dend = physics_variables.nd_fuel_ions * ( + 1.0 - current_drive_variables.f_beam_tritium + ) + dent = physics_variables.nd_fuel_ions * current_drive_variables.f_beam_tritium - # Wall plug power - heat_transport_variables.pinjwpfix = ( - current_drive_variables.p_hcd_secondary_injected_mw - / current_drive_variables.eta_lowhyb_injector_wall_plug - ) + # Power split to ions / electrons + f_p_beam_injected_ions = self.cfnbi( + physics_variables.m_beam_amu, + current_drive_variables.e_beam_kev, + physics_variables.ten, + physics_variables.dene, + dend, + dent, + physics_variables.zeffai, + physics_variables.dlamie, + ) - # Wall plug to injector efficiency - current_drive_variables.eta_hcd_secondary_injector_wall_plug = ( - current_drive_variables.eta_lowhyb_injector_wall_plug - ) + # Current drive efficiency + effnbss = current_drive_variables.frbeam * self.etanb( + physics_variables.m_beam_amu, + physics_variables.alphan, + physics_variables.alphat, + physics_variables.aspect, + physics_variables.dene, + current_drive_variables.e_beam_kev, + physics_variables.rmajor, + physics_variables.ten, + physics_variables.zeff, + ) - # Normalised current drive efficiency gamma - gamcdfix = effrfssfix * (dene20 * physics_variables.rmajor) + return effnbss, f_p_beam_injected_ions, fshine - # the fixed auxiliary current - auxiliary_cdfix = ( - effrfssfix - * ( - current_drive_variables.p_hcd_secondary_injected_mw - - current_drive_variables.p_hcd_secondary_extra_heat_mw - ) - * 1.0e6 - ) - aux_current_fraction_fix = ( - auxiliary_cdfix / physics_variables.plasma_current - ) - elif current_drive_variables.i_hcd_secondary in [3, 7, 10, 12, 13]: - # Injected power - pinjemwfix = current_drive_variables.p_hcd_secondary_injected_mw + def culnbi(self): + """Routine to calculate Neutral Beam current drive parameters + author: P J Knight, CCFE, Culham Science Centre + effnbss : output real : neutral beam current drive efficiency (A/W) + f_p_beam_injected_ions : output real : fraction of NB power given to ions + fshine : output real : shine-through fraction of beam + 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). + AEA FUS 172: Physics Assessment for the European Reactor Study + """ + if (1.0e0 + physics_variables.eps) < current_drive_variables.frbeam: + raise ProcessValueError( + "Imminent negative square root argument; NBI will miss plasma completely", + eps=physics_variables.eps, + frbeam=current_drive_variables.frbeam, + ) - # Wall plug power - heat_transport_variables.pinjwpfix = ( - current_drive_variables.p_hcd_secondary_injected_mw - / current_drive_variables.eta_ecrh_injector_wall_plug - ) + # Calculate beam path length to centre - # Wall plug to injector efficiency - current_drive_variables.eta_hcd_secondary_injector_wall_plug = ( - current_drive_variables.eta_ecrh_injector_wall_plug - ) + dpath = physics_variables.rmajor * np.sqrt( + (1.0e0 + physics_variables.eps) ** 2 - current_drive_variables.frbeam**2 + ) - # the fixed auxiliary current - auxiliary_cdfix = ( - effrfssfix - * ( - current_drive_variables.p_hcd_secondary_injected_mw - - current_drive_variables.p_hcd_secondary_extra_heat_mw - ) - * 1.0e6 - ) - aux_current_fraction_fix = ( - auxiliary_cdfix / physics_variables.plasma_current - ) - elif current_drive_variables.i_hcd_secondary in [5, 8]: - # Account for first orbit losses - # (power due to particles that are ionised but not thermalised) [MW]: - # This includes a second order term in shinethrough*(first orbit loss) - current_drive_variables.f_p_beam_orbit_loss = min( - 0.999, current_drive_variables.f_p_beam_orbit_loss - ) # Should never be needed + # Calculate beam stopping cross-section - pnbitotfix = current_drive_variables.p_hcd_secondary_injected_mw / ( - 1.0e0 - - current_drive_variables.f_p_beam_orbit_loss - + current_drive_variables.f_p_beam_orbit_loss - * current_drive_variables.f_p_beam_shine_through - ) + sigstop = self.sigbeam( + current_drive_variables.e_beam_kev / physics_variables.m_beam_amu, + physics_variables.te, + physics_variables.dene, + physics_variables.f_nd_alpha_electron, + physics_variables.rncne, + physics_variables.rnone, + physics_variables.rnfene, + ) - # Shinethrough power (atoms that are not ionised) [MW]: - nbshinemwfix = ( - pnbitotfix * current_drive_variables.f_p_beam_shine_through - ) + # Calculate number of decay lengths to centre - # First orbit loss - porbitlossmwfix = current_drive_variables.f_p_beam_orbit_loss * ( - pnbitotfix - nbshinemwfix - ) + current_drive_variables.n_beam_decay_lengths_core = ( + dpath * physics_variables.dnla * sigstop + ) - # Power deposited - pinjmwfix = pnbitotfix - nbshinemwfix - porbitlossmwfix - pinjimwfix = pinjmwfix * current_drive_variables.f_p_beam_injected_ions - pinjemwfix = pinjmwfix * ( - 1.0e0 - current_drive_variables.f_p_beam_injected_ions - ) + # Shine-through fraction of beam - current_drive_variables.pwpnb = ( - pnbitotfix / current_drive_variables.eta_beam_injector_wall_plug - ) # neutral beam wall plug power - heat_transport_variables.pinjwpfix = current_drive_variables.pwpnb - current_drive_variables.eta_hcd_secondary_injector_wall_plug = ( - current_drive_variables.eta_beam_injector_wall_plug - ) - gamnb = effnbssfix * (dene20 * physics_variables.rmajor) - gamcdfix = gamnb - beam_current_fix = ( - 1.0e-3 * (pnbitotfix * 1.0e6) / current_drive_variables.e_beam_kev - ) # Neutral beam current (A) - auxiliary_cdfix = ( - effnbssfix - * ( - current_drive_variables.p_hcd_secondary_injected_mw - - current_drive_variables.p_hcd_secondary_extra_heat_mw - ) - * 1.0e6 - ) - aux_current_fraction_fix = ( - auxiliary_cdfix / physics_variables.plasma_current - ) - - # Fenstermacher Lower Hybrid model - if current_drive_variables.i_hcd_primary == 1: - effrfss = ( - (0.36e0 * (1.0e0 + (physics_variables.te / 25.0e0) ** 1.16e0)) - / (physics_variables.rmajor * dene20) - * current_drive_variables.feffcd - ) - current_drive_variables.eta_cd_hcd_primary = effrfss - # Ion-Cyclotron current drive - elif current_drive_variables.i_hcd_primary == 2: - effrfss = ( - 0.63e0 - * 0.1e0 - * physics_variables.ten - / (2.0e0 + physics_variables.zeff) - / (physics_variables.rmajor * dene20) - * current_drive_variables.feffcd - ) - current_drive_variables.eta_cd_hcd_primary = effrfss - # Fenstermacher Electron Cyclotron Resonance model - elif current_drive_variables.i_hcd_primary == 3: - effrfss = ( - 0.21e0 - * physics_variables.ten - / (physics_variables.rmajor * dene20 * physics_variables.dlamee) - * current_drive_variables.feffcd - ) - current_drive_variables.eta_cd_hcd_primary = effrfss - # Ehst Lower Hybrid / Fast Wave current drive - elif current_drive_variables.i_hcd_primary == 4: - effrfss = ( - physics_variables.te**0.77e0 - * (0.034e0 + 0.196e0 * physics_variables.beta) - / (physics_variables.rmajor * dene20) - * ( - 32.0e0 / (5.0e0 + physics_variables.zeff) - + 2.0e0 - + (12.0e0 * (6.0e0 + physics_variables.zeff)) - / (5.0e0 + physics_variables.zeff) - / (3.0e0 + physics_variables.zeff) - + 3.76e0 / physics_variables.zeff - ) - / 12.507e0 - * current_drive_variables.feffcd - ) - current_drive_variables.eta_cd_hcd_primary = effrfss - # ITER Neutral Beam current drive - elif current_drive_variables.i_hcd_primary == 5: - ( - effnbss, - current_drive_variables.f_p_beam_injected_ions, - current_drive_variables.f_p_beam_shine_through, - ) = self.iternb() - effnbss = effnbss * current_drive_variables.feffcd - current_drive_variables.eta_cd_hcd_primary = effnbss - # Culham Lower Hybrid current drive model - elif current_drive_variables.i_hcd_primary == 6: - effrfss = self.cullhy() - effrfss = effrfss * current_drive_variables.feffcd - current_drive_variables.eta_cd_hcd_primary = effrfss - # Culham ECCD model - elif current_drive_variables.i_hcd_primary == 7: - effrfss = self.culecd() - effrfss = effrfss * current_drive_variables.feffcd - current_drive_variables.eta_cd_hcd_primary = effrfss - # Culham Neutral Beam model - elif current_drive_variables.i_hcd_primary == 8: - ( - effnbss, - current_drive_variables.f_p_beam_injected_ions, - current_drive_variables.f_p_beam_shine_through, - ) = self.culnbi() - effnbss = effnbss * current_drive_variables.feffcd - current_drive_variables.eta_cd_hcd_primary = effnbss - # ECRH user input gamma - elif current_drive_variables.i_hcd_primary == 10: - current_drive_variables.eta_cd_norm_hcd_primary = ( - current_drive_variables.eta_cd_norm_ecrh - ) - - # Absolute current drive efficiency - effrfss = current_drive_variables.eta_cd_norm_hcd_primary / ( - dene20 * physics_variables.rmajor - ) - current_drive_variables.eta_cd_hcd_primary = effrfss - # EBW scaling - elif current_drive_variables.i_hcd_primary == 12: - # Scaling author Simon Freethy - # Ref : PROCESS issue 1262 - - # Normalised current drive efficiency gamma - current_drive_variables.eta_cd_norm_hcd_primary = ( - current_drive_variables.xi_ebw / 32.7e0 - ) * physics_variables.te - - # Absolute current drive efficiency - effrfss = current_drive_variables.eta_cd_norm_hcd_primary / ( - dene20 * physics_variables.rmajor - ) - current_drive_variables.eta_cd_hcd_primary = effrfss - # EBWs can only couple to plasma if cyclotron harmonic is above plasma density cut-off; - # this behaviour is captured in the following function (ref issue #1262): - # current_drive_variables.n_ecrh_harmonic = cyclotron harmonic number (fundamental used as default) - # contant 'a' controls sharpness of transition - a = 0.1e0 - - fc = ( - 1.0e0 - / (2.0e0 * np.pi) - * current_drive_variables.n_ecrh_harmonic - * constants.electron_charge - * physics_variables.bt - / constants.electron_mass - ) - fp = ( - 1.0e0 - / (2.0e0 * np.pi) - * np.sqrt( - physics_variables.dene - * constants.electron_charge**2 - / (constants.electron_mass * constants.epsilon0) - ) - ) + fshine = np.exp(-2.0e0 * dpath * physics_variables.dnla * sigstop) + fshine = max(fshine, 1.0e-20) - density_factor = 0.5e0 * ( - 1.0e0 + np.tanh((2.0e0 / a) * ((fp - fc) / fp - a)) - ) + # Deuterium and tritium beam densities - current_drive_variables.eta_cd_hcd_primary = ( - current_drive_variables.eta_cd_hcd_primary * density_factor - ) - effrfss = effrfss * density_factor + dend = physics_variables.nd_fuel_ions * ( + 1.0e0 - current_drive_variables.f_beam_tritium + ) + dent = physics_variables.nd_fuel_ions * current_drive_variables.f_beam_tritium - elif current_drive_variables.i_hcd_primary == 13: - # ECCD model for O-mode cut-off with added Te and Zeff dependance - # Scaling author: Simon Freethy - # Ref : PROCESS issue #2994 + # Power split to ions / electrons - fc = ( - 1 - / (2 * np.pi) - * constants.electron_charge - * physics_variables.bt - / constants.electron_mass - ) - fp = ( - 1 - / (2 * np.pi) - * np.sqrt( - (physics_variables.dene / 1.0e19) - * constants.electron_charge**2 - / (constants.electron_mass * constants.epsilon0) - ) - ) + f_p_beam_injected_ions = self.cfnbi( + physics_variables.m_beam_amu, + current_drive_variables.e_beam_kev, + physics_variables.ten, + physics_variables.dene, + dend, + dent, + physics_variables.zeffai, + physics_variables.dlamie, + ) - xi_CD = 0.18e0 # This is tuned to the results of a GRAY study - xi_CD = xi_CD * ( - 4.8e0 / (2 + physics_variables.zeff) - ) # Zeff correction - effrfss = ( - xi_CD - * physics_variables.te - / ( - 3.27e0 - * physics_variables.rmajor - * (physics_variables.dene / 1.0e19) - ) - ) + # Current drive efficiency - # O-mode case - if current_drive_variables.i_ecrh_wave_mode == 0: - f_cutoff = fp + effnbss = self.etanb2( + physics_variables.m_beam_amu, + physics_variables.alphan, + physics_variables.alphat, + physics_variables.aspect, + physics_variables.dene, + physics_variables.dnla, + current_drive_variables.e_beam_kev, + current_drive_variables.frbeam, + fshine, + physics_variables.rmajor, + physics_variables.rminor, + physics_variables.ten, + physics_variables.zeff, + ) - # X-mode case - elif current_drive_variables.i_ecrh_wave_mode == 1: - f_cutoff = 0.5 * ( - fc - + np.sqrt( - current_drive_variables.n_ecrh_harmonic * fc**2 + 4 * fp**2 - ) - ) + return effnbss, f_p_beam_injected_ions, fshine - # 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 - cutoff_factor = 0.5 * ( - 1 - + np.tanh( - (2 / (a)) - * ( - (current_drive_variables.n_ecrh_harmonic * fc - f_cutoff) - / fp - - a - ) - ) - ) - current_drive_variables.eta_cd_hcd_primary = effrfss * cutoff_factor - else: - raise ProcessValueError( - f"Current drive switch is invalid: {current_drive_variables.i_hcd_primary = }" - ) + def etanb2( + self, + m_beam_amu, + alphan, + alphat, + aspect, + dene, + dnla, + e_beam_kev, + frbeam, + fshine, + rmajor, + rminor, + ten, + zeff, + ): + """Routine to find neutral beam current drive efficiency + using the ITER 1990 formulation, plus correction terms + outlined in Culham Report AEA FUS 172 + author: P J Knight, CCFE, Culham Science Centre + m_beam_amu : input real : beam ion mass (amu) + alphan : input real : density profile factor + alphat : input real : temperature profile factor + aspect : input real : aspect ratio + dene : input real : volume averaged electron density (m**-3) + dnla : input real : line averaged electron density (m**-3) + e_beam_kev : input real : neutral beam energy (keV) + frbeam : input real : R_tangent / R_major for neutral beam injection + fshine : input real : shine-through fraction of beam + rmajor : input real : plasma major radius (m) + rminor : input real : plasma minor radius (m) + ten : input real : density weighted average electron temperature (keV) + zeff : input real : plasma effective charge + 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. +
The formulae are from AEA FUS 172, unless denoted by IPDG89. + AEA FUS 172: Physics Assessment for the European Reactor Study + ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, + ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 + """ + # Charge of beam ions + zbeam = 1.0 - # Compute current drive wall plug and injected powers (MW) and efficiencies - auxiliary_cd = ( - physics_variables.aux_current_fraction - * physics_variables.plasma_current - ) + # Fitting factor (IPDG89) + bbd = 1.0 - # LHCD or ICCD - if current_drive_variables.i_hcd_primary in [1, 2, 4, 6]: - # Injected power - current_drive_variables.plhybd = ( - 1.0e-6 - * ( - physics_variables.aux_current_fraction - - aux_current_fraction_fix - ) - * physics_variables.plasma_current - / effrfss - + current_drive_variables.p_hcd_primary_extra_heat_mw - ) - pinjimw1 = 0.0e0 - pinjemw1 = current_drive_variables.plhybd + # Volume averaged electron density (10**20 m**-3) + dene20 = dene / 1e20 - # Wall plug power - current_drive_variables.pwplh = ( - current_drive_variables.plhybd - / current_drive_variables.eta_lowhyb_injector_wall_plug - ) - pinjwp1 = current_drive_variables.pwplh + # Line averaged electron density (10**20 m**-3) + dnla20 = dnla / 1e20 - # Wall plug to injector efficiency - current_drive_variables.eta_hcd_primary_injector_wall_plug = ( - current_drive_variables.eta_lowhyb_injector_wall_plug - ) + # Critical energy (MeV) (power to electrons = power to ions) (IPDG89) + # N.B. ten is in keV + ecrit = 0.01 * m_beam_amu * ten - # Normalised current drive efficiency gamma - gamrf = effrfss * (dene20 * physics_variables.rmajor) - current_drive_variables.eta_cd_norm_hcd_primary = gamrf - # ECCD - elif current_drive_variables.i_hcd_primary in [3, 7, 10, 12, 13]: - # Injected power (set to close to close the Steady-state current equilibrium) - current_drive_variables.p_ecrh_injected_mw = ( - 1.0e-6 - * ( - physics_variables.aux_current_fraction - - aux_current_fraction_fix - ) - * physics_variables.plasma_current - / effrfss - + current_drive_variables.p_hcd_primary_extra_heat_mw - ) - pinjemw1 = current_drive_variables.p_ecrh_injected_mw + # Beam energy in MeV + ebmev = e_beam_kev / 1e3 - # Wall plug power - current_drive_variables.echwpow = ( - current_drive_variables.p_ecrh_injected_mw - / current_drive_variables.eta_ecrh_injector_wall_plug - ) + # x and y coefficients of function J0(x,y) (IPDG89) + xjs = ebmev / (bbd * ecrit) + xj = np.sqrt(xjs) - # Wall plug to injector efficiency - pinjwp1 = current_drive_variables.echwpow - current_drive_variables.eta_hcd_primary_injector_wall_plug = ( - current_drive_variables.eta_ecrh_injector_wall_plug - ) - elif current_drive_variables.i_hcd_primary in [5, 8]: - # MDK. See Gitlab issue #248, and scanned note. - power1 = ( - 1.0e-6 - * ( - physics_variables.aux_current_fraction - - aux_current_fraction_fix - ) - * physics_variables.plasma_current - / effnbss - + current_drive_variables.p_hcd_primary_extra_heat_mw - ) + yj = 0.8 * zeff / m_beam_amu - # Account for first orbit losses - # (power due to particles that are ionised but not thermalised) [MW]: - # This includes a second order term in shinethrough*(first orbit loss) - current_drive_variables.f_p_beam_orbit_loss = min( - 0.999, current_drive_variables.f_p_beam_orbit_loss - ) # Should never be needed + # Fitting function J0(x,y) + j0 = xjs / (4.0 + 3.0 * yj + xjs * (xj + 1.39 + 0.61 * yj**0.7)) - current_drive_variables.p_beam_injected_mw = power1 / ( - 1.0e0 - - current_drive_variables.f_p_beam_orbit_loss - + current_drive_variables.f_p_beam_orbit_loss - * current_drive_variables.f_p_beam_shine_through - ) + # Effective inverse aspect ratio, with a limit on its maximum value + epseff = min(0.2, (0.5 / aspect)) - # Shinethrough power (atoms that are not ionised) [MW]: - current_drive_variables.p_beam_shine_through_mw = ( - current_drive_variables.p_beam_injected_mw - * current_drive_variables.f_p_beam_shine_through - ) + # Reduction in the reverse electron current + # due to neoclassical effects + gfac = (1.55 + 0.85 / zeff) * np.sqrt(epseff) - (0.2 + 1.55 / zeff) * epseff - # First orbit loss - current_drive_variables.p_beam_orbit_loss_mw = ( - current_drive_variables.f_p_beam_orbit_loss - * ( - current_drive_variables.p_beam_injected_mw - - current_drive_variables.p_beam_shine_through_mw - ) - ) + # Reduction in the net beam driven current + # due to the reverse electron current + ffac = 1.0 - (zbeam / zeff) * (1.0 - gfac) - # Power deposited - pinjmw1 = ( - current_drive_variables.p_beam_injected_mw - - current_drive_variables.p_beam_shine_through_mw - - current_drive_variables.p_beam_orbit_loss_mw - ) - pinjimw1 = pinjmw1 * current_drive_variables.f_p_beam_injected_ions - pinjemw1 = pinjmw1 * ( - 1.0e0 - current_drive_variables.f_p_beam_injected_ions - ) + # Normalisation to allow results to be valid for + # non-ITER plasma size and density: - current_drive_variables.pwpnb = ( - current_drive_variables.p_beam_injected_mw - / current_drive_variables.eta_beam_injector_wall_plug - ) # neutral beam wall plug power - pinjwp1 = current_drive_variables.pwpnb - current_drive_variables.eta_hcd_primary_injector_wall_plug = ( - current_drive_variables.eta_beam_injector_wall_plug - ) - gamnb = effnbss * (dene20 * physics_variables.rmajor) - current_drive_variables.eta_cd_norm_hcd_primary = gamnb - current_drive_variables.c_beam_total = ( - 1.0e-3 - * (current_drive_variables.p_beam_injected_mw * 1.0e6) - / current_drive_variables.e_beam_kev - ) # Neutral beam current (A) + # Line averaged electron density (10**20 m**-3) normalised to ITER + nnorm = 1.0 - # Total injected power - # sum contributions from primary and secondary systems - current_drive_variables.p_hcd_injected_total_mw = ( - pinjemw1 + pinjimw1 + pinjemwfix + pinjimwfix - ) - pinjmw1 = pinjemw1 + pinjimw1 - pinjmwfix = pinjemwfix + pinjimwfix - current_drive_variables.p_hcd_injected_electrons_mw = pinjemw1 + pinjemwfix - current_drive_variables.p_hcd_injected_ions_mw = pinjimw1 + pinjimwfix - heat_transport_variables.pinjwp = ( - pinjwp1 + heat_transport_variables.pinjwpfix + # Distance along beam to plasma centre + r = max(rmajor, rmajor * frbeam) + eps1 = rminor / r + + if (1.0 + eps1) < frbeam: + raise ProcessValueError( + "Imminent negative square root argument; NBI will miss plasma completely", + eps=eps1, + frbeam=frbeam, ) - # Reset injected power to zero for ignited plasma (fudge) - if physics_variables.ignite == 1: - heat_transport_variables.pinjwp = 0.0e0 + d = rmajor * np.sqrt((1.0 + eps1) ** 2 - frbeam**2) - # Ratio of fusion to input (injection+ohmic) power - if ( - abs( - current_drive_variables.p_hcd_injected_total_mw - + current_drive_variables.p_beam_orbit_loss_mw - + physics_variables.p_plasma_ohmic_mw - ) - < 1.0e-6 - ): - current_drive_variables.bigq = 1.0e18 - else: - current_drive_variables.bigq = physics_variables.fusion_power / ( - current_drive_variables.p_hcd_injected_total_mw - + current_drive_variables.p_beam_orbit_loss_mw - + physics_variables.p_plasma_ohmic_mw - ) + # Distance along beam to plasma centre for ITER + # assuming a tangency radius equal to the major radius + epsitr = 2.15 / 6.0 + dnorm = 6.0 * np.sqrt(2.0 * epsitr + epsitr**2) - if not output: - return + # Normalisation to beam energy (assumes a simplified formula for + # the beam stopping cross-section) + ebnorm = ebmev * ((nnorm * dnorm) / (dnla20 * d)) ** (1.0 / 0.78) - po.oheadr(self.outfile, "Current Drive System") + # A_bd fitting coefficient, after normalisation with ebnorm + abd = ( + 0.107 + * (1.0 - 0.35 * alphan + 0.14 * alphan**2) + * (1.0 - 0.21 * alphat) + * (1.0 - 0.2 * ebnorm + 0.09 * ebnorm**2) + ) - if current_drive_variables.i_hcd_calculations == 0: - po.ocmmnt(self.outfile, "No current drive used") - po.oblnkl(self.outfile) - return + # Normalised current drive efficiency (A/W m**-2) (IPDG89) + gamnb = 5.0 * abd * 0.1 * ten * (1.0 - fshine) * frbeam * j0 / 0.2 * ffac - if current_drive_variables.i_hcd_primary in [1, 4, 6]: - po.ocmmnt(self.outfile, "Lower Hybrid Current Drive") - elif current_drive_variables.i_hcd_primary == 2: - po.ocmmnt(self.outfile, "Ion Cyclotron Current Drive") - elif current_drive_variables.i_hcd_primary in [3, 7]: - po.ocmmnt(self.outfile, "Electron Cyclotron Current Drive") - elif current_drive_variables.i_hcd_primary in [5, 8]: - po.ocmmnt(self.outfile, "Neutral Beam Current Drive") - elif current_drive_variables.i_hcd_primary == 10: - po.ocmmnt( - self.outfile, "Electron Cyclotron Current Drive (user input gamma_CD)" - ) - elif current_drive_variables.i_hcd_primary == 12: - po.ocmmnt(self.outfile, "EBW current drive") - elif current_drive_variables.i_hcd_primary == 13: - po.ocmmnt( - self.outfile, - "Electron Cyclotron Current Drive (O-mode cutoff with Zeff & Te)", - ) + # Current drive efficiency (A/W) + return gamnb / (dene20 * rmajor) - po.ovarin( - self.outfile, - "Current drive efficiency model", - "(i_hcd_primary)", - current_drive_variables.i_hcd_primary, - ) + def etanb(self, m_beam_amu, alphan, alphat, aspect, dene, ebeam, rmajor, ten, zeff): + """Routine to find neutral beam current drive efficiency + using the ITER 1990 formulation + author: P J Knight, CCFE, Culham Science Centre + m_beam_amu : input real : beam ion mass (amu) + alphan : input real : density profile factor + alphat : input real : temperature profile factor + aspect : input real : aspect ratio + dene : input real : volume averaged electron density (m**-3) + 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 + This routine calculates the current drive efficiency of + a neutral beam system, based on the 1990 ITER model. + ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, + ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 + """ - if current_drive_variables.i_hcd_secondary in [1, 4, 6]: - po.ocmmnt(self.outfile, "Lower Hybrid Current Drive") - elif current_drive_variables.i_hcd_secondary == 2: - po.ocmmnt(self.outfile, "Ion Cyclotron Current Drive") - elif current_drive_variables.i_hcd_secondary in [3, 7]: - po.ocmmnt(self.outfile, "Electron Cyclotron Current Drive") - elif current_drive_variables.i_hcd_secondary in [5, 8]: - po.ocmmnt(self.outfile, "Neutral Beam Current Drive") - elif current_drive_variables.i_hcd_secondary == 10: - po.ocmmnt( - self.outfile, "Electron Cyclotron Current Drive (user input gamma_CD)" - ) - elif current_drive_variables.i_hcd_secondary == 12: - po.ocmmnt(self.outfile, "EBW current drive") - elif current_drive_variables.i_hcd_secondary == 13: - po.ocmmnt( - self.outfile, - "Electron Cyclotron Current Drive (O-mode cutoff with Zeff & Te)", - ) + zbeam = 1.0 + bbd = 1.0 - po.ovarin( - self.outfile, - "Secondary current drive efficiency model", - "(i_hcd_secondary)", - current_drive_variables.i_hcd_secondary, - ) + dene20 = 1e-20 * dene - if physics_variables.ignite == 1: - po.ocmmnt( - self.outfile, - "Ignited plasma; injected power only used for start-up phase", - ) + # Ratio of E_beam/E_crit + xjs = ebeam / (bbd * 10.0 * m_beam_amu * ten) + xj = np.sqrt(xjs) - po.oblnkl(self.outfile) + yj = 0.8 * zeff / m_beam_amu - if abs(physics_variables.inductive_current_fraction) > 1.0e-8: - po.ocmmnt(self.outfile, "Current is driven by both inductive") - po.ocmmnt(self.outfile, "and non-inductive means.") + rjfunc = xjs / (4.0 + 3.0 * yj + xjs * (xj + 1.39 + 0.61 * yj**0.7)) - po.ovarre( - self.outfile, - "Ratio of power for flat-top to start-up (MW)", - "(startupratio)", - cost_variables.startupratio, - ) - po.ovarre( - self.outfile, - "Auxiliary power used for plasma heating only (MW)", - "(p_hcd_primary_extra_heat_mw)", - current_drive_variables.p_hcd_primary_extra_heat_mw - + current_drive_variables.p_hcd_secondary_extra_heat_mw, - ) - po.ovarre( - self.outfile, - "Power injected for current drive (MW)", - "(pcurrentdrivemw)", - current_drive_variables.p_hcd_injected_total_mw - - current_drive_variables.p_hcd_primary_extra_heat_mw - - current_drive_variables.p_hcd_secondary_extra_heat_mw, - ) - po.ovarre( - self.outfile, - "Maximum Allowed Bootstrap current fraction", - "(f_c_plasma_bootstrap_max)", - current_drive_variables.f_c_plasma_bootstrap_max, - ) - if current_drive_variables.i_hcd_secondary != 0: - po.ovarre( - self.outfile, - "Power injected for main current drive (MW)", - "(pcurrentdrivemw1)", - pinjmw1 - current_drive_variables.p_hcd_primary_extra_heat_mw, - ) - po.ovarre( - self.outfile, - "Power injected for secondary current drive (MW)", - "(pcurrentdrivemw2)", - pinjmwfix - current_drive_variables.p_hcd_secondary_extra_heat_mw, - ) + epseff = 0.5 / aspect + gfac = (1.55 + 0.85 / zeff) * np.sqrt(epseff) - (0.2 + 1.55 / zeff) * epseff + ffac = 1.0 / zbeam - (1.0 - gfac) / zeff - po.ovarre( - self.outfile, - "Fusion gain factor Q", - "(bigq)", - current_drive_variables.bigq, - "OP ", + abd = ( + 0.107 + * (1.0 - 0.35 * alphan + 0.14 * alphan**2) + * (1.0 - 0.21 * alphat) + * (1.0 - 0.2e-3 * ebeam + 0.09e-6 * ebeam**2) ) - po.ovarre( - self.outfile, - "Auxiliary current drive (A)", - "(auxiliary_cd)", - auxiliary_cd, - "OP ", + + return abd * (5.0 / rmajor) * (0.1 * ten / dene20) * rjfunc / 0.2 * ffac + + def sigbeam(self, eb, te, ne, rnhe, rnc, rno, rnfe): + """Calculates the stopping cross-section for a hydrogen + beam in a fusion plasma + author: P J Knight, CCFE, Culham Science Centre + eb : input real : beam energy (kev/amu) + te : input real : electron temperature (keV) + ne : input real : electron density (10^20m-3) + rnhe : input real : alpha density / ne + rnc : input real : carbon density /ne + rno : input real : oxygen density /ne + rnfe : input real : iron density /ne + This function calculates the stopping cross-section (m^-2) + for a hydrogen beam in a fusion plasma. + Janev, Boley and Post, Nuclear Fusion 29 (1989) 2125 + """ + a = np.array([ + [ + [4.4, -2.49e-2], + [7.46e-2, 2.27e-3], + [3.16e-3, -2.78e-5], + ], + [ + [2.3e-1, -1.15e-2], + [-2.55e-3, -6.2e-4], + [1.32e-3, 3.38e-5], + ], + ]) + + b = np.array([ + [ + [[-2.36, -1.49, -1.41, -1.03], [0.185, -0.0154, -4.08e-4, 0.106]], + [ + [-0.25, -0.119, -0.108, -0.0558], + [-0.0381, -0.015, -0.0138, -3.72e-3], + ], + ], + [ + [ + [0.849, 0.518, 0.477, 0.322], + [-0.0478, 7.18e-3, 1.57e-3, -0.0375], + ], + [ + [0.0677, 0.0292, 0.0259, 0.0124], + [0.0105, 3.66e-3, 3.33e-3, 8.61e-4], + ], + ], + [ + [ + [-0.0588, -0.0336, -0.0305, -0.0187], + [4.34e-3, 3.41e-4, 7.35e-4, 3.53e-3], + ], + [ + [-4.48e-3, -1.79e-3, -1.57e-3, -7.43e-4], + [-6.76e-4, -2.04e-4, -1.86e-4, -5.12e-5], + ], + ], + ]) + + z = np.array([2.0, 6.0, 8.0, 26.0]) + nn = np.array([rnhe, rnc, rno, rnfe]) + + nen = ne * 1e-19 + + s1 = 0.0 + for k in range(2): + for j in range(3): + for i in range(2): + s1 += ( + a[i, j, k] + * (np.log(eb)) ** i + * (np.log(nen)) ** j + * (np.log(te)) ** k + ) + + sz = 0.0 + for l in range(4): # noqa: E741 + for k in range(2): + for j in range(2): + for i in range(3): + sz += ( + b[i, j, k, l] + * (np.log(eb)) ** i + * (np.log(nen)) ** j + * (np.log(te)) ** k + * nn[l] + * z[l] + * (z[l] - 1.0) + ) + + return max(1e-20 * (np.exp(s1) / eb * (1.0 + sz)), 1e-23) + + def cfnbi(self, afast, efast, te, ne, _nd, _nt, zeffai, xlmbda): + """Routine to calculate the fraction of the fast particle energy + coupled to the ions + author: P J Knight, CCFE, Culham Science Centre + afast : input real : mass of fast particle (units of proton mass) + efast : input real : energy of fast particle (keV) + te : input real : density weighted average electron temp. (keV) + ne : input real : volume averaged electron density (m**-3) + nd : input real : deuterium beam density (m**-3) + nt : input real : tritium beam density (m**-3) + zeffai : input real : mass weighted plasma effective charge + xlmbda : input real : ion-electron coulomb logarithm + f_p_beam_injected_ions : output real : fraction of fast particle energy coupled to ions + This routine calculates the fast particle energy coupled to + the ions in the neutral beam system. + """ + # atmd = 2.0 + atmdt = 2.5 + # atmt = 3.0 + c = 3.0e8 + me = constants.electron_mass + # zd = 1.0 + # zt = 1.0 + + # xlbd = self.xlmbdabi(afast, atmd, efast, te, ne) + # xlbt = self.xlmbdabi(afast, atmt, efast, te, ne) + + # sum = nd * zd * zd * xlbd / atmd + nt * zt * zt * xlbt / atmt + # ecritfix = 16.0e0 * te * afast * (sum / (ne * xlmbda)) ** (2.0e0 / 3.0e0) + + xlmbdai = self.xlmbdabi(afast, atmdt, efast, te, ne) + sumln = zeffai * xlmbdai / xlmbda + xlnrat = ( + 3.0e0 * np.sqrt(np.pi) / 4.0e0 * me / constants.proton_mass * sumln + ) ** (2.0e0 / 3.0e0) + ve = c * np.sqrt(2.0e0 * te / 511.0e0) + + ecritfi = ( + afast + * constants.proton_mass + * ve + * ve + * xlnrat + / (2.0e0 * constants.electron_charge * 1.0e3) ) - if current_drive_variables.i_hcd_secondary != 0: - po.ovarre( - self.outfile, - "Secondary auxiliary current drive (A)", - "(auxiliary_cdfix)", - auxiliary_cdfix, - "OP ", - ) - po.ovarre( - self.outfile, - "Current drive efficiency (A/W)", - "(eta_cd_hcd_primary)", - current_drive_variables.eta_cd_hcd_primary, - "OP ", + x = np.sqrt(efast / ecritfi) + t1 = np.log((x * x - x + 1.0e0) / ((x + 1.0e0) ** 2)) + thx = (2.0e0 * x - 1.0e0) / np.sqrt(3.0e0) + t2 = 2.0e0 * np.sqrt(3.0e0) * (np.arctan(thx) + np.pi / 6.0e0) + + return (t1 + t2) / (3.0e0 * x * x) + + def xlmbdabi(self, mb, mth, eb, t, nelec): + """Calculates the Coulomb logarithm for ion-ion collisions + author: P J Knight, CCFE, Culham Science Centre + mb : input real : mass of fast particle (units of proton mass) + mth : input real : mass of background ions (units of proton mass) + eb : input real : energy of fast particle (keV) + t : input real : density weighted average electron temp. (keV) + nelec : input real : volume averaged electron density (m**-3) + This function calculates the Coulomb logarithm for ion-ion + collisions where the relative velocity may be large compared + with the background ('mt') thermal velocity. + Mikkelson and Singer, Nuc Tech/Fus, 4, 237 (1983) + """ + + x1 = (t / 10.0) * (eb / 1000.0) * mb / (nelec / 1e20) + x2 = mth / (mth + mb) + + return 23.7 + np.log(x2 * np.sqrt(x1)) + + +class ElectronCyclotron: + def __init__(self, plasma_profile: PlasmaProfile): + self.outfile = constants.nout + self.plasma_profile = plasma_profile + + def culecd(self): + """Routine to calculate Electron Cyclotron current drive efficiency + author: M R O'Brien, CCFE, Culham Science Centre + author: P J Knight, CCFE, Culham Science Centre + effrfss : output real : electron cyclotron current drive efficiency (A/W) + This routine calculates the current drive parameters for a + electron cyclotron system, based on the AEA FUS 172 model. + AEA FUS 172: Physics Assessment for the European Reactor Study + """ + rrr = 1.0e0 / 3.0e0 + + # Temperature + tlocal = self.plasma_profile.teprofile.calculate_profile_y( + rrr, + physics_variables.rhopedt, + physics_variables.te0, + physics_variables.teped, + physics_variables.tesep, + physics_variables.alphat, + physics_variables.tbeta, ) - po.ovarre( - self.outfile, - "Normalised current drive efficiency of primary HCD system (10^20 A / W m^2)", - "(eta_cd_norm_hcd_primary)", - current_drive_variables.eta_cd_norm_hcd_primary, - "OP ", + + # Density (10**20 m**-3) + dlocal = 1.0e-20 * self.plasma_profile.neprofile.calculate_profile_y( + rrr, + physics_variables.rhopedn, + physics_variables.ne0, + physics_variables.neped, + physics_variables.nesep, + physics_variables.alphan, ) - po.ovarre( - self.outfile, - "Wall plug to injector efficiency", - "(eta_hcd_primary_injector_wall_plug)", - current_drive_variables.eta_hcd_primary_injector_wall_plug, + + # Inverse aspect ratio + epsloc = rrr * physics_variables.rminor / physics_variables.rmajor + + # Effective charge (use average value) + zlocal = physics_variables.zeff + + # Coulomb logarithm for ion-electron collisions + # (From J. A. Wesson, 'Tokamaks', Clarendon Press, Oxford, p.293) + coulog = 15.2e0 - 0.5e0 * np.log(dlocal) + np.log(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. + 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) + + # Normalised current drive efficiency (A/W m**-2) + ecgam = 0.25e0 * (ecgam1 + ecgam2 + ecgam3 + ecgam4) + + # Current drive efficiency (A/W) + return ecgam / (dlocal * physics_variables.rmajor) + + def eccdef(self, tlocal, epsloc, zlocal, cosang, coulog): + """Routine to calculate Electron Cyclotron current drive efficiency + author: M R O'Brien, CCFE, Culham Science Centre + author: P J Knight, CCFE, Culham Science Centre + tlocal : input real : local electron temperature (keV) + epsloc : input real : local inverse aspect ratio + zlocal : input real : local plasma effective charge + cosang : input real : cosine of the poloidal angle at which ECCD takes + place (+1 outside, -1 inside) + coulog : input real : local coulomb logarithm for ion-electron collisions + ecgam : output real : normalised current drive efficiency (A/W m**-2) + 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 nIR/P with n the
+ local density in 10**20 /m**3, I the driven current in MAmps,
+ R the major radius in metres, and P the absorbed power in MWatts.
+ AEA FUS 172: Physics Assessment for the European Reactor Study
+ ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al,
+ ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990
+ """
+ mcsq = (
+ constants.electron_mass * 2.9979e8**2 / (1.0e3 * constants.electron_volt)
+ ) # keV
+ f = 16.0e0 * (tlocal / mcsq) ** 2
+
+ # fp is the derivative of f with respect to gamma, the relativistic
+ # factor, taken equal to 1 + 2T/(m c**2)
+
+ fp = 16.0e0 * tlocal / mcsq
+
+ # 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) = -8/(1+zlocal). alpha is of the
+ # form (-1/2 + ix), with x a real number and i = sqrt(-1).
+
+ lam = 1.0e0
+ palpha, palphap = self.legend(zlocal, lam)
+
+ lams = np.sqrt(2.0e0 * epsloc / (1.0e0 + epsloc))
+ palphas, _ = self.legend(zlocal, lams)
+
+ # hp is the derivative of IPDG89's h function with respect to lam
+
+ h = -4.0e0 * lam / (zlocal + 5.0e0) * (1.0e0 - lams * palpha / (lam * palphas))
+ hp = -4.0e0 / (zlocal + 5.0e0) * (1.0e0 - lams * palphap / palphas)
+
+ # facm is IPDG89's momentum conserving factor
+
+ facm = 1.5e0
+ y = mcsq / (2.0e0 * tlocal) * (1.0e0 + epsloc * cosang)
+
+ # We take the negative of the IPDG89 expression to get a positive
+ # number
+
+ ecgam = (
+ -7.8e0
+ * facm
+ * np.sqrt((1.0e0 + epsloc) / (1.0e0 - epsloc))
+ / coulog
+ * (h * fp - 0.5e0 * y * f * hp)
)
- if current_drive_variables.i_hcd_primary == 10:
- po.ovarre(
- self.outfile,
- "ECRH plasma heating efficiency",
- "(eta_cd_norm_ecrh)",
- current_drive_variables.eta_cd_norm_ecrh,
- )
- if current_drive_variables.i_hcd_primary == 12:
- po.ovarre(
- self.outfile,
- "EBW plasma heating efficiency",
- "(xi_ebw)",
- current_drive_variables.xi_ebw,
- )
- if current_drive_variables.i_hcd_primary in [12, 13]:
- po.ovarre(
- self.outfile,
- "EC harmonic number",
- "(n_ecrh_harmonic)",
- current_drive_variables.n_ecrh_harmonic,
- )
- if current_drive_variables.i_hcd_primary == 13:
- po.ovarin(
- self.outfile,
- "EC cutoff wave mode switch",
- "(i_ecrh_wave_mode)",
- current_drive_variables.i_ecrh_wave_mode,
- )
+ if ecgam < 0.0e0:
+ raise ProcessValueError("Negative normalised current drive efficiency")
+ return ecgam
- if current_drive_variables.i_hcd_secondary != 0:
- po.ovarre(
- self.outfile,
- "Secondary current drive efficiency (A/W)",
- "(eta_cd_hcd_secondary)",
- eta_cd_hcd_secondary,
- "OP ",
- )
- po.ovarre(
- self.outfile,
- "Seconday wall plug to injector efficiency",
- "(eta_hcd_secondary_injector_wall_plug)",
- current_drive_variables.eta_hcd_secondary_injector_wall_plug,
- )
- po.ovarre(
- self.outfile,
- "Normalised secondary current drive efficiency, gamma (10^20 A/W-m2)",
- "(gamcdfix)",
- gamcdfix,
- "OP ",
- )
+ def electron_cyclotron_fenstermacher(
+ self,
+ ten: float,
+ rmajor: float,
+ dene20: float,
+ dlamee: float,
+ ) -> float:
+ """
+ Routine to calculate Fenstermacher Electron Cyclotron heating efficiency.
- po.osubhd(self.outfile, "Fractions of current drive :")
- po.ovarrf(
- self.outfile,
- "Bootstrap fraction",
- "(f_c_plasma_bootstrap)",
- current_drive_variables.f_c_plasma_bootstrap,
- "OP ",
- )
- po.ovarrf(
- self.outfile,
- "Diamagnetic fraction",
- "(f_c_plasma_diamagnetic)",
- current_drive_variables.f_c_plasma_diamagnetic,
- "OP ",
- )
- po.ovarrf(
- self.outfile,
- "Pfirsch-Schlueter fraction",
- "(f_c_plasma_pfirsch_schluter)",
- current_drive_variables.f_c_plasma_pfirsch_schluter,
- "OP ",
- )
- po.ovarrf(
- self.outfile,
- "Auxiliary current drive fraction",
- "(aux_current_fraction)",
- physics_variables.aux_current_fraction,
- "OP ",
- )
- po.ovarrf(
- self.outfile,
- "Inductive fraction",
- "(inductive_current_fraction)",
- physics_variables.inductive_current_fraction,
- "OP ",
- )
- # Add total error check.
- po.ovarrf(
- self.outfile,
- "Total",
- "(f_c_plasma_internal+aux_current_fraction+inductive_current_fraction)",
- current_drive_variables.f_c_plasma_internal
- + physics_variables.aux_current_fraction
- + physics_variables.inductive_current_fraction,
- )
- if (
- abs(
- current_drive_variables.f_c_plasma_internal
- + physics_variables.aux_current_fraction
- + physics_variables.inductive_current_fraction
- - 1.0e0
- )
- > 1.0e-8
- ):
- po.ocmmnt(self.outfile, "ERROR: current drive fractions do not add to 1")
+ :param ten: Density weighted average electron temperature keV.
+ :type ten: float
+ :param zeff: Plasma effective charge.
+ :type zeff: float
+ :param rmajor: Major radius of the plasma in meters.
+ :type rmajor: float
+ :param dene20: Volume averaged electron density in 1x10^20 m^-3.
+ :type dene20: float
+ :param dlamee: Electron collision frequency in 1/s.
+ :type dlamee: float
- # MDK Add physics_variables.fvsbrnni as it can be an iteration variable
- po.ovarrf(
- self.outfile,
- "Fraction of the plasma current produced by non-inductive means",
- "(fvsbrnni)",
- physics_variables.fvsbrnni,
- )
+ :return: The calculated electron cyclotron heating efficiency in A/W.
+ :rtype: float
- if (
- abs(
- current_drive_variables.f_c_plasma_bootstrap
- - current_drive_variables.f_c_plasma_bootstrap_max
- )
- < 1.0e-8
- ):
- po.ocmmnt(self.outfile, "Warning : bootstrap current fraction is at")
- po.ocmmnt(self.outfile, " its prescribed maximum.")
+ :notes:
- po.oblnkl(self.outfile)
+ :references:
+ - T.C. Hender et al., 'Physics Assessment of the European Reactor Study', AEA FUS 172, 1992.
- if abs(current_drive_variables.plhybd) > 1.0e-8:
- po.ovarre(self.outfile, "RF efficiency (A/W)", "(effrfss)", effrfss, "OP ")
- po.ovarre(self.outfile, "RF gamma (10^20 A/W-m2)", "(gamrf)", gamrf, "OP ")
- po.ovarre(
- self.outfile,
- "Lower hybrid injected power (MW)",
- "(plhybd)",
- current_drive_variables.plhybd,
- "OP ",
- )
- po.ovarre(
- self.outfile,
- "Lower hybrid wall plug efficiency",
- "(eta_lowhyb_injector_wall_plug)",
- current_drive_variables.eta_lowhyb_injector_wall_plug,
- )
- po.ovarre(
- self.outfile,
- "Lower hybrid wall plug power (MW)",
- "(pwplh)",
- current_drive_variables.pwplh,
- "OP ",
- )
+ """
- # MDK rearranged and added current_drive_variables.p_beam_shine_through_mw
- # if (abs(current_drive_variables.pnbeam) > 1.0e-8) :
- if (
- (current_drive_variables.i_hcd_primary == 5)
- or (current_drive_variables.i_hcd_primary == 8)
- or (current_drive_variables.i_hcd_secondary == 5)
- or (current_drive_variables.i_hcd_secondary == 8)
- ):
- po.ovarre(
- self.outfile,
- "Neutral beam energy (keV)",
- "(e_beam_kev)",
- current_drive_variables.e_beam_kev,
+ return (0.21e0 * ten) / (rmajor * dene20 * dlamee)
+
+ def electron_cyclotron_freethy(
+ self,
+ te: float,
+ zeff: float,
+ rmajor: float,
+ dene: float,
+ bt: float,
+ n_ecrh_harmonic: int,
+ i_ecrh_wave_mode: int,
+ ) -> float:
+ """
+ Calculate the Electron Cyclotron current drive efficiency using the Freethy model.
+
+ This function computes the ECCD efficiency based on the electron temperature,
+ effective charge, major radius, electron density, magnetic field, harmonic number,
+ and wave mode.
+
+ :param te: Volume averaged electron temperature in keV.
+ :type te: float
+ :param zeff: Plasma effective charge.
+ :type zeff: float
+ :param rmajor: Major radius of the plasma in meters.
+ :type rmajor: float
+ :param dene: Volume averaged electron density in m^-3.
+ :type dene: float
+ :param bt: Toroidal magnetic field in Tesla.
+ :type bt: float
+ :param n_ecrh_harmonic: Cyclotron harmonic number (fundamental used as default).
+ :type n_ecrh_harmonic: int
+ :param i_ecrh_wave_mode: Wave mode switch (0 for O-mode, 1 for X-mode).
+ :type i_ecrh_wave_mode: int
+
+ :return: The calculated absolute ECCD efficiency in A/W.
+ :rtype: float
+
+ :notes:
+ - Plasma coupling only occurs if the plasma cut-off is below the cyclotron harmonic.
+ - The density factor accounts for this behavior.
+
+ :references:
+ - Freethy, S., PROCESS issue #2994.
+ """
+
+ # Cyclotron frequency
+ fc = 1 / (2 * np.pi) * constants.electron_charge * bt / constants.electron_mass
+
+ # Plasma frequency
+ fp = (
+ 1
+ / (2 * np.pi)
+ * np.sqrt(
+ (dene / 1.0e19)
+ * constants.electron_charge**2
+ / (constants.electron_mass * constants.epsilon0)
)
- if (current_drive_variables.i_hcd_primary == 5) or (
- current_drive_variables.i_hcd_primary == 8
- ):
- po.ovarre(
- self.outfile,
- "Neutral beam current (A)",
- "(c_beam_total)",
- current_drive_variables.c_beam_total,
- "OP ",
- )
+ )
- if (current_drive_variables.i_hcd_secondary == 5) or (
- current_drive_variables.i_hcd_secondary == 8
- ):
- po.ovarre(
- self.outfile,
- "Secondary fixed neutral beam current (A)",
- "(beam_current_fix)",
- beam_current_fix,
- "OP ",
- )
+ # Scaling factor for ECCD efficiency
+ xi_CD = 0.18e0 # Tuned to the results of a GRAY study
+ xi_CD *= 4.8e0 / (2 + zeff) # Zeff correction
- if (current_drive_variables.i_hcd_primary == 5) or (
- current_drive_variables.i_hcd_primary == 8
- ):
- po.ovarre(
- self.outfile, "Beam efficiency (A/W)", "(effnbss)", effnbss, "OP "
- )
+ # ECCD efficiency
+ eta_cd = xi_CD * te / (3.27e0 * rmajor * (dene / 1.0e19))
- if (current_drive_variables.i_hcd_secondary == 5) or (
- current_drive_variables.i_hcd_secondary == 8
- ):
- po.ovarre(
- self.outfile,
- "Secondary fixed beam efficiency (A/W)",
- "(effnbssfix)",
- effnbssfix,
- "OP ",
- )
+ # Determine the cut-off frequency based on wave mode
+ if i_ecrh_wave_mode == 0: # O-mode case
+ f_cutoff = fp
+ elif i_ecrh_wave_mode == 1: # X-mode case
+ f_cutoff = 0.5 * (fc + np.sqrt(n_ecrh_harmonic * fc**2 + 4 * fp**2))
+ else:
+ raise ValueError("Invalid wave mode. Use 0 for O-mode or 1 for X-mode.")
- po.ovarre(
- self.outfile, "Beam gamma (10^20 A/W-m2)", "(gamnb)", gamnb, "OP "
- )
- po.ovarre(
- self.outfile,
- "Neutral beam wall plug efficiency",
- "(eta_beam_injector_wall_plug)",
- current_drive_variables.eta_beam_injector_wall_plug,
- )
- po.ovarre(
- self.outfile,
- "Beam decay lengths to centre",
- "(n_beam_decay_lengths_core)",
- current_drive_variables.n_beam_decay_lengths_core,
- "OP ",
- )
- po.ovarre(
- self.outfile,
- "Beam shine-through fraction",
- "(f_p_beam_shine_through)",
- current_drive_variables.f_p_beam_shine_through,
- "OP ",
- )
- po.ovarre(
- self.outfile,
- "Neutral beam wall plug power (MW)",
- "(pwpnb)",
- current_drive_variables.pwpnb,
- "OP ",
- )
+ # Plasma coupling factor
+ a = 0.1 # Controls sharpness of the transition
+ cutoff_factor = 0.5 * (
+ 1 + np.tanh((2 / a) * ((n_ecrh_harmonic * fc - f_cutoff) / fp - a))
+ )
- po.oblnkl(self.outfile)
- po.ocmmnt(self.outfile, "Neutral beam power balance :")
- po.ocmmnt(self.outfile, "----------------------------")
- if (current_drive_variables.i_hcd_primary == 5) or (
- current_drive_variables.i_hcd_primary == 8
- ):
- po.ovarrf(
- self.outfile,
- "Beam first orbit loss power (MW)",
- "(p_beam_orbit_loss_mw)",
- current_drive_variables.p_beam_orbit_loss_mw,
- "OP ",
- )
- po.ovarrf(
- self.outfile,
- "Beam shine-through power [MW]",
- "(p_beam_shine_through_mw)",
- current_drive_variables.p_beam_shine_through_mw,
- "OP ",
- )
- po.ovarrf(
- self.outfile,
- "Beam power deposited in plasma (MW)",
- "(p_hcd_injected_total_mw)",
- pinjmw1,
- "OP ",
- )
- po.ovarrf(
- self.outfile,
- "Maximum allowable beam power (MW)",
- "(p_hcd_injected_max)",
- current_drive_variables.p_hcd_injected_max,
- )
- po.ovarrf(
- self.outfile,
- "Total (MW)",
- "(current_drive_variables.p_beam_orbit_loss_mw+current_drive_variables.p_beam_shine_through_mw+current_drive_variables.p_hcd_injected_total_mw)",
- current_drive_variables.p_beam_orbit_loss_mw
- + current_drive_variables.p_beam_shine_through_mw
- + pinjmw1,
- )
- po.oblnkl(self.outfile)
- po.ovarrf(
- self.outfile,
- "Beam power entering vacuum vessel (MW)",
- "(p_beam_injected_mw)",
- current_drive_variables.p_beam_injected_mw,
- "OP ",
- )
+ # Final ECCD efficiency
+ return eta_cd * cutoff_factor
+
+ def legend(self, zlocal, arg):
+ """Routine to calculate Legendre function and its derivative
+ author: M R O'Brien, CCFE, Culham Science Centre
+ author: P J Knight, CCFE, Culham Science Centre
+ zlocal : input real : local plasma effective charge
+ arg : input real : argument of Legendre function
+ palpha : output real : value of Legendre function
+ palphap : output real : derivative of Legendre function
+ This routine calculates the Legendre function palpha
+ of argument arg and order
+ alpha = -0.5 + i sqrt(xisq),
+ and its derivative palphap.
+
This Legendre function is a conical function and we use the series
+ in xisq given in Abramowitz and Stegun. The
+ derivative is calculated from the derivative of this series.
+
The derivatives were checked by calculating palpha for
+ neighbouring arguments. The calculation of palpha for zero
+ argument was checked by comparison with the expression
+ palpha(0) = 1/sqrt(pi) * cos(pi*alpha/2) * gam1 / gam2
+ (Abramowitz and Stegun, eqn 8.6.1). Here gam1 and
+ gam2 are the Gamma functions of arguments
+ 0.5*(1+alpha) and 0.5*(2+alpha) respectively.
+ Abramowitz and Stegun, equation 8.12.1
+ """
+ if abs(arg) > (1.0e0 + 1.0e-10):
+ eh.fdiags[0] = arg
+ raise ProcessValueError("Invalid argument", arg=arg)
+
+ arg2 = min(arg, (1.0e0 - 1.0e-10))
+ sinsq = 0.5e0 * (1.0e0 - arg2)
+ xisq = 0.25e0 * (32.0e0 * zlocal / (zlocal + 1.0e0) - 1.0e0)
+ palpha = 1.0e0
+ pold = 1.0e0
+ pterm = 1.0e0
+ palphap = 0.0e0
+ poldp = 0.0e0
+
+ for n in range(10000):
+ # Check for convergence every 20 iterations
+
+ if (n > 1) and ((n % 20) == 1):
+ term1 = 1.0e-10 * max(abs(pold), abs(palpha))
+ term2 = 1.0e-10 * max(abs(poldp), abs(palphap))
- if (current_drive_variables.i_hcd_secondary == 5) or (
- current_drive_variables.i_hcd_secondary == 8
- ):
- po.oblnkl(self.outfile)
- po.ocmmnt(self.outfile, "Secondary fixed neutral beam power balance :")
- po.ocmmnt(self.outfile, "----------------------------")
- po.ovarrf(
- self.outfile,
- "Secondary fixed beam first orbit loss power (MW)",
- "(porbitlossmwfix)",
- porbitlossmwfix,
- "OP ",
- )
- po.ovarrf(
- self.outfile,
- "Secondary fixed beam shine-through power [MW]",
- "(nbshinemwfix)",
- nbshinemwfix,
- "OP ",
- )
- po.ovarrf(
- self.outfile,
- "Secondary fixed beam power deposited in plasma (MW)",
- "(pinjmwfix)",
- pinjmwfix,
- "OP ",
- )
- po.ovarrf(
- self.outfile,
- "Maximum allowable beam power (MW)",
- "(p_hcd_injected_max)",
- current_drive_variables.p_hcd_injected_max,
- )
- po.ovarrf(
- self.outfile,
- "Secondary fixed total (MW)",
- "(porbitlossmwfixed+nbshinemwfix+pinjmwfix)",
- porbitlossmwfix + nbshinemwfix + pinjmwfix,
- )
- po.oblnkl(self.outfile)
- po.ovarrf(
- self.outfile,
- "Secondary beam power entering vacuum vessel (MW)",
- "(pnbitotfix)",
- pnbitotfix,
- "OP ",
- )
+ if (abs(pold - palpha) < term1) and (abs(poldp - palphap) < term2):
+ return palpha, palphap
- po.oblnkl(self.outfile)
+ pold = palpha
+ poldp = palphap
- po.ovarre(
- self.outfile,
- "Fraction of beam energy to ions",
- "(f_p_beam_injected_ions)",
- current_drive_variables.f_p_beam_injected_ions,
- "OP ",
- )
- po.ovarre(
- self.outfile,
- "Beam duct shielding thickness (m)",
- "(dx_beam_shield)",
- current_drive_variables.dx_beam_shield,
- )
- po.ovarre(
- self.outfile,
- "Beam tangency radius / Plasma major radius",
- "(frbeam)",
- current_drive_variables.frbeam,
- )
- po.ovarre(
- self.outfile,
- "Beam centreline tangency radius (m)",
- "(rtanbeam)",
- current_drive_variables.rtanbeam,
- "OP ",
- )
- po.ovarre(
- self.outfile,
- "Maximum possible tangency radius (m)",
- "(rtanmax)",
- current_drive_variables.rtanmax,
- "OP ",
+ pterm = (
+ pterm
+ * (4.0e0 * xisq + (2.0e0 * n - 1.0e0) ** 2)
+ / (2.0e0 * n) ** 2
+ * sinsq
)
+ palpha = palpha + pterm
+ palphap = palphap - n * pterm / (1.0e0 - arg2)
+ else:
+ raise ProcessError("legend: Solution has not converged")
- if abs(current_drive_variables.p_ecrh_injected_mw) > 1.0e-8:
- po.ovarre(
- self.outfile,
- "Electron cyclotron injected power (MW)",
- "(p_ecrh_injected_mw)",
- current_drive_variables.p_ecrh_injected_mw,
- "OP ",
- )
- po.ovarrf(
- self.outfile,
- "Maximum allowable ECRH power (MW)",
- "(p_hcd_injected_max)",
- current_drive_variables.p_hcd_injected_max,
- )
- po.ovarre(
- self.outfile,
- "ECH wall plug efficiency",
- "(eta_ecrh_injector_wall_plug)",
- current_drive_variables.eta_ecrh_injector_wall_plug,
- )
- po.ovarre(
- self.outfile,
- "ECH wall plug power (MW)",
- "(echwpow)",
- current_drive_variables.echwpow,
- "OP ",
- )
- if abs(current_drive_variables.p_hcd_secondary_injected_mw) > 1.0e-8:
- po.ovarrf(
- self.outfile,
- "Fixed ECRH power (MW)",
- "(pinjmwfix)",
- current_drive_variables.pinjmwfix,
- )
- po.ovarre(
- self.outfile,
- "ECH wall plug efficiency",
- "(eta_ecrh_injector_wall_plug)",
- current_drive_variables.eta_ecrh_injector_wall_plug,
- )
- po.ovarre(
- self.outfile,
- "Secondary fixed ECH wall plug power (MW)",
- "(pinjwpfix)",
- current_drive_variables.pinjwpfix,
- "OP ",
- )
+class IonCyclotron:
+ def __init__(self, plasma_profile: PlasmaProfile):
+ self.outfile = constants.nout
+ self.plasma_profile = plasma_profile
- def iternb(self):
- """Routine to calculate ITER Neutral Beam current drive parameters
- author: P J Knight, CCFE, Culham Science Centre
- effnbss : output real : neutral beam current drive efficiency (A/W)
- f_p_beam_injected_ions : output real : fraction of NB power given to ions
- fshine : output real : shine-through fraction of beam
- This routine calculates the current drive parameters for a
- neutral beam system, based on the 1990 ITER model.
- ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al,
- ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990
+ def ion_cyclotron_ipdg89(
+ self, ten: float, zeff: float, rmajor: float, dene20: float
+ ) -> float:
"""
- # Check argument sanity
- if (1 + physics_variables.eps) < current_drive_variables.frbeam:
- raise ProcessValueError(
- "Imminent negative square root argument; NBI will miss plasma completely",
- eps=physics_variables.eps,
- frbeam=current_drive_variables.frbeam,
- )
+ Routine to calculate IPDG89 Ion Cyclotron heating efficiency.
- # Calculate beam path length to centre
- dpath = physics_variables.rmajor * np.sqrt(
- (1.0 + physics_variables.eps) ** 2 - current_drive_variables.frbeam**2
- )
+ This function computes the ion cyclotron heating efficiency based on
+ the electron temperature, effective charge, major radius, and electron density.
- # Calculate beam stopping cross-section
- sigstop = self.sigbeam(
- current_drive_variables.e_beam_kev / physics_variables.m_beam_amu,
- physics_variables.te,
- physics_variables.dene,
- physics_variables.f_nd_alpha_electron,
- physics_variables.rncne,
- physics_variables.rnone,
- physics_variables.rnfene,
- )
+ :param ten: Density weighted average electron temperature keV.
+ :type ten: float
+ :param zeff: Plasma effective charge.
+ :type zeff: float
+ :param rmajor: Major radius of the plasma in meters.
+ :type rmajor: float
+ :param dene: Volume averaged electron density in 1x10^20 m^-3.
+ :type dene: float
- # Calculate number of decay lengths to centre
- current_drive_variables.n_beam_decay_lengths_core = (
- dpath * physics_variables.dene * sigstop
- )
+ :return: The calculated ion cyclotron heating efficiency in A/W.
+ :rtype: float
- # Shine-through fraction of beam
- fshine = np.exp(-2.0 * dpath * physics_variables.dene * sigstop)
- fshine = max(fshine, 1e-20)
+ :notes:
+ - The 0.1 term is to convert the temperature into 10 keV units
+ - The original formula is for the normalised current drive efficnecy
+ hence the addition of the density and majro radius terms to get back to an absolute value
- # Deuterium and tritium beam densities
- dend = physics_variables.nd_fuel_ions * (
- 1.0 - current_drive_variables.f_beam_tritium
- )
- dent = physics_variables.nd_fuel_ions * current_drive_variables.f_beam_tritium
+ :references:
+ - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989',
+ https://inis.iaea.org/collection/NCLCollectionStore/_Public/21/068/21068960.pdf
- # Power split to ions / electrons
- f_p_beam_injected_ions = self.cfnbi(
- physics_variables.m_beam_amu,
- current_drive_variables.e_beam_kev,
- physics_variables.ten,
- physics_variables.dene,
- dend,
- dent,
- physics_variables.zeffai,
- physics_variables.dlamie,
+ - T.C. Hender et al., 'Physics Assessment of the European Reactor Study', AEA FUS 172, 1992.
+ """
+
+ return ((0.63e0 * 0.1e0 * ten) / (2.0e0 + zeff)) / (rmajor * dene20)
+
+
+class ElectronBernstein:
+ def __init__(self, plasma_profile: PlasmaProfile):
+ self.outfile = constants.nout
+ self.plasma_profile = plasma_profile
+
+ def electron_bernstein_freethy(
+ self,
+ te: float,
+ rmajor: float,
+ dene20: float,
+ bt: float,
+ n_ecrh_harmonic: int,
+ xi_ebw: float,
+ ) -> float:
+ """
+ Calculate the Electron Bernstein Wave (EBW) current drive efficiency using the Freethy model.
+
+ This function computes the EBW current drive efficiency based on the electron temperature,
+ major radius, electron density, magnetic field, harmonic number, and scaling factor.
+
+ :param te: Volume averaged electron temperature in keV.
+ :type te: float
+ :param rmajor: Major radius of the plasma in meters.
+ :type rmajor: float
+ :param dene20: Volume averaged electron density in units of 10^20 m^-3.
+ :type dene20: float
+ :param bt: Toroidal magnetic field in Tesla.
+ :type bt: float
+ :param n_ecrh_harmonic: Cyclotron harmonic number (fundamental used as default).
+ :type n_ecrh_harmonic: int
+ :param xi_ebw: Scaling factor for EBW efficiency.
+ :type xi_ebw: float
+
+ :return: The calculated absolute EBW current drive efficiency in A/W.
+ :rtype: float
+
+ :notes:
+ - EBWs can only couple to plasma if the cyclotron harmonic is above the plasma density cut-off.
+ - The density factor accounts for this behavior.
+
+ :references:
+ - Freethy, S., PROCESS issue #1262.
+ """
+
+ # Normalised current drive efficiency gamma
+ eta_cd_norm = (xi_ebw / 32.7e0) * te
+
+ # Absolute current drive efficiency
+ eta_cd = eta_cd_norm / (dene20 * rmajor)
+
+ # EBWs can only couple to plasma if cyclotron harmonic is above plasma density cut-off;
+ # this behavior is captured in the following function:
+ # constant 'a' controls sharpness of transition
+ a = 0.1e0
+
+ fc = (
+ 1.0e0
+ / (2.0e0 * np.pi)
+ * n_ecrh_harmonic
+ * constants.electron_charge
+ * bt
+ / constants.electron_mass
)
- # Current drive efficiency
- effnbss = current_drive_variables.frbeam * self.etanb(
- physics_variables.m_beam_amu,
- physics_variables.alphan,
- physics_variables.alphat,
- physics_variables.aspect,
- physics_variables.dene,
- current_drive_variables.e_beam_kev,
- physics_variables.rmajor,
- physics_variables.ten,
- physics_variables.zeff,
+ fp = (
+ 1.0e0
+ / (2.0e0 * np.pi)
+ * np.sqrt(
+ dene20
+ * constants.electron_charge**2
+ / (constants.electron_mass * constants.epsilon0)
+ )
)
- return effnbss, f_p_beam_injected_ions, fshine
+ density_factor = 0.5e0 * (1.0e0 + np.tanh((2.0e0 / a) * ((fp - fc) / fp - a)))
+
+ return eta_cd * density_factor
+
+
+class LowerHybrid:
+ def __init__(self, plasma_profile: PlasmaProfile):
+ self.outfile = constants.nout
+ self.plasma_profile = plasma_profile
def cullhy(self):
"""Routine to calculate Lower Hybrid current drive efficiency
@@ -1468,7 +1011,7 @@ def cullhy(self):
if term03 > term04:
raise ProcessValueError(
- "Normalised LH efficiency < 0; use a different value of iefrf",
+ "Normalised LH efficiency < 0; use a different value of i_hcd_primary",
term03=term03,
term04=term04,
)
@@ -1479,20 +1022,89 @@ def cullhy(self):
return gamlh / ((0.1e0 * dlocal) * physics_variables.rmajor)
- def culecd(self):
- """Routine to calculate Electron Cyclotron current drive efficiency
- author: M R O'Brien, CCFE, Culham Science Centre
+ def lhrad(self):
+ """Routine to calculate Lower Hybrid wave absorption radius
author: P J Knight, CCFE, Culham Science Centre
- effrfss : output real : electron cyclotron current drive efficiency (A/W)
- This routine calculates the current drive parameters for a
- electron cyclotron system, based on the AEA FUS 172 model.
+ rratio : output real : minor radius of penetration / rminor
+ This routine determines numerically the minor radius at which the
+ damping of Lower Hybrid waves occurs, using a Newton-Raphson method.
AEA FUS 172: Physics Assessment for the European Reactor Study
"""
- rrr = 1.0e0 / 3.0e0
+ # Correction to refractive index (kept within valid bounds)
+ drfind = min(0.7e0, max(0.1e0, 12.5e0 / physics_variables.te0))
+
+ # Use Newton-Raphson method to establish the correct minor radius
+ # ratio. g is calculated as a function of r / r_minor, where g is
+ # the difference between the results of the two formulae for the
+ # energy E given in AEA FUS 172, p.58. The required minor radius
+ # ratio has been found when g is sufficiently close to zero.
+
+ # Initial guess for the minor radius ratio
+
+ rat0 = 0.8e0
+
+ for _ in range(100):
+ # Minor radius ratios either side of the latest guess
+
+ r1 = rat0 - 1.0e-3 * rat0
+ r2 = rat0 + 1.0e-3 * rat0
+
+ # Evaluate g at rat0, r1, r2
+
+ g0 = self.lheval(drfind, rat0)
+ g1 = self.lheval(drfind, r1)
+ g2 = self.lheval(drfind, r2)
+
+ # Calculate gradient of g with respect to minor radius ratio
+
+ dgdr = (g2 - g1) / (r2 - r1)
+
+ # New approximation
+
+ rat1 = rat0 - g0 / dgdr
+
+ # Force this approximation to lie within bounds
+
+ rat1 = max(0.0001e0, rat1)
+ rat1 = min(0.9999e0, rat1)
+
+ if abs(g0) <= 0.01e0:
+ break
+ rat0 = rat1
+
+ else:
+ eh.report_error(16)
+ rat0 = 0.8e0
+
+ return rat0
+
+ def lheval(self, drfind, rratio):
+ """Routine to evaluate the difference between electron energy
+ expressions required to find the Lower Hybrid absorption radius
+ author: P J Knight, CCFE, Culham Science Centre
+ drfind : input real : correction to parallel refractive index
+ rratio : input real : guess for radius of penetration / rminor
+ ediff : output real : difference between the E values (keV)
+ This routine evaluates the difference between the values calculated
+ from the two equations for the electron energy E, given in
+ AEA FUS 172, p.58. This difference is used to locate the Lower Hybrid
+ wave absorption radius via a Newton-Raphson method, in calling
+ routine lhrad.
+ AEA FUS 172: Physics Assessment for the European Reactor Study
+ """
+ dlocal = 1.0e-19 * self.plasma_profile.neprofile.calculate_profile_y(
+ rratio,
+ physics_variables.rhopedn,
+ physics_variables.ne0,
+ physics_variables.neped,
+ physics_variables.nesep,
+ physics_variables.alphan,
+ )
+
+ # Local electron temperature
- # Temperature
tlocal = self.plasma_profile.teprofile.calculate_profile_y(
- rrr,
+ rratio,
physics_variables.rhopedt,
physics_variables.te0,
physics_variables.teped,
@@ -1501,714 +1113,1327 @@ def culecd(self):
physics_variables.tbeta,
)
- # Density (10**20 m**-3)
- dlocal = 1.0e-20 * self.plasma_profile.neprofile.calculate_profile_y(
- rrr,
- physics_variables.rhopedn,
- physics_variables.ne0,
- physics_variables.neped,
- physics_variables.nesep,
- physics_variables.alphan,
+ # Local toroidal field (evaluated at the inboard region of the flux surface)
+
+ blocal = (
+ physics_variables.bt
+ * physics_variables.rmajor
+ / (physics_variables.rmajor - rratio * physics_variables.rminor)
)
- # Inverse aspect ratio
- epsloc = rrr * physics_variables.rminor / physics_variables.rmajor
+ # Parallel refractive index needed for plasma access
- # Effective charge (use average value)
- zlocal = physics_variables.zeff
+ frac = np.sqrt(dlocal) / blocal
+ nplacc = frac + np.sqrt(1.0e0 + frac * frac)
- # Coulomb logarithm for ion-electron collisions
- # (From J. A. Wesson, 'Tokamaks', Clarendon Press, Oxford, p.293)
- coulog = 15.2e0 - 0.5e0 * np.log(dlocal) + np.log(tlocal)
+ # Total parallel refractive index
- # 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.
- 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)
+ refind = nplacc + drfind
- # Normalised current drive efficiency (A/W m**-2)
- ecgam = 0.25e0 * (ecgam1 + ecgam2 + ecgam3 + ecgam4)
+ # First equation for electron energy E
- # Current drive efficiency (A/W)
- return ecgam / (dlocal * physics_variables.rmajor)
+ e1 = 511.0e0 * (np.sqrt(1.0e0 + 1.0e0 / (refind * refind)) - 1.0e0)
- def eccdef(self, tlocal, epsloc, zlocal, cosang, coulog):
- """Routine to calculate Electron Cyclotron current drive efficiency
- author: M R O'Brien, CCFE, Culham Science Centre
- author: P J Knight, CCFE, Culham Science Centre
- tlocal : input real : local electron temperature (keV)
- epsloc : input real : local inverse aspect ratio
- zlocal : input real : local plasma effective charge
- cosang : input real : cosine of the poloidal angle at which ECCD takes
- place (+1 outside, -1 inside)
- coulog : input real : local coulomb logarithm for ion-electron collisions
- ecgam : output real : normalised current drive efficiency (A/W m**-2)
- 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 nIR/P with n the - local density in 10**20 /m**3, I the driven current in MAmps, - R the major radius in metres, and P the absorbed power in MWatts. - AEA FUS 172: Physics Assessment for the European Reactor Study - ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, - ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 - """ - mcsq = ( - constants.electron_mass * 2.9979e8**2 / (1.0e3 * constants.electron_volt) - ) # keV - f = 16.0e0 * (tlocal / mcsq) ** 2 + # Second equation for E - # fp is the derivative of f with respect to gamma, the relativistic - # factor, taken equal to 1 + 2T/(m c**2) + e2 = 7.0e0 * tlocal - fp = 16.0e0 * tlocal / mcsq + # Difference - # 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) = -8/(1+zlocal). alpha is of the - # form (-1/2 + ix), with x a real number and i = sqrt(-1). + return e1 - e2 - lam = 1.0e0 - palpha, palphap = self.legend(zlocal, lam) + def lower_hybrid_fenstermacher( + self, te: float, rmajor: float, dene20: float + ) -> float: + """ + Calculate the lower hybrid frequency using the Fenstermacher formula. + This function computes the lower hybrid frequency based on the electron + temperature, major radius, and electron density. - lams = np.sqrt(2.0e0 * epsloc / (1.0e0 + epsloc)) - palphas, _ = self.legend(zlocal, lams) + :param te: Volume averaged electron temperature in keV. + :type te: float + :param rmajor: Major radius of the plasma in meters. + :type rmajor: float + :param dene20: Volume averaged electron density in units of 10^20 m^-3. + :type dene20: float - # hp is the derivative of IPDG89's h function with respect to lam + :return: The calculated absolute current drive efficiency in A/W. + :rtype: float - h = -4.0e0 * lam / (zlocal + 5.0e0) * (1.0e0 - lams * palpha / (lam * palphas)) - hp = -4.0e0 / (zlocal + 5.0e0) * (1.0e0 - lams * palphap / palphas) + :notes: + - This forumla was originally in the Oak RidgeSystems Code, attributed to Fenstermacher + and is used in the AEA FUS 172 report. - # facm is IPDG89's momentum conserving factor + :references: + - T.C. Hender et al., 'Physics Assessment of the European Reactor Study', AEA FUS 172, 1992. - facm = 1.5e0 - y = mcsq / (2.0e0 * tlocal) * (1.0e0 + epsloc * cosang) + - R.L.Reid et al, Oak Ridge Report ORNL/FEDC-87-7, 1988 + """ - # We take the negative of the IPDG89 expression to get a positive - # number + return (0.36e0 * (1.0e0 + (te / 25.0e0) ** 1.16e0)) / (rmajor * dene20) - ecgam = ( - -7.8e0 - * facm - * np.sqrt((1.0e0 + epsloc) / (1.0e0 - epsloc)) - / coulog - * (h * fp - 0.5e0 * y * f * hp) + def lower_hybrid_ehst( + self, te: float, beta: float, rmajor: float, dene20: float, zeff: float + ) -> float: + """ + Calculate the Lower Hybrid current drive efficiency using the Ehst model. + + This function computes the current drive efficiency based on the electron + temperature, beta, major radius, electron density, and effective charge. + + :param te: Volume averaged electron temperature in keV. + :type te: float + :param beta: Plasma beta value (ratio of plasma pressure to magnetic pressure). + :type beta: float + :param rmajor: Major radius of the plasma in meters. + :type rmajor: float + :param dene20: Volume averaged electron density in units of 10^20 m^-3. + :type dene20: float + :param zeff: Plasma effective charge. + :type zeff: float + + :return: The calculated absolute current drive efficiency in A/W. + :rtype: float + + :notes: + + :references: + - Ehst, D.A., and Karney, C.F.F., "Lower Hybrid Current Drive in Tokamaks", + Nuclear Fusion, 31(10), 1933-1949, 1991. + """ + return ( + ((te**0.77 * (0.034 + 0.196 * beta)) / (rmajor * dene20)) + * ( + 32.0 / (5.0 + zeff) + + 2.0 + + (12.0 * (6.0 + zeff)) / (5.0 + zeff) / (3.0 + zeff) + + 3.76 / zeff + ) + / 12.507 ) - if ecgam < 0.0e0: - raise ProcessValueError("Negative normalised current drive efficiency") - return ecgam - def culnbi(self): - """Routine to calculate Neutral Beam current drive parameters - author: P J Knight, CCFE, Culham Science Centre - effnbss : output real : neutral beam current drive efficiency (A/W) - f_p_beam_injected_ions : output real : fraction of NB power given to ions - fshine : output real : shine-through fraction of beam - 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). - AEA FUS 172: Physics Assessment for the European Reactor Study +class CurrentDrive: + def __init__( + self, + plasma_profile: PlasmaProfile, + electron_cyclotron: ElectronCyclotron, + ion_cyclotron: IonCyclotron, + lower_hybrid: LowerHybrid, + neutral_beam: NeutralBeam, + electron_bernstein: ElectronBernstein, + ): + self.outfile = constants.nout + self.plasma_profile = plasma_profile + self.electron_cyclotron = electron_cyclotron + self.ion_cyclotron = ion_cyclotron + self.lower_hybrid = lower_hybrid + self.neutral_beam = neutral_beam + self.electron_bernstein = electron_bernstein + + def cudriv(self) -> None: """ - if (1.0e0 + physics_variables.eps) < current_drive_variables.frbeam: - raise ProcessValueError( - "Imminent negative square root argument; NBI will miss plasma completely", - eps=physics_variables.eps, - frbeam=current_drive_variables.frbeam, + Calculate the current drive power requirements. + + This method computes the power requirements of the current drive system + using a choice of models for the current drive efficiency. + + :param output: Flag indicating whether to write results to the output file. + :type output: bool + + :raises ProcessValueError: If an invalid current drive switch is encountered. + """ + + current_drive_variables.p_hcd_ecrh_injected_total_mw = 0.0e0 + current_drive_variables.p_hcd_beam_injected_total_mw = 0.0e0 + current_drive_variables.p_hcd_lowhyb_injected_total_mw = 0.0e0 + current_drive_variables.p_hcd_icrh_injected_total_mw = 0.0e0 + current_drive_variables.p_hcd_ebw_injected_total_mw = 0.0e0 + current_drive_variables.c_beam_total = 0.0e0 + current_drive_variables.p_beam_orbit_loss_mw = 0.0e0 + + pinjmw1 = 0.0 + p_hcd_primary_ions_mw = 0.0 + p_hcd_primary_electrons_mw = 0.0 + p_hcd_secondary_electrons_mw = 0.0 + p_hcd_secondary_ions_mw = 0.0 + + # To stop issues with input file we force + # zero secondary heating if no injection method + if current_drive_variables.i_hcd_secondary == 0: + current_drive_variables.p_hcd_secondary_extra_heat_mw = 0.0 + + # i_hcd_calculations | switch for current drive calculation + # = 0 | turned off + # = 1 | turned on + + if current_drive_variables.i_hcd_calculations != 0: + # Put electron density in desired units (10^-20 m-3) + dene20 = physics_variables.dene * 1.0e-20 + + # Calculate current drive efficiencies + # ============================================================== + + # Define a dictionary of lambda functions for current drive efficiency models + hcd_models = { + 1: lambda: self.lower_hybrid.lower_hybrid_fenstermacher( + physics_variables.te, physics_variables.rmajor, dene20 + ) + * current_drive_variables.feffcd, + 2: lambda: self.ion_cyclotron.ion_cyclotron_ipdg89( + ten=physics_variables.ten, + zeff=physics_variables.zeff, + rmajor=physics_variables.rmajor, + dene20=dene20, + ) + * current_drive_variables.feffcd, + 3: lambda: self.electron_cyclotron.electron_cyclotron_fenstermacher( + ten=physics_variables.ten, + rmajor=physics_variables.rmajor, + dene20=dene20, + dlamee=physics_variables.dlamee, + ) + * current_drive_variables.feffcd, + 4: lambda: self.lower_hybrid.lower_hybrid_ehst( + te=physics_variables.te, + beta=physics_variables.beta, + rmajor=physics_variables.rmajor, + dene20=dene20, + zeff=physics_variables.zeff, + ) + * current_drive_variables.feffcd, + 5: lambda: ( + self.neutral_beam.iternb()[0] * current_drive_variables.feffcd + ), + 6: lambda: self.lower_hybrid.cullhy() * current_drive_variables.feffcd, + 7: lambda: self.electron_cyclotron.culecd() + * current_drive_variables.feffcd, + 8: lambda: ( + self.neutral_beam.culnbi()[0] * current_drive_variables.feffcd + ), + 10: lambda: current_drive_variables.eta_cd_norm_ecrh + / (dene20 * physics_variables.rmajor), + 12: lambda: self.electron_bernstein.electron_bernstein_freethy( + te=physics_variables.te, + rmajor=physics_variables.rmajor, + dene20=dene20, + bt=physics_variables.bt, + n_ecrh_harmonic=current_drive_variables.n_ecrh_harmonic, + xi_ebw=current_drive_variables.xi_ebw, + ) + * current_drive_variables.feffcd, + 13: lambda: self.electron_cyclotron.electron_cyclotron_freethy( + te=physics_variables.te, + zeff=physics_variables.zeff, + rmajor=physics_variables.rmajor, + dene=physics_variables.dene, + bt=physics_variables.bt, + n_ecrh_harmonic=current_drive_variables.n_ecrh_harmonic, + i_ecrh_wave_mode=current_drive_variables.i_ecrh_wave_mode, + ) + * current_drive_variables.feffcd, + } + + # Assign outputs for models that return multiple values + if current_drive_variables.i_hcd_secondary in [5, 8]: + _, f_p_beam_injected_ions, f_p_beam_shine_through = ( + self.neutral_beam.iternb() + if current_drive_variables.i_hcd_secondary == 5 + else self.neutral_beam.culnbi() + ) + current_drive_variables.f_p_beam_injected_ions = f_p_beam_injected_ions + current_drive_variables.f_p_beam_shine_through = f_p_beam_shine_through + + # Calculate eta_cd_hcd_secondary based on the selected model + if current_drive_variables.i_hcd_secondary.item() in hcd_models: + current_drive_variables.eta_cd_hcd_secondary = hcd_models[ + current_drive_variables.i_hcd_secondary.item() + ]() + elif current_drive_variables.i_hcd_secondary != 0: + raise ProcessValueError( + f"Current drive switch is invalid: {current_drive_variables.i_hcd_secondary = }" + ) + + # Calculate eta_cd_hcd_primary based on the selected model + if current_drive_variables.i_hcd_primary.item() in hcd_models: + current_drive_variables.eta_cd_hcd_primary = hcd_models[ + current_drive_variables.i_hcd_primary.item() + ]() + else: + raise ProcessValueError( + f"Current drive switch is invalid: {current_drive_variables.i_hcd_primary = }" + ) + + # Calculate the normalised current drive efficieny for the primary heating method + current_drive_variables.eta_cd_norm_hcd_primary = ( + current_drive_variables.eta_cd_hcd_primary + * (dene20 * physics_variables.rmajor) + ) + # Calculate the normalised current drive efficieny for the secondary heating method + current_drive_variables.eta_cd_norm_hcd_secondary = ( + current_drive_variables.eta_cd_hcd_secondary + * (dene20 * physics_variables.rmajor) ) - # Calculate beam path length to centre + # Calculate the driven current for the secondary heating method + current_drive_variables.c_hcd_secondary_driven = ( + current_drive_variables.eta_cd_hcd_secondary + * current_drive_variables.p_hcd_secondary_injected_mw + * 1.0e6 + ) + # Calculate the fraction of the plasma current driven by the secondary heating method + current_drive_variables.f_c_plasma_hcd_secondary = ( + current_drive_variables.c_hcd_secondary_driven + / physics_variables.plasma_current + ) - dpath = physics_variables.rmajor * np.sqrt( - (1.0e0 + physics_variables.eps) ** 2 - current_drive_variables.frbeam**2 - ) + # Calculate the injected power for the primary heating method + current_drive_variables.p_hcd_primary_injected_mw = ( + 1.0e-6 + * ( + physics_variables.f_c_plasma_auxiliary + - current_drive_variables.f_c_plasma_hcd_secondary + ) + * physics_variables.plasma_current + / current_drive_variables.eta_cd_hcd_primary + ) - # Calculate beam stopping cross-section + # Calculate the driven current for the primary heating method + current_drive_variables.c_hcd_primary_driven = ( + current_drive_variables.eta_cd_hcd_primary + * current_drive_variables.p_hcd_primary_injected_mw + * 1.0e6 + ) + # Calculate the fraction of the plasma current driven by the primary heating method + current_drive_variables.f_c_plasma_hcd_primary = ( + current_drive_variables.c_hcd_primary_driven + / physics_variables.plasma_current + ) - sigstop = self.sigbeam( - current_drive_variables.e_beam_kev / physics_variables.m_beam_amu, - physics_variables.te, - physics_variables.dene, - physics_variables.f_nd_alpha_electron, - physics_variables.rncne, - physics_variables.rnone, - physics_variables.rnfene, - ) + # =========================================================== - # Calculate number of decay lengths to centre + # Calculate the wall plug power for the secondary heating method + # ============================================================== - current_drive_variables.n_beam_decay_lengths_core = ( - dpath * physics_variables.dnla * sigstop - ) + # Lower hybrid cases + if current_drive_variables.i_hcd_secondary in [1, 4, 6]: + # Injected power + p_hcd_secondary_electrons_mw = ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) - # Shine-through fraction of beam + # Wall plug power + heat_transport_variables.p_hcd_secondary_electric_mw = ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) / current_drive_variables.eta_lowhyb_injector_wall_plug - fshine = np.exp(-2.0e0 * dpath * physics_variables.dnla * sigstop) - fshine = max(fshine, 1.0e-20) + # Wall plug to injector efficiency + current_drive_variables.eta_hcd_secondary_injector_wall_plug = ( + current_drive_variables.eta_lowhyb_injector_wall_plug + ) - # Deuterium and tritium beam densities + current_drive_variables.p_hcd_lowhyb_injected_total_mw += ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) - dend = physics_variables.nd_fuel_ions * ( - 1.0e0 - current_drive_variables.f_beam_tritium - ) - dent = physics_variables.nd_fuel_ions * current_drive_variables.f_beam_tritium + # ========================================================== - # Power split to ions / electrons + # Ion cyclotron cases + if current_drive_variables.i_hcd_secondary in [2]: + # Injected power + p_hcd_secondary_ions_mw = ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) - f_p_beam_injected_ions = self.cfnbi( - physics_variables.m_beam_amu, - current_drive_variables.e_beam_kev, - physics_variables.ten, - physics_variables.dene, - dend, - dent, - physics_variables.zeffai, - physics_variables.dlamie, - ) + # Wall plug power + heat_transport_variables.p_hcd_secondary_electric_mw = ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) / current_drive_variables.eta_icrh_injector_wall_plug + + # Wall plug to injector efficiency + current_drive_variables.eta_hcd_secondary_injector_wall_plug = ( + current_drive_variables.eta_icrh_injector_wall_plug + ) + + current_drive_variables.p_hcd_icrh_injected_total_mw += ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) + + # ========================================================== + + # Electron cyclotron cases + if current_drive_variables.i_hcd_secondary in [3, 7, 10, 13]: + # Injected power + p_hcd_secondary_electrons_mw = ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) + + # Wall plug power + heat_transport_variables.p_hcd_secondary_electric_mw = ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) / current_drive_variables.eta_ecrh_injector_wall_plug + + # Wall plug to injector efficiency + current_drive_variables.eta_hcd_secondary_injector_wall_plug = ( + current_drive_variables.eta_ecrh_injector_wall_plug + ) + + current_drive_variables.p_hcd_ecrh_injected_total_mw += ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) + + # ========================================================== + + # Electron berstein cases + if current_drive_variables.i_hcd_secondary in [12]: + # Injected power + p_hcd_secondary_electrons_mw = ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) + + # Wall plug power + heat_transport_variables.p_hcd_secondary_electric_mw = ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) / current_drive_variables.eta_ebw_injector_wall_plug - # Current drive efficiency + # Wall plug to injector efficiency + current_drive_variables.eta_hcd_secondary_injector_wall_plug = ( + current_drive_variables.eta_ebw_injector_wall_plug + ) - effnbss = self.etanb2( - physics_variables.m_beam_amu, - physics_variables.alphan, - physics_variables.alphat, - physics_variables.aspect, - physics_variables.dene, - physics_variables.dnla, - current_drive_variables.e_beam_kev, - current_drive_variables.frbeam, - fshine, - physics_variables.rmajor, - physics_variables.rminor, - physics_variables.ten, - physics_variables.zeff, - ) + current_drive_variables.p_hcd_ebw_injected_total_mw += ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) - return effnbss, f_p_beam_injected_ions, fshine + # ========================================================== - def lhrad(self): - """Routine to calculate Lower Hybrid wave absorption radius - author: P J Knight, CCFE, Culham Science Centre - rratio : output real : minor radius of penetration / rminor - This routine determines numerically the minor radius at which the - damping of Lower Hybrid waves occurs, using a Newton-Raphson method. - AEA FUS 172: Physics Assessment for the European Reactor Study - """ - # Correction to refractive index (kept within valid bounds) - drfind = min(0.7e0, max(0.1e0, 12.5e0 / physics_variables.te0)) + # Neutral beam cases + elif current_drive_variables.i_hcd_secondary in [5, 8]: + # Account for first orbit losses + # (power due to particles that are ionised but not thermalised) [MW]: + # This includes a second order term in shinethrough*(first orbit loss) - # Use Newton-Raphson method to establish the correct minor radius - # ratio. g is calculated as a function of r / r_minor, where g is - # the difference between the results of the two formulae for the - # energy E given in AEA FUS 172, p.58. The required minor radius - # ratio has been found when g is sufficiently close to zero. + current_drive_variables.f_p_beam_orbit_loss = min( + 0.999, current_drive_variables.f_p_beam_orbit_loss + ) # Should never be needed - # Initial guess for the minor radius ratio + # Shinethrough power (atoms that are not ionised) [MW]: + current_drive_variables.p_beam_shine_through_mw = ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + * (1.0 - current_drive_variables.f_p_beam_shine_through) + ) - rat0 = 0.8e0 + # First orbit loss + current_drive_variables.p_beam_orbit_loss_mw = ( + current_drive_variables.f_p_beam_orbit_loss + * ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + - current_drive_variables.p_beam_shine_through_mw + ) + ) - for _ in range(100): - # Minor radius ratios either side of the latest guess + # Power deposited + current_drive_variables.p_beam_plasma_coupled_mw = ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + - current_drive_variables.p_beam_shine_through_mw + - current_drive_variables.p_beam_orbit_loss_mw + ) - r1 = rat0 - 1.0e-3 * rat0 - r2 = rat0 + 1.0e-3 * rat0 + p_hcd_secondary_ions_mw = ( + current_drive_variables.p_beam_plasma_coupled_mw + * current_drive_variables.f_p_beam_injected_ions + ) - # Evaluate g at rat0, r1, r2 + p_hcd_secondary_electrons_mw = ( + current_drive_variables.p_beam_plasma_coupled_mw + * (1.0e0 - current_drive_variables.f_p_beam_injected_ions) + ) - g0 = self.lheval(drfind, rat0) - g1 = self.lheval(drfind, r1) - g2 = self.lheval(drfind, r2) + current_drive_variables.pwpnb = ( + ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) + / current_drive_variables.eta_beam_injector_wall_plug + ) # neutral beam wall plug power - # Calculate gradient of g with respect to minor radius ratio + heat_transport_variables.p_hcd_secondary_electric_mw = ( + current_drive_variables.pwpnb + ) - dgdr = (g2 - g1) / (r2 - r1) + current_drive_variables.eta_hcd_secondary_injector_wall_plug = ( + current_drive_variables.eta_beam_injector_wall_plug + ) - # New approximation + current_drive_variables.c_beam_total = ( + 1.0e-3 + * ( + ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) + * 1.0e6 + ) + / current_drive_variables.e_beam_kev + ) # Neutral beam current (A) - rat1 = rat0 - g0 / dgdr + current_drive_variables.p_hcd_beam_injected_total_mw += ( + current_drive_variables.p_hcd_secondary_injected_mw + + current_drive_variables.p_hcd_secondary_extra_heat_mw + ) - # Force this approximation to lie within bounds + # ========================================================== - rat1 = max(0.0001e0, rat1) - rat1 = min(0.9999e0, rat1) + # Lower hybrid cases + if current_drive_variables.i_hcd_primary in [1, 4, 6]: + p_hcd_primary_electrons_mw = ( + current_drive_variables.p_hcd_primary_injected_mw + + current_drive_variables.p_hcd_primary_extra_heat_mw + ) - if abs(g0) <= 0.01e0: - break - rat0 = rat1 + current_drive_variables.p_hcd_lowhyb_injected_total_mw += ( + current_drive_variables.p_hcd_primary_injected_mw + + current_drive_variables.p_hcd_primary_extra_heat_mw + ) - else: - eh.report_error(16) - rat0 = 0.8e0 + # Wall plug power + heat_transport_variables.p_hcd_primary_electric_mw = ( + current_drive_variables.p_hcd_primary_injected_mw + + current_drive_variables.p_hcd_primary_extra_heat_mw + ) / current_drive_variables.eta_lowhyb_injector_wall_plug - return rat0 + # Wall plug to injector efficiency + current_drive_variables.eta_hcd_primary_injector_wall_plug = ( + current_drive_variables.eta_lowhyb_injector_wall_plug + ) - def etanb2( - self, - m_beam_amu, - alphan, - alphat, - aspect, - dene, - dnla, - e_beam_kev, - frbeam, - fshine, - rmajor, - rminor, - ten, - zeff, - ): - """Routine to find neutral beam current drive efficiency - using the ITER 1990 formulation, plus correction terms - outlined in Culham Report AEA FUS 172 - author: P J Knight, CCFE, Culham Science Centre - m_beam_amu : input real : beam ion mass (amu) - alphan : input real : density profile factor - alphat : input real : temperature profile factor - aspect : input real : aspect ratio - dene : input real : volume averaged electron density (m**-3) - dnla : input real : line averaged electron density (m**-3) - e_beam_kev : input real : neutral beam energy (keV) - frbeam : input real : R_tangent / R_major for neutral beam injection - fshine : input real : shine-through fraction of beam - rmajor : input real : plasma major radius (m) - rminor : input real : plasma minor radius (m) - ten : input real : density weighted average electron temperature (keV) - zeff : input real : plasma effective charge - 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. -
The formulae are from AEA FUS 172, unless denoted by IPDG89.
- AEA FUS 172: Physics Assessment for the European Reactor Study
- ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al,
- ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990
- """
- # Charge of beam ions
- zbeam = 1.0
+ # Wall plug power
+ current_drive_variables.p_hcd_lowhyb_electric_mw = (
+ current_drive_variables.p_hcd_lowhyb_injected_total_mw
+ / current_drive_variables.eta_lowhyb_injector_wall_plug
+ )
- # Fitting factor (IPDG89)
- bbd = 1.0
+ # ===========================================================
- # Volume averaged electron density (10**20 m**-3)
- dene20 = dene / 1e20
+ # Ion cyclotron cases
+ if current_drive_variables.i_hcd_primary in [2]:
+ p_hcd_primary_ions_mw = (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ )
- # Line averaged electron density (10**20 m**-3)
- dnla20 = dnla / 1e20
+ # Wall plug power
+ heat_transport_variables.p_hcd_primary_electric_mw = (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ ) / current_drive_variables.eta_icrh_injector_wall_plug
- # Critical energy (MeV) (power to electrons = power to ions) (IPDG89)
- # N.B. ten is in keV
- ecrit = 0.01 * m_beam_amu * ten
+ # Wall plug to injector efficiency
+ current_drive_variables.eta_hcd_primary_injector_wall_plug = (
+ current_drive_variables.eta_icrh_injector_wall_plug
+ )
- # Beam energy in MeV
- ebmev = e_beam_kev / 1e3
+ current_drive_variables.p_hcd_icrh_injected_total_mw += (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ )
- # x and y coefficients of function J0(x,y) (IPDG89)
- xjs = ebmev / (bbd * ecrit)
- xj = np.sqrt(xjs)
+ # Wall plug power
+ current_drive_variables.p_hcd_icrh_electric_mw = (
+ current_drive_variables.p_hcd_icrh_injected_total_mw
+ / current_drive_variables.eta_icrh_injector_wall_plug
+ )
- yj = 0.8 * zeff / m_beam_amu
+ # ===========================================================
- # Fitting function J0(x,y)
- j0 = xjs / (4.0 + 3.0 * yj + xjs * (xj + 1.39 + 0.61 * yj**0.7))
+ # Electron cyclotron cases
- # Effective inverse aspect ratio, with a limit on its maximum value
- epseff = min(0.2, (0.5 / aspect))
+ if current_drive_variables.i_hcd_primary in [3, 7, 10, 13]:
+ p_hcd_primary_electrons_mw = (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ )
- # Reduction in the reverse electron current
- # due to neoclassical effects
- gfac = (1.55 + 0.85 / zeff) * np.sqrt(epseff) - (0.2 + 1.55 / zeff) * epseff
+ # Wall plug to injector efficiency
+ heat_transport_variables.p_hcd_primary_electric_mw = (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ ) / current_drive_variables.eta_ecrh_injector_wall_plug
- # Reduction in the net beam driven current
- # due to the reverse electron current
- ffac = 1.0 - (zbeam / zeff) * (1.0 - gfac)
+ current_drive_variables.eta_hcd_primary_injector_wall_plug = (
+ current_drive_variables.eta_ecrh_injector_wall_plug
+ )
- # Normalisation to allow results to be valid for
- # non-ITER plasma size and density:
+ current_drive_variables.p_hcd_ecrh_injected_total_mw += (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ )
- # Line averaged electron density (10**20 m**-3) normalised to ITER
- nnorm = 1.0
+ # Wall plug power
+ current_drive_variables.p_hcd_ecrh_electric_mw = (
+ current_drive_variables.p_hcd_ecrh_injected_total_mw
+ / current_drive_variables.eta_ecrh_injector_wall_plug
+ )
- # Distance along beam to plasma centre
- r = max(rmajor, rmajor * frbeam)
- eps1 = rminor / r
+ # ===========================================================
- if (1.0 + eps1) < frbeam:
- raise ProcessValueError(
- "Imminent negative square root argument; NBI will miss plasma completely",
- eps=eps1,
- frbeam=frbeam,
- )
+ # Electron bernstein cases
- d = rmajor * np.sqrt((1.0 + eps1) ** 2 - frbeam**2)
+ if current_drive_variables.i_hcd_primary in [12]:
+ p_hcd_primary_electrons_mw = (
+ current_drive_variables.p_ebw_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ )
- # Distance along beam to plasma centre for ITER
- # assuming a tangency radius equal to the major radius
- epsitr = 2.15 / 6.0
- dnorm = 6.0 * np.sqrt(2.0 * epsitr + epsitr**2)
+ # Wall plug to injector efficiency
+ heat_transport_variables.p_hcd_primary_electric_mw = (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ ) / current_drive_variables.eta_ebw_injector_wall_plug
- # Normalisation to beam energy (assumes a simplified formula for
- # the beam stopping cross-section)
- ebnorm = ebmev * ((nnorm * dnorm) / (dnla20 * d)) ** (1.0 / 0.78)
+ current_drive_variables.eta_hcd_primary_injector_wall_plug = (
+ current_drive_variables.eta_ebw_injector_wall_plug
+ )
- # A_bd fitting coefficient, after normalisation with ebnorm
- abd = (
- 0.107
- * (1.0 - 0.35 * alphan + 0.14 * alphan**2)
- * (1.0 - 0.21 * alphat)
- * (1.0 - 0.2 * ebnorm + 0.09 * ebnorm**2)
- )
+ current_drive_variables.p_hcd_ebw_injected_total_mw += (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ )
+
+ # Wall plug power
+ current_drive_variables.p_hcd_ebw_electric_mw = (
+ current_drive_variables.p_ebw_injected_mw
+ / current_drive_variables.eta_ebw_injector_wall_plug
+ )
+
+ # ===========================================================
+
+ elif current_drive_variables.i_hcd_primary in [5, 8]:
+ # Account for first orbit losses
+ # (power due to particles that are ionised but not thermalised) [MW]:
+ # This includes a second order term in shinethrough*(first orbit loss)
+ current_drive_variables.f_p_beam_orbit_loss = min(
+ 0.999, current_drive_variables.f_p_beam_orbit_loss
+ ) # Should never be needed
- # Normalised current drive efficiency (A/W m**-2) (IPDG89)
- gamnb = 5.0 * abd * 0.1 * ten * (1.0 - fshine) * frbeam * j0 / 0.2 * ffac
+ # Shinethrough power (atoms that are not ionised) [MW]:
+ current_drive_variables.p_beam_shine_through_mw = (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ * (1.0 - current_drive_variables.f_p_beam_shine_through)
+ )
- # Current drive efficiency (A/W)
- return gamnb / (dene20 * rmajor)
+ # First orbit loss
+ current_drive_variables.p_beam_orbit_loss_mw = (
+ current_drive_variables.f_p_beam_orbit_loss
+ * (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ - current_drive_variables.p_beam_shine_through_mw
+ )
+ )
- def lheval(self, drfind, rratio):
- """Routine to evaluate the difference between electron energy
- expressions required to find the Lower Hybrid absorption radius
- author: P J Knight, CCFE, Culham Science Centre
- drfind : input real : correction to parallel refractive index
- rratio : input real : guess for radius of penetration / rminor
- ediff : output real : difference between the E values (keV)
- This routine evaluates the difference between the values calculated
- from the two equations for the electron energy E, given in
- AEA FUS 172, p.58. This difference is used to locate the Lower Hybrid
- wave absorption radius via a Newton-Raphson method, in calling
- routine lhrad.
- AEA FUS 172: Physics Assessment for the European Reactor Study
- """
- dlocal = 1.0e-19 * self.plasma_profile.neprofile.calculate_profile_y(
- rratio,
- physics_variables.rhopedn,
- physics_variables.ne0,
- physics_variables.neped,
- physics_variables.nesep,
- physics_variables.alphan,
- )
+ # Power deposited
+ current_drive_variables.p_beam_plasma_coupled_mw = (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ - current_drive_variables.p_beam_shine_through_mw
+ - current_drive_variables.p_beam_orbit_loss_mw
+ )
- # Local electron temperature
+ p_hcd_primary_ions_mw = (
+ pinjmw1 * current_drive_variables.f_p_beam_injected_ions
+ )
+ p_hcd_primary_electrons_mw = pinjmw1 * (
+ 1.0e0 - current_drive_variables.f_p_beam_injected_ions
+ )
- tlocal = self.plasma_profile.teprofile.calculate_profile_y(
- rratio,
- physics_variables.rhopedt,
- physics_variables.te0,
- physics_variables.teped,
- physics_variables.tesep,
- physics_variables.alphat,
- physics_variables.tbeta,
- )
+ current_drive_variables.pwpnb = (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ / current_drive_variables.eta_beam_injector_wall_plug
+ )
- # Local toroidal field (evaluated at the inboard region of the flux surface)
+ # Neutral beam wall plug power
+ heat_transport_variables.p_hcd_primary_electric_mw = (
+ current_drive_variables.pwpnb
+ )
+ current_drive_variables.eta_hcd_primary_injector_wall_plug = (
+ current_drive_variables.eta_beam_injector_wall_plug
+ )
- blocal = (
- physics_variables.bt
- * physics_variables.rmajor
- / (physics_variables.rmajor - rratio * physics_variables.rminor)
- )
+ current_drive_variables.c_beam_total = (
+ 1.0e-3
+ * (
+ (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ )
+ * 1.0e6
+ )
+ / current_drive_variables.e_beam_kev
+ ) # Neutral beam current (A)
- # Parallel refractive index needed for plasma access
+ current_drive_variables.p_hcd_beam_injected_total_mw += (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ )
- frac = np.sqrt(dlocal) / blocal
- nplacc = frac + np.sqrt(1.0e0 + frac * frac)
+ # ===========================================================
- # Total parallel refractive index
+ # Total injected power that contributed to heating
+ current_drive_variables.p_hcd_injected_total_mw = (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_primary_extra_heat_mw
+ + current_drive_variables.p_hcd_secondary_injected_mw
+ + current_drive_variables.p_hcd_secondary_extra_heat_mw
+ )
- refind = nplacc + drfind
+ # Total injected power that contributed to current drive
+ current_drive_variables.p_hcd_injected_current_total_mw = (
+ current_drive_variables.p_hcd_primary_injected_mw
+ + current_drive_variables.p_hcd_secondary_injected_mw
+ )
- # First equation for electron energy E
+ pinjmw1 = p_hcd_primary_electrons_mw + p_hcd_primary_ions_mw
- e1 = 511.0e0 * (np.sqrt(1.0e0 + 1.0e0 / (refind * refind)) - 1.0e0)
+ # Total injected power given to electrons
+ current_drive_variables.p_hcd_injected_electrons_mw = (
+ p_hcd_primary_electrons_mw + p_hcd_secondary_electrons_mw
+ )
- # Second equation for E
+ # Total injected power given to ions
+ current_drive_variables.p_hcd_injected_ions_mw = (
+ p_hcd_primary_ions_mw + p_hcd_secondary_ions_mw
+ )
- e2 = 7.0e0 * tlocal
+ # Total wall plug power for all heating systems
+ heat_transport_variables.p_hcd_electric_total_mw = (
+ heat_transport_variables.p_hcd_primary_electric_mw
+ + heat_transport_variables.p_hcd_secondary_electric_mw
+ )
- # Difference
+ # Reset injected power to zero for ignited plasma (fudge)
+ if physics_variables.ignite == 1:
+ heat_transport_variables.p_hcd_electric_total_mw = 0.0e0
- return e1 - e2
+ # Ratio of fusion to input (injection+ohmic) power
+ current_drive_variables.bigq = physics_variables.fusion_power / (
+ current_drive_variables.p_hcd_injected_total_mw
+ + current_drive_variables.p_beam_orbit_loss_mw
+ + physics_variables.p_plasma_ohmic_mw
+ )
- def etanb(self, m_beam_amu, alphan, alphat, aspect, dene, ebeam, rmajor, ten, zeff):
- """Routine to find neutral beam current drive efficiency
- using the ITER 1990 formulation
- author: P J Knight, CCFE, Culham Science Centre
- m_beam_amu : input real : beam ion mass (amu)
- alphan : input real : density profile factor
- alphat : input real : temperature profile factor
- aspect : input real : aspect ratio
- dene : input real : volume averaged electron density (m**-3)
- 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
- This routine calculates the current drive efficiency of
- a neutral beam system, based on the 1990 ITER model.
- ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al,
- ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990
+ def output_current_drive(self):
+ """
+ Output the current drive information to the output file.
+ This method writes the current drive information to the output file.
"""
- zbeam = 1.0
- bbd = 1.0
-
- dene20 = 1e-20 * dene
+ po.oheadr(self.outfile, "Heating & Current Drive System")
- # Ratio of E_beam/E_crit
- xjs = ebeam / (bbd * 10.0 * m_beam_amu * ten)
- xj = np.sqrt(xjs)
+ if physics_variables.ignite == 1:
+ po.ocmmnt(
+ self.outfile,
+ "Ignited plasma; injected power only used for start-up phase",
+ )
- yj = 0.8 * zeff / m_beam_amu
+ if abs(physics_variables.f_c_plasma_inductive) > 1.0e-8:
+ po.ocmmnt(
+ self.outfile,
+ "Current is driven by both inductive and non-inductive means.",
+ )
+ po.oblnkl(self.outfile)
- rjfunc = xjs / (4.0 + 3.0 * yj + xjs * (xj + 1.39 + 0.61 * yj**0.7))
+ po.ovarre(
+ self.outfile,
+ "Fusion gain factor Q",
+ "(bigq)",
+ current_drive_variables.bigq,
+ "OP ",
+ )
+ po.oblnkl(self.outfile)
- epseff = 0.5 / aspect
- gfac = (1.55 + 0.85 / zeff) * np.sqrt(epseff) - (0.2 + 1.55 / zeff) * epseff
- ffac = 1.0 / zbeam - (1.0 - gfac) / zeff
+ if current_drive_variables.i_hcd_calculations == 0:
+ po.ocmmnt(self.outfile, "No current drive used")
+ po.oblnkl(self.outfile)
+ return
- abd = (
- 0.107
- * (1.0 - 0.35 * alphan + 0.14 * alphan**2)
- * (1.0 - 0.21 * alphat)
- * (1.0 - 0.2e-3 * ebeam + 0.09e-6 * ebeam**2)
+ po.ovarin(
+ self.outfile,
+ "Primary current drive efficiency model",
+ "(i_hcd_primary)",
+ current_drive_variables.i_hcd_primary,
)
- return abd * (5.0 / rmajor) * (0.1 * ten / dene20) * rjfunc / 0.2 * ffac
+ if current_drive_variables.i_hcd_primary in [1, 4, 6]:
+ po.ocmmnt(self.outfile, "Lower Hybrid Current Drive")
+ elif current_drive_variables.i_hcd_primary == 2:
+ po.ocmmnt(self.outfile, "Ion Cyclotron Current Drive")
+ elif current_drive_variables.i_hcd_primary in [3, 7]:
+ po.ocmmnt(self.outfile, "Electron Cyclotron Current Drive")
+ elif current_drive_variables.i_hcd_primary in [5, 8]:
+ po.ocmmnt(self.outfile, "Neutral Beam Current Drive")
+ elif current_drive_variables.i_hcd_primary == 10:
+ po.ocmmnt(
+ self.outfile,
+ "Electron Cyclotron Current Drive (input normalised efficiency)",
+ )
+ elif current_drive_variables.i_hcd_primary == 12:
+ po.ocmmnt(self.outfile, "Electron Bernstein Wave Current Drive")
+ elif current_drive_variables.i_hcd_primary == 13:
+ po.ocmmnt(
+ self.outfile,
+ "Electron Cyclotron Current Drive (with Zeff & Te dependance)",
+ )
- def cfnbi(self, afast, efast, te, ne, _nd, _nt, zeffai, xlmbda):
- """Routine to calculate the fraction of the fast particle energy
- coupled to the ions
- author: P J Knight, CCFE, Culham Science Centre
- afast : input real : mass of fast particle (units of proton mass)
- efast : input real : energy of fast particle (keV)
- te : input real : density weighted average electron temp. (keV)
- ne : input real : volume averaged electron density (m**-3)
- nd : input real : deuterium beam density (m**-3)
- nt : input real : tritium beam density (m**-3)
- zeffai : input real : mass weighted plasma effective charge
- xlmbda : input real : ion-electron coulomb logarithm
- f_p_beam_injected_ions : output real : fraction of fast particle energy coupled to ions
- This routine calculates the fast particle energy coupled to
- the ions in the neutral beam system.
- """
- # atmd = 2.0
- atmdt = 2.5
- # atmt = 3.0
- c = 3.0e8
- me = constants.electron_mass
- # zd = 1.0
- # zt = 1.0
+ po.oblnkl(self.outfile)
+
+ po.ovarre(
+ self.outfile,
+ "Absolute current drive efficiency of primary system [A/W]",
+ "(eta_cd_hcd_primary)",
+ current_drive_variables.eta_cd_hcd_primary,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Normalised current drive efficiency of primary system [10^20 A / Wm^2]",
+ "(eta_cd_norm_hcd_primary)",
+ current_drive_variables.eta_cd_norm_hcd_primary,
+ "OP ",
+ )
+ if current_drive_variables.i_hcd_primary == 10:
+ po.ovarre(
+ self.outfile,
+ "ECRH plasma heating efficiency",
+ "(eta_cd_norm_ecrh)",
+ current_drive_variables.eta_cd_norm_ecrh,
+ )
+ po.ovarre(
+ self.outfile,
+ "Power injected into plasma by primary system for current drive (MW)",
+ "(p_hcd_primary_injected_mw)",
+ current_drive_variables.p_hcd_primary_injected_mw,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Extra power injected into plasma by primary system (MW)",
+ "(p_hcd_primary_extra_heat_mw)",
+ current_drive_variables.p_hcd_primary_extra_heat_mw,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Current driven in plasma by primary system (A)",
+ "(c_hcd_primary_driven)",
+ current_drive_variables.c_hcd_primary_driven,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Fraction of plasma current driven by primary system",
+ "(f_c_plasma_hcd_primary)",
+ current_drive_variables.f_c_plasma_hcd_primary,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Wall plug to injector efficiency of primary system",
+ "(eta_hcd_primary_injector_wall_plug)",
+ current_drive_variables.eta_hcd_primary_injector_wall_plug,
+ "IP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Wall plug electric power of primary system",
+ "(p_hcd_primary_electric_mw)",
+ heat_transport_variables.p_hcd_primary_electric_mw,
+ "OP ",
+ )
+
+ if current_drive_variables.i_hcd_primary in [12, 13]:
+ po.oblnkl(self.outfile)
+ po.ovarre(
+ self.outfile,
+ "ECRH / EBW harmonic number",
+ "(n_ecrh_harmonic)",
+ current_drive_variables.n_ecrh_harmonic,
+ )
+ if current_drive_variables.i_hcd_primary == 13:
+ po.ovarin(
+ self.outfile,
+ "Electron cyclotron cutoff wave mode switch",
+ "(i_ecrh_wave_mode)",
+ current_drive_variables.i_ecrh_wave_mode,
+ )
- # xlbd = self.xlmbdabi(afast, atmd, efast, te, ne)
- # xlbt = self.xlmbdabi(afast, atmt, efast, te, ne)
+ po.oblnkl(self.outfile)
- # sum = nd * zd * zd * xlbd / atmd + nt * zt * zt * xlbt / atmt
- # ecritfix = 16.0e0 * te * afast * (sum / (ne * xlmbda)) ** (2.0e0 / 3.0e0)
+ if current_drive_variables.i_hcd_primary in [5, 8]:
+ po.oblnkl(self.outfile)
+ po.ocmmnt(self.outfile, "Neutral beam power balance :")
+ po.ocmmnt(self.outfile, "----------------------------")
- xlmbdai = self.xlmbdabi(afast, atmdt, efast, te, ne)
- sumln = zeffai * xlmbdai / xlmbda
- xlnrat = (
- 3.0e0 * np.sqrt(np.pi) / 4.0e0 * me / constants.proton_mass * sumln
- ) ** (2.0e0 / 3.0e0)
- ve = c * np.sqrt(2.0e0 * te / 511.0e0)
+ po.ovarre(
+ self.outfile,
+ "Neutral beam energy (keV)",
+ "(e_beam_kev)",
+ current_drive_variables.e_beam_kev,
+ )
+ if (current_drive_variables.i_hcd_primary == 5) or (
+ current_drive_variables.i_hcd_primary == 8
+ ):
+ po.ovarre(
+ self.outfile,
+ "Neutral beam current (A)",
+ "(c_beam_total)",
+ current_drive_variables.c_beam_total,
+ "OP ",
+ )
- ecritfi = (
- afast
- * constants.proton_mass
- * ve
- * ve
- * xlnrat
- / (2.0e0 * constants.electron_charge * 1.0e3)
- )
+ po.ovarre(
+ self.outfile,
+ "Neutral beam wall plug efficiency",
+ "(eta_beam_injector_wall_plug)",
+ current_drive_variables.eta_beam_injector_wall_plug,
+ )
+ po.ovarre(
+ self.outfile,
+ "Beam decay lengths to centre",
+ "(n_beam_decay_lengths_core)",
+ current_drive_variables.n_beam_decay_lengths_core,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Beam shine-through fraction",
+ "(f_p_beam_shine_through)",
+ current_drive_variables.f_p_beam_shine_through,
+ "OP ",
+ )
- x = np.sqrt(efast / ecritfi)
- t1 = np.log((x * x - x + 1.0e0) / ((x + 1.0e0) ** 2))
- thx = (2.0e0 * x - 1.0e0) / np.sqrt(3.0e0)
- t2 = 2.0e0 * np.sqrt(3.0e0) * (np.arctan(thx) + np.pi / 6.0e0)
+ if (current_drive_variables.i_hcd_primary == 5) or (
+ current_drive_variables.i_hcd_primary == 8
+ ):
+ po.ovarrf(
+ self.outfile,
+ "Beam first orbit loss power (MW)",
+ "(p_beam_orbit_loss_mw)",
+ current_drive_variables.p_beam_orbit_loss_mw,
+ "OP ",
+ )
+ po.ovarrf(
+ self.outfile,
+ "Beam shine-through power [MW]",
+ "(p_beam_shine_through_mw)",
+ current_drive_variables.p_beam_shine_through_mw,
+ "OP ",
+ )
+ po.ovarrf(
+ self.outfile,
+ "Maximum allowable beam power (MW)",
+ "(p_hcd_injected_max)",
+ current_drive_variables.p_hcd_injected_max,
+ )
+ po.oblnkl(self.outfile)
+ po.ovarrf(
+ self.outfile,
+ "Beam power entering vacuum vessel (MW)",
+ "(p_beam_injected_mw)",
+ current_drive_variables.p_beam_injected_mw,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Fraction of beam energy to ions",
+ "(f_p_beam_injected_ions)",
+ current_drive_variables.f_p_beam_injected_ions,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Beam duct shielding thickness (m)",
+ "(dx_beam_shield)",
+ current_drive_variables.dx_beam_shield,
+ )
+ po.ovarre(
+ self.outfile,
+ "Beam tangency radius / Plasma major radius",
+ "(frbeam)",
+ current_drive_variables.frbeam,
+ )
+ po.ovarre(
+ self.outfile,
+ "Beam centreline tangency radius (m)",
+ "(rtanbeam)",
+ current_drive_variables.rtanbeam,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Maximum possible tangency radius (m)",
+ "(rtanmax)",
+ current_drive_variables.rtanmax,
+ "OP ",
+ )
- return (t1 + t2) / (3.0e0 * x * x)
+ po.ocmmnt(self.outfile, "----------------------------")
+ po.oblnkl(self.outfile)
+ po.ovarin(
+ self.outfile,
+ "Secondary current drive efficiency model",
+ "(i_hcd_secondary)",
+ current_drive_variables.i_hcd_secondary,
+ )
- def xlmbdabi(self, mb, mth, eb, t, nelec):
- """Calculates the Coulomb logarithm for ion-ion collisions
- author: P J Knight, CCFE, Culham Science Centre
- mb : input real : mass of fast particle (units of proton mass)
- mth : input real : mass of background ions (units of proton mass)
- eb : input real : energy of fast particle (keV)
- t : input real : density weighted average electron temp. (keV)
- nelec : input real : volume averaged electron density (m**-3)
- This function calculates the Coulomb logarithm for ion-ion
- collisions where the relative velocity may be large compared
- with the background ('mt') thermal velocity.
- Mikkelson and Singer, Nuc Tech/Fus, 4, 237 (1983)
- """
+ if current_drive_variables.i_hcd_secondary in [1, 4, 6]:
+ po.ocmmnt(self.outfile, "Lower Hybrid Current Drive")
+ elif current_drive_variables.i_hcd_secondary == 2:
+ po.ocmmnt(self.outfile, "Ion Cyclotron Current Drive")
+ elif current_drive_variables.i_hcd_secondary in [3, 7]:
+ po.ocmmnt(self.outfile, "Electron Cyclotron Current Drive")
+ elif current_drive_variables.i_hcd_secondary in [5, 8]:
+ po.ocmmnt(self.outfile, "Neutral Beam Current Drive")
+ elif current_drive_variables.i_hcd_secondary == 10:
+ po.ocmmnt(
+ self.outfile,
+ "Electron Cyclotron Current Drive (input normalised efficiency)",
+ )
+ elif current_drive_variables.i_hcd_secondary == 12:
+ po.ocmmnt(self.outfile, "Electron Bernstein Wave Current Drive")
+ elif current_drive_variables.i_hcd_secondary == 13:
+ po.ocmmnt(
+ self.outfile,
+ "Electron Cyclotron Current Drive (with Zeff & Te dependance)",
+ )
+ po.oblnkl(self.outfile)
- x1 = (t / 10.0) * (eb / 1000.0) * mb / (nelec / 1e20)
- x2 = mth / (mth + mb)
+ po.ovarre(
+ self.outfile,
+ "Absolute current drive efficiency of secondary system [A/W]",
+ "(eta_cd_hcd_secondary)",
+ current_drive_variables.eta_cd_hcd_secondary,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Normalised current drive efficiency of secondary system [10^20 A / Wm^2]",
+ "(eta_cd_norm_hcd_secondary)",
+ current_drive_variables.eta_cd_norm_hcd_secondary,
+ "OP ",
+ )
+ if current_drive_variables.i_hcd_secondary == 10:
+ po.ovarre(
+ self.outfile,
+ "ECRH plasma heating efficiency",
+ "(eta_cd_norm_ecrh)",
+ current_drive_variables.eta_cd_norm_ecrh,
+ )
- return 23.7 + np.log(x2 * np.sqrt(x1))
+ po.ovarre(
+ self.outfile,
+ "Power injected into plasma by secondary system (MW)",
+ "(p_hcd_secondary_injected_mw)",
+ current_drive_variables.p_hcd_secondary_injected_mw,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Extra power injected into plasma by secondary system (MW)",
+ "(p_hcd_secondary_extra_heat_mw)",
+ current_drive_variables.p_hcd_secondary_extra_heat_mw,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Current driven in plasma by secondary system (A)",
+ "(c_hcd_secondary_driven)",
+ current_drive_variables.c_hcd_secondary_driven,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Fraction of plasma current driven by secondary system",
+ "(f_c_plasma_hcd_secondary)",
+ current_drive_variables.f_c_plasma_hcd_secondary,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Wall plug to injector efficiency of secondary system",
+ "(eta_hcd_secondary_injector_wall_plug)",
+ current_drive_variables.eta_hcd_secondary_injector_wall_plug,
+ "IP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Wall plug electric power of secondary system",
+ "(p_hcd_secondary_electric_mw)",
+ heat_transport_variables.p_hcd_secondary_electric_mw,
+ "OP ",
+ )
- def legend(self, zlocal, arg):
- """Routine to calculate Legendre function and its derivative
- author: M R O'Brien, CCFE, Culham Science Centre
- author: P J Knight, CCFE, Culham Science Centre
- zlocal : input real : local plasma effective charge
- arg : input real : argument of Legendre function
- palpha : output real : value of Legendre function
- palphap : output real : derivative of Legendre function
- This routine calculates the Legendre function palpha
- of argument arg and order
- alpha = -0.5 + i sqrt(xisq),
- and its derivative palphap.
-
This Legendre function is a conical function and we use the series
- in xisq given in Abramowitz and Stegun. The
- derivative is calculated from the derivative of this series.
-
The derivatives were checked by calculating palpha for
- neighbouring arguments. The calculation of palpha for zero
- argument was checked by comparison with the expression
- palpha(0) = 1/sqrt(pi) * cos(pi*alpha/2) * gam1 / gam2
- (Abramowitz and Stegun, eqn 8.6.1). Here gam1 and
- gam2 are the Gamma functions of arguments
- 0.5*(1+alpha) and 0.5*(2+alpha) respectively.
- Abramowitz and Stegun, equation 8.12.1
- """
- if abs(arg) > (1.0e0 + 1.0e-10):
- eh.fdiags[0] = arg
- raise ProcessValueError("Invalid argument", arg=arg)
+ po.oblnkl(self.outfile)
- arg2 = min(arg, (1.0e0 - 1.0e-10))
- sinsq = 0.5e0 * (1.0e0 - arg2)
- xisq = 0.25e0 * (32.0e0 * zlocal / (zlocal + 1.0e0) - 1.0e0)
- palpha = 1.0e0
- pold = 1.0e0
- pterm = 1.0e0
- palphap = 0.0e0
- poldp = 0.0e0
+ if current_drive_variables.i_hcd_secondary in [5, 8]:
+ po.oblnkl(self.outfile)
+ po.ocmmnt(self.outfile, "Neutral beam power balance :")
+ po.ocmmnt(self.outfile, "----------------------------")
- for n in range(10000):
- # Check for convergence every 20 iterations
+ po.ovarre(
+ self.outfile,
+ "Neutral beam energy (keV)",
+ "(e_beam_kev)",
+ current_drive_variables.e_beam_kev,
+ )
+ if (current_drive_variables.i_hcd_primary == 5) or (
+ current_drive_variables.i_hcd_primary == 8
+ ):
+ po.ovarre(
+ self.outfile,
+ "Neutral beam current (A)",
+ "(c_beam_total)",
+ current_drive_variables.c_beam_total,
+ "OP ",
+ )
+
+ po.ovarre(
+ self.outfile,
+ "Neutral beam wall plug efficiency",
+ "(eta_beam_injector_wall_plug)",
+ current_drive_variables.eta_beam_injector_wall_plug,
+ )
+ po.ovarre(
+ self.outfile,
+ "Beam decay lengths to centre",
+ "(n_beam_decay_lengths_core)",
+ current_drive_variables.n_beam_decay_lengths_core,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Beam shine-through fraction",
+ "(f_p_beam_shine_through)",
+ current_drive_variables.f_p_beam_shine_through,
+ "OP ",
+ )
- if (n > 1) and ((n % 20) == 1):
- term1 = 1.0e-10 * max(abs(pold), abs(palpha))
- term2 = 1.0e-10 * max(abs(poldp), abs(palphap))
+ if (current_drive_variables.i_hcd_primary == 5) or (
+ current_drive_variables.i_hcd_primary == 8
+ ):
+ po.ovarrf(
+ self.outfile,
+ "Beam first orbit loss power (MW)",
+ "(p_beam_orbit_loss_mw)",
+ current_drive_variables.p_beam_orbit_loss_mw,
+ "OP ",
+ )
+ po.ovarrf(
+ self.outfile,
+ "Beam shine-through power [MW]",
+ "(p_beam_shine_through_mw)",
+ current_drive_variables.p_beam_shine_through_mw,
+ "OP ",
+ )
+ po.ovarrf(
+ self.outfile,
+ "Maximum allowable beam power (MW)",
+ "(p_hcd_injected_max)",
+ current_drive_variables.p_hcd_injected_max,
+ )
+ po.oblnkl(self.outfile)
+ po.ovarrf(
+ self.outfile,
+ "Beam power entering vacuum vessel (MW)",
+ "(p_beam_injected_mw)",
+ current_drive_variables.p_beam_injected_mw,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Fraction of beam energy to ions",
+ "(f_p_beam_injected_ions)",
+ current_drive_variables.f_p_beam_injected_ions,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Beam duct shielding thickness (m)",
+ "(dx_beam_shield)",
+ current_drive_variables.dx_beam_shield,
+ )
+ po.ovarre(
+ self.outfile,
+ "Beam tangency radius / Plasma major radius",
+ "(frbeam)",
+ current_drive_variables.frbeam,
+ )
+ po.ovarre(
+ self.outfile,
+ "Beam centreline tangency radius (m)",
+ "(rtanbeam)",
+ current_drive_variables.rtanbeam,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Maximum possible tangency radius (m)",
+ "(rtanmax)",
+ current_drive_variables.rtanmax,
+ "OP ",
+ )
- if (abs(pold - palpha) < term1) and (abs(poldp - palphap) < term2):
- return palpha, palphap
+ po.ocmmnt(self.outfile, "----------------------------")
- pold = palpha
- poldp = palphap
+ po.osubhd(self.outfile, "Totals :")
- pterm = (
- pterm
- * (4.0e0 * xisq + (2.0e0 * n - 1.0e0) ** 2)
- / (2.0e0 * n) ** 2
- * sinsq
- )
- palpha = palpha + pterm
- palphap = palphap - n * pterm / (1.0e0 - arg2)
- else:
- raise ProcessError("legend: Solution has not converged")
+ po.ovarre(
+ self.outfile,
+ "Total injected heating power that drove plasma current (MW)",
+ "(p_hcd_injected_current_total_mw)",
+ current_drive_variables.p_hcd_injected_current_total_mw,
+ )
+ po.ovarre(
+ self.outfile,
+ "Total injected heating power across all systems (MW)",
+ "(p_hcd_injected_total_mw)",
+ current_drive_variables.p_hcd_injected_total_mw,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Total injected heating power given to the electrons (MW)",
+ "(p_hcd_injected_electrons_mw)",
+ current_drive_variables.p_hcd_injected_electrons_mw,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Total injected heating power given to the ions (MW)",
+ "(p_hcd_injected_ions_mw)",
+ current_drive_variables.p_hcd_injected_ions_mw,
+ "OP ",
+ )
- def sigbeam(self, eb, te, ne, rnhe, rnc, rno, rnfe):
- """Calculates the stopping cross-section for a hydrogen
- beam in a fusion plasma
- author: P J Knight, CCFE, Culham Science Centre
- eb : input real : beam energy (kev/amu)
- te : input real : electron temperature (keV)
- ne : input real : electron density (10^20m-3)
- rnhe : input real : alpha density / ne
- rnc : input real : carbon density /ne
- rno : input real : oxygen density /ne
- rnfe : input real : iron density /ne
- This function calculates the stopping cross-section (m^-2)
- for a hydrogen beam in a fusion plasma.
- Janev, Boley and Post, Nuclear Fusion 29 (1989) 2125
- """
- a = np.array([
- [
- [4.4, -2.49e-2],
- [7.46e-2, 2.27e-3],
- [3.16e-3, -2.78e-5],
- ],
- [
- [2.3e-1, -1.15e-2],
- [-2.55e-3, -6.2e-4],
- [1.32e-3, 3.38e-5],
- ],
- ])
+ po.ovarre(
+ self.outfile,
+ "Upper limit on total plasma injected power (MW)",
+ "(p_hcd_injected_max)",
+ current_drive_variables.p_hcd_injected_max,
+ "OP ",
+ )
- b = np.array([
- [
- [[-2.36, -1.49, -1.41, -1.03], [0.185, -0.0154, -4.08e-4, 0.106]],
- [
- [-0.25, -0.119, -0.108, -0.0558],
- [-0.0381, -0.015, -0.0138, -3.72e-3],
- ],
- ],
- [
- [
- [0.849, 0.518, 0.477, 0.322],
- [-0.0478, 7.18e-3, 1.57e-3, -0.0375],
- ],
- [
- [0.0677, 0.0292, 0.0259, 0.0124],
- [0.0105, 3.66e-3, 3.33e-3, 8.61e-4],
- ],
- ],
- [
- [
- [-0.0588, -0.0336, -0.0305, -0.0187],
- [4.34e-3, 3.41e-4, 7.35e-4, 3.53e-3],
- ],
- [
- [-4.48e-3, -1.79e-3, -1.57e-3, -7.43e-4],
- [-6.76e-4, -2.04e-4, -1.86e-4, -5.12e-5],
- ],
- ],
- ])
+ po.osubhd(self.outfile, "Contributions:")
- z = np.array([2.0, 6.0, 8.0, 26.0])
- nn = np.array([rnhe, rnc, rno, rnfe])
+ po.ovarre(
+ self.outfile,
+ "Injected power into plasma from lower hybrid systems (MW)",
+ "(p_hcd_lowhyb_injected_total_mw)",
+ current_drive_variables.p_hcd_lowhyb_injected_total_mw,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Injected power into plasma from ion cyclotron systems (MW)",
+ "(p_hcd_icrh_injected_total_mw)",
+ current_drive_variables.p_hcd_icrh_injected_total_mw,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Injected power into plasma from electron cyclotron systems (MW)",
+ "(p_hcd_ecrh_injected_total_mw)",
+ current_drive_variables.p_hcd_ecrh_injected_total_mw,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Injected power into plasma from neutral beam systems (MW)",
+ "(p_hcd_beam_injected_total_mw)",
+ current_drive_variables.p_hcd_beam_injected_total_mw,
+ "OP ",
+ )
+ po.ovarre(
+ self.outfile,
+ "Injected power into plasma from lower hybrid systems (MW)",
+ "(p_hcd_ebw_injected_total_mw)",
+ current_drive_variables.p_hcd_ebw_injected_total_mw,
+ "OP ",
+ )
- nen = ne * 1e-19
+ po.osubhd(self.outfile, "Fractions of current drive :")
+ po.ovarrf(
+ self.outfile,
+ "Bootstrap fraction",
+ "(f_c_plasma_bootstrap)",
+ current_drive_variables.f_c_plasma_bootstrap,
+ "OP ",
+ )
+ po.ovarrf(
+ self.outfile,
+ "Diamagnetic fraction",
+ "(f_c_plasma_diamagnetic)",
+ current_drive_variables.f_c_plasma_diamagnetic,
+ "OP ",
+ )
+ po.ovarrf(
+ self.outfile,
+ "Pfirsch-Schlueter fraction",
+ "(f_c_plasma_pfirsch_schluter)",
+ current_drive_variables.f_c_plasma_pfirsch_schluter,
+ "OP ",
+ )
+ po.ovarrf(
+ self.outfile,
+ "Auxiliary current drive fraction",
+ "(f_c_plasma_auxiliary)",
+ physics_variables.f_c_plasma_auxiliary,
+ "OP ",
+ )
+ po.ovarrf(
+ self.outfile,
+ "Inductive fraction",
+ "(f_c_plasma_inductive)",
+ physics_variables.f_c_plasma_inductive,
+ "OP ",
+ )
- s1 = 0.0
- for k in range(2):
- for j in range(3):
- for i in range(2):
- s1 += (
- a[i, j, k]
- * (np.log(eb)) ** i
- * (np.log(nen)) ** j
- * (np.log(te)) ** k
- )
+ # MDK Add physics_variables.f_c_plasma_non_inductive as it can be an iteration variable
+ po.ovarrf(
+ self.outfile,
+ "Fraction of the plasma current produced by non-inductive means",
+ "(f_c_plasma_non_inductive)",
+ physics_variables.f_c_plasma_non_inductive,
+ )
- sz = 0.0
- for l in range(4): # noqa: E741
- for k in range(2):
- for j in range(2):
- for i in range(3):
- sz += (
- b[i, j, k, l]
- * (np.log(eb)) ** i
- * (np.log(nen)) ** j
- * (np.log(te)) ** k
- * nn[l]
- * z[l]
- * (z[l] - 1.0)
- )
+ if (
+ abs(
+ current_drive_variables.f_c_plasma_bootstrap
+ - current_drive_variables.f_c_plasma_bootstrap_max
+ )
+ < 1.0e-8
+ ):
+ po.ocmmnt(self.outfile, "Warning : bootstrap current fraction is at")
+ po.ocmmnt(self.outfile, " its prescribed maximum.")
- return max(1e-20 * (np.exp(s1) / eb * (1.0 + sz)), 1e-23)
+ po.oblnkl(self.outfile)
def init_current_drive_variables():
@@ -2233,7 +2458,7 @@ def init_current_drive_variables():
current_drive_variables.f_c_plasma_diamagnetic_hender = 0.0
current_drive_variables.f_c_plasma_diamagnetic_scene = 0.0
current_drive_variables.f_c_plasma_diamagnetic = 0.0
- current_drive_variables.p_ecrh_injected_mw = 0.0
+ current_drive_variables.p_hcd_ecrh_injected_total_mw = 0.0
current_drive_variables.echwpow = 0.0
current_drive_variables.eta_cd_hcd_primary = 0.0
current_drive_variables.n_ecrh_harmonic = 2.0
@@ -2266,10 +2491,11 @@ def init_current_drive_variables():
current_drive_variables.p_hcd_injected_electrons_mw = 0.0
current_drive_variables.p_hcd_injected_ions_mw = 0.0
current_drive_variables.p_hcd_injected_total_mw = 0.0
+ current_drive_variables.p_hcd_injected_current_total_mw = 0.0
current_drive_variables.p_hcd_secondary_injected_mw = 0.0
current_drive_variables.f_c_plasma_internal = 0.0
- current_drive_variables.plhybd = 0.0
- current_drive_variables.pnbeam = 0.0
+ current_drive_variables.p_hcs_lowhyb_injected_total_mw = 0.0
+ current_drive_variables.p_hcd_beam_injected_mw = 0.0
current_drive_variables.p_beam_orbit_loss_mw = 0.0
current_drive_variables.f_c_plasma_pfirsch_schluter = 0.0
current_drive_variables.pwplh = 0.0
diff --git a/process/ife.py b/process/ife.py
index 3391d558be..6daa2aff69 100644
--- a/process/ife.py
+++ b/process/ife.py
@@ -1954,11 +1954,13 @@ def ifepw1(self):
# Wall plug driver power (MW)
- heat_transport_variables.pinjwp = pdrvmw / ife_variables.etadrv
+ heat_transport_variables.p_hcd_electric_total_mw = pdrvmw / ife_variables.etadrv
# Waste driver power (MW)
- heat_transport_variables.pinjht = heat_transport_variables.pinjwp - pdrvmw
+ heat_transport_variables.pinjht = (
+ heat_transport_variables.p_hcd_electric_total_mw - pdrvmw
+ )
# Cryogenic power (MW)
# Cryogenic temperature is assumed to be 4.5K
@@ -2053,8 +2055,8 @@ def ifepw2(self, output: bool = False):
process_output.ovarre(
self.outfile,
"Driver wall plug power (MW)",
- "(pinjwp)",
- heat_transport_variables.pinjwp,
+ "(p_hcd_electric_total_mw)",
+ heat_transport_variables.p_hcd_electric_total_mw,
)
process_output.ovarre(
self.outfile,
@@ -2175,7 +2177,7 @@ def ifeacp(self, output: bool = False):
+ ife_variables.tfacmw
+ (ife_variables.htpmw_ife * ife_variables.reprat / 6.0)
+ heat_transport_variables.trithtmw
- + heat_transport_variables.pinjwp
+ + heat_transport_variables.p_hcd_electric_total_mw
+ basemw
+ (buildings_variables.efloor * pmwpm2)
+ ife_variables.lipmw
@@ -2216,8 +2218,8 @@ def ifeacp(self, output: bool = False):
process_output.ovarre(
self.outfile,
"Driver power supplies (MW)",
- "(pinjwp)",
- heat_transport_variables.pinjwp,
+ "(p_hcd_electric_total_mw)",
+ heat_transport_variables.p_hcd_electric_total_mw,
)
process_output.ovarre(
self.outfile,
diff --git a/process/input.py b/process/input.py
index be45298a1a..e52ffcce33 100644
--- a/process/input.py
+++ b/process/input.py
@@ -711,7 +711,9 @@ def __post_init__(self):
"fvolsi": InputVariable(fortran.fwbs_variables, float, range=(0.0, 10.0)),
"fvolso": InputVariable(fortran.fwbs_variables, float, range=(0.0, 10.0)),
"fvs": InputVariable(fortran.constraint_variables, float, range=(0.001, 10.0)),
- "fvsbrnni": InputVariable(fortran.physics_variables, float, range=(0.0, 1.0)),
+ "f_c_plasma_non_inductive": InputVariable(
+ fortran.physics_variables, float, range=(0.0, 1.0)
+ ),
"fvs_cs_pf_total_ramp": InputVariable(
fortran.pfcoil_variables, float, range=(0.001, 10.0)
),
diff --git a/process/io/mfile_comparison.py b/process/io/mfile_comparison.py
index 2a6c4658cc..d3cd5536b3 100644
--- a/process/io/mfile_comparison.py
+++ b/process/io/mfile_comparison.py
@@ -86,8 +86,8 @@
"p_plasma_separatrix_mw",
"p_hcd_primary_extra_heat_mw",
"f_c_plasma_bootstrap",
- "aux_current_fraction",
- "inductive_current_fraction",
+ "f_c_plasma_auxiliary",
+ "f_c_plasma_inductive",
"gamnb",
"e_beam_kev",
"p_plasma_loss_mw",
@@ -175,8 +175,8 @@
"p_hcd_injected_total_mw",
"p_hcd_primary_extra_heat_mw",
"f_c_plasma_bootstrap",
- "aux_current_fraction",
- "inductive_current_fraction",
+ "f_c_plasma_auxiliary",
+ "f_c_plasma_inductive",
"gamnb",
"e_beam_kev",
"p_plasma_loss_mw",
diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
index 4e8788ec8f..019623dcc2 100644
--- a/process/io/plot_proc.py
+++ b/process/io/plot_proc.py
@@ -3125,8 +3125,8 @@ def plot_current_drive_info(axis, mfile_data, scan):
(pinjie, "Steady state auxiliary power", "MW"),
("p_hcd_primary_extra_heat_mw", "Power for heating only", "MW"),
("f_c_plasma_bootstrap", "Bootstrap fraction", ""),
- ("aux_current_fraction", "Auxiliary fraction", ""),
- ("inductive_current_fraction", "Inductive fraction", ""),
+ ("f_c_plasma_auxiliary", "Auxiliary fraction", ""),
+ ("f_c_plasma_inductive", "Inductive fraction", ""),
("p_plasma_loss_mw", "Plasma heating used for H factor", "MW"),
(
"eta_cd_hcd_primary",
@@ -3155,8 +3155,8 @@ def plot_current_drive_info(axis, mfile_data, scan):
(pinjie, "Steady state auxiliary power", "MW"),
("p_hcd_primary_extra_heat_mw", "Power for heating only", "MW"),
("f_c_plasma_bootstrap", "Bootstrap fraction", ""),
- ("aux_current_fraction", "Auxiliary fraction", ""),
- ("inductive_current_fraction", "Inductive fraction", ""),
+ ("f_c_plasma_auxiliary", "Auxiliary fraction", ""),
+ ("f_c_plasma_inductive", "Inductive fraction", ""),
("gamnb", "NB gamma", "$10^{20}$ A W$^{-1}$ m$^{-2}$"),
("e_beam_kev", "NB energy", "keV"),
("p_plasma_loss_mw", "Plasma heating used for H factor", "MW"),
@@ -3181,8 +3181,8 @@ def plot_current_drive_info(axis, mfile_data, scan):
(pinjie, "Steady state auxiliary power", "MW"),
("p_hcd_primary_extra_heat_mw", "Power for heating only", "MW"),
("f_c_plasma_bootstrap", "Bootstrap fraction", ""),
- ("aux_current_fraction", "Auxiliary fraction", ""),
- ("inductive_current_fraction", "Inductive fraction", ""),
+ ("f_c_plasma_auxiliary", "Auxiliary fraction", ""),
+ ("f_c_plasma_inductive", "Inductive fraction", ""),
("p_plasma_loss_mw", "Plasma heating used for H factor", "MW"),
(
"eta_cd_norm_hcd_primary",
@@ -3210,8 +3210,8 @@ def plot_current_drive_info(axis, mfile_data, scan):
(pinjie, "Steady state auxiliary power", "MW"),
("p_hcd_primary_extra_heat_mw", "Power for heating only", "MW"),
("f_c_plasma_bootstrap", "Bootstrap fraction", ""),
- ("aux_current_fraction", "Auxiliary fraction", ""),
- ("inductive_current_fraction", "Inductive fraction", ""),
+ ("f_c_plasma_auxiliary", "Auxiliary fraction", ""),
+ ("f_c_plasma_inductive", "Inductive fraction", ""),
("p_plasma_loss_mw", "Plasma heating used for H factor", "MW"),
(
"eta_cd_norm_hcd_primary",
@@ -3239,8 +3239,8 @@ def plot_current_drive_info(axis, mfile_data, scan):
(pinjie, "Steady state auxiliary power", "MW"),
("p_hcd_primary_extra_heat_mw", "Power for heating only", "MW"),
("f_c_plasma_bootstrap", "Bootstrap fraction", ""),
- ("aux_current_fraction", "Auxiliary fraction", ""),
- ("inductive_current_fraction", "Inductive fraction", ""),
+ ("f_c_plasma_auxiliary", "Auxiliary fraction", ""),
+ ("f_c_plasma_inductive", "Inductive fraction", ""),
("p_plasma_loss_mw", "Plasma heating used for H factor", "MW"),
(
"eta_cd_norm_hcd_primary",
diff --git a/process/io/sankey_funcs.py b/process/io/sankey_funcs.py
index fef61d37b8..f66585cf55 100644
--- a/process/io/sankey_funcs.py
+++ b/process/io/sankey_funcs.py
@@ -635,7 +635,9 @@ def plot_sankey(mfilename="MFILE.DAT"): # Plot simplified power flow Sankey Dia
# Energy required for rest of power plant (MW)
pcoresystems = crypmw + fachtmw + tfacpd + trithtmw + vachtmw + pfwpmw + ppumpmw
- pinjwp = m_file.data["pinjwp"].get_scan(-1) # injector wall plug power (MW)
+ p_hcd_electric_total_mw = m_file.data["p_hcd_electric_total_mw"].get_scan(
+ -1
+ ) # injector wall plug power (MW)
htpmw = m_file.data["htpmw"].get_scan(
-1
) # heat transport system electrical pump power (MW)
@@ -762,13 +764,13 @@ def plot_sankey(mfilename="MFILE.DAT"): # Plot simplified power flow Sankey Dia
# -------------------------------- RECIRCULATING POWER - 5 --------------------------------
# Recirculated power, -Core Systems, -Heating System
- recirc = [precircmw, -pcoresystems - htpmw, -pinjwp + ppumpmw]
+ recirc = [precircmw, -pcoresystems - htpmw, -p_hcd_electric_total_mw + ppumpmw]
# Check if difference >2 between recirculated power and the output sum
if sum(recirc) ** 2 > 2:
print(
"Recirc. Power Balance",
precircmw,
- -pcoresystems + ppumpmw - pinjwp - htpmw,
+ -pcoresystems + ppumpmw - p_hcd_electric_total_mw - htpmw,
)
sankey.add(
flows=recirc,
@@ -783,9 +785,9 @@ def plot_sankey(mfilename="MFILE.DAT"): # Plot simplified power flow Sankey Dia
# HCD: Heating system, -Plasma heating, -losses
hcd = [
- pinjwp - ppumpmw,
+ p_hcd_electric_total_mw - ppumpmw,
-p_hcd_injected_total_mw,
- -pinjwp + p_hcd_injected_total_mw + ppumpmw,
+ -p_hcd_electric_total_mw + p_hcd_injected_total_mw + ppumpmw,
]
sankey.add(
flows=hcd,
@@ -860,12 +862,12 @@ def plot_sankey(mfilename="MFILE.DAT"): # Plot simplified power flow Sankey Dia
if pnetelmw >= 1:
t.set_position((
pos[0] + 0.15,
- pos[1] + 0.5 * (pinjwp / totalplasma) + 0.2,
+ pos[1] + 0.5 * (p_hcd_electric_total_mw / totalplasma) + 0.2,
))
if pnetelmw < 1:
t.set_position((
pos[0] + 0.15,
- pos[1] + 0.5 * (pinjwp / totalplasma) + 0.2,
+ pos[1] + 0.5 * (p_hcd_electric_total_mw / totalplasma) + 0.2,
))
if t == diagrams[6].texts[1]: # Plasma Heating
t.set_horizontalalignment("left")
@@ -878,6 +880,10 @@ def plot_sankey(mfilename="MFILE.DAT"): # Plot simplified power flow Sankey Dia
t.set_position((
pos[0] + 0.15,
pos[1]
- - 0.5 * ((pinjwp - p_hcd_injected_total_mw) / totalplasma)
+ - 0.5
+ * (
+ (p_hcd_electric_total_mw - p_hcd_injected_total_mw)
+ / totalplasma
+ )
- 0.2,
))
diff --git a/process/io/variable_metadata.py b/process/io/variable_metadata.py
index 1220eccadc..9aa213f9c1 100644
--- a/process/io/variable_metadata.py
+++ b/process/io/variable_metadata.py
@@ -232,10 +232,10 @@ class VariableMetadata:
latex=r"$\eta_{\mathrm{CD}}$[$A/W$]", description="CD efficiency", units="A/W"
),
"bigq": VariableMetadata(latex=r"$Q$", description="Plasma Q value", units=""),
- "aux_current_fraction": VariableMetadata(
+ "f_c_plasma_auxiliary": VariableMetadata(
latex=r"$f_{\mathrm{CD}}$", description="CD factor", units=""
),
- "inductive_current_fraction": VariableMetadata(
+ "f_c_plasma_inductive": VariableMetadata(
latex=r"$f_{\mathrm{CD,ind}}$", description="Inductive CD factor", units=""
),
"f_c_plasma_bootstrap": VariableMetadata(
diff --git a/process/iteration_variables.py b/process/iteration_variables.py
index 1bec00495c..411dd6efda 100644
--- a/process/iteration_variables.py
+++ b/process/iteration_variables.py
@@ -83,7 +83,9 @@ class IterationVariable:
"f_j_cs_start_pulse_end_flat_top", fortran.pfcoil_variables, 0.001, 1.0
),
42: IterationVariable("dr_cs_tf_gap", fortran.build_variables, 0.001, 10.00),
- 44: IterationVariable("fvsbrnni", fortran.physics_variables, 0.001, 1.0),
+ 44: IterationVariable(
+ "f_c_plasma_non_inductive", fortran.physics_variables, 0.001, 1.0
+ ),
45: IterationVariable("fqval", fortran.constraint_variables, 0.001, 1.0),
46: IterationVariable("fpinj", fortran.constraint_variables, 0.001, 1.0),
47: IterationVariable("feffcd", fortran.current_drive_variables, 0.001, 1.0),
diff --git a/process/main.py b/process/main.py
index c1e022ba8e..8b3c8c5312 100644
--- a/process/main.py
+++ b/process/main.py
@@ -59,7 +59,14 @@
from process.costs_2015 import Costs2015
from process.cryostat import Cryostat
from process.cs_fatigue import CsFatigue
-from process.current_drive import CurrentDrive
+from process.current_drive import (
+ CurrentDrive,
+ ElectronBernstein,
+ ElectronCyclotron,
+ IonCyclotron,
+ LowerHybrid,
+ NeutralBeam,
+)
from process.dcll import DCLL
from process.divertor import Divertor
from process.fw import Fw
@@ -675,7 +682,14 @@ def __init__(self):
self.fw = Fw()
self.blanket_library = BlanketLibrary(fw=self.fw)
self.ccfe_hcpb = CCFE_HCPB(blanket_library=self.blanket_library)
- self.current_drive = CurrentDrive(plasma_profile=self.plasma_profile)
+ self.current_drive = CurrentDrive(
+ plasma_profile=self.plasma_profile,
+ electron_cyclotron=ElectronCyclotron(plasma_profile=self.plasma_profile),
+ ion_cyclotron=IonCyclotron(plasma_profile=self.plasma_profile),
+ lower_hybrid=LowerHybrid(plasma_profile=self.plasma_profile),
+ neutral_beam=NeutralBeam(plasma_profile=self.plasma_profile),
+ electron_bernstein=ElectronBernstein(plasma_profile=self.plasma_profile),
+ )
self.physics = Physics(
plasma_profile=self.plasma_profile, current_drive=self.current_drive
)
diff --git a/process/output.py b/process/output.py
index 16cc4f67d9..98c95c4060 100644
--- a/process/output.py
+++ b/process/output.py
@@ -43,7 +43,7 @@ def write(models, _outfile):
models.physics.outplas()
# TODO what is this? Not in caller.f90?
- models.current_drive.cudriv(output=True)
+ models.current_drive.output_current_drive()
# Pulsed reactor model
models.pulse.run(output=True)
diff --git a/process/pfcoil.py b/process/pfcoil.py
index 2b62b0de63..7b13760b45 100644
--- a/process/pfcoil.py
+++ b/process/pfcoil.py
@@ -1216,7 +1216,7 @@ def ohcalc(self):
# Calculation of CS fatigue
# this is only valid for pulsed reactor design
- if pv.inductive_current_fraction > 0.0e-4:
+ if pv.f_c_plasma_inductive > 0.0e-4:
csfv.n_cycle, csfv.t_crack_radial = self.cs_fatigue.ncycle(
pf.sig_hoop,
csfv.residual_sig_hoop,
@@ -2330,7 +2330,7 @@ def outpf(self):
tfv.tmargmin_cs,
)
# only output CS fatigue model for pulsed reactor design
- if pv.inductive_current_fraction > 0.0e-4:
+ if pv.f_c_plasma_inductive > 0.0e-4:
op.ovarre(
self.outfile,
"Residual hoop stress in CS Steel (Pa)",
diff --git a/process/physics.py b/process/physics.py
index 3b42a64da8..ea172b1602 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -63,7 +63,7 @@ def rether(alphan, alphat, dene, dlamie, te, ti, zeffai):
def calculate_volt_second_requirements(
csawth: float,
eps: float,
- inductive_current_fraction: float,
+ f_c_plasma_inductive: float,
ejima_coeff: float,
kappa: float,
rmajor: float,
@@ -79,8 +79,8 @@ def calculate_volt_second_requirements(
:type csawth: float
:param eps: Inverse aspect ratio
:type eps: float
- :param inductive_current_fraction: Fraction of plasma current produced inductively
- :type inductive_current_fraction: float
+ :param f_c_plasma_inductive: Fraction of plasma current produced inductively
+ :type f_c_plasma_inductive: float
:param ejima_coeff: Ejima coefficient for resistive start-up V-s component
:type ejima_coeff: float
:param kappa: Plasma elongation
@@ -163,7 +163,7 @@ def calculate_volt_second_requirements(
# Include enhancement factor in flattop V-s requirement
# to account for MHD sawtooth effects.
- v_plasma_loop_burn = plasma_current * res_plasma * inductive_current_fraction
+ v_plasma_loop_burn = plasma_current * res_plasma * f_c_plasma_inductive
v_burn_resistive = v_plasma_loop_burn * csawth
@@ -2021,26 +2021,30 @@ def physics(self):
# produced by non-inductive means (which also includes
# the current drive proportion)
physics_module.err243 = 0
- if current_drive_variables.f_c_plasma_internal > physics_variables.fvsbrnni:
+ if (
+ current_drive_variables.f_c_plasma_internal
+ > physics_variables.f_c_plasma_non_inductive
+ ):
current_drive_variables.f_c_plasma_internal = min(
current_drive_variables.f_c_plasma_internal,
- physics_variables.fvsbrnni,
+ physics_variables.f_c_plasma_non_inductive,
)
physics_module.err243 = 1
# Fraction of plasma current produced by inductive means
- physics_variables.inductive_current_fraction = max(
- 1.0e-10, (1.0e0 - physics_variables.fvsbrnni)
+ physics_variables.f_c_plasma_inductive = max(
+ 1.0e-10, (1.0e0 - physics_variables.f_c_plasma_non_inductive)
)
# Fraction of plasma current produced by auxiliary current drive
- physics_variables.aux_current_fraction = (
- physics_variables.fvsbrnni - current_drive_variables.f_c_plasma_internal
+ physics_variables.f_c_plasma_auxiliary = (
+ physics_variables.f_c_plasma_non_inductive
+ - current_drive_variables.f_c_plasma_internal
)
# Auxiliary current drive power calculations
if current_drive_variables.i_hcd_calculations != 0:
- self.current_drive.cudriv(False)
+ self.current_drive.cudriv()
# ***************************** #
# FUSION REACTIONS #
@@ -2074,7 +2078,7 @@ def physics(self):
physics_variables.beta_beam,
physics_variables.beam_density_out,
physics_variables.alpha_power_beams,
- ) = physics_funcs.beam_fusion(
+ ) = reactions.beam_fusion(
physics_variables.beamfus0,
physics_variables.betbm0,
physics_variables.bp,
@@ -2243,7 +2247,7 @@ def physics(self):
physics_variables.f_res_plasma_neo,
physics_variables.res_plasma,
) = self.plasma_ohmic_heating(
- physics_variables.inductive_current_fraction,
+ physics_variables.f_c_plasma_inductive,
physics_variables.kappa95,
physics_variables.plasma_current,
physics_variables.rmajor,
@@ -2400,7 +2404,7 @@ def physics(self):
) = calculate_volt_second_requirements(
physics_variables.csawth,
physics_variables.eps,
- physics_variables.inductive_current_fraction,
+ physics_variables.f_c_plasma_inductive,
physics_variables.ejima_coeff,
physics_variables.kappa,
physics_variables.rmajor,
@@ -3173,7 +3177,7 @@ def phyaux(
@staticmethod
def plasma_ohmic_heating(
- inductive_current_fraction: float,
+ f_c_plasma_inductive: float,
kappa95: float,
plasma_current: float,
rmajor: float,
@@ -3186,7 +3190,7 @@ def plasma_ohmic_heating(
Calculate the ohmic heating power and related parameters.
Args:
- inductive_current_fraction (float): Fraction of plasma current driven inductively.
+ f_c_plasma_inductive (float): Fraction of plasma current driven inductively.
kappa95 (float): Plasma elongation at 95% surface.
plasma_current (float): Plasma current (A).
rmajor (float): Major radius (m).
@@ -3238,14 +3242,10 @@ def plasma_ohmic_heating(
error_handling.report_error(83)
# Ohmic heating power per unit volume
- # Corrected from: pden_plasma_ohmic_mw = (inductive_current_fraction*plasma_current)**2 * ...
+ # Corrected from: pden_plasma_ohmic_mw = (f_c_plasma_inductive*plasma_current)**2 * ...
pden_plasma_ohmic_mw = (
- inductive_current_fraction
- * plasma_current**2
- * res_plasma
- * 1.0e-6
- / vol_plasma
+ f_c_plasma_inductive * plasma_current**2 * res_plasma * 1.0e-6 / vol_plasma
)
# Total ohmic heating power
@@ -8063,8 +8063,8 @@ def init_physics_variables():
physics_variables.nd_impurities = 0.0
physics_variables.beta_poloidal_eps_max = 1.38
physics_variables.eps = 0.34399724802
- physics_variables.aux_current_fraction = 0.0
- physics_variables.inductive_current_fraction = 0.0
+ physics_variables.f_c_plasma_auxiliary = 0.0
+ physics_variables.f_c_plasma_inductive = 0.0
physics_variables.f_alpha_electron = 0.0
physics_variables.f_alpha_plasma = 0.95
physics_variables.f_alpha_ion = 0.0
@@ -8082,7 +8082,7 @@ def init_physics_variables():
physics_variables.f_tritium = 0.5
physics_variables.fusion_rate_density_total = 0.0
physics_variables.fusion_rate_density_plasma = 0.0
- physics_variables.fvsbrnni = 1.0
+ physics_variables.f_c_plasma_non_inductive = 1.0
physics_variables.ejima_coeff = 0.4
physics_variables.f_beta_alpha_beam_thermal = 0.0
physics_variables.hfac[:] = 0.0
diff --git a/process/power.py b/process/power.py
index 58d8d5f41a..8a9e6f2e87 100644
--- a/process/power.py
+++ b/process/power.py
@@ -441,7 +441,9 @@ def acpow(self, output: bool):
ppfmw = ppfmw + heat_transport_variables.peakmva
# Power to plasma heating supplies, MW
- pheatingmw = heat_transport_variables.pinjwp # Should be zero if ignite==1
+ pheatingmw = (
+ heat_transport_variables.p_hcd_electric_total_mw
+ ) # Should be zero if ignite==1
# Power to cryogenic comp. motors, MW
crymw = heat_transport_variables.crypmw
@@ -790,15 +792,15 @@ def power1(self):
# Secondary heat (some of it... rest calculated in POWER2)
# Wall plug injection power
# MDK
- # heat_transport_variables.pinjwp = (current_drive_variables.p_hcd_injected_total_mw + current_drive_variables.p_beam_orbit_loss_mw + physics_variables.p_fw_alpha_mw)/eta_hcd_primary_injector_wall_plug
- # heat_transport_variables.pinjwp calculated in current_drive.f90
+ # heat_transport_variables.p_hcd_electric_total_mw = (current_drive_variables.p_hcd_injected_total_mw + current_drive_variables.p_beam_orbit_loss_mw + physics_variables.p_fw_alpha_mw)/eta_hcd_primary_injector_wall_plug
+ # heat_transport_variables.p_hcd_electric_total_mw calculated in current_drive.f90
# Waste injection power
if physics_variables.ignite == 0:
# MDK
- # pinjht = heat_transport_variables.pinjwp - current_drive_variables.p_hcd_injected_total_mw - current_drive_variables.p_beam_orbit_loss_mw - physics_variables.p_fw_alpha_mw
+ # pinjht = heat_transport_variables.p_hcd_electric_total_mw - current_drive_variables.p_hcd_injected_total_mw - current_drive_variables.p_beam_orbit_loss_mw - physics_variables.p_fw_alpha_mw
heat_transport_variables.pinjht = (
- heat_transport_variables.pinjwp
+ heat_transport_variables.p_hcd_electric_total_mw
- current_drive_variables.p_hcd_injected_total_mw
)
else:
@@ -944,7 +946,7 @@ def power2(self, output: bool):
# Total recirculating power
heat_transport_variables.precircmw = (
self.pcoresystems
- + heat_transport_variables.pinjwp
+ + heat_transport_variables.p_hcd_electric_total_mw
+ heat_transport_variables.htpmw
)
@@ -1918,8 +1920,8 @@ def power2(self, output: bool):
po.ovarrf(
self.outfile,
"Electric power for heating and current drive (MW)",
- "(pinjwp)",
- heat_transport_variables.pinjwp,
+ "(p_hcd_electric_total_mw)",
+ heat_transport_variables.p_hcd_electric_total_mw,
"OP ",
)
po.ovarrf(
@@ -1971,7 +1973,7 @@ def power2(self, output: bool):
)
total_plant_power = (
heat_transport_variables.pnetelmw
- + heat_transport_variables.pinjwp
+ + heat_transport_variables.p_hcd_electric_total_mw
+ heat_transport_variables.htpmw
+ heat_transport_variables.vachtmw
+ heat_transport_variables.trithtmw
@@ -2179,7 +2181,7 @@ def power3(self, output: bool):
heat_transport_variables.pinjmax
/ current_drive_variables.eta_hcd_primary_injector_wall_plug
)
- p_hcd[3] = heat_transport_variables.pinjwp
+ p_hcd[3] = heat_transport_variables.p_hcd_electric_total_mw
p_hcd[4] = (
heat_transport_variables.pinjmax
/ current_drive_variables.eta_hcd_primary_injector_wall_plug
@@ -3103,8 +3105,8 @@ def init_heat_transport_variables():
heat_transport_variables.pgrossmw = 0.0
heat_transport_variables.pinjht = 0.0
heat_transport_variables.pinjmax = 120.0
- heat_transport_variables.pinjwp = 0.0
- heat_transport_variables.pinjwpfix = 0.0
+ heat_transport_variables.p_hcd_electric_total_mw = 0.0
+ heat_transport_variables.p_hcd_secondary_electric_mw = 0.0
heat_transport_variables.pnetelmw = 0.0
heat_transport_variables.precircmw = 0.0
heat_transport_variables.priheat = 0.0
diff --git a/process/pulse.py b/process/pulse.py
index a5dee23026..669dac2069 100644
--- a/process/pulse.py
+++ b/process/pulse.py
@@ -177,7 +177,7 @@ def burn(self, output: bool):
v_plasma_loop_burn = (
physics_variables.plasma_current
* physics_variables.res_plasma
- * physics_variables.inductive_current_fraction
+ * physics_variables.f_c_plasma_inductive
* physics_variables.csawth
)
diff --git a/process/stellarator.py b/process/stellarator.py
index 9f2a1ae178..d675e66655 100644
--- a/process/stellarator.py
+++ b/process/stellarator.py
@@ -4215,7 +4215,7 @@ def stphys(self, output):
# Calculate neutral beam slowing down effects
# If ignited, then ignore beam fusion effects
- if (current_drive_variables.pnbeam != 0.0e0) and (
+ if (current_drive_variables.p_hcd_beam_injected_total_mw != 0.0e0) and (
physics_variables.ignite == 0
):
(
@@ -4943,32 +4943,32 @@ def stheat(self, output: bool):
AEA FUS 172: Physics Assessment for the European Reactor Study
"""
if stellarator_variables.isthtr == 1:
- current_drive_variables.p_ecrh_injected_mw = (
+ current_drive_variables.p_hcd_ecrh_injected_total_mw = (
current_drive_variables.p_hcd_primary_extra_heat_mw
)
current_drive_variables.p_hcd_injected_ions_mw = 0
current_drive_variables.p_hcd_injected_electrons_mw = (
- current_drive_variables.p_ecrh_injected_mw
+ current_drive_variables.p_hcd_ecrh_injected_total_mw
)
current_drive_variables.eta_hcd_primary_injector_wall_plug = (
current_drive_variables.eta_ecrh_injector_wall_plug
)
- current_drive_variables.pinjwp = (
+ current_drive_variables.p_hcd_electric_total_mw = (
current_drive_variables.p_hcd_injected_ions_mw
+ current_drive_variables.p_hcd_injected_electrons_mw
) / current_drive_variables.eta_hcd_primary_injector_wall_plug
elif stellarator_variables.isthtr == 2:
- current_drive_variables.plhybd = (
+ current_drive_variables.p_hcd_lowhyb_injected_total_mw = (
current_drive_variables.p_hcd_primary_extra_heat_mw
)
current_drive_variables.p_hcd_injected_ions_mw = 0
current_drive_variables.p_hcd_injected_electrons_mw = (
- current_drive_variables.plhybd
+ current_drive_variables.p_hcd_lowhyb_injected_total_mw
)
current_drive_variables.eta_hcd_primary_injector_wall_plug = (
current_drive_variables.eta_lowhyb_injector_wall_plug
)
- current_drive_variables.pinjwp = (
+ current_drive_variables.p_hcd_electric_total_mw = (
current_drive_variables.p_hcd_injected_ions_mw
+ current_drive_variables.p_hcd_injected_electrons_mw
) / current_drive_variables.eta_hcd_primary_injector_wall_plug
@@ -4978,7 +4978,7 @@ def stheat(self, output: bool):
f_p_beam_injected_ions,
current_drive_variables.f_p_beam_shine_through,
) = self.current_drive.culnbi()
- current_drive_variables.pnbeam = (
+ current_drive_variables.p_hcd_beam_injected_total_mw = (
current_drive_variables.p_hcd_primary_extra_heat_mw
* (1 - current_drive_variables.f_p_beam_orbit_loss)
)
@@ -4987,15 +4987,17 @@ def stheat(self, output: bool):
* current_drive_variables.f_p_beam_orbit_loss
)
current_drive_variables.p_hcd_injected_ions_mw = (
- current_drive_variables.pnbeam * f_p_beam_injected_ions
+ current_drive_variables.p_hcd_beam_injected_total_mw
+ * f_p_beam_injected_ions
)
current_drive_variables.p_hcd_injected_electrons_mw = (
- current_drive_variables.pnbeam * (1 - f_p_beam_injected_ions)
+ current_drive_variables.p_hcd_beam_injected_total_mw
+ * (1 - f_p_beam_injected_ions)
)
current_drive_variables.eta_hcd_primary_injector_wall_plug = (
current_drive_variables.eta_beam_injector_wall_plug
)
- current_drive_variables.pinjwp = (
+ current_drive_variables.p_hcd_electric_total_mw = (
current_drive_variables.p_hcd_injected_ions_mw
+ current_drive_variables.p_hcd_injected_electrons_mw
) / current_drive_variables.eta_hcd_primary_injector_wall_plug
@@ -5013,10 +5015,10 @@ def stheat(self, output: bool):
# Calculate neutral beam current
- if abs(current_drive_variables.pnbeam) > 1e-8:
+ if abs(current_drive_variables.p_hcd_beam_injected_total_mw) > 1e-8:
current_drive_variables.c_beam_total = (
1e-3
- * (current_drive_variables.pnbeam * 1e6)
+ * (current_drive_variables.p_hcd_beam_injected_total_mw * 1e6)
/ current_drive_variables.e_beam_kev
)
else:
@@ -5071,7 +5073,7 @@ def stheat(self, output: bool):
current_drive_variables.bigq,
)
- if abs(current_drive_variables.pnbeam) > 1e-8:
+ if abs(current_drive_variables.p_hcd_beam_injected_total_mw) > 1e-8:
po.ovarre(
self.outfile,
"Neutral beam energy (KEV)",
diff --git a/scripts/create_dicts_config.py b/scripts/create_dicts_config.py
index 7bf42251ba..f602e942cd 100644
--- a/scripts/create_dicts_config.py
+++ b/scripts/create_dicts_config.py
@@ -1,2 +1,7 @@
# parameters that start with f, but are not f-values
-NON_F_VALUES = ["f_j_cs_start_pulse_end_flat_top", "fvsbrnni", "feffcd", "fcutfsu"]
+NON_F_VALUES = [
+ "f_j_cs_start_pulse_end_flat_top",
+ "f_c_plasma_non_inductive",
+ "feffcd",
+ "fcutfsu",
+]
diff --git a/source/fortran/current_drive_variables.f90 b/source/fortran/current_drive_variables.f90
index f1ff175c11..8b1f3b5c1a 100644
--- a/source/fortran/current_drive_variables.f90
+++ b/source/fortran/current_drive_variables.f90
@@ -78,14 +78,35 @@ module current_drive_variables
real(dp) :: f_c_plasma_diamagnetic
!! diamagnetic current fraction
- real(dp) :: p_ecrh_injected_mw
+ real(dp) :: p_hcd_ecrh_injected_total_mw
!! ECH power (MW)
- real(dp) :: echwpow
+ real(dp) :: p_ebw_injected_mw
+ !! Electron bernstein power (MW)
+
+ real(dp) :: p_hcd_ecrh_electric_mw
!! ECH wall plug power (MW)
+ real(dp) :: p_hcd_ebw_electric_mw
+ !! Electron bernstein wall plug power (MW)
+
real(dp) :: eta_cd_hcd_primary
- !! current drive efficiency (A/W)
+ !! Current drive efficiency of primary HCD system (A/W)
+
+ real(dp) :: eta_cd_hcd_secondary
+ !! Current drive efficiency of secondary HCD system (A/W)
+
+ real(dp) :: c_hcd_primary_driven
+ !! Current in plasma driven by primary HCD system (A)
+
+ real(dp) :: c_hcd_secondary_driven
+ !! Current in plasma driven by secondary HCD system (A)
+
+ real(dp) :: f_c_plasma_hcd_primary
+ !! Fraction of plasma current driven by primary HCD system
+
+ real(dp) :: f_c_plasma_hcd_secondary
+ !! Fraction of plasma current driven by secondary HCD system
real(dp) :: n_ecrh_harmonic
!! cyclotron harmonic frequency number, used in cut-off function
@@ -111,6 +132,12 @@ module current_drive_variables
real(dp) :: eta_lowhyb_injector_wall_plug
!! lower hybrid wall plug to injector efficiency
+ real(dp) :: eta_icrh_injector_wall_plug
+ !! Ion cyclotron wall plug to injector efficiency
+
+ real(dp) :: eta_ebw_injector_wall_plug
+ !! Electron bernstein wave wall plug to injector efficiency
+
real(dp) :: eta_beam_injector_wall_plug
!! neutral beam wall plug to injector efficiency
@@ -214,16 +241,31 @@ module current_drive_variables
real(dp) :: p_hcd_injected_total_mw
!! total auxiliary injected power (MW)
+ real(dp) :: p_hcd_injected_current_total_mw
+ !! total auxiliary injected power (MW)
+
real(dp) :: p_hcd_secondary_injected_mw
!! secondary total fixed auxiliary injected power (MW)
+ real(dp) :: p_hcd_primary_injected_mw
+ !! primary auxiliary injected power (MW)
+
real(dp) :: f_c_plasma_internal
!! plasma current fraction driven internally (Bootstrap + Diamagnetic + PS)
- real(dp) :: plhybd
- !! lower hybrid injection power (MW)
+ real(dp) :: p_hcd_lowhyb_injected_total_mw
+ !! Total lower hybrid injection power (MW)
+
+ real(dp) :: p_hcd_icrh_injected_total_mw
+ !! Total ion cyclotron injection power (MW)
+
+ real(dp) :: p_hcd_ebw_injected_total_mw
+ !! Total electron bernstein wave injection power (MW)
+
+ real(dp) :: p_beam_plasma_coupled_mw
+ !! Total neutral beam power that is coupled to plasma after losses (MW)
- real(dp) :: pnbeam
+ real(dp) :: p_hcd_beam_injected_total_mw
!! neutral beam injection power (MW)
real(dp) :: p_beam_orbit_loss_mw
@@ -232,7 +274,7 @@ module current_drive_variables
real(dp) :: f_c_plasma_pfirsch_schluter
!! Pfirsch-Schlüter current fraction
- real(dp) :: pwplh
+ real(dp) :: p_hcd_lowhyb_electric_mw
!! lower hybrid wall plug power (MW)
real(dp) :: pwpnb
diff --git a/source/fortran/heat_transport_variables.f90 b/source/fortran/heat_transport_variables.f90
index 4406089039..9666ad385b 100644
--- a/source/fortran/heat_transport_variables.f90
+++ b/source/fortran/heat_transport_variables.f90
@@ -133,11 +133,14 @@ module heat_transport_variables
real(dp) :: pinjmax
!! maximum injector power during pulse (heating and ramp-up/down phase) (MW)
- real(dp) :: pinjwp
+ real(dp) :: p_hcd_electric_total_mw
!! injector wall plug power (MW)
- real(dp) :: pinjwpfix
- !! secondary injector wall plug power (MW)
+ real(dp) :: p_hcd_secondary_electric_mw
+ !! Secondary HCD system injector wall plug power (MW)
+
+ real(dp) :: p_hcd_primary_electric_mw
+ !! Primary HCD system injector wall plug power (MW)
real(dp) :: pnetelmw
!! net electric power (MW)
diff --git a/source/fortran/numerics.f90 b/source/fortran/numerics.f90
index b99323775e..d940e5fe39 100755
--- a/source/fortran/numerics.f90
+++ b/source/fortran/numerics.f90
@@ -249,7 +249,7 @@ module numerics
!!