Offer an alternative calculation for the normalised internal inductance, $`l_i`$ #1878

Open jonmaddock opened 1 year ago

In GitLab by @mcoleman on Jun 20, 2023, 15:29

Summary

At present, PROCESS exclusively uses an empirical fit for $l_i$, found in e.g. ("Tokamaks" 4th Edition, Wesson, page 116), and ("Fusion" 2nd edition, Stacey, page 53). In the latter, one can see that the data come from DIII-D in the 1990's. I have no idea what kind of plasmas they were running back then, but it's fair to say that this fit gives high values for $l_i$ compared to what we typically see in larger, higher power, future tokamak reactors.

For EU-DEMO machines, PROCESS typically gives $l_i > 1.15$, compared to values from transport solvers of 0.7-0.9.

Obviously, PROCESS' profiles are not particularly to be trusted, so any $\rho$-based calculation will fall over pretty quickly. Instead, I propose an energy-based calculation, which is akin to ITER's $l_{i}(3)$ definition:

W_{p} = \int_{V}\dfrac{B_{\theta}^2}{2\mu_{0}}dV = \dfrac{1}{2}L_{i}^{W}I_{p}^2 \\
L_{i}^{W} = \dfrac{2W_{p}}{I_{p}^2} = \dfrac{l_{i}\mu_{0}R_{0}}{2} \\
l_{i} = \dfrac{4\int_{V}\dfrac{B_{\theta}^2}{2\mu_{0}}dV}{\mu_{0}R_{0}I_{p}^2} = \dfrac{2\langle B_{\theta}\rangle^2 V}{\mu_{0}^2R_{0}I_{p}^2}

Noting that:

\langle B_{\theta}\rangle^2 \textrm{not technically but for all intents and purposes (?) equal} \langle B_{\theta}^2\rangle

rli = 2 * bp**2 * vol / (rmu0**2 * rmajor * plascur**2)

Of course, if $\langle B_{\theta}\rangle = \mu_{0}I_{p}/P$, where $P$ is the perimeter of the LCFS, as appears to be the case most of the time in PROCESS, then this simplifies down to a purely geometric relation:

l_i = \dfrac{2V}{R_{0}P^2}

For a typical EU-DEMO, using this I get $l_i \approx$ 0.8.

Given that this approach to calculating $l_i$ is more-or-less what ITER is using, and captures at least some of the important terms (with the notable exception of $q$), I think it should be at least offered as an alternative to the Wesson/Stacey empirical fit.

Checklist

After implementing issue do the following

[ ] Run test suite
[ ] Does it change the output? If so, give justifications.
[ ] Does it change the baseline 2019 run? (you may need to run on Freia)
[ ] Merge request created

In GitLab by @mcoleman on Jun 20, 2023, 15:37

I've just realised I'm not allowed to push code to this repository.

In GitLab by @mcoleman on Jun 21, 2023, 08:14

In GitLab by @mkovari on Jun 21, 2023, 14:57

I wouldn't rush into this ideally. @ajpearcey @mcoleman

There is already an option for the user to specify $l_i$ (providing iprofile=0).

More generally, $l_i$ will depend on the current distribution that is required for the desired scenario.

The PROCESS formula relates $l_i$ to the q profile via alphaj, the current density profile index:

    if (iprofile == 1) then
       alphaj = qstar/q0 - 1.0D0
       rli = log(1.65D0 + 0.89D0*alphaj)  !  Tokamaks 4th Edition, Wesson, page 116
    end if

(As far as I know alphaj is not used anywhere else.) The value of q0 is an input (although it is usually set to 1.0) and the lower bound for q95 (boundl(18)) is usually input, so we sort of want these to be consistent with $l_i$.

I certainly couldn't tell you whether the formulae in the code are any good.

In GitLab by @mcoleman on Jun 22, 2023, 08:03

I mean the above calculation is totally flawed.

\langle B_{\theta}\rangle^2 \textrm{is definitely not equal to} \langle B_{\theta}^2\rangle

