Add option to use PLASMOD plasma resistance calculation #1616

Closed jonmaddock closed 1 month ago

In GitLab by @ajpearcey on Mar 23, 2022, 17:04

add switch to allow for use of PLASMOD plasma resistance calculation instead of Uckan scaling relation.

In GitLab by @ajpearcey on Mar 24, 2022, 14:44

To include this we can use the temperature profile to make the profile for the inductive current

\sigma(\rho) = 601.2 \frac{T^{3/2}(\rho) (0.76 + Z_{\mathrm{eff}})}{Z_{\mathrm{eff}}(1.18+0.58Z_{\mathrm{eff}})\ln\Lambda(\rho)}

Using this we can compute the parallel electric field as

E_{\parallel} = \frac{2\pi I_{\mathrm{ind}}}{F(\psi)B_T \int \frac{\sigma(\rho)}{F(\psi)^2} dV}

and the plasma resistance as

R_{\mathrm{plas}} = \frac{E_{\parallel} 2 \pi B_T}{I_{\mathrm{ind}}F(0)G_3(0)}

Here we are being imprecise and mixing our $\rho$'s and $\psi$'s so we will make the assumption that we can let $F(\psi)$ be treated as only the constant vacuum field contribution $F=B_{T}R_{\mathrm{major}}$. PLASMOD also includes an initial guess for $G_3(\psi)$ and if we take only the leading contribution it reduces too $G_3(0) = R_{\mathrm{major}}^{-2}$.

I have been able to find solution close to the DEMO 2018 baseline using both these new plasma resistivity calculation and implementing the same impurity profiles as used in PLASMOD (see issue #1486). The new solution found by PROCESS with these models gives a much lower loop voltage, lower even than the full PLASMOD calculation. The PROCESS 2018 baseline had a V_loop = 46mV, while the 2019 baseline that uses PLASMOD found V_loop = 16mV, while the new PROCESS solution using these new calculation find V_loop = 6.5mV. For example see below I have a PROCESS solution which scans the required burn time from 2 to 4 hours.

I have spent time trying to understand the source of the discrepancies in the calculations, and I believe it is due to the different profiles assumed in each code. PROCESS assumes some parameterised profile shapes with temperature at the pedestal and separatrix fixed and allows the volume averaged temperature to vary to find a feasible solution and also assumes T_e = T_i. I have made a plot to compare the temperature profiles created for both PROCESS and PLASMOD, shown below.

In GitLab by @ajpearcey on Mar 24, 2022, 15:34

Here are the input and output files used to make the scan plotted above: baseline_2018_IN.DAT baseline_2018_OUT.DAT

In GitLab by @mkovari on Mar 25, 2022, 11:56

The dependence on Zeff is different from the original PROCESS dependence, which is very primitive. Do you have a reference?

In GitLab by @mkovari on Mar 25, 2022, 12:06

Also - are you taking into account the neo-classical resistivity enhancement caused by particle trapping (the factor in R/a)? This is pretty significant.

In GitLab by @mkovari on Mar 25, 2022, 12:54

Could I simplify your algebra a little to make things easier for me, using your assumptions for F and G3.

R_{plas}=\frac{4\pi^2 R_{major}^2}{\int{\sigma}{dV}}

V_{loop}=I_{ind}R_{plas}

In GitLab by @mkovari on Mar 25, 2022, 13:56

A more complicated formula for the neoclassical resistivity is given in Wesson.
Neoclassical_resistivity_Wesson.pdf

Note that this is resistivity not resistance and presumably the correct aspect ratio to use is the aspect ratio of the individual flux surface, which tends to infinity at the plasma centre. (The inverse aspect ratio $\epsilon$ tends to zero.)

In GitLab by @mkovari on Mar 25, 2022, 14:06

Do you have any published values for loop voltage?

In GitLab by @ajpearcey on Apr 12, 2022, 11:04

Note this change needs investigating with the PROCESS implementation.

Dear All,
I just discovered a bug in the calculation of the loop voltage in PLASMOD. There was a missing dimensionless geometrical factor in the calculation, 
which makes the ouput basically a factor 2 times larger. that is, the actual loop voltage was underestimated by roughly a factor 2.
The incriminated lines are:

in equil.f90
the lines above are the corrected ones, while in the git version those lines have the missing geometrical terms at line 277 and 280.
Please change the equil.f90 according to these lines and push it (I am at the moment out of the possibility to do that).
Regards and sorry for the inconvenience,
Emiliano

We have now removed PLASMOD from PROCESS therefore I do not think this is so relevant anymore. I do think looking into more detailed plasma physics models is a good idea, but that is beyond the scope of this issue.

ukaea / PROCESS

Add option to use PLASMOD plasma resistance calculation #1616