Skip to content

Commit

Permalink
Merge pull request #999 from qualand/etes_sco2_dev
Browse files Browse the repository at this point in the history 
Improves the ETES model to handle user-defined HTF properties that include phase change specific heat curves.
  • Loading branch information
@qualand
qualand committed Mar 29, 2023
2 parents cf6a427 + 1049a41 commit a7606d1
Show file tree
Hide file tree
Showing 14 changed files with 215 additions and 544 deletions.
10 changes: 8 additions & 2 deletions ssc/cmod_sco2_csp_system.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,14 +62,20 @@ static var_info _cm_vtab_sco2_csp_system[] = {
"4) f_N_mc (=1 use design, <0, frac_des = abs(input),"
"5) f_N_pc (=1 use design, <0, frac_des = abs(input),"
"6) PHX_f_dP (=1 use design, <0 = abs(input)", "", "", "", "", "", "" },
{ SSC_INPUT, SSC_MATRIX, "od_max_htf_m_dot", "Columns: T_htf_C, T_amb_C, f_N_rc (=1 use design, <0, frac_des = abs(input), f_N_mc (=1 use design, <0, frac_des = abs(input), PHX_f_dP (=1 use design, <0 = abs(input), Rows: cases", "", "", "", "", "", "" },
{ SSC_INPUT, SSC_MATRIX, "od_max_htf_m_dot", "Columns: 0) T_htf_C, 1) T_amb_C,"
"2) f_N_rc (=1 use design, <0, frac_des = abs(input),"
"3) f_N_mc (=1 use design, <0, frac_des = abs(input),"
"4) PHX_f_dP (=1 use design, <0 = abs(input), Rows: cases", "", "", "", "", "", "" },
{ SSC_INPUT, SSC_MATRIX, "od_set_control", "Columns: 0) T_htf_C, 1) m_dot_htf_ND, 2) T_amb_C,"
" 3) P_LP_in_MPa, 4) T_mc_in_C, 5) T_pc_in_C,"
" 6) f_N_rc (=1 use design, <0, frac_des = abs(input),"
" 7) f_N_mc (=1 use design, <0, frac_des = abs(input),"
" 8) f_N_pc (=1 use design, =0 optimize, <0, frac_des = abs(input)),"
" 9) PHX_f_dP (=1 use design, <0 = abs(input), Rows: cases", "", "", "", "", "", "" },
{ SSC_INPUT, SSC_ARRAY, "od_generate_udpc", "True/False, f_N_rc (=1 use design, =0 optimize, <0, frac_des = abs(input), f_N_mc (=1 use design, =0 optimize, <0, frac_des = abs(input), PHX_f_dP (=1 use design, <0 = abs(input)", "", "", "", "", "", "" },
{ SSC_INPUT, SSC_ARRAY, "od_generate_udpc", "Columns 0) True/False,"
"1) f_N_rc (=1 use design, =0 optimize, <0, frac_des = abs(input),"
"2) f_N_mc (=1 use design, =0 optimize, <0, frac_des = abs(input),"
"3) PHX_f_dP (=1 use design, <0 = abs(input)", "", "", "", "", "", "" },
{ SSC_INPUT, SSC_NUMBER, "is_gen_od_polynomials","Generate off-design polynomials for Generic CSP models? 1 = Yes, 0 = No", "", "", "", "?=0", "", "" },

