Refactor fwbs subroutine #180

Closed jonmaddock closed 1 year ago

In GitLab by @jmorris-uk on Jan 23, 2015, 09:48

Refactor the fwbs subroutine in PROCESS (as it is currently exceeding 1000 lines for a single subroutine) and have a think about the easiest layout for the blanket model. Do we need lblnkt, smstr and blnktmodel?

Discuss after ongoing work on the module is finished.

In GitLab by @mkovari on Feb 24, 2015, 11:45

This code in initial.f90 also needs attention as part of this task.

  !  But... set coolant to water if blktmodel > 0
  !  Although the *blanket* is by definition helium-cooled in this case,
  !  the shield etc. are assumed to be water-cooled, and since water is
  !  heavier (and the unit cost of pumping it is higher), the calculation
  !  for coolmass is better done with coolwh=2 if blktmodel > 0 to give
  !  slightly pessimistic results.

  if (blktmodel > 0) then
     blktcycle = 0
     blkttype = 3  !  HCPB
     coolwh = 2
  end if

  !  Ensure that blanket material fractions add up to 1.0

  if (blkttype < 3) then
     fsum = fblli2o + fblbe + vfblkt + fblss + fblvd
     if (abs(fsum-1.0D0) > 1.0D-4) then
        idiags(1) = blkttype
        fdiags(1) = fblli2o
        fdiags(2) = fblbe
        fdiags(3) = vfblkt
        fdiags(4) = fblss
        fdiags(5) = fblvd
        fdiags(6) = fsum
        call report_error(165)
     end if
  else
     fsum = fbllipb + fblli + vfblkt + fblss + fblvd
     if (abs(fsum-1.0D0) > 1.0D-4) then
        idiags(1) = blkttype
        fdiags(1) = fbllipb
        fdiags(2) = fblli
        fdiags(3) = vfblkt
        fdiags(4) = fblss
        fdiags(5) = fblvd
        fdiags(6) = fsum
        call report_error(165)
     end if
  end if

In GitLab by @mkovari on Mar 16, 2015, 09:05

 MDK queries up to line 1333 of hcpb.f90

A:  Important and/or urgent
---------------------------
!  No power lost. Assume all power is absorbed
    pnucloss = 0.0D0
This doesn't seem consistent with the calculation of pnucloss further up.

    !  First wall mass (kg)
    fwmass = denstl * volfw
The tungsten seems to have been omitted.

B:  Less so
-----------

Is the first wall volume calculated anywhere?
Can you explain these defaults?-
fvoldw /1.4/ : area coverage factor for vacuum vessel volume
fvolsi /0.64/ : area coverage factor for inboard shield volume
fvolso /0.64/ : area coverage factor for outboard shield volume 
 fhole /0.15/ : area fraction taken up by other holes 

Why is external_cryo_geometry in hcpb.f90?
The build diagrams in the User Guide need to be updated to show the 
cryostat lid correctly.  The lid now consists of a a thick support 
structure (hcryopf), plus a not very important skin (ddwex).
I should really have looked at the ITER mass as well, but of course 
this won't affect the cost in the new cost model.

The first wall void fractions do not really belong in nuclear_heating_magnets

 "Nuclear heating in TF coil" is actually :
"nuclear heating in TF+PF coils (CS is negligible)"

"Unit heating of FW armour" should read "Unit heating of FW and armour"

vardes:
fhole /0.15/ : area fraction taken up by other holes 
This default is much too large.  I suggest 0 as the default.

!+ad_summ  Calculations for powerflow
This doesn't seem to be enough to explain what this routine does, and whether it is different from the "detailed_powerflow_model".

    !  Inboard and outboard FW coolant void fraction
    vffwi = afwi*afwi/(bfwi*bfwi)
    vffwo = afwo*afwo/(bfwo*bfwo)
Haven't we already had this earlier?

vardes:
The following are used in the new thermodynamic blanket model (secondary_cycle > 0):
secondary_cycle /0/ : Switch for thermodynamic model of power conversion (secon dary) cy 
These two lines seem to be in the wrong order!

