Design Point Module Equations

Revision: 2/2012 by L. Carlson

Code Version:ASC v.35

Code Modules:DesignPoint.cpp, DesignPoint.h

History: Z. Dragojlovic received from R. Raffray ~ 2/2007

Module note: This is the main engineering module that reads in multiple input files from powerflow.data such as power fractions, efficiencies, etc. It calculates CD powers, divertor and blanket pumping powers, plasma surface area, volumes, neutron wall loading and gross cycle efficiency. It also compiles and displays the System.data output file.

Reactor Building Volume

Calculates the exterior reactor building volume in m3 from the cryodome part.

[already in code, no reference]

vol_reactor_bldg = (pow((ksihi + 9.0), 2.0)*6.0*etahi + 1.55e5)*cATfit;
cATfit = 138698.0/181655.0 (investigation found that 138698 is Aries-AT volume from published literature and 181655 is Prometheus-HI volume, no reference)
ksihi = outboard max width of the cryodome
etahi = vertical max height of the cryodome

Support Structure Building Volume

Approximated as 20% of the fusion power core volume.

[already in code, no reference]

vol_str = 0.20*FPC_vol; // FPC_vol calculated below

Plasma Surface Area

Calculates the surface area of the plasma. The average neutron wall flux uses this by dividing the neutron power by the plasma surface area.

Elliptical equations: Good for estimation but not accurate when plasma elongation, squareness, and triangularity push the plasma shape beyond an ellipse.

[C. Kessel 5/2011 “Plasma Contour2 - Kessel plasma geometry discretized formulas.doc”]

A_plasma_surface = 4.0*pi*pi*R * (R / A + SOL) * sqrt((1.0 + pow(k, 2.0))/ 2.0);
SOL = SOL thickness from input file
A = aspect ratio from physics module output
k = kappa, plasma elongation
R = plasma major radius

Perimeter method: Uses the perimeter function in Part.cpp to find the plasma’s perimeter and then revolve it about 2R.

[from ASC Part.cpp module]

A_plasma_surface = (*Plasma).surface_area("plasma_perimeter")*2.0*pi*plasma_major_R;

Discretized equations: These take into account all the plasma factors that affect its size. This is the most accurate formula and is currently used in the code. The contour is multiplied by two to take advantage of symmetry.

[C. Kessel 5/2011 “Plasma Contour2 - Kessel plasma geometry discretized formulas.doc”]

R = Ro + a cos(ϑ + sin-1(δ) sin(ϑ))
Z = a κ sin(ϑ + ςo sin(2ϑ)) for ϑ < π/2
Z = a κ sin(ϑ + ςi sin(2ϑ)) for ϑ > π/2
Ro = plasma major radius
a = plasma minor radius (add SOL thickness here to find surface area incl. SOL)
κ = plasma elongation
δ = plasma triangularity
ςo = plasma squareness on outboard side (~0)
ςi = plasma squareness on inboard side (~0)
ϑ = poloidal angle

NWL Peaking Factor

Two tables, one for SCLL and one for DCLL blanket, give peaking factors depending on plasma aspect ratio. If the aspect ratio is not found in the table, then the code interpolates.

[L. El-Guebaly updated 5/2011 at project meeting.]

DCLL: at A = 1.6, PF = 1.94; at A = 4.0, PF = 1.5
DCLL: at A = 1.7, PF = 2.0; at A = 4.0, PF = 1.55

Blanket Pumping Powers

This calculates the pumping power requirement for the blanket in use, SCLL or DCLL.

SCLL Blanket: (no helium, liquid metal-cooled only). The old method simply divided the pumping power 50:50 between the blanket and divertor. The correlation below only calculates the IB first wall and blanket pumping power now. The OB will be much lower due to lower field strength and more space available. Therefore, the code only uses P_rad_first_wall_IB and NWL_ave_IB to calculate P_blanket_IB. The output is filtered out if the average neutron wall flux exceeds the data range given in the spreadsheet, which is 2 – 4 MW/m2.

