plasmod.f90 #644

Closed jonmaddock closed 1 year ago

In GitLab by @efable on Mar 23, 2018, 07:27

please consider this file and now it is consistent with PLASMOD @hlux @kellis plasmod.f90

In GitLab by @hlux on Mar 23, 2018, 09:34

Can you please clarify the impurity concentrations in PLASMOD again. I somehow thought you are not calculating impurity transport and are therefore also using fixed fractions. You say now that you are only fixing the fraction at the pedestal and then calculate impurity transport in the core??

In GitLab by @efable on Mar 23, 2018, 09:36

Well, in reality later on I will put transport of everything: He, H, Xe, Ar, etc. which we know

for now, H, He have fixed conc

for Ar and Xe, one can choose fixed conc, or fixed density with conc given at the pedestal top

there is no transport now

In GitLab by @hlux on Mar 23, 2018, 09:38

Sorry, just for me again, with fixed density you mean that you do not vary the Ar or Xe density and you match it to the concentration at the pedestal top? What describes the density distribution then, if this is fixed?

In GitLab by @efable on Mar 23, 2018, 09:38

fixed means flat. ne(0) = ne(pedestal top) (for Xe and Ar only)

In GitLab by @efable on Mar 23, 2018, 09:39

not ne in the sense of electrons, it is only for Xe and Ar. H and He are fixed concentration

In GitLab by @efable on Mar 23, 2018, 09:40

the reason is because transport would give that, but I will put it later

In GitLab by @hlux on Mar 23, 2018, 10:02

Okay, we would like to add the impurity densities to our plots, to have a visual aid. So now that you have clarified that the impurity densities for Ar/Xe in case of impmodel=1 are constant at ne(pedestal top)* f_imp , I assumes this is only true inside the pedestal top? What do the profiles look like in the actual pedestal region?

In GitLab by @efable on Mar 23, 2018, 10:02

in the pedestal they are also fixed density, so basically is flat density everywhere.

In GitLab by @hlux on Mar 23, 2018, 10:03

Ah, okay. Thanks! Hanni

In GitLab by @hlux on Mar 23, 2018, 10:11

Just to clarify before I push

if (plasmod_i_impmodel == 0 ) then ! Total impurity density (/m3) dnz = 0.0D0 do imp = 1,nimp if (impurity_arr(imp)%Z > 2) then dnz = dnz + impurity_arr(imp)%frac*dene end if end do

   if ( (dnz - radp%av_nz*1.d19)/dnz > 1e-6) then
      fdiags(1) = dnz; fdiags(2) = radp%av_nz*1.d19
      call report_error(194)
   endif
else
   dnz = radp%av_nz*1.d19
endif

This would then be the right check?

Implemented in Commit c7f03d0f

In GitLab by @efable on Mar 23, 2018, 10:14

yes

In GitLab by @efable on Mar 23, 2018, 10:32

ok pulled, running it now, lets see

In GitLab by @efable on Mar 23, 2018, 10:34

seems to run, lets see if I get a feasible solution, I attach here the IN.DAT:

-------------------------------------------------- verbose = 1

---------------Constraint Equations---------------

icc = 8 Neutron wall load upper limit icc = 10 Toroidal field 1 icc = 11 Radial build icc = 13 Burn time lower limit icc = 16 Net electric power lower limit icc = 25 Peak toroidal field upper limit icc = 26 Central solenoid EOF current density upper limit icc = 27 Central solenoid BOP current density upper limit icc = 31 TF coil case stress upper limit icc = 32 TF coil conduit stress upper limit icc = 33 I_op icc = 34 Dump voltage upper limit icc = 35 J_winding pack icc = 36 TF coil temperature margin lower limit icc = 60 Central solenoid temperature margin lower limit icc = 65 Dump time set by VV loads icc = 72 central solenoid Tresca stress limit icc = 68 Psep

---------------Iteration Variables----------------

boundl(103) = 1.0

ixc = 145 fgwped boundl(145)=0.5 boundu(145)=0.7

ixc = 36 fbetatry

ixc = 117 * fpsepbqar

ixc = 44 fvsbrnni boundl(44)=0.1 boundu(44)=0.6

ixc = 1 aspect boundu(1) = 5. boundl(1) = 1.5 ixc = 2 bt boundu(2) = 20.0 ixc = 3 rmajor boundu(3) = 13 ixc = 18 * q boundl(18) = 3.5

