Central Solenoid current calculation #98

Closed jonmaddock closed 1 year ago

In GitLab by @mkovari on Aug 15, 2014, 13:51

The central solenoid routine pfjalw is obsolete and gives unphysical and erratic results over 14.5 Tesla.
(The file K:\Power Plant Physics and Technology\PROCESS\Solver studies\OH coil.pdf shows this.) This routine should be replaced by the standard ITER NB3Sn routine, which should be checked to make sure that it does not give unphysical results outside its range of validity.

In GitLab by @pknight on Oct 13, 2014, 13:50

This should be straightforward. I've spotted that the cross-sectional area of supporting steel for the CS is calculated as

areaspf = 0.666D0 * forcepf / 5.0D8

instead of the newer method we introduced for the other PF coils:

areaspf = sigpfcf * forcepf / (sigpfcalw*1.0D6)

Should I make this consistent?

In GitLab by @pknight on Oct 14, 2014, 10:58

Following on from our discussion yesterday, the supporting steel issue appears to be more difficult to solve than I thought.

For the non-CS PF coils I think we can take the cross-sectional dimensions of the 'coil' as those of the winding pack. Then everything can be calculated in the existing order okay. The steel may be assumed to surround the winding pack, as its thickness (calculable) does not affect any build calculations.

However, the CS may have to be treated differently, as the coil cross-sectional dimensions are effectively inputs; thus the winding pack dimensions need to take into account the steel area fraction. The forcepf calculation above uses the magnetic field at the inboard and outboard edges, and forcepf determines areaspf, which in turn affects the winding pack dimensions via the steel area fraction, so there is a circular calculation needed. I propose that we introduce a new consistency equation to solve this.

In GitLab by @mkovari on Oct 15, 2014, 10:20

Peter,

We have the following constraints, which are based on the superconductor critical current (I think):
• (26) OH coil EOF current density upper limit
• (27) * OH coil BOP current density upper limit

I think these should be sufficient, without introducing a new constraint. The calculation sequence would be as follows, which is much as it is now:

Inputs:
ohcth /0.811/ : central solenoid thickness (m) (iteration variable 16)
bore /1.42/ : central solenoid inboard radius (m) (iteration variable 29)
coheof /1.85e7/ : central solenoid overall current density at end of flat-top ( A/m2) (iteration variable 37),
csutf /1.32e9/ : ultimate strength of case (Pa) (default value from DDD11-2 v2 2 (2009)) ALSO for CS
csytf /1.0005e9/ : yield strength of case (Pa) (default value from DDD11-2 v2 2 (2009)) ALSO for CS
tftmp /4.5/ : peak TF coil helium coolant temperature (K) ALSO for CS

Calculated:
total coil cross-section
alstrtf : allowable stress in TF coil (Pa) ALSO in CS.
Allowable hoop stress (=alstrtf?)
current = total coil cross-section coheof = I
magnetic field
forcepf = hoop force
areaspf = steel cross-section = forcepf / alstrtf
helium cross-section
conductor cross-section = total coil cross-section - helium cross-section – areaspf
conductor current density j = current / conductor cross-section
Superconductor critical current density (function of field, temperature, strain=0)
The conductor current density is then limited by the constraint.

The stress in the case would always be limiting (as now). The strain in the strand would be ignored, which makes sense as it should be possible to wind the coil to eliminate strain.

*The remaining issue is what happens when the conductor cross-section comes out negative. I suggest:

δ = 0.0001  
If conductor cross-section a > δ  
     conductor current density j = I/a  
Else If conductor cross-section a < δ  
     conductor current density j = I(2δ-a)/ δ^2  
End if

This creates a continuous and smooth function that increases monotonically as conductor cross-section decreases.

Have I missed something?

In GitLab by @pknight on Oct 16, 2014, 15:14

The PF and CS superconductor critical current densities are now calculated using the same routines as for the TF coils; commit 3592c2a6c (r350). I have introduced new inputs isumatoh, and strand copper fractions fcupfsu, fcuohsu (in fact this is just fcuoh renamed for consistency). However, I'm finding it very difficult (actually I haven't succeeded at all yet) to get a successful run with isumatoh = 3 (i.e. NbTi), as the field at the coil is higher than bcrit for NbTi, so the default value is isumatoh=1 (Nb3Sn) for now.

