Skip to content

Commit

Permalink
Updated clearsky calculations, added debugging outputs for now
Browse files Browse the repository at this point in the history
  • Loading branch information
@mjprilliman
mjprilliman committed Oct 18, 2023
1 parent c4971d6 commit d774469
Show file tree
Hide file tree
Showing 5 changed files with 47 additions and 11 deletions.
30 changes: 23 additions & 7 deletions shared/lib_irradproc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1692,7 +1692,7 @@ perez(double, double dn, double df, double alb, double inc, double tilt, double
}
}

void ineichen(double clearsky_results[3], double apparent_zenith, double absolute_airmass, double linke_turbidity = 1.0, double altitude = 0.0, double dni_extra = 1364.0, bool perez_enhancement = false) {
void ineichen(double clearsky_results[3], double apparent_zenith, int month, int day, double pressure = 101325.0, double linke_turbidity = 1.0, double altitude = 0.0, double dni_extra = 1364.0, bool perez_enhancement = false) {
double cos_zenith = Max(cosd(apparent_zenith), 0);
double tl = linke_turbidity;

Expand All @@ -1701,15 +1701,25 @@ void ineichen(double clearsky_results[3], double apparent_zenith, double absolut
double cg1 = 5.09e-5 * altitude + 0.868;
double cg2 = 3.92e-5 * altitude + 0.0387;

double ghi = exp(-cg2 * absolute_airmass * (fh1 + fh2 * (tl - 1)));
if (perez_enhancement) ghi *= exp(0.01 * pow(absolute_airmass, 1.8));
//double am = Min(15.25, 1.0 / (cosd(apparent_zenith) + 0.15 * (pow(93.9 - apparent_zenith, -1.253)))); // air mass
double am = Min(15.25, 1.0 / (cosd(apparent_zenith) + 0.50572 * (pow(6.07995 + (90 - apparent_zenith), -1.6364)))); // air mass kastenyoung1989 pvlib

ghi = cg1 * dni_extra * cos_zenith * tl / tl * Max(ghi, 0);
double abs_am = am * pressure / 101325.0;

double ghi = exp(-cg2 * abs_am * (fh1 + fh2 * (tl - 1)));
if (perez_enhancement) ghi *= exp(0.01 * pow(abs_am, 1.8));

double Gon = 1367 * (1 + 0.033 * cos(360.0 / 365.0 * day_of_year(month, day) * M_PI /
180));

if (dni_extra != 0) Gon = dni_extra;

ghi = cg1 * Gon * cos_zenith * tl / tl * Max(ghi, 0);

double b = 0.664 + 0.163 / fh1;

double bnci = b * exp(-0.09 * absolute_airmass * (tl - 1));
bnci = dni_extra * Max(bnci, 0);
double bnci = b * exp(-0.09 * abs_am * (tl - 1));
bnci = Gon * Max(bnci, 0);

double bnci_2 = ((1 - (0.1 - 0.2 * exp(-tl)) / (0.1 + 0.882 / fh1)) / cos_zenith);
bnci_2 = ghi * Min(Max(bnci_2, 0), 1e20);
Expand Down Expand Up @@ -1958,6 +1968,12 @@ void irrad::get_irrad(double *ghi, double *dni, double *dhi) {
*dhi = diffuseHorizontal;
}

void irrad::get_clearsky_irrad(double* ghi_cs, double* dni_cs, double* dhi_cs) {
*ghi_cs = clearskyIrradiance[0];
*dni_cs = clearskyIrradiance[1];
*dhi_cs = clearskyIrradiance[2];
}

void irrad::set_time(int y, int m, int d, int h, double min, double delt_hr) {
this->year = y;
this->month = m;
Expand Down Expand Up @@ -2160,7 +2176,7 @@ int irrad::calc() {

//clearsky
if (enableSubhourlyClipping) {
ineichen(clearskyIrradiance, RTOD * sunAnglesRadians[1], 1.5, 1.0, elevation);
ineichen(clearskyIrradiance, RTOD * sunAnglesRadians[1], month, day, pressure * 100.0, 1.0, elevation, 0, true);
}


Expand Down
4 changes: 3 additions & 1 deletion shared/lib_irradproc.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -778,7 +778,7 @@ void perez(double hextra, double dn, double df, double alb, double inc, double t
* \param[out] clearsky_results[1] clear sky dni
* \param[out] clearsky_results[2] clear sky dhi
*/
void ineichen(double clearsky_results[3], double apparent_zenith, double absolute_airmass, double linke_turbidity, double altitude, double dni_extra, bool perez_enhancement);
void ineichen(double clearsky_results[3], double apparent_zenith, int month, int day, double pressure, double linke_turbidity, double altitude, double dni_extra, bool perez_enhancement);
/**
* Isotropic sky model for diffuse irradiance on a tilted surface, see also perez(), hdkr().
*
Expand Down Expand Up @@ -1168,6 +1168,8 @@ class irrad
/// Return the front-side irradiance components
void get_irrad(double* ghi, double* dni, double* dhi);

void get_clearsky_irrad(double* ghi_cs, double* dni_cs, double* dhi_cs);

/// Return the effective hour and fraction used for the sun position calculation
double get_sunpos_calc_hour();

Expand Down
4 changes: 4 additions & 0 deletions shared/lib_pv_io_manager.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -933,6 +933,10 @@ void PVSystem_IO::AllocateOutputs(compute_module* cm)
p_inverterACOutputPreLoss = cm->allocate("ac_gross", numberOfWeatherFileRecords);
p_acWiringLoss = cm->allocate("ac_wiring_loss", numberOfWeatherFileRecords);
p_ClippingPotential = cm->allocate("clipping_potential", numberOfWeatherFileRecords);
p_CPBin = cm->allocate("clipping_potential_bin", numberOfWeatherFileRecords);
//p_DNIIndex = cm->allocate("dni_index", numberOfWeatherFileRecords);
p_DNIIndexBin = cm->allocate("dni_index_bin", numberOfWeatherFileRecords);

p_transmissionLoss = cm->allocate("ac_transmission_loss", numberOfWeatherFileRecords);
p_acPerfAdjLoss = cm->allocate("ac_perf_adj_loss", numberOfWeatherFileRecords);
p_acLifetimeLoss = cm->allocate("ac_lifetime_loss", numberOfWeatherFileRecords);
Expand Down
3 changes: 3 additions & 0 deletions shared/lib_pv_io_manager.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -422,6 +422,9 @@ struct PVSystem_IO
ssc_number_t *p_subhourlyClippingLoss;
ssc_number_t* p_subhourlyClippingLossFactor;
ssc_number_t* p_ClippingPotential;
//ssc_number_t* p_DNIIndex;
ssc_number_t* p_CPBin;
ssc_number_t* p_DNIIndexBin;
};

/**
Expand Down
17 changes: 14 additions & 3 deletions ssc/cmod_pvsamv1.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -652,6 +652,7 @@ static var_info _cm_vtab_pvsamv1[] = {
{ SSC_OUTPUT, SSC_ARRAY, "subarray1_dc_gross", "Subarray 1 DC power gross", "kW", "", "Time Series (Subarray 1)", "*", "", "" },
{ SSC_OUTPUT, SSC_ARRAY, "subarray1_voc", "Subarray 1 Open circuit DC voltage", "V", "", "Time Series (Subarray 1)", "", "", "" },
{ SSC_OUTPUT, SSC_ARRAY, "subarray1_isc", "Subarray 1 String short circuit DC current", "A", "", "Time Series (Subarray 1)", "", "", "" },
{ SSC_OUTPUT, SSC_ARRAY, "subarray1_dni_index", "Subarray 1 DNI Index", "W/m2", "", "Time Series (Subarray 1)", "", "", "" },

{ SSC_OUTPUT, SSC_ARRAY, "subarray2_surf_tilt", "Subarray 2 Surface tilt", "degrees", "", "Time Series (Subarray 2)", "", "", "" },
{ SSC_OUTPUT, SSC_ARRAY, "subarray2_surf_azi", "Subarray 2 Surface azimuth", "degrees", "", "Time Series (Subarray 2)", "", "", "" },
Expand Down Expand Up @@ -805,6 +806,8 @@ static var_info _cm_vtab_pvsamv1[] = {
{ SSC_OUTPUT, SSC_ARRAY, "ac_gross", "Inverter AC output power", "kW", "", "Time Series (Array)", "*", "", "" },
{ SSC_OUTPUT, SSC_ARRAY, "clipping_potential", "Clipping potential", "", "", "Time Series (Inverter)", "", "", "" },
{ SSC_OUTPUT, SSC_ARRAY, "subhourly_clipping_loss", "Subhourly clipping correction loss", "kW", "", "Time Series (Inverter)", "", "", "" },
{ SSC_OUTPUT, SSC_ARRAY, "clipping_potential_bin", "Clipping potential bin", "", "", "Time Series (Inverter)", "", "", "" },
{ SSC_OUTPUT, SSC_ARRAY, "dni_index_bin", "DNI Index bin", "", "", "Time Series (Inverter)", "", "", "" },

// transformer model outputs
{ SSC_OUTPUT, SSC_ARRAY, "xfmr_nll_ts", "Transformer no load loss", "kW", "", "Time Series (Transformer)", "", "", "" },
Expand Down Expand Up @@ -960,6 +963,7 @@ static var_info _cm_vtab_pvsamv1[] = {
{ SSC_OUTPUT, SSC_NUMBER, "annual_ac_inv_eff_loss_percent", "AC inverter efficiency loss", "%", "", "Loss", "", "", "" },
{ SSC_OUTPUT, SSC_NUMBER, "annual_ac_wiring_loss_percent", "AC wiring loss", "%", "", "Loss", "", "", "" },
{ SSC_OUTPUT, SSC_NUMBER, "annual_subhourly_clipping_loss_percent", "Subhourly clipping correction loss percent", "%", "", "Loss", "", "", "" },
{ SSC_OUTPUT, SSC_NUMBER, "nominal_annual_clipping_output", "Test output for nominal annual clipping output", "kWh", "", "", "", "", "" },

{ SSC_OUTPUT, SSC_NUMBER, "annual_transmission_loss_percent", "AC transmission loss", "%", "", "Loss", "", "", "" },
// { SSC_OUTPUT, SSC_NUMBER, "annual_ac_transformer_loss_percent", "AC step-up transformer loss", "%", "", "Loss", "", "", "" },
Expand Down Expand Up @@ -1532,8 +1536,11 @@ void cm_pvsamv1::exec()
// beam, skydiff, and grounddiff IN THE PLANE OF ARRAY (W/m2)
double ibeam, iskydiff, ignddiff;
double ibeam_csky, iskydiff_csky, ignddiff_csky;
double ghi_cs, dni_cs, dhi_cs;
double aoi, stilt, sazi, rot, btd;

irr.get_clearsky_irrad(&ghi_cs, &dni_cs, &dhi_cs);

// Ensure that the usePOAFromWF flag is false unless a reference cell has been used.
// This will later get forced to false if any shading has been applied (in any scenario)
// also this will also be forced to false if using the cec mcsp thermal model OR if using the spe module model with a diffuse util. factor < 1.0
Expand Down Expand Up @@ -1918,8 +1925,8 @@ void cm_pvsamv1::exec()
PVSystem->p_poaBeamFrontCS[nn][idx] = (ssc_number_t)ibeam_csky;
PVSystem->p_poaDiffuseFrontCS[nn][idx] = (ssc_number_t)(iskydiff_csky);
PVSystem->p_poaDiffuseFrontCS[nn][idx] = (ssc_number_t)(ignddiff_csky);
if (ibeam_csky != 0) {
PVSystem->p_DNIIndex[nn][idx] = (ssc_number_t)(ibeam / ibeam_csky);
if (dni_cs != 0) {
PVSystem->p_DNIIndex[nn][idx] = (ssc_number_t)(Irradiance->p_weatherFileDNI[idx] / dni_cs);
}
else {
PVSystem->p_DNIIndex[nn][idx] = 0;
Expand Down Expand Up @@ -2591,6 +2598,7 @@ void cm_pvsamv1::exec()
annual_subhourly_clipping_loss += sub_clipping_matrix.at(dni_row, clip_pot_col);
*/
}
assign("nominal_annual_clipping_output", nominal_annual_clipping_output);
}

for (size_t iyear = 0; iyear < nyears; iyear++)
Expand Down Expand Up @@ -2714,6 +2722,7 @@ void cm_pvsamv1::exec()
double clip_pot = (dcPower_kW_csky - paco ) / (paco);

PVSystem->p_ClippingPotential[idx] = clip_pot;
//PVSystem->p_DNIIndex[idx] = dni_clearness_index;

//Lookup matrix for percentage effect based on DNI index, Clipping potential

Expand All @@ -2731,6 +2740,7 @@ void cm_pvsamv1::exec()
}
}
}
PVSystem->p_DNIIndexBin[idx] = dni_row;

//Clipping potential indexing
if (clip_pot < sub_clipping_matrix.at(0, 1)) clip_pot_col = 1;
Expand All @@ -2742,9 +2752,10 @@ void cm_pvsamv1::exec()
}
}
}
PVSystem->p_CPBin[idx] = clip_pot_col;

//acpwr_gross *= (1 - sub_clipping_matrix.at(dni_row, clip_pot_col));
if (dcPower_kW_csky > 0.0) {
if (dcPower_kW > 0.0) {
ac_subhourlyclipping_loss = sub_clipping_matrix.at(dni_row, clip_pot_col) * nominal_annual_clipping_output;
//ac_subhourlyclipping_loss = sub_clipping_matrix.at(dni_row, clip_pot_col) * sharedInverter->powerAC_kW_clipping;
}
Expand Down

0 comments on commit d774469

Please sign in to comment.