vardes:
 secondary_cycle /0/ : Switch for thermodynamic model of power conversion (secon dary) cy
= 0 simple model: a set efficiency for the chosen blanket design is use
= 1 simple model: a set efficiency for the chosen blanket design is use
> 1 detailed thermo-hydraulic and balance-of-plant model -
    = 2 use input thermal-electric efficiency (etath);
    = 3 steam Rankine cycle;
    = 4 supercritical CO2 cycle

This definition has been truncated so it doesn't make any sense.

!  If not superconducting coils then power to the TF is zero
Why?

subroutine detailed_powerflow_model
    !+ad_name  detailed_powerflow_model
    !+ad_summ  Calculations for detailed powerflow
I would describe this a first wall/blanket thermo-hydraulic model, not a powerflow model.

    !  Length of coolant pipes is assumed to be 80% of total radial space
    !  available for blanket, allowing for connections, manifolds etc.
This needs a diagram (sent to Chris H)

It would seem to make sense to include the lost alpha power (palpfwmw) and the first orbit loss power (porbitlossmw) in these equations, rather than piecemeal elsewhere:
    !  Surface heat flux on first wall (MW) (sum = pradfw)
    psurffwi = pradfw * fwareaib/fwarea
    psurffwo = pradfw * fwareaob/fwarea
For example, why are they not included here?-
    !  Total mass flow rate to remove inboard first wall power (kg/s)
    mffwi = 1.0D6*(pnucfwi + psurffwi) / (cf*(outlet_temp-inlet_temp))

    !  Repeat thermal hydraulic calculations for outboard side
    !  Calculation of max FW temp. Include NBI orbit loss power..
Are we assuming that all the first orbit losses are on the outboard side?  If so, should be stated explicitly.

In the long run, we should consider rewriting the first wall thermo-hydraulic calculation using constraints rather than a local iteration.

    !  Assumes up/down flow, two 90 deg bends per length
I think this is alreday stated earlier.

vardes:
The following are used in the (old) full thermodynamic blanket model (lblnkt=1):
What does this mean?  If this is Panos's model, I think we are removing it.

vardes:
 pdivt : power to divertor (MW) 
Should read 
 pdivt : power conducted to the divertor region (MW) 

    divplt /0.035/ : divertor plate thickness (m) (from Spears, Sept 1990) 
This is used to calculate the divertor volume and mass, but clearly this will give a very small volume.  The entire divertor assembly will be much thicker.

   !  Mass of He coolant = volume * density at typical coolant temperatures and pressures (kg)
I thought we were ignoring the mass of helium?

subroutine component_masses
It would be useful to give a reference for the densities used in this section.

In GitLab by @jmorris-uk on Mar 16, 2015, 11:13

Michael,

Urgent:

I have created a new variable for the power going to the vacuum vessel and beyond called pnucvvplus. Pnucloss is now only attributed to fhole like before.
I have added the variables fw_w_mass and fw_w_vol for the first wall armour mass and volume

Less so:

fvoldw has remained unchanged, fvolsi has been changed 1.0 as the inboard shield probably won't have penetrations and fvolso has also remained unchanged.
the default for fhole has been changed to 0.0
Changed "nuclear heating in TF coil" description in output
Changed "unit heating of FW armour" to "Unit heating of FW and armour"
removed second instance of vffwi and vffwo calculation
sorted "secondary_cycle" entry in vardes
removed the setting of ptfnuc to zero for resistive TF coils
changed detailed_powerflow_model to thermo_hydraulic_model
changed references to "first wall tungsten coating" to "first wall armour"
explicitly stated that iterate_fw for the outboard side considers orbit losses and that this is only for the outboard side.

In GitLab by @mkovari on Mar 16, 2015, 13:48

vardes: fwerlim /0.005/ : erosion thickness allowance for first wall (m) 
Should read:
fwerlim /0.005/ : erosion thickness allowance for first wall armour (m) 