The coil steel case areas (and equivalent thicknesses) are now calculated as requested. For the non-CS coils the cross-sectional dimensions of the coils given in the output files are for the winding pack - the cases are assumed to extend outwards (which may very occasionally lead to overlaps in the build if cases are too thick, especially for ipfloc(i)=1 coils), but for the CS coil the winding pack is smaller than the coil dimensions given, and the case fills the rest of the space. I believe that everything is calculated self-consistently.

In GitLab by @pknight on Oct 16, 2014, 15:32

New output file with default input values OUT.DAT_git350

In GitLab by @pknight on Oct 16, 2014, 15:32

New DEMO1 output file OUT.DAT_git350

In GitLab by @pknight on Oct 16, 2014, 15:32

New DEMO2 output file OUT.DAT_git350

In GitLab by @pknight on Oct 16, 2014, 15:33

New inductive ITER output file OUT.DAT_git350

In GitLab by @pknight on Oct 16, 2014, 15:33

New ITER98 case with minimal constraints output file OUT.DAT_git350

In GitLab by @mkovari on Oct 17, 2014, 08:41

Peter, Many thanks for this. I haven't looked at the code yet, but I have a few minor comments on the output.

Global_variables.f90 has

!+ad_vars cohbop : central solenoid overall current density at beginning of pulse (A/m2)
!+ad_vars coheof /1.85e7/ : central solenoid overall current density at end of flat-top (A/m2)

This doesn't seem consistent with the description in OUT.DAT:

Actual winding pack current density at EOF (A/m2) (coheof) 7.092E+06
Actual winding pack current density at BOP (A/m2) (cohbop) 7.092E+06

I would prefer if "OH coil" was replaced with "Central Solenoid" throughout.
The header in OUT.DAT would be more accurate if it reads:

***** Central Solenoid and PF Coils **
The CS doesn't really have a "winding pack", so we shouldn't really use this term. Nor does it have a "case." It has steel, and superconductor strand.
The OUT.DAT should give the hoop stress in the steel of the CS (which will equal the maximum permissible stress, I think).
The header

OH Coil Stress Calculations

isn't really accurate and could probably be omitted.

In GitLab by @pknight on Oct 20, 2014, 09:13

Okay, all outputs changed as request in commit aac768911 (r351).

In GitLab by @mkovari on Oct 21, 2014, 13:03

strncon /-0.005/ : strain in superconductor material (used in ITER Nb3Sn critical surface model)
Please use this same variable for CS.

Steel thickness: on sides only (not top and bottom). Please output a diagnostic eg strand current density with key parameters adjacent: Cu fraction, temp, B, strain.

In GitLab by @mkovari on Nov 6, 2014, 11:02

Peter - Have you incorporated a minimum temperature margin or maximum ratio (current density/critical current density) for CS and PF as you have for the TF? I think this is what @rkemp means by "a lack of margin for error." Clearly we need one or both of these.

In GitLab by @pknight on Nov 6, 2014, 11:17

For the CS one could simply reduce the upper limit of fjohc and fjohc0 from 1.0 to 0.9, say, to give a 10% safety margin.

For the PF coils there are no relevant constraint equations, for good reason. We could write a constraint to prevent any ratio of J_actual(ipf) / J_critical(ipf) from exceeding 1.0 (or 0.9) for all PF coils ipf, but this may well break the solution process; during the course of iteration one coil's J might suddenly become okay, only for another coil's limit to be reached, which will muck up the requirement of smooth gradients.

In GitLab by @mkovari on Nov 6, 2014, 11:21

For the TF coils don't we usually set the current density to 50% of the critical value? If so we should do the same for CS - setting the upper limits for fjohc and fjohc0 to 0.5, I suppose.
I am a bit confused: when NbTi is selected for the PF coils then the code is exactly the same as before - is that correct? @rkemp

In GitLab by @rkemp on Nov 6, 2014, 12:18

