Files

612 lines
22 KiB
Python

# -*- coding: utf-8 -*-
"""Component models used by Brayton cycle simulations."""
import math
import ctREFPROP.ctREFPROP as ct
class ComponentMassMixin:
"""Shared mass-calculation interface for cycle components."""
component_type = "component"
def _init_mass_interface(self):
self.mass = None
self.mass_variables = None
def calculate_mass(self, mass_model=None, **kwargs):
"""Calculate and store component mass.
The future mass model can be provided as either a callable accepting
``(component, **kwargs)`` or an object exposing ``calculate_mass`` or
``calculate_<component_type>_mass``.
"""
if self.variables is None:
raise ValueError("Run component calculator before mass calculation")
if mass_model is None:
result = self._calculate_mass(**kwargs)
else:
result = self._run_external_mass_model(mass_model, **kwargs)
return self._store_mass_result(result)
def mass_calculator(self, mass_model=None, **kwargs):
return self.calculate_mass(mass_model=mass_model, **kwargs)
def _calculate_mass(self, **kwargs):
raise NotImplementedError(
f"{self.__class__.__name__} mass model is not implemented yet. "
"Pass a mass_model or override _calculate_mass()."
)
def _run_external_mass_model(self, mass_model, **kwargs):
method_name = f"calculate_{self.component_type}_mass"
if hasattr(mass_model, method_name):
return getattr(mass_model, method_name)(self, **kwargs)
if hasattr(mass_model, "calculate_mass"):
return mass_model.calculate_mass(self, **kwargs)
if callable(mass_model):
return mass_model(self, **kwargs)
raise TypeError("mass_model must be callable or expose a supported method")
def _store_mass_result(self, result):
if isinstance(result, dict):
if "mass" in result:
mass = result["mass"]
elif "total_mass" in result:
mass = result["total_mass"]
else:
raise ValueError("Mass result dictionary must include 'mass'")
mass_variables = dict(result)
else:
mass = result
mass_variables = {"mass": result}
self.mass = mass
self.mass_variables = mass_variables
self.variables["mass"] = self.mass
self.variables["mass_variables"] = self.mass_variables
return self.mass
class Compressor():
"""压缩机类"""
def __init__(self, name, eff):
"""
初始化参数
name: 名称
eff: 等熵效率
Wc: 压缩功
inlet_state: 入口参数
outlet_state: 出口参数
outlet_state_is: 等熵状态下出口参数
"""
self.name = name
self.eff = eff
self.variables = None
def calculator(self, p_in, T_in, p_out, property_calculator):
# 先计算熵值
inlet_state = property_calculator.calculate_properties(T=T_in, P=p_in)
mw = property_calculator.mw
s = inlet_state['s']
h_in = inlet_state['h']
outlet_state_is = property_calculator.calculate_properties(P=p_out, s=s)
h_out_is = outlet_state_is['h']
h_out = h_in + (h_out_is - h_in) / self.eff
outlet_state = property_calculator.calculate_properties(P=p_out, h=h_out)
Wc = h_out - h_in
T_out = outlet_state['T']
# 计算结果
self.variables = {
'name': self.name,
'inlet_state':{
'P': p_in, # 压强(kPa)
'T': T_in,
'h': h_in/mw, # 比焓(J/mol)->(J/kg)
's': s/mw, # 比熵(J/mol.K)->(J/kg.K)
},
'outlet_state':{
'P': p_out, # 压强(kPa)
'T': T_out,
'h': h_out/mw, # 比焓(kJ/mol)->(kJ/kg)
's': s/mw, # 比熵(kJ/mol.K)->(kJ/kg.K)
},
'eff': self.eff,
'Wc': Wc/mw, # 压缩功(kJ/mol)->(kJ/kg)
'pi': p_out / p_in
}
class Turbine(ComponentMassMixin):
"""透平类"""
component_type = "turbine"
def __init__(self, name, eff):
"""
初始化参数
name: 名称
eff: 透平效率
"""
self.name = name
self.eff = eff
self.variables = None
self._init_mass_interface()
def calculator(self, p_in, T_in, p_out, property_calculator):
"""涡轮参数计算"""
inlet_state = property_calculator.calculate_properties(P=p_in, T=T_in)
mw = property_calculator.mw
s = inlet_state['s']
h_in = inlet_state['h']
outlet_state_is = property_calculator.calculate_properties(P=p_out, s=s)
h_out_is = outlet_state_is['h']
h_out = h_in - (h_in - h_out_is) * self.eff
outlet_state = property_calculator.calculate_properties(P=p_out, h=h_out)
T_out = outlet_state['T']
Wt = h_in - h_out
self.variables = {
'name': self.name,
'inlet_state':{
'P': p_in, # 压强(kPa)
'T': T_in,
'h': h_in/mw, # 比焓(J/mol)->(J/kg)
's': s/mw, # 比熵(J/mol.K)->(J/kg.K)
},
'outlet_state':{
'P': p_out, # 压强(kPa)
'T': T_out,
'h': h_out/mw, # 比焓(J/mol)->(J/kg)
's': s/mw, # 比熵(J/mol.K)->(J/kg.K)
},
'eff': self.eff,
'Wt': Wt/mw, # 透平做功(J/mol)->(J/kg)
'pi': p_in / p_out
}
def _calculate_mass(self, Pe=None, A=0.5, **kwargs):
"""Estimate turbine/TAC mass with the empirical TAC formula.
Formula from the provided reference:
M = A * sqrt(0.5*pi*(30.522*ln(Pe) - 5.7178))
If Pe is not supplied, the current turbine specific work ``Wt`` is used
directly as Pe without unit conversion.
"""
if not 0.4 <= A <= 0.8:
raise ValueError("A should be within the recommended range 0.4-0.8")
if Pe is None:
Pe = abs(self.variables['Wt'])
if Pe <= 0:
raise ValueError("Pe must be positive")
fit_term = 30.522 * math.log(Pe) - 5.7178
if fit_term <= 0:
raise ValueError(
"Pe is outside the valid logarithmic domain for this "
"empirical mass formula"
)
mass = A * math.sqrt(0.5 * math.pi * fit_term)
return {
'mass': mass,
'A': A,
'Pe': Pe,
'fit_term': fit_term,
'formula': 'A*sqrt(0.5*pi*(30.522*ln(Pe)-5.7178))',
}
class Recuperator(ComponentMassMixin):
"""换热器类"""
component_type = "recuperator"
def __init__(self, name, eff, x=0):
"""
初始化参数
name: 名称
eff: 换热效率(基于焓的计算方法)
"""
self.name = name
self.eff = eff
self.Q_ex = None
self.variables = None
self.x = x
self._init_mass_interface()
def calculator(self, cold_inlet_state, hot_inlet_state, bypass_info, ploss=0.0, property_calculator=None):
"""
计算换热器两侧参数,默认逆流
bypass_info: 是否存在分流,0为不存在,1为存在
下标含义
----------
1: 换热器冷端入口
2: 换热器冷端出口
3: 换热器热端入口
4: 换热器热端出口
ass: 迭代中间变量,假设值
"""
# 计算两入口参数, 这里单位是kg
if property_calculator is None:
if hasattr(ploss, "calculate_properties"):
property_calculator = ploss
ploss = 0.0
else:
raise ValueError("property_calculator is required")
p1 = cold_inlet_state['P']
T1 = cold_inlet_state['T']
h1 = cold_inlet_state['h']
p2 = p1 * (1 - ploss)
p3 = hot_inlet_state['P']
T3 = hot_inlet_state['T']
h3 = hot_inlet_state['h']
p4 = p3 * (1 - ploss)
mw = property_calculator.mw
# 设定质量流量
m_cold = (1-self.x) if bypass_info == 1 else 1.0
m_hot = 1.0
# =============================================================================
# # 采用焓差效能的方式来计算换热器进出口参数
# # 假设最大温差发生在冷端, 计算冷端出口温度和比焓
# cold_outlet_state = property_calculator.calculate_properties(T=T3, P=p2)
# h2 = cold_outlet_state['h'] / mw
# Q_ass_cold = m_cold * abs(h2 - h1)
#
# # 假设最大温差发生在热端, 计算热端出口温度和比焓
# hot_outlet_state = property_calculator.calculate_properties(T=T1, P=p4)
# h4 = hot_outlet_state['h'] / mw
# Q_ass_hot = m_hot * abs(h3 - h4)
#
# # 比较两个可能的Q,取最小值与焓差效能的乘积作为实际换热量
# self.Q_ex = min(Q_ass_hot, Q_ass_cold) * self.eff
# # 由实际换热量计算出口焓和出口状态
# h2 = h1 + self.Q_ex / m_cold
# h4 = h3 - self.Q_ex / m_hot
# cold_outlet_state = property_calculator.calculate_properties(P=p2, h=h2*mw)
# hot_outlet_state = property_calculator.calculate_properties(P=p4, h=h4*mw)
# T2 = cold_outlet_state['T']
# T4 = hot_outlet_state['T']
# =============================================================================
# 采用温差效能的方式来计算换热器进出口参数
T_ass_max = abs(T1 - T3)
# 假设最大温差发生在冷端, 计算冷端出口温度和比焓
T2 = T1 + T_ass_max * self.eff
cold_outlet_state = property_calculator.calculate_properties(T=T2, P=p2)
h2 = cold_outlet_state['h'] / mw
Q_ass_cold = m_cold * abs(h2 - h1)
# 假设最大温差发生在热端, 计算热端出口温度和比焓
T4 = T3 - T_ass_max * self.eff
hot_outlet_state = property_calculator.calculate_properties(T=T4, P=p4)
h4 = hot_outlet_state['h'] / mw
Q_ass_hot = m_hot * abs(h3 - h4)
# 比较两个可能的Q,取最小值与焓差效能的乘积作为实际换热量
self.Q_ex = min(Q_ass_hot, Q_ass_cold)
# 由实际换热量计算出口焓和出口状态
h2 = h1 + self.Q_ex / m_cold
h4 = h3 - self.Q_ex / m_hot
cold_outlet_state = property_calculator.calculate_properties(P=p2, h=h2*mw)
hot_outlet_state = property_calculator.calculate_properties(P=p4, h=h4*mw)
T2 = cold_outlet_state['T']
T4 = hot_outlet_state['T']
# 拼装变量
self.variables = {
'name': self.name,
'cold_inlet_state':{
'P': p1,
'T': T1,
'h': h1,
's': cold_inlet_state['s']
},
'cold_outlet_state':{
'P': p2,
'T': cold_outlet_state['T'],
'h': h2,
's': cold_outlet_state['s']/mw
},
'hot_inlet_state':{
'P': p3,
'T': T3,
'h': h3,
's': hot_inlet_state['s']
},
'hot_outlet_state':{
'P': p4,
'T': T4,
'h': h4,
's': hot_outlet_state['s']/mw
},
'eff': self.eff,
'Q_exchange': self.Q_ex
}
def check_pinch_point(self, property_calculator, num_segments=20):
"""
换热器内部夹点校验
将换热量均分为 num_segments 段,检查内部每个微元的冷热流体温度
"""
h_cold_in = self.variables['cold_inlet_state']['h']
h_hot_in = self.variables['hot_inlet_state']['h']
p_cold = self.variables['cold_inlet_state']['P']
p_hot = self.variables['hot_inlet_state']['P']
m_cold = (1 - self.x) if self.name == "Low Temperature recuprerator" else 1.0
m_hot = 1.0
dQ = self.Q_ex / num_segments
# 沿冷流体流动方向步进检查
for i in range(num_segments + 1):
q_current = i * dQ
# 当前微元截面的焓值
h_cold_local = h_cold_in + q_current / m_cold
h_hot_local = (h_hot_in - self.Q_ex / m_hot) + q_current / m_hot
# 查温度
T_cold_local = property_calculator.calculate_properties(P=p_cold, h=h_cold_local * property_calculator.mw)['T']
T_hot_local = property_calculator.calculate_properties(P=p_hot, h=h_hot_local * property_calculator.mw)['T']
# 如果热流体温度低于等于冷流体温度 (设定一个 0.1K 的最小逼近温差容差)
if T_hot_local - T_cold_local < 0.1:
return False # 发生温度交叉,物理不可行!
return True
class Heater(ComponentMassMixin):
"""加热器类"""
component_type = "heater"
def __init__(self, name):
self.name = name
self.variables = None
self._init_mass_interface()
def calculator(self, inlet_state, outlet_state):
h_in = inlet_state['h']
h_out = outlet_state['h']
Q_input = h_out - h_in
self.variables = {
'name': self.name,
'inlet_state':{
'P': inlet_state['P'], # 压强(kPa)
'T': inlet_state['T'],
'h': inlet_state['h'],
's': inlet_state['s'],
},
'outlet_state':{
'P': outlet_state['P'], # 压强(kPa)
'T': outlet_state['T'],
'h': outlet_state['h'],
's': outlet_state['s'],
},
'Q_in': Q_input,
}
def _calculate_mass(
self,
P_heat_mwt=None,
mass_flow_rate=None,
shielding_mass_ton=2.8,
include_shielding=True,
**kwargs,
):
"""Estimate reactor and shielding mass for the heater module.
Reactor empirical formula from the provided reference:
M_reactor = 0.2195 * P_heat + 0.09836
P_heat is the reactor thermal power in MWt, and masses are in tons.
The shielding mass is added as a constant 2.8 ton by default.
If P_heat_mwt is not supplied, it is estimated from Q_in and
mass_flow_rate, assuming Q_in is kJ/kg and mass_flow_rate is kg/s:
P_heat_mwt = abs(Q_in) * mass_flow_rate / 1000
"""
if P_heat_mwt is None:
if mass_flow_rate is None:
raise ValueError(
"Heater mass calculation requires P_heat_mwt, or "
"mass_flow_rate to estimate P_heat_mwt from Q_in."
)
P_heat_mwt = abs(self.variables['Q_in']) * mass_flow_rate / 1000
if P_heat_mwt <= 0:
raise ValueError("P_heat_mwt must be positive")
if shielding_mass_ton < 0:
raise ValueError("shielding_mass_ton must be non-negative")
reactor_mass_ton = 0.2195 * P_heat_mwt + 0.09836
shielding_mass = shielding_mass_ton if include_shielding else 0.0
total_mass_ton = reactor_mass_ton + shielding_mass
return {
'mass': total_mass_ton,
'reactor_mass_ton': reactor_mass_ton,
'shielding_mass_ton': shielding_mass,
'P_heat_mwt': P_heat_mwt,
'include_shielding': include_shielding,
'formula': '0.2195*P_heat+0.09836+shielding_mass',
}
class Condenser(ComponentMassMixin):
"""冷凝器类"""
component_type = "condenser"
def __init__(self, name):
self.name = name
self.variables = None
self._init_mass_interface()
def calculator(self, inlet_state, outlet_state):
h_in = inlet_state['h']
h_out = outlet_state['h']
Q_output = h_in - h_out
self.variables = {
'name': self.name,
'inlet_state':{
'P': inlet_state['P'], # 压强(kPa)
'T': inlet_state['T'],
'h': inlet_state['h'],
's': inlet_state['s'],
},
'outlet_state':{
'P': outlet_state['P'], # 压强(kPa)
'T': outlet_state['T'],
'h': outlet_state['h'],
's': outlet_state['s'],
},
'Q_out': Q_output,
}
def _calculate_mass(
self,
Qc_kw=None,
co2_mass_flow_rate=None,
coolant_outlet_T=None,
water_inlet_T=None,
water_mass_flow_rate=None,
water_pressure_kpa=101.325,
refprop_path="C:/Program Files (x86)/REFPROP 10.0+/REFPROP",
emissivity=0.92,
surface_temperature=210.15,
area_density=6.75,
**kwargs,
):
"""Estimate radiator mass for the condenser module.
Radiator heat rejection model:
Qc = phi * sigma * A_rad * (T**4 - T0**4)
M_rad = kappa * A_rad
The NaK coolant in the reference is represented here by water. If
coolant_outlet_T is not supplied, water outlet temperature is evaluated
with REFPROP from water_inlet_T, water_mass_flow_rate, and Qc_kw.
Args:
Qc_kw: Total heat rejection in kW.
co2_mass_flow_rate: Optional kg/s. If Qc_kw is not supplied, Qc_kw
is estimated as Q_out * co2_mass_flow_rate, assuming Q_out is
kJ/kg.
coolant_outlet_T: Radiator coolant outlet temperature T in K.
water_inlet_T: Water inlet temperature in K for REFPROP calculation.
water_mass_flow_rate: Water mass flow rate in kg/s.
water_pressure_kpa: Water pressure for REFPROP calculation.
refprop_path: REFPROP root path.
emissivity: Radiator surface emissivity phi.
surface_temperature: Ambient/surface temperature T0 in K. Default is
-63 degC for Mars, 210.15 K.
area_density: Radiator face density kappa in kg/m2.
"""
if Qc_kw is None:
if co2_mass_flow_rate is None:
raise ValueError(
"Condenser mass calculation requires Qc_kw, or "
"co2_mass_flow_rate to estimate Qc_kw from Q_out."
)
Qc_kw = abs(self.variables['Q_out']) * co2_mass_flow_rate
if Qc_kw <= 0:
raise ValueError("Qc_kw must be positive")
if emissivity <= 0:
raise ValueError("emissivity must be positive")
if area_density <= 0:
raise ValueError("area_density must be positive")
if coolant_outlet_T is None:
coolant_outlet_T = self._calculate_water_outlet_temperature(
Qc_kw=Qc_kw,
water_inlet_T=water_inlet_T,
water_mass_flow_rate=water_mass_flow_rate,
water_pressure_kpa=water_pressure_kpa,
refprop_path=refprop_path,
)
temperature_term = coolant_outlet_T**4 - surface_temperature**4
if temperature_term <= 0:
raise ValueError(
"coolant_outlet_T must be higher than surface_temperature for "
"radiative heat rejection"
)
stefan_boltzmann = 5.670374419e-8
Qc_w = Qc_kw * 1000
area_m2 = Qc_w / (emissivity * stefan_boltzmann * temperature_term)
mass_kg = area_density * area_m2
return {
'mass': mass_kg,
'radiator_area_m2': area_m2,
'Qc_kw': Qc_kw,
'coolant_outlet_T': coolant_outlet_T,
'surface_temperature': surface_temperature,
'emissivity': emissivity,
'area_density': area_density,
'mass_unit': 'kg',
'formula': 'M_rad=kappa*Qc/(phi*sigma*(T^4-T0^4))',
}
def _calculate_water_outlet_temperature(
self,
Qc_kw,
water_inlet_T,
water_mass_flow_rate,
water_pressure_kpa,
refprop_path,
):
if water_inlet_T is None or water_mass_flow_rate is None:
raise ValueError(
"Provide coolant_outlet_T directly, or provide water_inlet_T "
"and water_mass_flow_rate for REFPROP water calculation."
)
if water_mass_flow_rate <= 0:
raise ValueError("water_mass_flow_rate must be positive")
water = ct.REFPROPFunctionLibrary(refprop_path)
water.SETUPdll(1, 'WATER.FLD', 'HMX.BNC', 'DEF')
water.SETUPdll(2, 'SI', 'SI', 'DEF')
z = [1.0]
mw = water.WMOLdll(z)
inlet = water.TPFLSHdll(water_inlet_T, water_pressure_kpa, z)
if inlet.ierr > 0:
raise ValueError(f"REFPROP water inlet calculation error: {inlet.ierr}")
h_in_mass = inlet.h / mw
h_out_mass = h_in_mass + Qc_kw / water_mass_flow_rate
outlet = water.PHFLSHdll(water_pressure_kpa, h_out_mass * mw, z)
if outlet.ierr > 0:
raise ValueError(f"REFPROP water outlet calculation error: {outlet.ierr}")
return outlet.T
class Concentrator(ComponentMassMixin):
"""汇流组件"""
component_type = "concentrator"
def __init__(self, name):
self.name = name
self.variables = None
self._init_mass_interface()
def calculator(self, inlet_state_bypass, inlet_state_mroad, x, property_calculator):
h_in_bypass = inlet_state_bypass['h']
h_in_mroad = inlet_state_mroad['h']
p_in = inlet_state_bypass['P']
mw = property_calculator.mw
h_out = x * h_in_bypass + (1-x) * h_in_mroad
outlet_state = property_calculator.calculate_properties(P=p_in, h=h_out*mw)
self.variables = {
'name': self.name,
'inlet_state_bypass': inlet_state_bypass,
'inlet_state_mroad': inlet_state_mroad,
'outlet_state':{
'P': p_in,
'T': outlet_state['T'],
'h': h_out,
's': outlet_state['s']/mw
}
}