subroutine iterate_fw
---------------------
Overall this routine is hard to understand, with key questions about the distinctions between bulk temperature and wall   
temperature, fluence and wall loading, still obscure to me.  Probably we should replace it with something based on a   
flat first wall. (see issue #236.)

The quantities 'tav' and 'tmpdif' are not explained.

       !  Coolant properties at average coolant temperature
Should read
       !  Coolant properties at average coolant temperature (tb)

vardes:       
outlet_temp /598.0/ : outlet temperature of coolant for blanket and first wall (K) (sec input if coolwh=1 (helium),  
calculated if coolwh=2 (water)   
What is "sec"?
How does the different treatment of water and helium affect the iterate_fw loop?

       !  tb is not known at this point,
Why is tb not known?  A few lines back (line 1407), both the inlet and outlet temperatures are known.         

        but since we only need the cf output (which does not depend on tb)  
Why doesn't cf depend on tb?

        the temperature on the inner
       !  wall (in contact with the coolant) and therefore equal to the bulk
       !  coolant temperature,
The wall temperature and bulk coolant temperature are not the same.

The loop to determine the average first wall temperature is a very crude approach.

          !  Temperature in Kelvin and deg C
Once function tk(t) is removed (see blow), there will be no need to calculate the temperature in C.

          !  Heat transfer coefficient calculated using Sieder-Tate correlation, 
          !  valid for Re > 1.0e3, 0.7 < Pr < 16700, L/D > 10
The limits of applicability should be checked on the last run through before the final output.
I suggest calculating Re, Pr and L/D separately, then using the standard formula for the Sieder-Tate correlation using   
them.
Note that Wikipedia says Re>10000, not 1000 - needs to be checked.         

       !  Fluence
       flnce = 1.0D-6*qppp * fwvol/area * fwlife
Please state units of 'flnce'.       

             write(*,*) 'Warning in routine ITERATE_FW:'
             write(*,*) 'Swelling limit exceeded, and'
             write(*,*) 'optimisation is failing to find a'
             write(*,*) 'suitable first wall thickness...'
             write(*,*) 'PROCESS continuing.'
It would be nice to get this warning into a single line of output, as it can appear many times,   
pushing other warnings off the screen.

    The lower limit on the first wall thickness is derived from the
    constraint that the first wall must possess the ability to withstand
    the internal coolant pressure.
Please confirm in this comment that this is the correct thick wall formula (see issue #174).

function smt(t, fwlifs)
-----------------------
Calculates the maximum stress intensity for the first wall, from fits via the ASME code.

Can we discuss if we want to keep this?  It will be opaque to anyone looking at the final PROCESS output, but it does   
have the smooth dependence on engineering parameters that Rich likes.

It was formerly called from:
       maxstress = 1.0D9 !smt(tpeakfw_c, fwlifs)

At the very least the maximum stress needs to be a user input, together with a maximum temperature at that stress.       

A more flexible approach, similar to that used in the new availability code, allows the user to set a chisel-shaped stress curve:
stress v temperature:
________________
                \
                 \
                  \
                   \
____________________\___________________________  

If we remove this routine (which is currently not called), then I think first wall lifetime no longer affects anything in  
 the loop - we can remove it from the loop and calculate it only once.

First wall lifetime
-------------------
I can't understand the apparently circular reasoning here:
line 214:    fwlife = min(abktflnc/wallmw, tlife)
line 1517:   flnce = 1.0D-6*qppp * fwvol/area * fwlife
Furthermore, the term fluence is seriously ambiguous here.  In line 1517 it is calculated by dividing the neutron power  
deposited in the first wall by the area, whereas  line 214 uses 'wallmw : average neutron wall load', a completely    
different quantity, which takes the total neutron power divided by wall area - not just the power absorbed in the first   
wall, which is much less.

Obviously, we need to look at the original data relating to fluence limit - if we knew where was.

function tk(t)
----------------
This function doesn't make a lot of sense at present.
I suggest we remove the unreferenced stainless steel formula, and use the  Eurofer formula.  
  Note that the formula uses T in KELVIN, so it is not correct as coded.  I have checked the reference, and the formula   
is correct except for one error in the last decimal place:
0.00023863D0  should be 0.00023862D0.
The function is just one line of code - I suggest removing it and placing the formula in the body of the calling routine.

subroutine cprops(tb, cf, rhof, viscf, viscfs, kf)
----------------------------------------------------

    x = 0.5D0*(outlet_temp + inlet_temp)
    y = x + 0.5D0*( (tb - outlet_temp) + (tb - inlet_temp) )
Remarkable code!  It seems that y=tb  !

After we have done a thorough comparison, I suggest we remove the Panos correlations.

This whole subroutine can then disappear, with fluid_properties being called directly instead.     

Note that this subroutine is in any case coded in a misleading way, using both arguments and module-level   
variables.      

function tsat(p)
----------------
This can also be removed, as it will just be a single call to tsat_refprop.

function pumppower
------------------
The term "power" should be replaced by "mechanical power" throughout this routine, as there is also a   
separate electrical-to-mechanical conversion efficiency:
etahtp /0.95/ : electrical efficiency of primary coolant pumps 

    For a water coolant the pumping power is
    calculated as the volume flow rate multiplied by the total pressure
    drop, ..
I suggest using the exact same method for water as for helium.  This will make the results accurate even if   
the water is in a compressible regime (near to or above the critical temperature).

In GitLab by @mkovari on Mar 16, 2015, 15:45

plant_power.f90

Important/Urgent

subroutine power1

    if (secondary_cycle == 0) then
        !  Primary thermal power (MW)
        pthermmw = pthermfw + pthermblkt + iprimshld*pthermshld 
        !  Secondary thermal power deposited in divertor (MW)
        psecdiv = pthermdiv     
        ! Divertor primary/secondary power switch value
        iprimdiv = 0

This doesn't seem to be consistent with Table 5 in the engineering paper, which allows the divertor heat to be used (or not). Maybe I have missed something.

Less So

vardes:
tmpcry /4.5/ : cryostat temperature for cryogenic plant power calculation (K)

This is the coil temperature, not the cryostat temperature.

subroutine power2(outfile,iprint)
---------------------------------

    !  Centrepost coolant pump power (ST)

What's this line doing here?

vardes:  tfcmw : peak power per TF power supply (MW)

This can't be right (I hope), since this is used as total electrical power required for the TF coil supply. Also, it is constant, and doesn't have a peak.

!  Power consumed by fusion power core systems

Should read

!  Electrical power consumed by fusion power core systems

vardes:

htpsecmw : waste power lost from heat transport system (MW)

Should read

htpsecmw : Waste power lost from primary coolant pumps (MW)

In GitLab by @mkovari on Mar 17, 2015, 15:02

fw_armour_thickness : first wall armour thickness (m)
Should this be an input with a default?

ipowerflow /1/ : switch for power flow model:

    = 0 pre-2014 version;
    = 1 comprehensive 2014 model

Should this still be here?

In GitLab by @jmorris-uk on Mar 18, 2015, 09:44

Changes

For the urgent issue I think what is there is already what we want. secondary_cycle = 0 the divertor power doesn't go to the primary loop and for secondary_cycle >= 1 the divertor power goes to the primary loop.
fw_armour_thickness had a default of 0.005 but I hadn't put it in the comment line. It's now there.
ipowerflow will be removed
changed waste power comment
changed electrical power consumed comment
have now put the ppumpmw (centre post pump power (ST)) into a if statement for itart to make it clearer.
cprops has been changed to just call fluid_properties. Panos' correlations have been removed. irefprop can be removed from the code
have removed tsat and replaced with tsat_refprop for when it is called. (Which is currently commented out anyway as none of the options have water coolant at the moment)
units of fluence are MW-yr/m^2
...

In GitLab by @jmorris-uk on Mar 18, 2015, 14:35

mentioned in commit 0d3e93b87140809178b6d25608f475f606a5d33d

In GitLab by @jmorris-uk on Apr 28, 2015, 14:04

Mostly done in r389 and r390. Outstanding items have their own issues.

ukaea / PROCESS

Refactor fwbs subroutine #180