change implementation of optional user-defined thrust coefficient cur…

…ve to belong to wind turbine class so that it will be available to all wake models

shared/lib_windwakemodel.cpp

@@ -49,6 +49,18 @@ bool windTurbine::setPowerCurve(std::vector<double> windSpeeds, std::vector<doub
 	return 1;
 }
 
+bool windTurbine::setCtCurve(std::vector<double> thrustCoeffCurve)
+{
+    if (thrustCoeffCurve.size() != powerCurveWS.size())
+    {
+        errDetails = "Coefficient of thrust curve must have the same number of values as the power curve wind speeds";
+        return 0;
+    }
+    ctCurve.resize(thrustCoeffCurve.size());
+    ctCurve = thrustCoeffCurve;
+    return 1;
+}
+
 double windTurbine::tipSpeedRatio(double windSpeed)
 {
 	if (powerCurveRPM[0] == -1) return 7.0;
@@ -132,6 +144,14 @@ void windTurbine::turbinePower(double windVelocity, double airDensity, double *t
 		*turbineOutput = out_pwr;
 		if (fPowerCoefficient >= 0.0)
 			*thrustCoefficient = max_of(0.0, -1.453989e-2 + 1.473506*fPowerCoefficient - 2.330823*pow(fPowerCoefficient, 2) + 3.885123*pow(fPowerCoefficient, 3));
+
+        // overwrite the coefficient of thrust if it has been specified by the user
+        // if it has not been specified by the user, the thrust curve vector is {0.}
+        if (ctCurve.size() != 1)
+        {
+            // do something here
+        }
+
 	} // out_pwr > (rated power * 0.001)
 
 	return;
@@ -214,13 +234,6 @@ double parkWakeModel::delta_V_Park(double Uo, double Ui, double distCrosswind, d
 	// bound the coeff of thrust
 	double Ct = max_of(min_of(0.999, dThrustCoeff), minThrustCoeff);
 
-    // overwrite the coefficient of thrust if it has been specified by the user
-    // if it has not been specified by the user, the thrust curve vector is {0.}
-    if (ctCurve.size() != 1)
-    {
-        // do something here
-    }
-
 	double k = wakeDecayConstant;
 
 	double dRadiusOfWake = dRadiusUpstream + (k * distDownwind); // radius of circle formed by wake from upwind rotor

shared/lib_windwakemodel.h

@@ -47,10 +47,13 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 class windTurbine
 {
 private:
-	std::vector<double> powerCurveWS,			// windspeed: x-axis on turbine power curve
-						powerCurveKW,			// power output: y-axis
-						densityCorrectedWS,
-						powerCurveRPM;
+    std::vector<double> powerCurveWS,			// windspeed: x-axis on turbine power curve
+                        powerCurveKW,			// power output: y-axis
+                        densityCorrectedWS,
+                        powerCurveRPM;
+    // vector that stores the optional coefficient of thrust curve input
+    // set to a default value of length 1 to mean that it's not assigned, and check for that length before using it
+    std::vector<double> ctCurve = { 0. };
 	double cutInSpeed;
 	double previousAirDensity;
 public:
@@ -70,6 +73,7 @@ class windTurbine
         previousAirDensity = physics::AIR_DENSITY_SEA_LEVEL;
 	}
 	bool setPowerCurve(std::vector<double> windSpeeds, std::vector<double> powerOutput);
+    bool setCtCurve(std::vector<double> thrustCoeffCurve);
 
 	double tipSpeedRatio(double windSpeed);
 
@@ -153,15 +157,16 @@ class parkWakeModel : public wakeModelBase{
 	double rotorDiameter;
 	double wakeDecayConstant = 0.07,
 		   minThrustCoeff = 0.02;
-    std::vector<double> ctCurve; //vector that stores the optional coefficient of thrust curve input
 	double delta_V_Park(double dVelFreeStream, double dVelUpwind, double dDistCrossWind, double dDistDownWind, double dRadiusUpstream, double dRadiusDownstream, double dThrustCoeff);
 	double circle_overlap(double dist_center_to_center, double rad1, double rad2);
 
 public:
 	parkWakeModel(){ nTurbines = 0; }
-    parkWakeModel(size_t numberOfTurbinesInFarm, windTurbine* wt, double wdc, std::vector<double> ctc)
+    parkWakeModel(size_t numberOfTurbinesInFarm, windTurbine* wt, double wdc)
     {
-        nTurbines = numberOfTurbinesInFarm; wTurbine = wt; setWakeDecayConstant(wdc); ctCurve = ctc;
+        nTurbines = numberOfTurbinesInFarm;
+        wTurbine = wt;
+        setWakeDecayConstant(wdc);
     }
 	virtual ~parkWakeModel() {};
 	std::string getModelName() override { return "Park"; }

ssc/cmod_windpower.cpp

@@ -74,7 +74,7 @@ static var_info _cm_vtab_windpower[] = {
 
     // optional SDK only wake loss inputs
     { SSC_INPUT  , SSC_NUMBER , "wake_loss_multiplier"               , "Multiplier for the calculated wake loss"  , ""        ,">1 increases loss, <1 decreases loss", "Farm"                                 , ""                                                , "MIN=0"                                           , "" } ,
-    { SSC_INOUT  , SSC_ARRAY  , "wind_turbine_ct_curve"              , "Park model coeff of thrust curve vs WS"   , ""        ,"uses same wind speeds as power curve", "Turbine"                              , ""                                                , "LENGTH_EQUAL=wind_turbine_powercurve_windspeeds" , "GROUP=WTPCD" } ,
+    { SSC_INOUT  , SSC_ARRAY  , "wind_turbine_ct_curve"              , "User-defined Ct curve vs WS for wake models", ""      ,"uses same wind speeds as power curve", "Turbine"                              , ""                                                , "LENGTH_EQUAL=wind_turbine_powercurve_windspeeds" , "GROUP=WTPCD" } ,
 
     { SSC_INPUT  , SSC_NUMBER , "wake_int_loss"                      , "Constant Wake Model, internal wake loss"  , "%"       ,""                                    , "Losses"                               , "wind_farm_wake_model=3"                          , "MIN=0,MAX=100"                                   , "" } ,
     { SSC_INPUT  , SSC_NUMBER , "wake_ext_loss"                      , "External Wake loss"                       , "%"       ,""                                    , "Losses"                               , "?=0"                                             , "MIN=0,MAX=100"                                   , "" } ,
@@ -302,6 +302,17 @@ void cm_windpower::exec()
 	}
 	wt.setPowerCurve(windSpeeds, powerOutput);
 
+    // get optional thrust curve for wind turbine
+    if (is_assigned("wind_turbine_ct_curve"))
+    {
+        size_t* ctCurveLength = 0;
+        ssc_number_t* ctc = as_array("wind_turbine_ct_curve", ctCurveLength);
+        std::vector<double> ct_curve(*ctCurveLength);
+        for (size_t i = 0; i < *ctCurveLength; i++)
+            ct_curve[i] = ctc[i];
+        wt.setCtCurve(ct_curve);
+    }
+
 	// create windPowerCalculator using windTurbine
 	windPowerCalculator wpc;
 	wpc.windTurb = &wt;
@@ -406,26 +417,13 @@ void cm_windpower::exec()
         wakeModel = std::make_shared<simpleWakeModel>(simpleWakeModel(wpc.nTurbines, &wt));
     else if (wakeModelChoice == PARK)
     {
-        // get optional thrust curve
-        std::vector<double> ct_curve = { 0. }; //initialize this to length 1: a power curve length of 1 throws an exec error above & the length of ct curve is tied to the power curve, so this is a safe null value
-        if (is_assigned("wind_turbine_ct_curve"))
-        {
-            size_t *ctCurveLength = 0;
-            ssc_number_t *ctc = as_array("wind_turbine_ct_curve", ctCurveLength);
-            if (*ctCurveLength != wt.powerCurveArrayLength)
-                throw exec_error("windpower", "Thrust curve must have same number of values as the power curve");
-            ct_curve.resize(*ctCurveLength);
-            for (size_t i = 0; i < *ctCurveLength; i++)
-                ct_curve[i] = ctc[i];
-        }
-
         // get optional wake decay constant
         double wdc = 0.07; //this is the default value
         if (is_assigned("park_wake_decay_constant"))
             wdc = as_double("park_wake_decay_constant");
 
         // create the wake model
-        wakeModel = std::make_shared<parkWakeModel>(parkWakeModel(wpc.nTurbines, &wt, wdc, ct_curve));
+        wakeModel = std::make_shared<parkWakeModel>(parkWakeModel(wpc.nTurbines, &wt, wdc));
     }
     else if (wakeModelChoice == EDDYVISCOSITY)
     {

test/shared_test/lib_windwakemodel_test.h

@@ -154,7 +154,7 @@ class parkWakeModelTest : public ::testing::Test{
 		windSpeed.resize(numberTurbines);
 		turbIntensity.resize(numberTurbines, 0.1);
 		createDefaultTurbine(&wt);
-        pm = parkWakeModel(numberTurbines, &wt, 0.07, { 0. }); //use default values for the last two inputs
+        pm = parkWakeModel(numberTurbines, &wt, 0.07);
 		for (int i = 0; i < numberTurbines; i++){
 			windSpeed[i] = 10.;
 		}