[M. Tillack 3/2011 “SiC Blanket Pumping Power.xls"]

P_pump_blanket = (0.43840 * (flux_neutron_wall_Ave/1.0e6) + 0.00036) * P_blanket_sic / 100.0;
P_rad_first_wall_IB = (f_rad_blanket_ib) * P_rad_chamb + (1.0 - f_rad_edge) * P_rad_edge * f_rad_edge_blanket_ib;
A_FW_IB = 16.0/2.0*( (*FirstWallInb).surface_area("outer") );
flux_neutron_wall_IB = P_neutron * f_n_blanket_ib / A_FW_IB;
P_blanket_sic = (1.0 - f_alpha_div) * P_alpha_loss + f_n_blanket_ib * P_neutron * f_nmul_fw_inboard + P_rad_first_wall_IB;
P_blanket = P_blanket_sic;
P_pump_blanket_loss = (1.0 - eta_pump_lm_mech) * P_pump_blanket / eta_pump_lm_mech;

DCLL Blanket: (liquid-metal and helium cooling). The correlation below combines both the liquid-metal and helium cooling, as well as the inboard and outboard pumping powers. The output is filtered out if the average neutron wall flux exceeds the data range given in the spreadsheet, which is 2 – 4 MW/m2.

[M. Tillack 1/2011 “DCLL Blanket Pumping Power.xls"]

P_pump_blanket = (0.03602*pow(flux_neutron_wall_Ave/1.0e6, 2) +

0.03056 * (flux_neutron_wall_Ave/1.0e6) +

0.38359) * P_blanket_dcll/100.0;

P_blanket_dcll = (1.0 - f_alpha_div) * P_alpha_loss

+ (f_n_blanket_ib) * P_neutron * f_nmul_fw_inboard

+ (f_n_blanket_ob) * P_neutron * f_nmul_fw_outboard

+ P_rad_first_wall; // both IB+OB

// Need coolant fractional power split for DCLL blanket pump losses (LM pump vs He-blower); Mark's spreadsheet has PbLi 62.3%, He 37.7%
double f_pump_PbLi = 0.623;
double f_pump_He = 0.377;
P_pump_blanket_loss = ((1.0 - eta_pump_lm_mech) * f_pump_PbLi +

(1.0 - eta_pump_he) * f_pump_He) *

P_pump_blanket / ((eta_pump_he + eta_pump_lm_mech)/2);

Max Heat Flux on Divertor

This section of the code stems from C. Kessel’s sysengr.f systems code. The output is filtered if the maximum heat flux on the divertor is > 8 MW/m2.

[C. Kessel, existing code, no specific reference. Edited “adivo” and “adivi” with CK 1/2011]

Definitions:

sym = 1.0 for up-down symmetric, 0.0 for up-down asymmetric
facdiv = factor for splitting power if up-down symmetric or not
adivo = area in outboard divertor slot for absorbing radiated power
adivi = area in inboard divertor slot for absorbing radiated power
fdivrad = total radiated power fraction in divertor
fdivcondo = fraction of conducted power to divertor going to the outboard
fdivcondi = fraction of conducted power to divertor going to the inboard
fdivrado = fraction of radiated power to divertor going to the outboard
fdivradi= fraction of radiated power to divertor going to the inboard
fluxexpo = poloidal mag. flux exp. on outboard from midplane to divertor target
fluxexpi = poloidal mag. flux exp. on inboard from midplane to divertor target
pdiv= total power going to divertors
pdivrad= total power radiated in divertors
pdiv2= total power conducted to divertors
pdiv2o= power conducted to outboard divertor
pdiv2i= power conducted to inboard divertor
pdivrado= power radiated in outboard divertor
pdivradi= power radiated in inboard divertor
pdivtoto= conducted + radiated power in outboard divertor
pdivtoti= conducted + radiated power in inboard divertor
den= density factor in power scrape-off formula
xlam= power scrape-off layer width at plasma midplane
qpeakcondo= peak conducted power heat flux to outboard divertor
qpeakcondi= peak conducted power heat flux to inboard divertor
qpeakrado= peak radiated power heat flux to outboard divertor
qpeakradi= peak radiated power heat flux to inboard divertor
f_greenwald = greenwald fraction
Bt= toroidal magnetic field at plasma major radius
Aip = plasma current
Qmhd = qmhd

if(sym == 1.0) { facdiv = 0.65; } else { facdiv = 1.00; }
Rmaj = plasma_major_radius();
Asp = plasma_major_radius()/plasma_minor_radius();
adivo = 2.*pi*(Rmaj -((Rmaj/Asp)/2.0))*2.0*(Rmaj/(Asp*2.0));
adivi = 2.*pi*(Rmaj-(Rmaj/Asp))*2.0*(Rmaj/(Asp*4.0));
P_div = (P_plasma - P_rad_first_wall)/1.e6;
P_div_rad = f_div_rad*P_div;
P_div_2 = P_div - P_div_rad;
P_div_2_o = facdiv*fdivcondo*P_div_2;
P_div_2_i = facdiv*fdivcondi*P_div_2;
P_div_rad_o = facdiv*fdivrado*P_div_rad;
P_div_rad_i = facdiv*fdivradi*P_div_rad;
den = f_greenwald*(Aip/pi*pow(Rmaj/Asp, 2.0))*0.35;
xlam = 2.9e-2*2.5*pow(Qmhd, 0.75)*pow(den, 0.15)/(pow(P_div_2_o, 0.4)*Bt); // power scrape off width
Q_peak_cond_o = P_div_2_o/(2.0*pi*(Rmaj-((Rmaj/Asp)/2.0))*fluxexpo*xlam); Q_peak_cond_i = P_div_2_i/(2.0*pi*(Rmaj-(Rmaj/Asp)) *fluxexpi*xlam);
Q_peak_rad_o = P_div_rad_o/adivo;
Q_peak_rad_i = P_div_rad_i/adivi;
Q_peak_tot_o = Q_peak_cond_o + Q_peak_rad_o;
Q_peak_tot_i = Q_peak_cond_i + Q_peak_rad_i;
get_divertor_loads(Q_peak_tot_i*1.e6, Q_peak_tot_o*1.e6); // Get's peak divertor loads from a simple comparison
flux_divertor_Max = max(Q_peak_tot_o, Q_peak_tot_i)*1.e6;

Divertor Pumping Powers

This calculates the pumping power requirement for the divertor in use. The divertor in use is selected from the input file Setup.data. Note that the P_thermal_divertor is used, not the overall P_thermal. The output is filtered if the maximum heat flux on the divertor is less than 5 or greater than 15 MW/m2, which is the range of the tables.

[M. Tillack, X. Wang, "Divertor Pumping Powers - ppcompare2.xls" Dec. 2010]

P_thermal_divertor = P_alpha_loss *f_alpha_div

+ P_neutron *f_n_div *f_nmul_divertor

+ P_rad_edge *f_rad_edge

+ P_rad_chamb *f_rad_div

+ P_cond_div;

P_pump_divertor = pumping_power_divertor(flux_divertor_Max) *P_thermal_divertor/100.0;
P_pump_divertor_loss = (1.0 - eta_pump_he)*P_pump_divertor/eta_pump_he;

T-tube design: J. Burke’s T-Tube 2010 divertor design pumping power curve. Linear-taper with fixed slot width of 0.45 mm, inlet temp 600C, outlet 677C, P = 10 MPa; surface heat flux (MW/m^2) vs P_pump/P_thermal (%).

ppttube = 0.05594 * pow(flux_divertor_Max/1.0e6,2) + 0.04931 * (flux_divertor_Max/1.0e6) - 0.08702;

Plate design: X. Wang’s plate divertor design (HCFP). Helium-cooled flat plate, inlet temp 600C, outlet 677C, P = 10 MPa.

ppplate = 0.0392 * pow(flux_divertor_Max/1.0e6,3) - 0.743 * pow(flux_divertor_Max/1.0e6,2) + 5.44 * (flux_divertor_Max/1.0e6) - 11.1;

Combined plate-finger design: X. Wang’s finger divertor design(HCPF). Helium-cooled combined plate and finger, inlet temp 600C, outlet 700C, P = 10 MPa.

ppfinger = 0.0335 * pow(flux_divertor_Max/1.0e6,2) + 0.22 * (flux_divertor_Max/1.0e6) - 0.892;

Powers

These equations describe the net electric power and power losses. See also the Power Flow Diagram under Documentation. There is a check to ensure that the power generated is equal to the net electric power plus the unrecoverable power lost through heat, friction, or inefficiencies.

[This already existed in the code; L. Carlson added accounting for losses 8/2011.]

P_net_electric = P_gross_electric –

P_plasma_heating / eta_plasma_heating –

P_cd_generic / eta_cd_generic –

P_aux_func –

P_cryo –

P_pump_blanket / eta_pump_lm_mech –

P_pump_divertor / eta_pump_he;

P_generation = P_fusion*(0.8*1.1 + 0.2)/1.0e6;
P_losses = P_plasma_heating_loss +

P_cd_loss +

P_brayton_loss +

P_pump_blanket_loss +

P_pump_divertor_loss +

P_aux_func +

P_cryo; //+ LT_fwbs_loss, + LT_div_loss to be added

P_leftover = P_generation - (P_net_electric/1.0e6) - (P_losses/1.0e6);

Fusion Power Core (FPC)

The volume of the fusion power core is calculated for use in other areas of the code. The SiC and DCLL FPCs are slightly different. The volume is multiplied by 16x at the end since each component is built as a 1/16th sector.

[No specific reference.]

DCLL:

fpc_vol = ((*FirstWallInb).volume() +

(*FirstWallOutb).volume() +

(*BackWallInboard).volume() +

(*BackWallOutboard).volume() +

(*BlanketInboard).volume() +

(*BlanketOutboard).volume() +

(*RepHTShield)[0].volume() +

(*RepHTShield)[1].volume() +

(1.0+pen_shielding)*((*HTShieldDCLL).volume()) +

(*SkeletonRing).volume() +

(*VacuumVesselDCLL)[0].volume() +

(*VacuumVesselDCLL)[1].volume() +

(*VacuumVesselDCLL)[2].volume() +

(*VacuumVesselDCLL)[3].volume() +

DivVol + (*TFCoil).volume() +

PoloFCVolume + PoloFCVolumeSpares)*16.0;

SCLL:

fpc_vol = ((*FirstWallInb).volume() +

(*FirstWallOutb).volume() +

(*BlanketInboard).volume() +

(*BlanketOutboard).volume() +

(*BlanketII).volume() +

(*RepHTShield)[0].volume() +

(*RepHTShield)[1].volume() +

(1.0+pen_shielding)*((*HTShield)[0].volume() + (*HTShield)[1].volume()) +

(*VacuumVessel)[0].volume() +

(*VacuumVessel)[1].volume() +

DivVol + (*TFCoil).volume() +

PoloFCVolume + PoloFCVolumeSpares)*16.0;

Gross Cycle Efficiency of Brayton Cycle

This calculates the gross cycle efficiency of the Brayton power cycle. It takes the maximum neutron wall loading in MW/m2 at the first wall for different surface heat fluxes (q’’max) in MW/m2 and returns the gross cycle efficiency Q.

[R. Raffray, 8/2010 "Raffray - Pumping Power systems code input.xls"; polynomial fits added by L. Carlson ~5/2011.]

[M. Tillack concluded that we can continue to use R. Raffray's NWL vs cycle efficiency for the near term since our NWL_max and q" are low (~4.2 MW/m^2, ~0.25 MW/m^2 for ~AT). 11/2010.]

DCLL: choose polynomial fit depending on q’’

if (q2 < 0.250) { eta = 0; cout < "Maximum surface heat flux q'' is too low for the Brayton cycle tabulated interval " < q2 < "\n"; }
if (q2 >= 0.250 & q2 <= 0.500) { eta = 0.00013*pow(nwl2,3) - 0.00190*pow(nwl2,2) + 0.00087*nwl2 + 0.44578;
if (nwl2 < 2.0 & nwl2 > 7.5) cout < "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " < nwl2 < "\n"; }
if (q2 > 0.500 & q2 <= 0.750) { eta = 0.00002*pow(nwl2,3) - 0.00009*pow(nwl2,2) - 0.00814*nwl2 + 0.45770;
if (nwl2 < 2.4 & nwl2 > 7.5) cout < "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " < nwl2 < "\n"; }
if (q2 > 0.750 & q2 <= 1.000) { eta = 0.00018*pow(nwl2,3) - 0.00275*pow(nwl2,2) + 0.00599*nwl2 + 0.43195;
if (nwl2 < 4.5 & nwl2 > 7.5) cout < "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " < nwl2 < "\n"; }
if (q2 > 1.000 & q2 <= 1.250) { eta = -0.00031*pow(nwl2,2) + 0.00290*nwl2 + 0.38929;
if (nwl2 < 4.8 & nwl2 > 7.5) cout < "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " < nwl2 < "\n"; }
if (q2 > 0.750) { eta = 0; cout < "Maximum surface heat flux q'' (" < q2 < ") is too high for the Brayton cycle tabulated interval.\n";}

SCLL: choose polynomial fit depending on q’’

if (q2 < 0.250) { eta = 0; cout < "Maximum surface heat flux q'' is too low for the Brayton cycle tabulated interval " < q2 < "\n"; }
if (q2 >= 0.250 & q2 <= 0.375) { eta = -0.00347*nwl2 + 0.58859;
if (nwl2 < 1.5 & nwl2 > 6.0) cout < "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " < nwl2 < "\n"; }
if (q2 > 0.375 & q2 <= 0.500) { eta = -0.00361*nwl2 + 0.56934;
if (nwl2 < 1.8 & nwl2 > 6.0) cout < "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " < nwl2 < "\n"; }
if (q2 > 0.500 & q2 <= 0.625) { eta = -0.00362*nwl2 + 0.54828;
if (nwl2 < 2.4 & nwl2 > 6.0) cout < "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " < nwl2 < "\n"; }
if (q2 > 0.625 & q2 <= 0.750) { eta = -0.00363*nwl2 + 0.52605;
if (nwl2 < 3.0 & nwl2 > 6.0) cout < "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " < nwl2 < "\n"; }
if (q2 > 0.750) { eta = 0; cout < "Maximum surface heat flux q'' (" < q2 < ") is too high for the Brayton cycle tabulated interval.\n";}

Current Drive Parameters

The current drive parameters are called by reading part of the powerflow.data input file. This defines parameters such as efficiencies, frequencies, and unit costs of potential CD units such as lower hybrid (LH), neutral beam (NB), and fast wave (FW). The CD efficiencies are interpolated if the value is not directly found in the table. The frequency is set in the input file as well.

[R. Raffray, "Power Flow.doc” 5/2008, implemented into code by Z. Dragojlovic]

frequency[6] = {80.0e6, 158.0e6, 250.0e6, 800.0e6, 2500.0e6, 8000.0e6};
efficiency[6] = {0.84, 0.72, 0.65, 0.63, 0.48, 0.44 };
unit_cost[6] = {1.15, 1.15, 2.30, 2.30, 2.30, 2.87 };