ixc = 12 oacdcp ixc = 13 tfcth ixc = 14 fwalld ixc = 16 ohcth ixc = 29 bore boundl(29) = 0.1 ixc = 37 coheof ixc = 38 fjohc boundu(38) = 1.0 ixc = 39 fjohc0 boundu(39) = 1.0 ixc = 41 fcohbop ixc = 42 gapoh boundl(42) = 0.05 boundu(42) = 0.1 ixc = 48 fstrcase ixc = 49 fstrcond ixc = 50 fiooic boundu(50) = 1.0 ixc = 51 fvdump ixc = 52 vdalw boundu(52) = 10.0 ixc = 53 fjprot ixc = 54 ftmargtf ixc = 56 tdmptf ixc = 57 thkcas ixc = 58 thwcndut boundl(58) = 8.0d-3 ixc = 59 fcutfsu boundl(59) = 0.50 boundu(59) = 0.94 ixc = 60 cpttf boundl(60) = 6.0d4 boundu(60) = 9.0d4 ixc = 61 gapds boundl(61) = 0.02 ixc = 106 ftmargoh ixc = 113 ftaucq ixc = 122 oh_steel_frac ixc = 123 * foh_stress

-----------------Build Variables------------------

blnkith = 0.755 Inboard blanket thickness (m); blnkoth = 0.982 Outboard blanket thickness (m); bore = 2.3261 Central solenoid inboard radius (m) ddwex = 0.15 Cryostat thickness (m) ddwi = 0.30 Vacuum vessel thickness (tf coil / shield) (m) gapds = 0.02 Gap between inboard vacuum vessel and tf coil (m) gapoh = 0.05 Gap between central solenoid and tf coil (m) gapomin = 0.20 Minimum gap between outboard vacuum vessel and tf coil (m) iohcl = 1 Switch for existence of central solenoid; ohcth = 0.54761 Central solenoid thickness (m) scrapli = 0.225 Gap between plasma and first wall; inboard side (m) scraplo = 0.225 Gap between plasma and first wall; outboard side (m) shldith = 0.30 Inboard shield thickness (m) shldoth = 0.80 Outboard shield thickness (m) shldtth = 0.30 Upper/lower shield thickness (m); tfcth = 1.2931 Inboard tf coil thickness; (centrepost for st) (m) tftsgap = 0.05 Manufacturing/thermal expansion gap between tf and thermal shield (m) vgap2 = 0.05 Vertical gap between vacuum vessel and tf coil (m) vvblgap = 0.02 * Gap between vacuum vessel and blanket (m)

---------------Buildings Variables----------------

---------------Constraint Variables---------------

bmxlim = 12.5 Maximum peak toroidal field (t) fbetatry = 0.41378 F-value for beta limit fdene = 1.2 F-value for density limit ffuspow = 1 F-value for maximum fusion power fiooic = 0.60263 F-value for tf coil operating current / critical fjohc = 0.58039 F-value for central solenoid current at end-of-flattop fjohc0 = 0.54091 F-value for central solenoid current at beginning of pulse fjprot = 1.0 F-value for tf coil winding pack current density flhthresh = 1.4114 F-value for l-h power threshold foh_stress = 1.0 F-value for tresca stress in oh coil fpeakb = 1.0 F-value for maximum toroidal field fpinj = 1.0 F-value for injection power fpnetel = 1.0 F-value for net electric power fpsepbqar = 1.0 F-value for maximum psepbt/qar limit fstrcase = 1.0 F-value for tf coil case stress fstrcond = 0.84474 F-value for tf coil conduit stress ftaucq = 1.0 F-value for calculated minimum tf quench time ftburn = 1.00e+00 F-value for minimum burn time ftmargoh = 1.0 F-value for central solenoid temperature margin ftmargtf = 1.0 F-value for tf coil temperature margin fvdump = 0.97322 F-value for dump voltage fwalld = 0.12856 F-value for maximum wall load pnetelin = 500. Required net electric power (mw) psepbqarmax = 9.2 Maximum ratio of psepbt/qar (mwt/m) tbrnmn = 7.2e3 Minimum burn time (s) walalw = 8.0 Allowable wall-load (mw/m2) ftaulimit = 1.0 F-value for lower limit on taup/taueff the ratio of alpha particle to energy confinement times ------------------Cost Variables------------------*

