Skip to content

Commit

Permalink
Adding flux profile interpolation, receiver curtain discretization in…
Browse files Browse the repository at this point in the history 
…puts, and receiver orientation
  • Loading branch information
@jannamartinek
jannamartinek committed May 10, 2024
1 parent 57d4fda commit e7a5a4f
Show file tree
Hide file tree
Showing 3 changed files with 137 additions and 36 deletions.
13 changes: 9 additions & 4 deletions ssc/cmod_csp_tower_particle.cpp
View file
Expand Up @@ -186,6 +186,8 @@ static var_info _cm_vtab_csp_tower_particle[] = {
{ SSC_INPUT, SSC_NUMBER, "rec_adv_model_type", "Receiver advection model type (0 = Fixed advective loss coefficient, 1 = Correlations from Sandia)", "", "", "Tower and Receiver", "?=1", "", ""},
{ SSC_INPUT, SSC_NUMBER, "rec_eta_fixed", "User-provided fixed receiver efficiency (fraction)", "", "", "Tower and Receiver", "rec_model_type=0", "", ""},
{ SSC_INPUT, SSC_NUMBER, "rec_hadv_fixed", "User-provided fixed receiver advective loss coefficient", "W/m2/K", "", "Tower and Receiver", "rec_adv_model_type=0", "", ""},
{ SSC_INPUT, SSC_NUMBER, "rec_ny", "Number of curtain elements in the vertical dimension", "", "", "Tower and Receiver", "?=12", "", "" },
{ SSC_INPUT, SSC_NUMBER, "rec_nx", "Number of curtain elements in the width dimension", "", "", "Tower and Receiver", "?=1", "", "" },
{ SSC_INPUT, SSC_NUMBER, "rec_rad_ny", "Number of curtain element groups in vertical dimension for radiative exchange", "", "", "Tower and Receiver", "?=5", "", ""},
{ SSC_INPUT, SSC_NUMBER, "rec_rad_nx", "Number of curtain element groups in width dimension for radiative exchange", "", "", "Tower and Receiver", "?=3", "", ""},
{ SSC_INPUT, SSC_NUMBER, "particle_dp", "Particle diameter", "m", "", "Tower and Receiver", "?=350e-6", "", ""},
Expand Down Expand Up @@ -1154,7 +1156,7 @@ class cm_csp_tower_particle : public compute_module
double user_hadv = as_integer("rec_adv_model_type") == 0 ? as_double("rec_hadv_fixed") : 0.0;

int input_idx;
double ap_height, ap_width, curtain_height, curtain_width, ap_curtain_depth_ratio;
double ap_height, ap_width, curtain_height, curtain_width, ap_curtain_depth_ratio, rec_orientation;
for (int i = 0; i < num_recs; i++) { // loop through receivers
input_idx = duplicate_recs ? 0 : i; // Use the first element if duplicate receivers
ap_height = rec_height[input_idx];
Expand All @@ -1164,11 +1166,14 @@ class cm_csp_tower_particle : public compute_module
A_rec_curtain = curtain_height * curtain_width; // This receiver area is used to define the flux distribution. Particle receiver model assumes that the flux distribution is defined based on the curtain area.
A_rec_aperture = ap_height * ap_width; // The aperture area should be used in cost calculations
ap_curtain_depth_ratio = max_curtain_depth[input_idx] / ap_height;
rec_orientation = rec_azimuth[input_idx];
if (rec_orientation < 0)
rec_orientation += 360; // Relative wind directions in receiver code require receiver orientation witin [0,360]

// TODO (Janna): update to initialize individual receivers
}
int n_x = mt_flux_maps.ncols() / num_recs;
int n_y = (mt_flux_maps.nrows() / mt_eta_map.nrows()) + 1;
int n_x = as_integer("rec_nx"); //mt_flux_maps.ncols() / num_recs;
int n_y = as_integer("rec_ny"); //(mt_flux_maps.nrows() / mt_eta_map.nrows()) + 1;


if (as_integer("curtain_type") == 1) {
Expand All @@ -1183,7 +1188,7 @@ class cm_csp_tower_particle : public compute_module
as_double("csp.pt.rec.max_oper_frac"), as_double("eta_lift"),
as_integer("rec_htf"), as_matrix("field_fl_props"),
as_integer("rec_model_type"), user_efficiency, as_integer("rec_rad_model_type"), as_integer("rec_adv_model_type"), user_hadv,
ap_height, ap_width, norm_curtain_height[input_idx], norm_curtain_width[input_idx], ap_curtain_depth_ratio,
ap_height, ap_width, norm_curtain_height[input_idx], norm_curtain_width[input_idx], ap_curtain_depth_ratio, rec_orientation,
as_double("particle_dp"), as_double("particle_abs"), as_double("curtain_emis"), as_double("curtain_dthdy"),
as_double("cav_abs"), as_double("cav_twall"), as_double("cav_kwall"), as_double("cav_hext"),
as_double("transport_deltaT_cold"), as_double("transport_deltaT_hot"),
Expand Down
151 changes: 120 additions & 31 deletions tcs/csp_solver_falling_particle_receiver.cpp
View file
Expand Up @@ -47,7 +47,7 @@ C_falling_particle_receiver::C_falling_particle_receiver(double h_tower /*m*/,
double m_dot_htf_max_frac /*-*/, double eta_pump /*-*/,
int field_fl, util::matrix_t<double> field_fl_props,
int model_type /*-*/, double fixed_efficiency /*-*/, int rad_model_type /*-*/, int hadv_model_type /*-*/, double hadv_user /*-*/,
double ap_height /*m*/, double ap_width /*m*/, double ap_height_ratio /*-*/, double ap_width_ratio /*-*/, double ap_curtain_depth_ratio /*-*/,
double ap_height /*m*/, double ap_width /*m*/, double ap_height_ratio /*-*/, double ap_width_ratio /*-*/, double ap_curtain_depth_ratio /*-*/, double rec_orientation /*deg*/,
double particle_dp /*m*/, double particle_abs /*-*/, double curtain_emis /*-*/, double dthdy /*-*/,
double cav_emis /*-*/, double cav_twall /*m*/, double cav_kwall /*m*/, double cav_hext /*W/m2/K*/,
double deltaT_transport_cold /*K*/, double deltaT_transport_hot /*K*/,
Expand Down Expand Up @@ -80,6 +80,7 @@ C_falling_particle_receiver::C_falling_particle_receiver(double h_tower /*m*/,
m_ap_height_ratio = ap_height_ratio;
m_ap_width_ratio = ap_width_ratio;
m_ap_curtain_depth_ratio = ap_curtain_depth_ratio;
m_rec_orientation = rec_orientation;

m_is_ap_at_bot = true; // Hard-coded to true to match SolarPILOT
m_is_curtain_flat = true; // Curved curtain is not implemented yet
Expand Down Expand Up @@ -124,6 +125,7 @@ C_falling_particle_receiver::C_falling_particle_receiver(double h_tower /*m*/,
m_eta_therm_des_est = 0.9; //[-] Estimated and used to calculate min incident flux
m_include_back_wall_convection = false;
m_include_wall_axial_conduction = false;


m_invert_matrices = true; // Invert coefficient matrix for radiative exchange? Used multiple times during temperature iterations - test case solution time was faster when the coefficient matrix was inverted up front

Expand Down Expand Up @@ -312,16 +314,9 @@ void C_falling_particle_receiver::call(const C_csp_weatherreader::S_outputs& wea
throw(C_csp_exception("Wind direction is not defined", "Falling particle receiver timestep performance call"));
}

// Check resolution of flux map input compared to resolution of the particle curtain discretization
// TODO: Generalize to allow different resolution
// Calculate total flux input prior to interpolation
int n_flux_y = (int)flux_map_input->nrows();
int n_flux_x = (int)flux_map_input->ncols();
if (n_flux_y + 1 != m_n_y || n_flux_x != m_n_x) {
error_msg = "The falling particle receiver model flux map input must match the specified discretization resolution of the particle curtain";
throw(C_csp_exception(error_msg, "Falling particle receiver timestep performance call"));
}


double flux_sum = 0.0;
for (int i = 0; i < n_flux_y; i++) {
for (int j = 0; j < n_flux_x; j++) {
Expand Down Expand Up @@ -841,39 +836,125 @@ int C_falling_particle_receiver::use_previous_solution(const s_steady_state_soln

}

util::matrix_t<double> C_falling_particle_receiver::interpolate_flux_1d(util::matrix_t<double>& flux_input, int n, bool is_rows)
{
int n1, n2, nflux, nc, i0, i1, nelem;
double dx, dxflux, d1, d2, x;
util::matrix_t<double> flux_input_loc, flux_interp, flux;

flux_input_loc = flux_input;
if (!is_rows) // Transpose input array
{
n1 = flux_input.nrows();
n2 = flux_input.ncols();
flux_input_loc.resize_fill(n2, n1, 0.0);
for (int i = 0; i < n2; i++)
{
for (int j = 0; j < n1; j++)
flux_input_loc.at(i, j) = flux_input.at(j, i);
}
}

nflux = flux_input_loc.nrows();
nc = flux_input_loc.ncols();
flux_interp.resize_fill(n, nc, 0.0);
dx = 1.0 / n;
dxflux = 1.0 / nflux;
for (int i = 0; i < n; i++)
{
if (n <= nflux) // Curtain resolution < flux profile resolution
{
i0 = fmin(floor(i * dx / dxflux), nflux - 1); // Flux profile element containing lower bound of element j
i1 = fmin(floor((i + 1) * dx / dxflux), nflux - 1); // Flux profile element containing upper bound of element j
d1 = ((i0 + 1) * dxflux - (i * dx)); // Distance from upper bound of flux element j0 to lower bound of element j
d2 = ((i + 1) * dx - i1 * dxflux); // Distance from upper bound of element j to lower bound of element j1
nelem = i1 - i0;
for (int j = 0; j < nc; j++)
{
if (nelem < 1)
flux_interp.at(i, j) = flux_input_loc.at(i0, j);
else
{
flux_interp.at(i, j) = d1 * flux_input_loc.at(i0, j) + d2 * flux_input_loc.at(i1, j);
for (int k = 0; k < nelem - 1; k++)
flux_interp.at(i, j) += dxflux * flux_input_loc.at(i0+1+k, j);
flux_interp.at(i, j) /= dx;
flux_interp.at(i, j) = fmax(flux_interp.at(i, j), 0.0);
}
}
}
else // Curtain resolution > flux profile resolution
{
x = (i + 0.5) * dx;
i0 = fmax(fmin(floor(x / dxflux - 0.5), nflux - 2), 0);
for (int j = 0; j < nc; j++)
{
if (x < 0.5 * dxflux)
flux_interp.at(i, j) = flux_input_loc.at(0, j);
else if (x > 1 - 0.5 * dxflux)
flux_interp.at(i, j) = flux_input_loc.at(nflux - 1, j);
else
flux_interp.at(i, j) = flux_input_loc.at(i0, j) + (flux_input_loc.at(i0+1, j) - flux_input_loc.at(i0, j)) / dxflux * (x - (i0 + 0.5) * dxflux);
flux_interp.at(i, j) = fmax(flux_interp.at(i, j), 0.0);
}
}
}

// Calculate flux profiles (interpolated to curtain elements at specified DNI)
// TODO: Update with interpolation to allow the curtain discretization resolution to differ from the provided flux profile resolution
// Transpose interpolated array
flux = flux_interp;
if (!is_rows) // Tranpose interpolated array
{
n1 = flux_interp.nrows();
n2 = flux_interp.ncols();
flux.resize_fill(n2, n1, 0.0);
for (int i = 0; i < n2; i++)
{
for (int j = 0; j < n1; j++)
flux.at(i, j) = flux_interp.at(j, i);
}
}
return flux;
}



// Calculate flux profiles (interpolated to curtain elements)
util::matrix_t<double> C_falling_particle_receiver::calculate_flux_profiles(double flux_sum /*W/m2*/, double flux_scale /*-*/, double plant_defocus /*-*/,
double od_control, const util::matrix_t<double>* flux_map_input)
{
util::matrix_t<double> q_dot_inc;
int n_flux_y = (int)flux_map_input->nrows();
int n_flux_x = (int)flux_map_input->ncols();
q_dot_inc.resize_fill(m_n_y, m_n_x, 0.0);

double total_defocus = plant_defocus * od_control;
q_dot_inc.resize_fill(m_n_y, m_n_x, 0.0);
double total_defocus = plant_defocus * od_control; // Combine controller-specified defocus (plant_defocus) and receiver defocus (od_control)

if (flux_sum > 1.0)
{
for (int i = 0; i < n_flux_x; i++)
//--- Scale input flux
int n_flux_y = (int)flux_map_input->nrows();
int n_flux_x = (int)flux_map_input->ncols();
util::matrix_t<double>flux_scaled(n_flux_y, n_flux_x);
for (int j = 0; j < n_flux_y; j++)
{
for (int j = 0; j < n_flux_y; j++)
for (int i = 0; i < n_flux_x; i++)
{
q_dot_inc.at(j, i) = (*flux_map_input)(j, i) * total_defocus * flux_scale * 1000; // Flux incident on curtain (assuming curtain discretization resolution and flux profile resolution match) [W/m^2];

flux_scaled.at(j,i) = (*flux_map_input)(j, i) * total_defocus * flux_scale * 1000; // Incident flux adjusted for DNI and defocus [W/m2]
}
}

// Fill in interpolated value at last axial node in fall direction. This doesn't count toward energy balance, value here is just to define reasonable wall temperature
int j = m_n_y - 1;
q_dot_inc.at(j, i) = q_dot_inc.at(j-1, i) + (q_dot_inc.at(j-1,i) - q_dot_inc.at(j-2, i));

//--- Interpolate
util::matrix_t<double> flux_x = (m_n_y-1 == n_flux_y) ? flux_scaled : interpolate_flux_1d(flux_scaled, m_n_y - 1, true); // Interpolate to curtain y-resolution
util::matrix_t<double> flux = (m_n_x == n_flux_x) ? flux_x : interpolate_flux_1d(flux_x, m_n_x, false); // Interpolate to curtain x-resolution
for (int j = 0; j < m_n_y; j++)
{
for (int i = 0; i < m_n_x; i++)
{
if (j < m_n_y - 1)
q_dot_inc.at(j, i) = flux.at(j, i);
else
q_dot_inc.at(j, i) = flux.at(j - 1, i); // Note the flux at the last height node doesn't contribute to the energy balance, this value will only be used to calculation the wall temperature solution at the last node (which also doesn't contribute to the energy balance)
}
}
}
else
{
q_dot_inc.fill(0.0);
}

return q_dot_inc;
}

Expand Down Expand Up @@ -1670,8 +1751,13 @@ double C_falling_particle_receiver::sandia_efficiency_correlation(bool is_multis
H = 5000.;
}

// Adjust wind direction relative to receiver orientation
double wdir_eff = wind_direc - m_rec_orientation;
if (wdir_eff < 0)
wdir_eff = 360 + wdir_eff;

q = exp(-Q_inc * 1.e-6 / m_ap_area);
val = 180 - fabs(180 - wind_direc);
val = 180 - fabs(180 - wdir_eff);
theta = pow(val, F) * exp(-val / G) / H;
eta = A + B * q + C * pow(q, 2) + D * q * v_wind * theta + E * pow(v_wind, 2) * theta;
eta = fmax(eta, 0.0);
Expand Down Expand Up @@ -1765,12 +1851,15 @@ void C_falling_particle_receiver::calculate_advection_coeff_sandia(double vel, d
hadv = Nu * k_air / m_curtain_height; // Advective loss coefficient without wind (W/m2/K)

double wind_factor = 1.0;
double phi_wind = exp(-pow((std::abs(wdir - F) - G) / H, 2));
double wdir_eff = wdir - m_rec_orientation; // Adjust wind direction relative to receiver orientation
if (wdir_eff < 0)
wdir_eff = 360 + wdir_eff;
double phi_wind = exp(-pow((std::abs(wdir_eff - F) - G) / H, 2));
fwind = 1.0 + (D - E * m_ap_height) * wspd * phi_wind;

// Scale loss coefficient to account for differences in aperture and curtain area
// Loss coefficient from Gonzalez et al. assumes equal aperture area and curtain area and will overestimate convective losses if the loss coefficient is applied for cases where curtain area > aperture area
hadv *= m_ap_height / m_curtain_height;
hadv *= m_ap_area / m_curtain_area;

return;
}
Expand Down
9 changes: 8 additions & 1 deletion tcs/csp_solver_falling_particle_receiver.h
View file
Expand Up @@ -155,6 +155,9 @@ class C_falling_particle_receiver : public C_pt_receiver
double m_initial_fall_height_fixed; // Fixed initial fall height [m]
double m_initial_fall_height_fraction; // Initial fall height of particle curtain before reaching top of cavity / cavity height [-]

// Receiver orientation
double m_rec_orientation; // Receiver orientation (0 = north, 90 = east, 180 = south, 270 = west


// Particle and curtain properties
double m_particle_dp; // Particle diameter [m]
Expand Down Expand Up @@ -262,8 +265,12 @@ class C_falling_particle_receiver : public C_pt_receiver
double scale_wind_speed(double v_wind_10);

int use_previous_solution(const s_steady_state_soln& soln, const std::vector<s_steady_state_soln>& stored_solns);

util::matrix_t<double> interpolate_flux_1d(util::matrix_t<double>& flux_input, int n, bool is_rows);

util::matrix_t<double> calculate_flux_profiles(double flux_sum /*W/m2*/, double dni_scale /*-*/, double plant_defocus /*-*/,
double od_control /*-*/, const util::matrix_t<double>* flux_map_input);

void calculate_steady_state_soln(s_steady_state_soln& soln, double tol = 1.0e-4, bool init_from_existing = false, int max_iter = 50);
void solve_for_mass_flow(s_steady_state_soln& soln);
void solve_for_mass_flow_and_defocus(s_steady_state_soln& soln, double m_dot_htf_max, const util::matrix_t<double>* flux_map_input);
Expand Down Expand Up @@ -322,7 +329,7 @@ class C_falling_particle_receiver : public C_pt_receiver
double m_dot_htf_max_frac /*-*/, double eta_pump /*-*/,
int field_fl, util::matrix_t<double> field_fl_props,
int model_type /*-*/, double fixed_efficiency /*-*/, int rad_model_type /*-*/, int hadv_model_type /*-*/, double hadv_user /*-*/,
double ap_height /*m*/, double ap_width /*m*/, double ap_height_ratio /*-*/, double ap_width_ratio /*-*/, double ap_curtain_depth_ratio /*-*/,
double ap_height /*m*/, double ap_width /*m*/, double ap_height_ratio /*-*/, double ap_width_ratio /*-*/, double ap_curtain_depth_ratio /*-*/, double rec_orientation,
double particle_dp /*m*/, double particle_abs /*-*/, double curtain_emis /*-*/, double dthdy /*-*/,
double cav_emis /*-*/, double cav_twall /*m*/, double cav_kwall /*m*/, double cav_hext /*W/m2/K*/,
double deltaT_transport_cold /*K*/, double deltaT_transport_hot /*K*/,
Expand Down

0 comments on commit e7a5a4f

Please sign in to comment.