From 5dec06fd83f1fa02b7e2170c65353449f0fba7bf Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Thu, 30 Jul 2026 16:13:12 +0100 Subject: [PATCH 1/3] Simplify water use --- process/models/water_use.py | 217 ++++++++++++++++++++---------------- 1 file changed, 120 insertions(+), 97 deletions(-) diff --git a/process/models/water_use.py b/process/models/water_use.py index bc603609a9..8aa551bd47 100644 --- a/process/models/water_use.py +++ b/process/models/water_use.py @@ -1,5 +1,7 @@ """Water use model for calculating water usage during plant operation.""" +from dataclasses import dataclass + import numpy as np from process.core import constants @@ -139,122 +141,90 @@ def cooling_water_body(self, wastetherm: float, output: bool): """ evapsum = 0.0e0 - for icool in range(1, 4): - if icool == 1: + icools = [ + # coefficients as per Brady et al. 1969: + CoolingWaterBodyCoeffs( + # wind function coefficients + a=2.47e0, + b=0e0, + c=0.12e0, + # fitted coefficients of heat loading + d=3061.331e0, + e=-48.810e0, + f=-78.559e0, + g=-291.820e0, + h=0.267e0, + i=-0.610e0, + j=33.497e0, # small pond as a cooling body # heat loading, MW/acre, based on estimations from US power plants - heatload = 0.35e0 - # coefficients as per Brady et al. 1969: + heatload=0.35e0, + ), + # coefficients as per Webster et al. 1995: + CoolingWaterBodyCoeffs( # wind function coefficients - a = 2.47e0 - b = 0e0 - c = 0.12e0 + a=1.04e0, + b=1.05e0, + c=0.0e0, # fitted coefficients of heat loading - d = 3061.331e0 - e = -48.810e0 - f = -78.559e0 - g = -291.820e0 - h = 0.267e0 - i = -0.610e0 - j = 33.497e0 - - elif icool == 2: + d=3876.843e0, + e=-49.071e0, + f=-295.246e0, + g=-327.935e0, + h=0.260e0, + i=10.528e0, + j=40.188e0, # large lake or reservoir as a cooling body # heat loading, MW/acre, based on estimations from US power plants - heatload = 0.10e0 - # coefficients as per Webster et al. 1995: + heatload=0.10e0, + ), + # coefficients as per Gulliver et al. 1986: + CoolingWaterBodyCoeffs( # wind function coefficients - a = 1.04e0 - b = 1.05e0 - c = 0.0e0 + a=2.96e0, + b=0.64e0, + c=0.0e0, # fitted coefficients of heat loading - d = 3876.843e0 - e = -49.071e0 - f = -295.246e0 - g = -327.935e0 - h = 0.260e0 - i = 10.528e0 - j = 40.188e0 - - elif icool == 3: + d=2565.009e0, + e=-43.636e0, + f=-93.834e0, + g=-203.767e0, + h=0.257e0, + i=2.408e0, + j=20.596e0, # stream or river as a cooling body # heat loading, MW/acre, based on estimations from US power plants - heatload = 0.20e0 - # coefficients as per Gulliver et al. 1986: - # wind function coefficients - a = 2.96e0 - b = 0.64e0 - c = 0.0e0 - # fitted coefficients of heat loading - d = 2565.009e0 - e = -43.636e0 - f = -93.834e0 - g = -203.767e0 - h = 0.257e0 - i = 2.408e0 - j = 20.596e0 + heatload=0.20e0, + ), + ] + for icool in icools: # Unfortunately, the source spreadsheet was from the US, so the fits for # water body heating due to heat loading and the cooling wind functions # are in non-metric units, hence the conversions required here. # Limitations: maximum wind speed of ~5 m/s; initial # self.data.water_use.watertemp < 25 degC - # convert self.data.water_use.windspeed to mph - windspeedmph = self.data.water_use.windspeed * 2.237e0 - - # convert heat loading into cal/(cm2.sec) - heatloadimp = heatload * 1000000.0e0 * 0.239e0 / 40469000.0e0 - - # estimate how heat loading will raise temperature, for this water body - heatratio = ( - d - + (e * self.data.water_use.watertemp) - + (f * windspeedmph) - + (g * heatload) - + (h * self.data.water_use.watertemp**2) - + (i * windspeedmph**2) - + (j * heatload**2) - ) - - # estimate resultant heated water temperature - watertempheated = self.data.water_use.watertemp + (heatloadimp * heatratio) - - # find wind function, m/(day.kPa), applicable to this water body: - windfunction = ( - a - + (b * self.data.water_use.windspeed) - + (c * self.data.water_use.windspeed**2) - ) / 1000.0e0 - - # difference in saturation vapour pressure (Clausius-Clapeyron approximation) - satvapdelta = ( - 0.611e0 - * np.exp((17.27e0 * watertempheated) / (237.3e0 + watertempheated)) - ) - ( - 0.611e0 - * np.exp( - (17.27e0 * self.data.water_use.watertemp) - / (237.3e0 + self.data.water_use.watertemp) + # find evaporation ratio: ratio of the heat used to evaporate water + # to the total heat discharged through the tower + self.data.water_use.evapratio = ( + icool.delta_e( + self.data.water_use.windspeed, + # estimate resultant heated water temperature + self.data.water_use.watertemp + + ( + icool.imp_heatload_conversion() + * icool.heat_ratio( + self.data.water_use.watertemp, self.data.water_use.windspeed + ) + ), + self.data.water_use.watertemp, + self.data.water_use.waterdens, + self.data.water_use.latentheat, ) + / icool.met_heatload_conversion() ) - - # find 'forced evaporation' driven by heat inserted into system - deltae = ( - self.data.water_use.waterdens - * self.data.water_use.latentheat - * windfunction - * satvapdelta - ) - - # convert heat loading to J/(m2.day) - heatloadmet = heatload * 1000000.0e0 / 4046.85642e0 * SECDAY - - # find evaporation ratio: ratio of the heat used to evaporate water - # to the total heat discharged through the tower - self.data.water_use.evapratio = deltae / heatloadmet # Diehl et al. USGS Report 2013-5188, http://dx.doi.org/10.3133/sir20135188 - self.data.water_use.volheat = ( self.data.water_use.waterdens * self.data.water_use.latentheat ) @@ -274,7 +244,7 @@ def cooling_water_body(self, wastetherm: float, output: bool): # stage of calculation evapsum += self.data.water_use.evapvol - evapsum /= icool + evapsum /= len(icools) # water volume withdrawn from external source depends on recirculation or # 'once-through' system choice. Estimated as a ratio to evaporated water @@ -305,3 +275,56 @@ def cooling_water_body(self, wastetherm: float, output: bool): self.data.water_use.wateruseonethru, "OP ", ) + + +@dataclass +class CoolingWaterBodyCoeffs: + """Cooling water body coefficients""" + + a: float + b: float + c: float + d: float + e: float + f: float + g: float + h: float + i: float + j: float + heatload: float + + def heat_ratio(self, watertemp, windspeed): + """Estimate how heat loading will raise temperature, for this water body""" + # convert self.data.water_use.windspeed to mph + windspeedmph = windspeed * 2.237e0 + return ( + self.d + + (self.e * watertemp) + + (self.f * windspeedmph) + + (self.g * self.heatload) + + (self.h * watertemp**2) + + (self.i * windspeedmph**2) + + (self.j * self.heatload**2) + ) + + def imp_heatload_conversion(self): + """Convert heat loading into cal/(cm2.sec)""" + return self.heatload * 1000000.0e0 * 0.239e0 / 40469000.0e0 + + def met_heatload_conversion(self): + """Convert heat loading to J/(m2.day)""" + return self.heatload * 1000000.0e0 / 4046.85642e0 * SECDAY + + def delta_e(self, windspeed, watertempheated, watertemp, waterdens, latentheat): + """Find 'forced evaporation' driven by heat inserted into system""" + # find wind function, m/(day.kPa), applicable to this water body: + windfunction = ( + self.a + (self.b * windspeed) + (self.c * windspeed**2) + ) / 1000.0e0 + + # difference in saturation vapour pressure (Clausius-Clapeyron approximation) + satvapdelta = ( + 0.611e0 * np.exp((17.27e0 * watertempheated) / (237.3e0 + watertempheated)) + ) - (0.611e0 * np.exp((17.27e0 * watertemp) / (237.3e0 + watertemp))) + + return waterdens * latentheat * windfunction * satvapdelta From 9b748f08919681004727c25b8e42cf89c1c2155b Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Thu, 30 Jul 2026 16:33:04 +0100 Subject: [PATCH 2/3] Some vacuum improvements --- process/models/vacuum.py | 559 ++++++++++++++++++++------------------- 1 file changed, 285 insertions(+), 274 deletions(-) diff --git a/process/models/vacuum.py b/process/models/vacuum.py index 5fc3c998e5..feb8936e2a 100644 --- a/process/models/vacuum.py +++ b/process/models/vacuum.py @@ -282,34 +282,28 @@ def vacuum( # speed of 2.0 m^3/s (1.95 for N2, 1.8 for He, 1.8 for DT) # i_vacuum_pump_type = 1 for compound cryopump with nominal speed of 10 m^3/s # (9.0 for N2, 5.0 for He and 25. for DT) - pfus = pfusmw * 1.0e6 # Fusion power (W) ntf = int(n_tf_coils) # Feed rate (gas load) of D-T into chamber (pellets + gas puffing + # NBI + ...) = load from fueller + load from NBI # frate (kg/s) = gasld (kg/s) + qtorus (D2/s) * 6.64e-27 (kg/D2) - frate = gasld + qtorus * 6.64e-27 # Set duct shield thickness to zero for no biological shielding # instead of thshldo/3.0e0 - thdsh = 0.0e0 # Shielding (m) between duct and TF coils is scaled from inboard shield # thickness - thcsh = thshldi / 3.0e0 # Multiplier to convert conductance from gas species i to nitrogen xmult = [1.0e0, 0.423e0, 0.378e0, 0.423e0] # nitrogen, D-T, helium, D-T again - nduct = ntf * ndiv # Speed of high-vacuum pumps (m^3/s) - # nitrogen, DT, helium, DT again sp = ( [1.95, 1.8, 1.8, 1.8] @@ -319,7 +313,6 @@ def vacuum( ) # Calculate required pumping speeds - s = [] # Initial pumpdown based on outgassing @@ -331,9 +324,7 @@ def vacuum( # Old method: area = 4.0e0 * pi*pi * r0 * aw # * sqrt(0.5e0*(1.0e0 + kappa*kappa)) - area = plasma_sarea * (aw + dsol) / aw - ogas = self.data.vacuum.outgrat_fw * area * 10.0e0 # Outgassing rate (Pa-m^3/s) s.append(ogas / self.data.vacuum.pres_vv_chamber_base) @@ -342,15 +333,13 @@ def vacuum( # temp_vv_chamber_gas_burn_end = temperature of neutral gas in chamber (K) # t_plant_pulse_dwell = dwell time between burns (s) - pend = ( - 0.5e0 * nplasma * k * self.data.vacuum.temp_vv_chamber_gas_burn_end - ) # pressure in plasma chamber after burn (Pa) - pstart = 0.01e0 * pend # pressure in chamber before start of burn (Pa) + # pressure in plasma chamber after burn (Pa) + pend = 0.5e0 * nplasma * k * self.data.vacuum.temp_vv_chamber_gas_burn_end + # pressure in chamber before start of burn (Pa) + pstart = 0.01e0 * pend # Chamber volume (m^3) - # Old method: volume = 2.0e0 * pi*pi * r0 * aw*aw * kappa - volume = plasma_vol * (aw + dsol) * (aw + dsol) / (aw * aw) # dwell pumping options @@ -368,7 +357,6 @@ def vacuum( # source = alpha production rate (pa - m^3/s) # fhe = fraction of neutral gas in divertor chamber that is helium # pres_div_chamber_burn = pressure in divertor chamber during burn (Pa) - source = pfus * 1.47e-09 fhe = source / (frate * 4.985e5) s.extend( @@ -384,13 +372,9 @@ def vacuum( ) # Calculate conductance of a single duct - imax = 1 cmax = 0.01e0 pumpn = 1.0e0 - nflag = 0 # Control option if ducts are too small in x-sectional area - # = 1 if problem is identified in output, but run continues - # = 0 otherwise l1 = thshldo + thtf # Length of passage from divertor to ducts (m) l2 = thshldo + 4.0e0 # Length of ducts from divertor passage to elbow (m) @@ -415,87 +399,15 @@ def vacuum( pumpn = max(pumpn, pumpn1, pumpn2) ceff[i] = 1.0e0 / (nduct / s[i] - 1.0e0 / (sp[i] * pumpn)) - # Newton's method solution for duct diameter - while True: - d[i] = 1.0e0 - - for _ in range(100): - a1 = ( - 0.25e0 * math.pi * d[i] * d[i] - ) # Area of aperture and duct (m^2) - a2 = 1.44e0 * a1 - a3 = a2 - k1 = 4.0e0 / 3.0e0 * d[i] / (l1 + 4.0e0 / 3.0e0 * d[i]) - k2 = ( - 4.0e0 - / 3.0e0 - * d[i] - * 1.2e0 - / (l2 + 4.0e0 / 3.0e0 * d[i] * 1.2e0) - ) - k3 = ( - 4.0e0 - / 3.0e0 - * d[i] - * 1.2e0 - / (l3 + 4.0e0 / 3.0e0 * d[i] * 1.2e0) - ) - cap = 119.0e0 * a1 / xmult[i] - dcap = 2.0e0 * cap / d[i] - c1 = 119.0e0 * a1 * k1 / xmult[i] - dc1 = c1 / d[i] * (3.0e0 - k1) - c2 = 119.0e0 * a2 * k2 / xmult[i] - dc2 = c2 / d[i] / 1.2e0 * (3.0e0 - k2) - c3 = 119.0e0 * a3 * k3 / xmult[i] - dc3 = c3 / d[i] / 1.2e0 * (3.0e0 - k3) - cnew = 1.0e0 / (1.0e0 / cap + 1.0e0 / c1 + 1.0e0 / c2 + 1.0e0 / c3) - y = -ceff[i] + cnew - dy = ( - cnew - * cnew - * ( - dcap / cap / cap - + dc1 / c1 / c1 - + dc2 / c2 / c2 - + dc3 / c3 / c3 - ) - ) - dnew = d[i] - y / dy - dd = abs((d[i] - dnew) / d[i]) - d[i] = dnew - if dd <= 0.01e0: - break - - else: - logger.error( - f"Newton's method not converging; check fusion power, te " - f"{self.data.physics.p_fusion_total_mw=} " - f"{self.data.physics.temp_plasma_electron_vol_avg_kev=}" - ) - - theta = math.pi / ntf - - # Area between adjacent TF coils available for pump ducts - # ritf = outer radius of inboard leg of TF coil (m) - - a1max = (r0 + aw - ritf - thcsh / math.tan(theta)) ** 2 * math.tan(theta) - d1max = math.sqrt(4.0e0 * a1max / math.pi) # Equivalent diameter - if a1 < a1max: - break - - ceff[i] *= 0.9e0 - if ceff[i] <= (1.1e0 * s[i]): - # Ducts are not big enough. Flag and continue. - nflag = 1 - break - + ceff, nflag, d1max = self._newton_method_duct_diameter( + d, i, s, xmult, l1, l2, l3, ntf, r0, aw, ritf, thcsh, ceff + ) cmax = ceff[i] pumpn *= nduct # d[imax]= diameter of passage from divertor to pumping ducts (m) # dout = diameter of ducts from passage to hi-vac pumps (m) - dout = d[imax] * 1.2e0 # Net pumping speeds provided by vacuum pumping system @@ -503,17 +415,14 @@ def vacuum( # snet(2) - net pump speed (D-T) provided (m^3/s) # snet(3) - net pump speed (He) provided (m^3/s) # snet(4) - snet(2) - snet = [] - for i in range(4): - ceff1 = ceff[imax] * nduct - snet.append( - 1.0e0 - / (1.0e0 / (ceff1 * xmult[imax] / xmult[i]) + 1.0e0 / sp[i] / pumpn) - ) + ceff1 = ceff[imax] * nduct + snet = [ + 1 / (1 / (ceff1 * xmult[imax] / xmult[i]) + 1 / sp[i] / pumpn) + for i in range(4) + ] # If cryopumps are used then an additional pump is required # for continuous operation with regeneration. - if ( VacuumPumpType(self.data.vacuum.i_vacuum_pump_type) == VacuumPumpType.COMPOUND_CRYOPUMP @@ -521,11 +430,9 @@ def vacuum( pumpn *= 2.0e0 # Information for costing routine - dlscalc = l1 * d[imax] ** 1.4e0 + (ltot - l1) * (d[imax] * 1.2e0) ** 1.4e0 # Mass of duct shielding - arsh = ( 0.25e0 * math.pi * ((d[imax] * 1.2e0 + thdsh) ** 2 - (d[imax] * 1.2e0) ** 2) ) @@ -534,195 +441,299 @@ def vacuum( dimax = d[imax] if output: - # Output section - - process_output.oheadr(self.outfile, "Vacuum System") - - process_output.ocmmnt(self.outfile, "Pumpdown to Base Pressure :") - process_output.oblnkl(self.outfile) - process_output.ovarre( - self.outfile, - "First wall outgassing rate (Pa m/s)", - "(outgrat_fw)", - self.data.vacuum.outgrat_fw, - ) - process_output.ovarre( - self.outfile, "Total outgassing load (Pa m3/s)", "(ogas)", ogas, "OP " - ) - process_output.ovarre( - self.outfile, - "Base pressure required (Pa)", - "(pres_vv_chamber_base)", - self.data.vacuum.pres_vv_chamber_base, - ) - process_output.ovarre( - self.outfile, "Required N2 pump speed (m3/s)", "(s(1))", s[0], "OP " - ) - process_output.ovarre( - self.outfile, - "N2 pump speed provided (m3/s)", - "(snet(1))", - snet[0], - "OP ", - ) - - process_output.osubhd(self.outfile, "Pumpdown between Burns :") - process_output.ovarre( - self.outfile, "Plasma chamber volume (m3)", "(volume)", volume, "OP " - ) - process_output.ovarre( - self.outfile, "Chamber pressure after burn (Pa)", "(pend)", pend, "OP " - ) - process_output.ovarre( - self.outfile, "Chamber pressure before burn (Pa)", "(pstart)", pstart - ) - process_output.ovarre( - self.outfile, - "Allowable pumping time switch", - "(i_vac_pump_dwell)", - self.data.vacuum.i_vac_pump_dwell, - ) - process_output.ovarre( - self.outfile, - "Dwell time between burns (s)", - "(t_plant_pulse_dwell.)", + self._write_to_outfile( + ogas, + s, + snet, + volume, + pend, + pstart, t_plant_pulse_dwell, - ) - process_output.ovarre( - self.outfile, - "CS ramp-up time burns (s)", - "(t_plant_pulse_coil_precharge.)", - self.data.times.t_plant_pulse_coil_precharge, - ) - process_output.ovarre( - self.outfile, - "Allowable pumping time between burns (s)", - "(tpump)", tpump, - ) - process_output.ovarre( - self.outfile, "Required D-T pump speed (m3/s)", "(s(2))", s[1], "OP " - ) - process_output.ovarre( - self.outfile, - "D-T pump speed provided (m3/s)", - "(snet(2))", - snet[1], - "OP ", - ) - - process_output.osubhd(self.outfile, "Helium Ash Removal :") - process_output.ovarre( - self.outfile, - "Divertor chamber gas pressure (Pa)", - "(pres_div_chamber_burn)", - self.data.vacuum.pres_div_chamber_burn, - ) - process_output.ovarre( - self.outfile, - "Helium gas fraction in divertor chamber", - "(fhe)", fhe, - "OP ", - ) - process_output.ovarre( - self.outfile, "Required helium pump speed (m3/s)", "(s(3))", s[2], "OP " - ) - process_output.ovarre( - self.outfile, - "Helium pump speed provided (m3/s)", - "(snet(3))", - snet[2], - "OP ", + frate, + nflag, + d1max, + nduct, + dimax, + imax, + l1, + l2, + l3, + dout, + pumpn, ) - process_output.osubhd(self.outfile, "D-T Removal at Fuelling Rate :") - process_output.ovarre( - self.outfile, "D-T fuelling rate (kg/s)", "(frate)", frate, "OP " - ) - process_output.ovarre( - self.outfile, "Required D-T pump speed (m3/s)", "(s(4))", s[3], "OP " - ) - process_output.ovarre( - self.outfile, - "D-T pump speed provided (m3/s)", - "(snet(4))", - snet[3], - "OP ", - ) + return pumpn, nduct, dlscalc, mvdsh, dimax - if nflag == 1: - process_output.oblnkl(self.outfile) - process_output.ocmmnt( - self.outfile, "Vacuum pumping ducts are space limited." - ) - process_output.ocmmnt( - self.outfile, f"Maximum duct diameter is only {d1max} m" + def _newton_method_duct_diameter( + self, d, i, s, xmult, l1, l2, l3, ntf, r0, aw, ritf, thcsh, ceff + ): + nflag = 0 # Control option if ducts are too small in x-sectional area + # = 1 if problem is identified in output, but run continues + # = 0 otherwise + # Newton's method solution for duct diameter + while True: + d[i] = 1.0e0 + + for _ in range(100): + a1 = 0.25e0 * math.pi * d[i] * d[i] # Area of aperture and duct (m^2) + a2 = 1.44e0 * a1 + a3 = a2 + k1 = 4.0e0 / 3.0e0 * d[i] / (l1 + 4.0e0 / 3.0e0 * d[i]) + k2 = 4.0e0 / 3.0e0 * d[i] * 1.2e0 / (l2 + 4.0e0 / 3.0e0 * d[i] * 1.2e0) + k3 = 4.0e0 / 3.0e0 * d[i] * 1.2e0 / (l3 + 4.0e0 / 3.0e0 * d[i] * 1.2e0) + cap = 119.0e0 * a1 / xmult[i] + dcap = 2.0e0 * cap / d[i] + c1 = 119.0e0 * a1 * k1 / xmult[i] + dc1 = c1 / d[i] * (3.0e0 - k1) + c2 = 119.0e0 * a2 * k2 / xmult[i] + dc2 = c2 / d[i] / 1.2e0 * (3.0e0 - k2) + c3 = 119.0e0 * a3 * k3 / xmult[i] + dc3 = c3 / d[i] / 1.2e0 * (3.0e0 - k3) + cnew = 1.0e0 / (1.0e0 / cap + 1.0e0 / c1 + 1.0e0 / c2 + 1.0e0 / c3) + y = -ceff[i] + cnew + dy = ( + cnew + * cnew + * (dcap / cap / cap + dc1 / c1 / c1 + dc2 / c2 / c2 + dc3 / c3 / c3) ) - process_output.ocmmnt(self.outfile, "Conductance is inadequate.") - process_output.oblnkl(self.outfile) - - i_fw_blkt_shared_coolant = ( - "cryo " - if VacuumPumpType(self.data.vacuum.i_vacuum_pump_type) - == VacuumPumpType.COMPOUND_CRYOPUMP - else "turbo" - ) - - process_output.oblnkl(self.outfile) - process_output.ocmmnt( - self.outfile, "The vacuum pumping system size is governed by the" - ) + dnew = d[i] - y / dy + dd = abs((d[i] - dnew) / d[i]) + d[i] = dnew + if dd <= 0.01e0: + break - if imax == 1: - process_output.ocmmnt( - self.outfile, "requirements for pumpdown to base pressure." - ) - elif imax == 2: - process_output.ocmmnt( - self.outfile, "requirements for pumpdown between burns." - ) - elif imax == 3: - process_output.ocmmnt( - self.outfile, "requirements for helium ash removal." - ) else: - process_output.ocmmnt( - self.outfile, "requirements for D-T removal at fuelling rate." + logger.error( + f"Newton's method not converging; check fusion power, te " + f"{self.data.physics.p_fusion_total_mw=} " + f"{self.data.physics.temp_plasma_electron_vol_avg_kev=}" ) + theta = math.pi / ntf + + # Area between adjacent TF coils available for pump ducts + # ritf = outer radius of inboard leg of TF coil (m) + + a1max = (r0 + aw - ritf - thcsh / math.tan(theta)) ** 2 * math.tan(theta) + d1max = math.sqrt(4.0e0 * a1max / math.pi) # Equivalent diameter + if a1 < a1max: + break + + ceff[i] *= 0.9e0 + if ceff[i] <= (1.1e0 * s[i]): + # Ducts are not big enough. Flag and continue. + nflag = 1 + break + return ceff, nflag, d1max + + def _write_to_outfile( + self, + ogas, + s, + snet, + volume, + pend, + pstart, + t_plant_pulse_dwell, + tpump, + fhe, + frate, + nflag, + d1max, + nduct, + dimax, + imax, + l1, + l2, + l3, + dout, + pumpn, + ): + # Output section + + process_output.oheadr(self.outfile, "Vacuum System") + + process_output.ocmmnt(self.outfile, "Pumpdown to Base Pressure :") + process_output.oblnkl(self.outfile) + process_output.ovarre( + self.outfile, + "First wall outgassing rate (Pa m/s)", + "(outgrat_fw)", + self.data.vacuum.outgrat_fw, + ) + process_output.ovarre( + self.outfile, "Total outgassing load (Pa m3/s)", "(ogas)", ogas, "OP " + ) + process_output.ovarre( + self.outfile, + "Base pressure required (Pa)", + "(pres_vv_chamber_base)", + self.data.vacuum.pres_vv_chamber_base, + ) + process_output.ovarre( + self.outfile, "Required N2 pump speed (m3/s)", "(s(1))", s[0], "OP " + ) + process_output.ovarre( + self.outfile, + "N2 pump speed provided (m3/s)", + "(snet(1))", + snet[0], + "OP ", + ) + + process_output.osubhd(self.outfile, "Pumpdown between Burns :") + process_output.ovarre( + self.outfile, "Plasma chamber volume (m3)", "(volume)", volume, "OP " + ) + process_output.ovarre( + self.outfile, "Chamber pressure after burn (Pa)", "(pend)", pend, "OP " + ) + process_output.ovarre( + self.outfile, "Chamber pressure before burn (Pa)", "(pstart)", pstart + ) + process_output.ovarre( + self.outfile, + "Allowable pumping time switch", + "(i_vac_pump_dwell)", + self.data.vacuum.i_vac_pump_dwell, + ) + process_output.ovarre( + self.outfile, + "Dwell time between burns (s)", + "(t_plant_pulse_dwell.)", + t_plant_pulse_dwell, + ) + process_output.ovarre( + self.outfile, + "CS ramp-up time burns (s)", + "(t_plant_pulse_coil_precharge.)", + self.data.times.t_plant_pulse_coil_precharge, + ) + process_output.ovarre( + self.outfile, + "Allowable pumping time between burns (s)", + "(tpump)", + tpump, + ) + process_output.ovarre( + self.outfile, "Required D-T pump speed (m3/s)", "(s(2))", s[1], "OP " + ) + process_output.ovarre( + self.outfile, + "D-T pump speed provided (m3/s)", + "(snet(2))", + snet[1], + "OP ", + ) + + process_output.osubhd(self.outfile, "Helium Ash Removal :") + process_output.ovarre( + self.outfile, + "Divertor chamber gas pressure (Pa)", + "(pres_div_chamber_burn)", + self.data.vacuum.pres_div_chamber_burn, + ) + process_output.ovarre( + self.outfile, + "Helium gas fraction in divertor chamber", + "(fhe)", + fhe, + "OP ", + ) + process_output.ovarre( + self.outfile, "Required helium pump speed (m3/s)", "(s(3))", s[2], "OP " + ) + process_output.ovarre( + self.outfile, + "Helium pump speed provided (m3/s)", + "(snet(3))", + snet[2], + "OP ", + ) + + process_output.osubhd(self.outfile, "D-T Removal at Fuelling Rate :") + process_output.ovarre( + self.outfile, "D-T fuelling rate (kg/s)", "(frate)", frate, "OP " + ) + process_output.ovarre( + self.outfile, "Required D-T pump speed (m3/s)", "(s(4))", s[3], "OP " + ) + process_output.ovarre( + self.outfile, + "D-T pump speed provided (m3/s)", + "(snet(4))", + snet[3], + "OP ", + ) + + if nflag == 1: process_output.oblnkl(self.outfile) - process_output.ovarre( - self.outfile, "Number of large pump ducts", "(nduct)", nduct - ) - process_output.ovarre( - self.outfile, - "Passage diameter, divertor to ducts (m)", - "(d(imax))", - d[imax], - "OP ", + process_output.ocmmnt( + self.outfile, "Vacuum pumping ducts are space limited." ) - process_output.ovarre(self.outfile, "Passage length (m)", "(l1)", l1, "OP ") - process_output.ovarre( - self.outfile, "Diameter of ducts (m)", "(dout)", dout, "OP " + process_output.ocmmnt( + self.outfile, f"Maximum duct diameter is only {d1max} m" ) + process_output.ocmmnt(self.outfile, "Conductance is inadequate.") + process_output.oblnkl(self.outfile) - process_output.ovarre( - self.outfile, "Duct length, divertor to elbow (m)", "(l2)", l2, "OP " - ) - process_output.ovarre( - self.outfile, "Duct length, elbow to pumps (m)", "(l3)", l3 + i_fw_blkt_shared_coolant = ( + "cryo " + if VacuumPumpType(self.data.vacuum.i_vacuum_pump_type) + == VacuumPumpType.COMPOUND_CRYOPUMP + else "turbo" + ) + + process_output.oblnkl(self.outfile) + process_output.ocmmnt( + self.outfile, "The vacuum pumping system size is governed by the" + ) + + if imax == 1: + process_output.ocmmnt( + self.outfile, "requirements for pumpdown to base pressure." ) - process_output.ovarre( - self.outfile, "Number of pumps", "(pumpn)", pumpn, "OP " + elif imax == 2: + process_output.ocmmnt( + self.outfile, "requirements for pumpdown between burns." ) - process_output.oblnkl(self.outfile) + elif imax == 3: + process_output.ocmmnt(self.outfile, "requirements for helium ash removal.") + else: process_output.ocmmnt( - self.outfile, - f"The vacuum system uses {i_fw_blkt_shared_coolant} pumps.", + self.outfile, "requirements for D-T removal at fuelling rate." ) - return pumpn, nduct, dlscalc, mvdsh, dimax + process_output.oblnkl(self.outfile) + process_output.ovarre( + self.outfile, "Number of large pump ducts", "(nduct)", nduct + ) + process_output.ovarre( + self.outfile, + "Passage diameter, divertor to ducts (m)", + "(d(imax))", + dimax, + "OP ", + ) + process_output.ovarre(self.outfile, "Passage length (m)", "(l1)", l1, "OP ") + process_output.ovarre( + self.outfile, "Diameter of ducts (m)", "(dout)", dout, "OP " + ) + + process_output.ovarre( + self.outfile, "Duct length, divertor to elbow (m)", "(l2)", l2, "OP " + ) + process_output.ovarre( + self.outfile, "Duct length, elbow to pumps (m)", "(l3)", l3 + ) + process_output.ovarre(self.outfile, "Number of pumps", "(pumpn)", pumpn, "OP ") + process_output.oblnkl(self.outfile) + process_output.ocmmnt( + self.outfile, + f"The vacuum system uses {i_fw_blkt_shared_coolant} pumps.", + ) class VacuumVessel(Model): From 91682c75ac7c1e1cf3751b79a707a203c98f9e0c Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Thu, 30 Jul 2026 16:44:30 +0100 Subject: [PATCH 3/3] factor out calculation in newton --- process/models/vacuum.py | 58 +++++++++++++++++++++------------------- 1 file changed, 31 insertions(+), 27 deletions(-) diff --git a/process/models/vacuum.py b/process/models/vacuum.py index feb8936e2a..8d5bfed60a 100644 --- a/process/models/vacuum.py +++ b/process/models/vacuum.py @@ -473,35 +473,14 @@ def _newton_method_duct_diameter( # = 1 if problem is identified in output, but run continues # = 0 otherwise # Newton's method solution for duct diameter + while True: d[i] = 1.0e0 - for _ in range(100): - a1 = 0.25e0 * math.pi * d[i] * d[i] # Area of aperture and duct (m^2) - a2 = 1.44e0 * a1 - a3 = a2 - k1 = 4.0e0 / 3.0e0 * d[i] / (l1 + 4.0e0 / 3.0e0 * d[i]) - k2 = 4.0e0 / 3.0e0 * d[i] * 1.2e0 / (l2 + 4.0e0 / 3.0e0 * d[i] * 1.2e0) - k3 = 4.0e0 / 3.0e0 * d[i] * 1.2e0 / (l3 + 4.0e0 / 3.0e0 * d[i] * 1.2e0) - cap = 119.0e0 * a1 / xmult[i] - dcap = 2.0e0 * cap / d[i] - c1 = 119.0e0 * a1 * k1 / xmult[i] - dc1 = c1 / d[i] * (3.0e0 - k1) - c2 = 119.0e0 * a2 * k2 / xmult[i] - dc2 = c2 / d[i] / 1.2e0 * (3.0e0 - k2) - c3 = 119.0e0 * a3 * k3 / xmult[i] - dc3 = c3 / d[i] / 1.2e0 * (3.0e0 - k3) - cnew = 1.0e0 / (1.0e0 / cap + 1.0e0 / c1 + 1.0e0 / c2 + 1.0e0 / c3) - y = -ceff[i] + cnew - dy = ( - cnew - * cnew - * (dcap / cap / cap + dc1 / c1 / c1 + dc2 / c2 / c2 + dc3 / c3 / c3) - ) - dnew = d[i] - y / dy + dnew, a1 = self._newton_function(d[i], l1, l2, l3, xmult, i, ceff) dd = abs((d[i] - dnew) / d[i]) d[i] = dnew - if dd <= 0.01e0: + if dd <= 0.01: break else: @@ -517,17 +496,42 @@ def _newton_method_duct_diameter( # ritf = outer radius of inboard leg of TF coil (m) a1max = (r0 + aw - ritf - thcsh / math.tan(theta)) ** 2 * math.tan(theta) - d1max = math.sqrt(4.0e0 * a1max / math.pi) # Equivalent diameter + d1max = math.sqrt(4.0 * a1max / math.pi) # Equivalent diameter if a1 < a1max: break - ceff[i] *= 0.9e0 - if ceff[i] <= (1.1e0 * s[i]): + ceff[i] *= 0.9 + if ceff[i] <= (1.1 * s[i]): # Ducts are not big enough. Flag and continue. nflag = 1 break return ceff, nflag, d1max + @staticmethod + def _newton_function(d_i, l1, l2, l3, xmult, i, ceff): + a1 = 0.25 * math.pi * d_i * d_i # Area of aperture and duct (m^2) + a2 = 1.44 * a1 + a3 = a2 + k1 = 4 / 3 * d_i / (l1 + 4 / 3 * d_i) + k2 = 4 / 3 * d_i * 1.2 / (l2 + 4 / 3 * d_i * 1.2) + k3 = 4 / 3 * d_i * 1.2 / (l3 + 4 / 3 * d_i * 1.2) + cap = 119 * a1 / xmult[i] + dcap = 2 * cap / d_i + c1 = 119 * a1 * k1 / xmult[i] + dc1 = c1 / d_i * (3 - k1) + c2 = 119 * a2 * k2 / xmult[i] + dc2 = c2 / d_i / 1.2 * (3 - k2) + c3 = 119 * a3 * k3 / xmult[i] + dc3 = c3 / d_i / 1.2 * (3 - k3) + cnew = 1 / (1 / cap + 1 / c1 + 1 / c2 + 1 / c3) + y = -ceff[i] + cnew + dy = ( + cnew + * cnew + * (dcap / cap / cap + dc1 / c1 / c1 + dc2 / c2 / c2 + dc3 / c3 / c3) + ) + return d_i - y / dy, a1 + def _write_to_outfile( self, ogas,