abktflnc = 15 Allowable first wall/blanket neutron adivflnc = 20.0 Allowable divertor heat fluence (mw-yr/m2) cfactr = 0.75 Total plant availability fraction; cost_model = 0 Switch for cost model; dintrt = 0.00 Diff between borrowing and saving interest rates fcap0 = 1.15 Average cost of money for construction of plant fcap0cp = 1.06 Average cost of money for replaceable components fcontng = 0.15 Project contingency factor fcr0 = 0.065 Fixed charge rate during construction fkind = 1.0 Multiplier for nth of a kind costs iavail = 0 Switch for plant availability model; ifueltyp = 1 Switch; lsa = 2 Level of safety assurance switch (generally; use 3 or 4); output_costs = 0 Switch for costs output; ratecdol = 0.06 Effective cost of money in constant dollars tlife = 40 Plant life (years) ucblvd = 280.0 Unit cost for blanket vanadium ($/kg) ucdiv = 5.0d5 Cost of divertor blade ($) ucme = 3.0d8 * Unit cost of maintenance equipment ($/w**0;3)

-------------Current Drive Variables--------------

bscfmax = 0.99 Maximum fraction of plasma current from bootstrap; etanbi = 0.4 Ech wall plug to injector efficiency gamma_ecrh = 0.30 User input ecrh gamma iefrf = 8 Switch for current drive efficiency model; pheat = 50. Heating power not used for current drive (mw) pinjalw = 50. Maximum allowable value for injected power (mw)

----------Divertor Kallenbach Variables-----------

-------------------Divertor Ode-------------------

----------------Divertor Variables----------------

divdum = 1 Switch for divertor zeff model; 0=calc; 1=input divfix = 0.621 Divertor structure vertical thickness (m) hldivlim = 10 Heat load limit (mw/m2) ksic = 1.4 Power fraction for outboard double-null scrape-off plasma prn1 = 0.4 N-scrape-off / n-average plasma; zeffdiv = 3.5 Zeff in the divertor region (if divdum /= 0)

------------------Fwbs Variables------------------

inuclear = 1 Switch for nuclear heating in the coils; qnuc = 1.292e4 Nuclear heating in the coils (w) (inuclear=1) primary_pumping = 3 Switch for pumping power for primary coolant (06/01/2016); secondary_cycle = 2 Switch for power conversion cycle; vfshld = 0.60 Coolant void fraction in shield etaiso = 0.9 Isentropic efficiency of fw and blanket coolant pumps etahtp = 0.87 * Electrical efficiency of primary coolant pumps

-----------------Global Variables-----------------

-------------Heat Transport Variables-------------

etath = 0.375d0 Thermal to electric conversion efficiency ipowerflow = 0 Switch for power flow model; iprimshld = 1 * Switch for shield thermal power destiny;

------------------Ife Variables-------------------

------------Impurity Radiation Module-------------

imprad_model = 1 Switch for impurity radiation model; coreradius = 0.75 Normalised radius defining the 'core' region coreradiationfraction = 0.6 Fraction of radiation from 'core' region that is subtracted from the loss power fimp(1) = 1.0 fimp(2) = 0.1 fimp(3) = 0.0 fimp(4) = 0.0 fimp(5) = 0.0 fimp(6) = 0.0 fimp(7) = 0.0 fimp(8) = 0.0 fimp(9) = 0.0 fimp(10) = 0.0 fimp(11) = 0.0 fimp(12) = 0.0 fimp(13) = 0.00044 fimp(14) = 0.0 fimpvar = 0.00035148 Impurity fraction to be used as fimp(impvar) impvar = 13 Impurity to be iterated (deprecated) iradloss = 0 ---------------------Numerics---------------------*

ioptimz = 1 for optimisation VMCON only minmax = 1 Switch for figure-of-merit (see lablmm for descriptions) epsvmc = 1.0e-4 Error tolerance for vmcon plasmod_tol= 1.0e-7 tolerance to be reached in % variation at each time step plasmod_dtmin=0.01d0 min time step plasmod_dtmax=0.05d0 max time step plasmod_eopt=0.d0 exponent of jipperdo plasmod_dtmaxmin=0.d0 exponent of jipperdo2 plasmod_test=1000. max iteration number maxcal = 300 ----------------Pf Power Variables----------------*

-----------------Pfcoil Variables-----------------