The q dependency is hidden by that (in)equality.

My problem is that to calculate $l_i$ one needs to have good profile information, which PROCESS does not really have.

I'm happy to sit on this until I can think of a better way of coming up with a calculation of $\langle B_{\theta}^2\rangle$, which I can imagine being possible from the $q$ profile, which is perhaps the only profile we could imagine trusting.

All I can say is that the present empirical fit is not fit for EU-DEMO purposes.

I've also got fits to ASTRA data for $q_{95}$ and plasma volume, but I don't have enough data to make a serious one.

In GitLab by @mcoleman on Jun 28, 2023, 13:49

Alright, sorry for the delay, I've been meaning to write this up for a few days.

My faith in PROCESS profiles is basically zero, so I have been trying to get $l_i$ exclusively in 0-D. Also, the PROCESS plasma geometry parameterisations (with circle / ellipse(?) arcs) are a little inconvenient, so in terms of internal inductance there is not much solid stuff to work with.

Note that if you're running PLASMOD, there are plenty of easy and convenient 1-D integrals that could be added to calculate $l_i$ self-consistently there. My assumption here is that when running PROCESS without PLASMOD, one will want a decent guess as to the internal inductance of the plasma.

If my above flawed calculation didn't convince you because of a lack of explicit $q$-dependence, let me offer another:

l_i \equiv \dfrac{\big\langle B_{\theta}^2\big\rangle}{\bigg(\dfrac{\mu_{0}I_{p}}{\oint dl}\bigg)^2}

I found a cute approximation to $\big\langle B_{\theta}^2\big\rangle$ in Freidberg's "Plasma Physics and Fusion Energy" (13.181, pp 417), which I think is valid for large aspect ratio:

\big\langle B_{\theta}^2\big\rangle = \dfrac{\pi^2}{32}\dfrac{(\kappa^2+1)^2}{\kappa^2+4/3}\dfrac{B_{T,0}^2}{A^2q_{*}^2}

If you take (Freidberg, 13.159, pp 406):

q_{*} = \dfrac{2\pi a^2 \kappa B_{T,0}}{\mu_{0} R_{0} I_{p}}\dfrac{1+\dfrac{4}{\pi^2}(\kappa^2-1)}{\kappa}

Noting that $q_{*}$ not equal $q_{cyl}$ which seems to be assumed in PROCESS.

Which leaves me with:

l_i = \dfrac{1}{128}\dfrac{(\kappa^2+1)^2}{(\kappa^2+4/3)(1+\dfrac{4}{\pi^2}(\kappa^2-1))^2}\dfrac{L_{p}^2}{a^2}

which is a little inconvenient as it depends on minor radius, as opposed to aspect ratio. You can of course take different approximations for $q_{*}$, and play around with $\kappa$ or $\kappa_{95}$

Shown as li proposal(2) here (where I have chosen to use $\kappa_{95}$):

Keen to hear thoughts on this. Is it possible there is some definition issue with the Wesson fit? It's about a factor 2/3 off, 1/2 if you're being optimistic.

Matti, Sorry, I didn't see your comments as I wasn't tagged.

We seem to be going round in circles. The formula for $q_{*}$ (whatever that is - I have to admit I did think it was the same as $q_{cyl}$) still doesn't say anything about the current profile - which is the whole point of $l_i$, as far as I can see.

If most of the current flows near the edge, then $B_{\theta}$ will be small over much of the volume and $\langle B_{\theta}^2\big\rangle$ will be small. This won't affect $q_{95}$, and presumably $q_{*}$, so it won't be included in your formula.

I am prepared to open Freidberg if you insist, but why change the habit of a lifetime?

I have to say to I preferred the Gitlab rendering - Githib shows the equations very black.

Hi @mkovari, @mcoleman here.

The primary purpose of $l_i$ in PROCESS from an EU-DEMO perspective is its use in the calculation of the pulse duration. We do not use the PROCESS profiles beyond the initial design point.

