Olivier Sauter's formulae including negative triangularity

In GitLab by @Hlux on Jun 2, 2016, 10:42

Please include Oliviers new formula in PROCESS, so that we can run comparisons with M. Kikuchi's suggested reactor version.

In GitLab by @rkemp on Jun 2, 2016, 11:41

Added using icurr=8 and tested. Negative triangularity for other current scalings will give an error and stop PROCESS. Note that at this stage I have only implemented the current scaling and not the other geometrical factors for poloidal length and plasma volume, for example, as these were not very wrong. There may be other issues in the code, e.g. TF inductance, with negative triangularity so use with care...

In GitLab by @Hlux on Jun 2, 2016, 12:14

Thanks,
Hanni

In GitLab by @Hlux on May 16, 2018, 09:53

Documentation does not include this model. Possibly it should be checked which other models are affected by negative triangularity and whether this is valid.

In GitLab by @Hlux on May 16, 2018, 09:53

reopened

In GitLab by @Hlux on May 16, 2018, 09:57

Check that the new shape fits inside the TF coils!

In GitLab by @rkemp on May 16, 2018, 09:58

Reassigning due to no longer working on this task. Other things to check:

Plasma chamber / TF coil geometry (I think I did implement this -- it should be sufficient to run PROCESS and plot cross-section).
Any post hoc calculations in the summary sheet that may assume particular shaping.
Relationships between k95/kx and d95/dx (should be described in the paper), although I think these were okay.

Paper reference: https://www.sciencedirect.com/science/article/pii/S0920379616303234
Sauter, FED Volume 112, 2016, Pages 633-645

In GitLab by @Hlux on May 16, 2018, 10:17

Attaching Kikuchi presentation for complete reference. NegativeTriangularityKikuchi.pptx

In GitLab by @jmorris-uk on Jul 17, 2019, 18:03

removed milestone

The reactor parameters for Kikuchi's negative triangularity configuration, using Sauter's new formulae valid for negative delta, are given on pp 40-41 of the link above.

icurr=8 Sauter scaling allowing negative triangularity is listed in Variable Descriptions.

Calculation of the plasma current using the safety factor requires the factor fq:

        elif (
            icurr == 8
        ):  # Sauter scaling allowing negative triangularity [FED May 2016]
            # Assumes zero squareness, note takes kappa, delta at separatrix not _95

            w07 = 1.0  # zero squareness - can be modified later if required

            fq = (
                (1.0 + 1.2 * (kappa - 1.0) + 0.56 * (kappa - 1.0) ** 2)
                * (1.0 + 0.09 * triang + 0.16 * triang**2)
                * (1.0 + 0.45 * triang * eps)
                / (1.0 - 0.74 * eps)
                * (1.0 + 0.55 * (w07 - 1.0))
            )
and later
        if icurr == 8:
            curhat = 4.1e6 * rminor**2 / (rmajor * qpsi) * fq

Sauter:

where the geometric parameters are defined at the last closed flux surface.

The new parameter $w_{07}$ is related to the radial width of the plasma shape at 70% of the maximum height, and takes into account the squareness present in most plasma shapes, in particular single-null ones. It is defined at the bottom of page 4 and top of page 5. The problem is that we don't know its value most of the time. Richard's assumption that $w_{07}=1$ is reasonable.

The PROCESS docs give:

but Richard has failed to use this definition of fq, since

$$\frac{2\pi}{\mu_0}=5\times10^6 \text{A T}^{-1} \text{m}^{-1}$$

qpsi : input real : plasma edge safety factor (= q-bar for icurr=2)
In the calling function we find physics_variables.q. This is passed to culcur where it is named qpsi for some reason. This is consistent with the documentation:

q	Input	real	3.0	Safety factor 'near' plasma edge (iteration variable 18) equal to q95 (unless icurr=2 (ST current scaling), in which case q = mean edge safety factor qbar)

Update docstring in culcur.
Rename qpsi in culcur as q95.
Eliminate curhat since it's not necessary.
Test for a range of values of triangularity $\delta$.

I have now included all the geometric formulae.

When the triangularity is negative, the vacuum vessel may clash with the TF coil, not surprisingly:

I have not attempted to fix this.

Unfortunately, Kikuchi's presentation gives ambiguous values of elongation or triangularity, so I can't use it for benchmarking. I won't be able to do a quantitative check on the geometric formulae.

I will try a qualitative check, using two runs with triangularity of opposite sign.

The poloidal perimeter is the same, because the formula contains only $\delta^2$.
The cross-section is also the same, because the formula doesn't include $\delta$ at all.
The surface area and volume are bigger for negative $\delta$, because the mean radius is bigger - the plasma is further from the axis on average.

This seems OK.

		Negative	Positive
Major radius (m)	(rmajor)	8	8	ITV
Minor radius (m)	(rminor)	2.667	2.667	OP
Aspect ratio	(aspect)	3	3
Elongation, X-point (input value used)	(kappa)	1.85	1.85	IP
Elongation, 95% surface (calculated from kappa)	(kappa95)	1.652	1.652	OP
Elongation, area ratio calc.	(kappaa)	1.85	1.85	OP
Triangularity, X-point (input value used)	(triang)	-0.5	0.5	IP
Triangularity, 95% surface (calculated from triang)	(triang95)	-0.333	0.333	OP
Plasma poloidal perimeter (m)	(pperim)	25.08	25.08	OP
Plasma cross-sectional area (m2)	(xarea)	41.329	41.329	OP
Plasma surface area (m2)	(sarea)	1.33E+03	1.19E+03	OP
Plasma volume (m3)	(vol)	2.16E+03	1.99E+03	OP

triang_minus0.5_icurr8_IN.DAT.txt
triang_minus0.5_icurr8_MFILE.DAT.txt
triang_minus0.5_icurr8_OUT.DAT.txt
triang_plus0.5_icurr8_IN.DAT.txt
triang_plus0.5_icurr8_MFILE.DAT.txt
triang_plus0.5_icurr8_OUT.DAT.txt

Sauter also gives a formula for the trapped fraction, using ref 16 . Although it looks simple, the reference doesn't seem to make sense. The point corresponding to the plasma edge ($\epsilon=0.32$) in Fig 6(a) gives values of trapped fraction which don't agree with the fitting formula (0.746 for $\delta=0.5$ and 0.824 for $\delta=-0.5$).

I have not implemented the formula for trapped fraction.

I have used Sauter's formulae only for the edge parameters. I think they should work for every flux surface, but one would need to know the values of $\kappa, \delta$ and $\epsilon$ at each flux surface, which we don't have in Process.