alstroh = 6.6d8 Allowable hoop stress in central solenoid structural material (pa) coheof = 20837000.0 Central solenoid overall current density at end of flat-top (a/m2) cptdin = 4.22d4, 4.22d4, 4.22d4, 4.22d4, 4.3d4, 4.3d4, 4.3d4, 4.3d4, Peak current per turn input for pf coil i (a) fcohbop = 0.93313 Ratio of central solenoid overall current density at fcuohsu = 0.70 Copper fraction of strand in central solenoid cable ipfloc = 2,2,3,3 Switch for locating scheme of pf coil group i; isumatoh = 5 Switch for superconductor material in central solenoid; isumatpf = 3 Switch for superconductor material in pf coils; ncls = 1,1,2,2, Number of pf coils in group j ngrp = 4 Number of groups of pf coils; ohhghf = 0.9 Central solenoid height / tf coil internal height oh_steel_frac = 0.57951 Central solenoid steel fraction (iteration variable 122) rjconpf = 1.1d7, 1.1d7, 6.d6, 6.d6, 8.d6, 8.0d6, 8.0d6, 8.0d6, Average winding pack current density of pf coil i (a/m2) rpf2 = -1.825 Offset (m) of radial position of ipfloc=2 pf coils zref(1) = 3.6 zref(2) = 1.2 zref(3) = 1.0 zref(4) = 2.8 zref(5) = 1.0 zref(6) = 1.0 zref(7) = 1.0 zref(8) = 1.0

----------------Physics Variables-----------------

alphan = 1.00 Density profile index alphat = 1.45 Temperature profile index aspect = 3.1 Aspect ratio (iteration variable 1) beta = 0.025927 Total plasma beta (iteration variable 5) bt = 5.8547 Toroidal field on axis (t) (iteration variable 2) dene = 7.2607e+19 Electron density (/m3) (iteration variable 6) dnbeta = 3.0 (troyon-like) coefficient for beta scaling; fgwped = 0.85 Fraction of greenwald density to set as pedestal-top density fgwsep = 0.5 Fraction of greenwald density to set as separatrix density fkzohm = 1.0245 Zohm elongation scaling adjustment factor (ishape=2; 3) fvsbrnni = 0.38686 Fraction of the plasma current produced by gamma = 0.3 Ejima coefficient for resistive startup v-s formula hfact = 0.95 H factor on energy confinement times (iteration variable 10) ibss = 4 Switch for bootstrap current scaling; iculbl = 1 Switch for beta limit scaling (constraint equation 24); icurr = 4 Switch for plasma current scaling to use; idensl = 7 Switch for density limit to enforce (constraint equation 5); ifalphap = 1 Switch for fast alpha pressure calculation; ifispact = 0 Switch for neutronics calculations; iinvqd = 1 Switch for inverse quadrature in l-mode scaling laws 5 and 9; impc = 0. Carbon impurity multiplier (imprad_model=0 only) impo = 0. Oxygen impurity multiplier (imprad_model=0 only) ipedestal = 3 Switch for pedestal profiles; ieped = 1 Switch for scaling pedestal-top temperature with plasma parameters; neped = 0.678e20 Electron density of pedestal [m-3] (ipedestal=1) nesep = 0.2e20 Electron density at separatrix [m-3] (ipedestal=1) rhopedn = 0.94 R/a of density pedestal (ipedestal=1) rhopedt = 0.94 R/a of temperature pedestal (ipedestal=1) tbeta = 2.0 Temperature profile index beta (ipedestal=1) teped = 5.5 Electron temperature of pedestal (kev) (ipedestal=1; ieped=0) tesep = 0.1 Electron temperature at separatrix (kev) (ipedestal=1) iprofile = 1 Switch for current profile consistency; isc = 34 Switch for energy confinement time scaling law ishape = 0 Switch for plasma cross-sectional shape calculation; kappa = 1.848 Plasma separatrix elongation (calculated if ishape > 0) q = 3.8871 Safety factor 'near' plasma edge (iteration variable 18); q0 = 1.0 Safety factor on axis ralpne = 0.070993 Thermal alpha density / electron density (iteration variable 109) rmajor = 9.0019 Plasma major radius (m) (iteration variable 3) snull = 1 Switch for single null / double null plasma; ssync = 0.6 Synchrotron wall reflectivity factor te = 12.485 Volume averaged electron temperature (kev) triang = 0.5 Plasma separatrix triangularity (calculated if ishape=1; 3 or 4) zfear = 1 High-z impurity switch; 0=iron; 1=argon -----------------Pulse Variables------------------

lpulse = 1 * Switch for reactor model;

------------------Rfp Variables-------------------

-------------------Scan Module--------------------

--------------Stellarator Variables---------------

-----------------Tfcoil Variables-----------------