We do usually use 50% of jcrit as a limit, but there is the additional constraint applied by the stress on the coil, which helps prevent the field going too high. Do we have such a limit in the CS case?

I don't think we changed the NbTi model, but this still runs to 20T... I will need to check some references to see how realistic this is.

In GitLab by @rkemp on Nov 6, 2014, 14:49

I think we should have a think about the structural content, and what the winding pack actually consists of, to make sure we're being consistent in the use of the critical current parameterisation. This is the only limit, I think, on the size of the CS so it's quite critical. We could also consider adding a limit equation to set a maximum field on the coil. I think it should also be pretty easy to estimate the bursting stress and therefore the stresses in the structure as a sanity check. Can @mkovari @jmorris-uk @pknight and I get together sometime early next week to discuss?

In GitLab by @pknight on Nov 6, 2014, 16:32

Okay - added requested diagnostic outputs, and now the CS and PF calculations use strncon instead of 0.0 for the strain on the superconductor. This in itself makes a difference by a factor of two to the critical J. (Commit 49489522f - r355).

In GitLab by @mkovari on Nov 13, 2014, 13:59

I believe everything is sorted except for fixing what happens when the conductor cross-section comes out negative. I suggest:

δ = 0.0001
If conductor cross-section a > δ
   conductor current density j = I/a
Else If conductor cross-section a < δ
   conductor current density j = I(2δ-a)/ δ^2
End if

This creates a continuous and smooth function that increases monotonically as conductor cross-section decreases.

In GitLab by @rkemp on Nov 13, 2014, 14:29

That has to be worth a try: I'm in favour of limits which avoid singularities (provided there is a check that the area is positive at the end of the iteration).

In GitLab by @pknight on Nov 13, 2014, 16:22

Done in commit b718cdff6 (r365). Thankfully doesn't change any test results.

In GitLab by @mkovari on Nov 14, 2014, 08:49

I am a bit confused by the code.
! Non-steel cross-sectional area
awpoh = areaoh - areaspf
But then later this parameter is divided by the number of turns to become
aturn : input real : Area per turn (i.e. entire jacketed cable)
and then, oddly,
! Conduit and insulation around each turn is neglected, i.e. ! cable space area is set equal to the total area per turn
acs = aturn
and then
! Critical current density in winding pack
jcritwp = jcritstr * (1.0D0-fhe) * acs/aturn

I don't think that subroutine superconpf needs to know the cross-section, as it just calculates current density.

In GitLab by @rkemp on Nov 14, 2014, 09:36

The other subroutines (itersc, etc.) calculate the current density and don't need to know the cross-section; I think superconpf does need to know it to convert between the superconductor and the winding pack. Unfortunately this seems to be done in a number of different places:

In the main code:

awpoh = areaoh - areaspf

Within each superconpf for each isumat option:

jcritstr = jcritsc * (1.0D0-fcu)

(This is different for the Bi-2212 HTS, which is why it is separate.)

And then in the other section you quote:

jcritwp = jcritstr * (1.0D0-fhe) * acs/aturn

Overall this means that the total conductor current density is:

jcritwp = jcritstr * (1.0D0-fhe) * (1.0d0-fcu)

with the steel accounted for elsewhere. I think it's correct, albeit confusing. I assume there's a historical reason for the acs/aturn oddness -- Pete might know -- but it doesn't affect anything. However the definition of aturn at the start of superconpf is wrong: it is the conductor pack only, not including jacket.

In GitLab by @pknight on Nov 17, 2014, 09:49

Perhaps I will just remove acs and aturn completely. They are present in superconpf as this routine is derived from supercon in sctfcoil.f90, and this does differentiate between the two areas - each turn has a conducting part and a surrounding jacket as in Fig.3.5 of the User Guide.

Just for clarity, Rich's second jcritwp line above should read

jcritwp = jcritsc * (1.0D0-fhe) * (1.0d0-fcu)

In GitLab by @pknight on Nov 17, 2014, 10:34

Okay, done in commit bf499604c (r366).

ukaea / PROCESS

Central Solenoid current calculation #98