The primary motivation to change the calculation of $l_i$ is because it is consistently a factor 1.4--1.6 higher than more detailed estimates. The secondary motivation is that the existing formula is a fit to data from the late 80s and early 90s (IIRC) - of questionable validity given the regimes we now consider for power reactors.

This fit, incidentally, makes use of $q^{*}$ as the sole fitting parameter rather than the current profile. The fact that neither the present method, nor any of the ones I proposed above make use of the current profile is rather moot, is it not?

All this said, there is a way to run PROCESS with a fixed internal inductance, so it's a take it or leave it job as you suggest. It's just a shame because when we perform scans where the pulse length is of interest, we're ignoring any influence on internal inductance that would vary with some of the scanned parameters.

I have my github on in dark mode and I see the equations as white on black, as with the rest of the text.

Hi @mcoleman, I am beginning to understand your point now. I didn't appreciate that the main use of $l_i$ in PROCESS is to calculate the pulse length. In EU-DEMO (according to PROCESS) the internal inductance is about 45% of the total inductance, so it will have a signficant effect on the pulse length.

Referring to Freidberg's formula 13.181 (p. 417). It is based on the large aspect ratio (as you say). It also assumes a particular current profile:

A simple model in which the current density increases with major radius, as it does for realistic profiles...

and the equation for the vector potential is indeed expressed in terms of the major radius. This is odd, since I thought the current density was a flux function, so it depends on the flux (or, loosely, the minor radius) not the major radius.

The formula also assumes a particular plasma shape - possibly the half-ellipse shown on p. 415.

Freidberg introduces a new definition of $q_{*}$ in 13.171 (p.414). Since the equation for $l_i$ contains $q_{*}^2$, this is not a trivial matter, so we would need to make sure we use the right one.

The second term in $q^*$ is a shape factor, which is usually close to 1 for normal elongations, and similar in nature to a lot of the approximations for various integrals used in 0-D plasma physics. Or at least they are ubiquitous in Freidberg...

In the first definition of $q^*$ (13.160), there is the approximation that $G(\kappa)\approx1$, which the second term in 13.171 attempts to correct without resorting to an elliptical integral.

I vaguely remember trying out this second term, 1, and the actual term with elliptical integral. The results varied, but not that significantly. This was a while ago now, so I'm afraid I have lost my workings.

Not really my field, but I don't think plasma current density is really a flux function, or better said a function solely of flux. It might be a matter of terminology, but the current density does depend on the major radius.

@mcoleman. I never claimed to know much about this, but it is embarassing to discover my knowledge is actually negative. So the plots that people show of current density versus flux are actually $j/R$ or something?

The formula we use for the external inductance seems to be fairly soundly based*, given an elliptical cross-section. To prevent integrating over portions of the same volume twice it would be nice to use a formula for internal inductance that is based on the same shape - not the case for Freidberg if he has indeed used the D-shape in Figure 13.37.

However, I agree that we shouldn't be using a formula that overestimates by 50%. I also agree that it is good to use a calculation that depends in a sensible way on aspect ratio and elongation, and that the current profile is less important at this stage.

The next step (!?) is to explain why there are three or more different definitions of $l_i$ . The definitions on p. 281 of Freidberg look unique. It would be nice to make sure we are using the same definition as Hirshman and Neilson do in their calculation of external inductance.

What I am saying is that possibly we should use the Freidberg formula, but ideally we should check our values of both internal and external inductance against those from detailed calculations. Are you able to do this?

S.P. Hirshman, G.H. Neilson, External inductance of an axisymmetric plasma, doi:10.1063/1.865934, eq 14a, 15a and 15b.

Hi @mkovari ,

All I can say re current density "profile" is that, as far as they are in G-S solvers, generally $j_{tor} = 2\pi r p^{'} + \dfrac{FF^{'}}{\mu_{0}r}$ or some similar formulation, but with a r + 1/r relation, with r being in cylindrical coordinates. Meaning that current density is not solely a function of flux. I don't know much more than that re terminology etc.