alstrtf = 6.6e8 Allowable von mises stress in tf coil structural material (pa) casthi = 0.06 Inboard tf coil case inner (plasma side) thickness (m) casths = 0.05 Inboard tf coil sidewall case thickness (m) cpttf = 90000.0 Tf coil current per turn (a); dhecoil = 0.010 Diameter of he coil in tf winding (m) fcutfsu = 0.77798 Copper fraction of cable conductor (tf coils) isumattf = 5 Switch for superconductor material in tf coils; oacdcp = 8935500.0 Overall current density in tf coil inboard legs (a/m2) ripmax = 0.6 Maximum allowable toroidal field ripple amplitude tdmptf = 29.265 Dump time for tf coil (s) tfno = 16 Number of tf coils (default = 50 for stellarators) tftmp = 4.750 Peak helium coolant temperature in tf coils and pf coils (k) thicndut = 1.5d-3 Conduit insulation thickness (m) thkcas = 0.60557 Inboard tf coil case outer (non-plasma side) thickness (m) thwcndut = 0.008 Tf coil conduit case thickness (m) tinstf = 0.008 Ground insulation thickness surrounding winding pack (m) tmargmin = 1.500 Minimum allowable temperature margin (cs and tf coils) (k) vdalw = 9.3753 Max voltage across tf coil during quench (kv) vftf = 0.300 * Coolant fraction of tfc 'cable' (itfsup=1); or of tfc leg (itfsup=0)

-----------------Times Variables------------------

pulsetimings = 0 Switch for pulse timings (if lpulse=1); tburn = 1.0d4 Burn time (s) (calculated if lpulse=1) tdwell = 0 Time between pulses in a pulsed reactor (s) tramp = 500.0 Initial pf coil charge time (s); if pulsed; = tohs

-----------------Vacuum Variables-----------------

----------PLASMOD 1D Transport Variables---------- plasmod_dt=0.01d0 time step plasmod_dtinc=2.d0 decrease of dt plasmod_ainc=1.1d0 increase of dt plasmod_tolmin=10.1d0 multiplier of etolm that should not be overcome plasmod_capa=0.1d0 first radial grid point plasmod_maxa=0.d0 diagz 0 or 1 plasmod_dgy=1.d-5 Newton differential plasmod_i_modeltype=1 1 - simple gyrobohm scaling plasmod_i_equiltype=1 1 - EMEQ solve equilibrium with given q95 with sawteeth 2- EMEQ solve with given Ip with sawteeth plasmod_nx=41 number of interpolated grid points plasmod_nxt=7 number of reduced grid points plasmod_nchannels=3 leave this at 3 plasmod_i_impmodel=1 * impurity model: 0 - fixed concentration 1 - concentration fixed at pedestal top then fixed density plasmod_globtau(1)=4.0d0 plasmod_globtau(2)=4.0d0 plasmod_globtau(3)=4.0d0 plasmod_globtau(4)=4.0d0 plasmod_globtau(5)=1.0d0 impurity_enrichment(9)=20.d0 plasmod_qdivt=0.0d0 plasmod_v_loop=-1.d4 plasmod_qnbi_psepfac = 100.0d0 plasmod_cxe_psepfac = 1.0d-3 plasmod_car_qdivt = 1.0d-4 plasmod_x_heat(1)=0.0d0 plasmod_x_heat(2)=0.0d0 plasmod_dx_heat(1)=0.2d0 plasmod_dx_heat(2)=0.03d0 plasmod_x_cd(1)=0.0d0 plasmod_x_cd(2)=0.0d0 plasmod_dx_cd(1)=0.2d0 plasmod_dx_cd(2)=0.03d0 plasmod_x_fus(1)=0.0d0 plasmod_x_fus(2)=0.0d0 plasmod_dx_fus(1)=0.2d0 plasmod_dx_fus(2)=0.03d0 plasmod_x_control(1)=0.0d0 plasmod_x_control(2)=0.0d0 plasmod_dx_control(1)=0.2d0 plasmod_dx_control(2)=0.03d0 plasmod_nbi_energy = 1000.0d0 csawth = 1.4 plasmod_pedscal=1.d0

In GitLab by @efable on Mar 23, 2018, 10:34

obviously it doesnt give a duck about the enter button...

In GitLab by @efable on Mar 23, 2018, 10:34

IN_plasmod_converges_1D.DAT

In GitLab by @efable on Mar 26, 2018, 07:13

closed

ukaea / PROCESS

plasmod.f90 #644