// ** Off-Design Outputs **
Expand Down
2 changes: 2 additions & 0 deletions ssc/csp_common.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -934,6 +934,7 @@ var_info vtab_sco2_design[] = {
{ SSC_OUTPUT, SSC_NUMBER, "PHX_co2_deltaP_des", "PHX co2 side design pressure drop", "-", "PHX Design Solution", "", "*", "", "" },
{ SSC_OUTPUT, SSC_NUMBER, "PHX_cost_equipment", "PHX cost equipment", "M$", "PHX Design Solution", "", "*", "", "" },
{ SSC_OUTPUT, SSC_NUMBER, "PHX_cost_bare_erected","PHX cost equipment and install", "M$", "PHX Design Solution", "", "*", "", "" },
{ SSC_OUTPUT, SSC_NUMBER, "PHX_min_dT", "PHX min temperature difference", "C", "PHX Design Solution", "", "*", "", "" },
// main compressor cooler
{ SSC_OUTPUT, SSC_NUMBER, "mc_cooler_T_in", "Low pressure cross flow cooler inlet temperature", "C", "Low Pressure Cooler", "", "*", "", "" },
{ SSC_OUTPUT, SSC_NUMBER, "mc_cooler_P_in", "Low pressure cross flow cooler inlet pressure", "MPa", "Low Pressure Cooler", "", "*", "", "" },
Expand Down Expand Up @@ -1686,6 +1687,7 @@ int sco2_design_cmod_common(compute_module *cm, C_sco2_phx_air_cooler & c_sco2_c
cm->assign("P_co2_PHX_in", (ssc_number_t)(c_sco2_cycle.get_design_solved()->ms_rc_cycle_solved.m_pres[C_sco2_cycle_core::HTR_HP_OUT] * 1.E-3)); //[MPa] convert from kPa
cm->assign("deltaT_HTF_PHX", (ssc_number_t)s_sco2_des_par.m_T_htf_hot_in - T_htf_cold_calc); //[K]
cm->assign("q_dot_PHX", (ssc_number_t)(c_sco2_cycle.get_design_solved()->ms_phx_des_solved.m_Q_dot_design*1.E-3)); //[MWt] convert from kWt
cm->assign("PHX_min_dT", (ssc_number_t)(c_sco2_cycle.get_design_solved()->ms_phx_des_solved.m_min_DT_design)); //[C/K]
double PHX_deltaP_frac = 1.0 - c_sco2_cycle.get_design_solved()->ms_rc_cycle_solved.m_pres[C_sco2_cycle_core::TURB_IN] /
c_sco2_cycle.get_design_solved()->ms_rc_cycle_solved.m_pres[C_sco2_cycle_core::HTR_HP_OUT]; //[-] Fractional pressure drop through co2 side of PHX
cm->assign("PHX_co2_deltaP_des", (ssc_number_t)PHX_deltaP_frac);
Expand Down
16 changes: 10 additions & 6 deletions tcs/csp_solver_cr_electric_resistance.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#include "csp_solver_core.h"

#include "htf_props.h"
#include "lib_physics.h"

static C_csp_reported_outputs::S_output_info S_cr_electric_resistance_output_info[] =
{
Expand Down Expand Up @@ -127,7 +128,7 @@ void C_csp_cr_electric_resistance::init(const C_csp_collector_receiver::S_csp_cr
m_dP_htf = 0.0;

// Calculate the design point HTF mass flow rate
m_cp_htf_des = mc_pc_htfProps.Cp_ave(m_T_htf_cold_des + 273.15, m_T_htf_hot_des + 273.15, 5); //[kJ/kg-K]
m_cp_htf_des = mc_pc_htfProps.Cp_ave(physics::CelciusToKelvin(m_T_htf_cold_des), physics::CelciusToKelvin(m_T_htf_hot_des)); //[kJ/kg-K]
m_m_dot_htf_des = m_q_dot_heater_des*1.E3 / (m_cp_htf_des*(m_T_htf_hot_des - m_T_htf_cold_des)); //[kg/s]
m_W_dot_heater_des = m_q_dot_heater_des / m_heater_efficiency; //[MWe]

Expand All @@ -140,10 +141,10 @@ void C_csp_cr_electric_resistance::init(const C_csp_collector_receiver::S_csp_cr
m_E_su_des = m_q_dot_su_max*m_hrs_startup_at_max_rate; //[MWt-hr]
m_t_su_des = m_E_su_des / m_q_dot_su_max; //[hr]

solved_params.m_T_htf_cold_des = m_T_htf_cold_des + 273.15; //[K]
solved_params.m_T_htf_cold_des = physics::CelciusToKelvin(m_T_htf_cold_des); //[K]
solved_params.m_P_cold_des = std::numeric_limits<double>::quiet_NaN(); //[kPa]
solved_params.m_x_cold_des = std::numeric_limits<double>::quiet_NaN(); //[-]
solved_params.m_T_htf_hot_des = m_T_htf_hot_des + 273.15; //[K]
solved_params.m_T_htf_hot_des = physics::CelciusToKelvin(m_T_htf_hot_des); //[K]
solved_params.m_q_dot_rec_des = m_q_dot_heater_des; //[MWt]
solved_params.m_A_aper_total = 0.0; //[m2]
solved_params.m_dP_sf = m_dP_htf; //[bar]
Expand Down Expand Up @@ -323,7 +324,10 @@ void C_csp_cr_electric_resistance::on(const C_csp_weatherreader::S_outputs& weat
m_E_su_calculated = 0.0; //[MWt-hr]
}

double m_dot_htf = q_dot_elec * 1.E3 / (m_cp_htf_des*(m_T_htf_hot_des - htf_state_in.m_temp)); //[kg/s]
double W_dot_heater = q_dot_elec; //[MWe]

double cp_ave = mc_pc_htfProps.Cp_ave(physics::CelciusToKelvin(htf_state_in.m_temp), physics::CelciusToKelvin(m_T_htf_hot_des));
double m_dot_htf = q_dot_elec * 1.E3 / (cp_ave*(m_T_htf_hot_des - htf_state_in.m_temp)); //[kg/s]

double q_startup = 0.0; //[MWt-hr]
double q_dot_startup = 0.0; //[MWt-hr]
Expand Down Expand Up @@ -370,8 +374,8 @@ void C_csp_cr_electric_resistance::estimates(const C_csp_weatherreader::S_output
// 2) heater is always capable of design output
// 3) no mass flow rate bounds (for now)
// 4) heater is controlled to always return HTF at design hot temperature

double m_dot_htf = m_q_dot_heater_des * 1.E3 / (m_cp_htf_des * (m_T_htf_hot_des - htf_state_in.m_temp)); //[kg/s]
double cp_ave = mc_pc_htfProps.Cp_ave(physics::CelciusToKelvin(htf_state_in.m_temp), physics::CelciusToKelvin(m_T_htf_hot_des));
double m_dot_htf = m_q_dot_heater_des * 1.E3 / (cp_ave * (m_T_htf_hot_des - htf_state_in.m_temp)); //[kg/s]

E_csp_cr_modes mode = get_operating_state();

Expand Down
54 changes: 37 additions & 17 deletions tcs/csp_solver_pc_Rankine_indirect_224.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -154,7 +154,7 @@ void C_pc_Rankine_indirect_224::init(C_csp_power_cycle::S_solved_params &solved_
// ***********************************************************************
// Design-point calculations before initializing Rankine or UDPC model
// Calculate design point HTF mass flow rate
m_cp_htf_design = mc_pc_htfProps.Cp(physics::CelciusToKelvin((ms_params.m_T_htf_hot_ref + ms_params.m_T_htf_cold_ref) / 2.0)); //[kJ/kg-K]
m_cp_htf_design = mc_pc_htfProps.Cp_ave(physics::CelciusToKelvin(ms_params.m_T_htf_hot_ref), physics::CelciusToKelvin(ms_params.m_T_htf_cold_ref)); //[kJ/kg-K]

ms_params.m_P_ref *= 1000.0; //[kW] convert from MW
m_q_dot_design = ms_params.m_P_ref / 1000.0 / ms_params.m_eta_ref; //[MWt]
Expand Down Expand Up @@ -1004,7 +1004,7 @@ double C_pc_Rankine_indirect_224::get_efficiency_at_load(double load_frac, doubl

if( !ms_params.m_is_user_defined_pc )
{
double cp = mc_pc_htfProps.Cp( (ms_params.m_T_htf_cold_ref + ms_params.m_T_htf_hot_ref)/2. +273.15); //kJ/kg-K
double cp = mc_pc_htfProps.Cp_ave(physics::CelciusToKelvin(ms_params.m_T_htf_cold_ref), physics::CelciusToKelvin(ms_params.m_T_htf_hot_ref)); //kJ/kg-K

//calculate mass flow [kg/hr]
double mdot = ms_params.m_P_ref /* kW */ / ( ms_params.m_eta_ref * cp * (ms_params.m_T_htf_hot_ref - ms_params.m_T_htf_cold_ref) ) *3600.;
Expand Down Expand Up @@ -1166,8 +1166,7 @@ void C_pc_Rankine_indirect_224::call(const C_csp_weatherreader::S_outputs &weath
{
case STARTUP:
{
double c_htf = mc_pc_htfProps.Cp(physics::CelciusToKelvin((T_htf_hot + ms_params.m_T_htf_cold_ref) / 2.0)); //[kJ/kg-K]

double c_htf = mc_pc_htfProps.Cp_ave(physics::CelciusToKelvin(ms_params.m_T_htf_cold_ref), physics::CelciusToKelvin(T_htf_hot)); //[kJ/kg-K]
double time_required_su_energy = m_startup_energy_remain_prev / (m_dot_htf*c_htf*(T_htf_hot - ms_params.m_T_htf_cold_ref)/3600.0); //[hr]
double time_required_su_ramping = m_startup_time_remain_prev; //[hr]

Expand Down Expand Up @@ -1399,8 +1398,7 @@ void C_pc_Rankine_indirect_224::call(const C_csp_weatherreader::S_outputs &weath
// Calculate other important metrics
eta = P_cycle / 1.E3 / q_dot_htf; //[-]

// Want to iterate to fine more accurate cp_htf?
T_htf_cold = T_htf_hot - q_dot_htf / (m_dot_htf / 3600.0*m_cp_htf_design / 1.E3); //[MJ/s * hr/kg * s/hr * kg-K/kJ * MJ/kJ] = C/K
T_htf_cold = Calculate_T_htf_cold_Converge_Cp(q_dot_htf * 1.E3, physics::CelciusToKelvin(T_htf_hot), m_dot_htf / 3600.0) - 273.15; //[K] -> [C]

was_method_successful = true;
}
Expand Down Expand Up @@ -1428,8 +1426,7 @@ void C_pc_Rankine_indirect_224::call(const C_csp_weatherreader::S_outputs &weath

case STANDBY:
{
double c_htf = mc_pc_htfProps.Cp(physics::CelciusToKelvin((T_htf_hot + ms_params.m_T_htf_cold_ref) / 2.0)); //[kJ/kg-K]
// double c_htf = specheat(m_pbp.HTF, physics::CelciusToKelvin((m_pbi.T_htf_hot + m_pbp.T_htf_cold_ref)/2.0), 1.0);
double c_htf = mc_pc_htfProps.Cp_ave(physics::CelciusToKelvin(ms_params.m_T_htf_cold_ref), physics::CelciusToKelvin(T_htf_hot)); //[kJ/kg-K]
double q_tot = ms_params.m_P_ref / ms_params.m_eta_ref;

// Calculate the actual q_sby_needed from the reference flows
Expand Down Expand Up @@ -1574,7 +1571,7 @@ void C_pc_Rankine_indirect_224::call(const C_csp_weatherreader::S_outputs &weath
// simultaneously with the required startup time. If the timestep is less than the required startup time
// scale the mass flow rate appropriately

double c_htf = mc_pc_htfProps.Cp(physics::CelciusToKelvin((T_htf_hot + ms_params.m_T_htf_cold_ref) / 2.0)); //[kJ/kg-K]
double c_htf = mc_pc_htfProps.Cp_ave(physics::CelciusToKelvin(ms_params.m_T_htf_cold_ref), physics::CelciusToKelvin(T_htf_hot)); //[kJ/kg-K]

// Maximum thermal power to power cycle based on heat input constraint parameters:
double q_dot_to_pc_max_q_constraint = ms_params.m_cycle_max_frac * ms_params.m_P_ref / ms_params.m_eta_ref; //[kWt]
Expand Down Expand Up @@ -1731,7 +1728,7 @@ void C_pc_Rankine_indirect_224::call(const C_csp_weatherreader::S_outputs &weath
m_startup_time_remain_calc = fmax(m_startup_time_remain_prev - step_hrs, 0.0);
m_startup_energy_remain_calc = fmax(m_startup_energy_remain_prev - startup_e_used, 0.0);

double c_htf = mc_pc_htfProps.Cp(physics::CelciusToKelvin((T_htf_hot + ms_params.m_T_htf_cold_ref) / 2.0)); //[kJ/kg-K]
double c_htf = mc_pc_htfProps.Cp_ave(physics::CelciusToKelvin(ms_params.m_T_htf_cold_ref), physics::CelciusToKelvin(T_htf_hot)); //[kJ/kg-K]
// If still starting up, then all of energy input going to startup
if(m_startup_time_remain_calc + m_startup_energy_remain_calc > 0.0)
{
Expand Down Expand Up @@ -2034,11 +2031,6 @@ void C_pc_Rankine_indirect_224::RankineCycle_V2(double T_db /*K*/, double T_wb /
double P_cond_min = ms_params.m_P_cond_min;

water_state wp;

// Calculate the specific heat before converting to Kelvin
double c_htf_ref = mc_pc_htfProps.Cp(physics::CelciusToKelvin((T_htf_hot_ref + T_htf_cold_ref) / 2.0));
double c_htf = mc_pc_htfProps.Cp(physics::CelciusToKelvin((T_htf_hot + T_htf_cold_ref) / 2.0)); //[kJ/kg-k]

// Convert units
// **Temperatures from Celcius to Kelvin
T_htf_hot = physics::CelciusToKelvin(T_htf_hot); //[K]
Expand Down Expand Up @@ -2139,8 +2131,7 @@ void C_pc_Rankine_indirect_224::RankineCycle_V2(double T_db /*K*/, double T_wb /

// Final performance calcs
double q_dot_cycle = P_cycle / eta;

T_htf_cold = T_htf_hot - q_dot_cycle / (m_dot_htf * c_htf); //[K]
T_htf_cold = Calculate_T_htf_cold_Converge_Cp(q_dot_cycle, T_htf_hot, m_dot_htf); //[K]
m_dot_demand = fmax(m_dot_htf_ND * m_dot_htf_ref, 0.00001); //[kg/s]

// Finally, convert to output units/names
Expand Down Expand Up @@ -2286,3 +2277,32 @@ double C_pc_Rankine_indirect_224::Interpolate(int YT, int XT, double X, double Z

} // Interpolate

double C_pc_Rankine_indirect_224::Calculate_T_htf_cold_Converge_Cp(double q_dot_htf /*kWt*/, double T_htf_hot /*K*/, double m_dot_htf /*kg/s*/)
{
double alpha = 0.3;
int iter = 0;
double T_htf_cold = physics::CelciusToKelvin(ms_params.m_T_htf_cold_ref);
double c_htf, T_error = 1.0, T_htf_cold_prev;
while (fabs(T_error) > 1e-4 && iter < 30) {
// TODO: set up C_monotonic_equation? and C_monotonic_eq_solver?
T_htf_cold_prev = T_htf_cold;
try {
c_htf = mc_pc_htfProps.Cp_ave(T_htf_cold, T_htf_hot); //[kJ/kg-K]
}
catch (C_csp_exception) {
// Recover T_htf_cold or T_htf_hot is < 0
T_error = 1;
break;
}
T_htf_cold = T_htf_hot - q_dot_htf / (m_dot_htf * c_htf); //[kJ/s * s/kg * kg-K/kJ] = C/K
T_htf_cold = T_htf_cold * alpha + T_htf_cold_prev * (1 - alpha);
T_error = (T_htf_cold - T_htf_cold_prev) / T_htf_cold_prev;
iter++;
}
if (fabs(T_error) > 1e-4) {
// Divergent - Use the initial iteration and deal with error
c_htf = mc_pc_htfProps.Cp_ave(physics::CelciusToKelvin(ms_params.m_T_htf_cold_ref), T_htf_hot); //[kJ/kg-K]
T_htf_cold = T_htf_hot - q_dot_htf / (m_dot_htf * c_htf); //[kJ/s * s/kg * kg-K/kJ] = C/K
}
return T_htf_cold; //[K]
}
2 changes: 2 additions & 0 deletions tcs/csp_solver_pc_Rankine_indirect_224.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,6 +168,8 @@ class C_pc_Rankine_indirect_224 : public C_csp_power_cycle

double Interpolate(int YT, int XT, double X, double Z = std::numeric_limits<double>::quiet_NaN());

double Calculate_T_htf_cold_Converge_Cp(double q_dot_htf /*kWt*/, double T_htf_hot /*K*/, double m_dot_htf /*kg/s*/);

// Isopentane
double T_sat4(double P/*Bar*/)
{
Expand Down
4 changes: 2 additions & 2 deletions tcs/csp_solver_pc_heat_sink.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,7 +112,7 @@ void C_pc_heat_sink::init(C_csp_power_cycle::S_solved_params &solved_params)
}

// Calculate the design point HTF mass flow rate
double cp_htf_des = mc_pc_htfProps.Cp_ave(ms_params.m_T_htf_cold_des+273.15, ms_params.m_T_htf_hot_des+273.15, 5); //[kJ/kg-K]
double cp_htf_des = mc_pc_htfProps.Cp_ave(ms_params.m_T_htf_cold_des+273.15, ms_params.m_T_htf_hot_des+273.15); //[kJ/kg-K]

m_m_dot_htf_des = ms_params.m_q_dot_des*1.E3 / (cp_htf_des*(ms_params.m_T_htf_hot_des - ms_params.m_T_htf_cold_des)); //[kg/s]

Expand Down Expand Up @@ -223,7 +223,7 @@ void C_pc_heat_sink::call(const C_csp_weatherreader::S_outputs &weather,
double T_htf_hot = htf_state_in.m_temp; //[C]
double m_dot_htf = inputs.m_m_dot/3600.0; //[kg/s]

double cp_htf = mc_pc_htfProps.Cp_ave(ms_params.m_T_htf_cold_des+273.15, T_htf_hot+273.15, 5); //[kJ/kg-K]
double cp_htf = mc_pc_htfProps.Cp_ave(ms_params.m_T_htf_cold_des+273.15, T_htf_hot+273.15); //[kJ/kg-K]

// For now, let's assume the Heat Sink can always return the HTF at the design cold temperature
double q_dot_htf = m_dot_htf*cp_htf*(T_htf_hot - ms_params.m_T_htf_cold_des)/1.E3; //[MWt]
Expand Down
2 changes: 1 addition & 1 deletion tcs/csp_solver_trough_collector_receiver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1519,7 +1519,7 @@ int C_csp_trough_collector_receiver::loop_energy_balance_T_t_int(const C_csp_wea
m_E_dot_HR_hot_subts = E_HR_hot / sim_info.ms_ts.m_step; //[MWt]

// HTF out of system
m_c_htf_ave_ts_ave_temp = m_htfProps.Cp_ave(T_htf_cold_in, m_T_sys_h_t_int, 5)*1000.0; //[J/kg-K]
m_c_htf_ave_ts_ave_temp = m_htfProps.Cp_ave(T_htf_cold_in, m_T_sys_h_t_int)*1000.0; //[J/kg-K]
m_q_dot_htf_to_sink_subts = m_m_dot_htf_tot*m_c_htf_ave_ts_ave_temp*(m_T_sys_h_t_int - T_htf_cold_in)*1.E-6;

double Q_dot_balance_subts = m_q_dot_sca_abs_summed_subts - m_q_dot_xover_loss_summed_subts -
Expand Down
Loading

0 comments on commit a7606d1

Please sign in to comment.