Some time ago we had a discussion with CREATE at EUROfusion regarding external inductance. The Hirshman and Neilson calculation is totally fine. There are two approaches of calculating achievable pulse length, one that uses external inductance, and one that does not. In both these approaches, the definitions of internal inductance are the same, so there are no issues.

CREATE believe, and I agree, that it is simpler and easier to calculate the achievable pulse length without going down the external inductance calculation route. That said, when I looked at this in PROCESS, and given the time constraints we were operating under, it became clear that there was no point in changing anything (the results were basically identical).

The multiple definitions of dimensionless internal inductance, $l_i$ are just different ways of approximating the same thing, with corrections for things like the plasma circumference differing slightly. All I know is that ITER use li(3), and that seems to be the accepted approach. I don't think they differ heavily, unless of course you always assume you plasma has a circular cross-section when it doesn't.

Surely what matters for the inductive flux consumption during ramp-up is the total plasma inductance, and this is calculated as internal + external?

@mkovari, long story short, I think you can ignore this. As far as I can tell, the PROCESS calculations for values associated with the pulse length are remarkably good, even if the internal inductance does not match up with what more detailed codes say. There may be some difference in definition of "normalised internal".

Apologies for probably having created entropy.

OK, I will close the issue. It's a pity we can't sort this out in a vaguely self-consistent way, but hey.

The PhD thesis of Tobias Hartmann has helpful info. and then we get a

Since we very much do not have a circular cross-section, we really want to get rid of the Process formula: alphaj = qstar / q0 - 1.0

Also, we usually use a pedestal, so the plasma current density profile does not have the "parabolic" form. Indeed, there may be a large amount of bootstrap current flowing in the pedestal, reducing $l_i$. @CoronelBuendia

There is another formula in Fusion Plasma Physics (Stacey, 2nd ed, eq 19.7, p. 589), although no derivation or justification is given:

l_i=ln(1.65+0.89(q_{95}-1))

This contrasts oddly with the PROCESS formula: rli = np.log(1.65 + 0.89 * alphaj). They can't both be right.

The Stacey formula is not correct. For a very reasonable value of $q_{95}=5$, and using $g=4 l_i$ (p. 584 of Fusion Plasma Physics, and PROCESS), we get a limit on normalised beta $g=6.6$ - unreasonably large for a conventional aspect ratio.

Hi @mkovari,

On leave, so I don't have access to papers etc., but as far as I can tell, they are effectively equivalent.

$\alpha_j = q^{\textrm{*}}/q_0 - 1$

$q^{\textrm*}$, $q{cyl}$, and $q{95}$ are all basically approximate measures of the safety factor near the edge of the plasma, pretending there isn't an X-point. As in some of the discussion above, different choices of metric here do make a difference.

If you take the safety factor on axis, $q_0$, to be 1.0, which is the case for convential plasmas in conventional aspect ratio machines (and probably a default value in PROCESS), then $q0$ disappears, and the only discrepancy is in the choice of measure for the safety factor near the edge; $q{95}$ or $q^{\textrm}$. $q_{95}$ is typically taken as > 3.0 in conventional aspect ratio reactors, and in ITER (not 1.0 as you write above). $q^{\textrm}$ is likely a bit lower than $q_{95}$, but still above 2.5 or so I'd guess. Does that help?

Sorry - I meant to write 5, not 1. As I understand it, $q_{95}$ is strongly influenced by the current density near the edge, and can be very high. Obviously the 'parabolic' density profile does not capture this.

My point is that since everybody always quotes $q_{95}$, and never $\alphaj$ which is not a real quantity but a ropey parameter, Stacey's formula looked useful. However, it is clearly wrong since $q{95}$ can be very high and normalised beta cannot be.

ukaea / PROCESS

Offer an alternative calculation for the normalised internal inductance, $`l_i`$ #1878

Summary

Checklist