Merge branch 'develop' into iph_fin_cmods

NREL · Oct 16, 2024 · 36d08c9 · 36d08c9
2 parents d39923e + 8e8c06f
commit 36d08c9
Show file tree

Hide file tree

Showing 5 changed files with 261 additions and 56 deletions.
diff --git a/shared/lib_windwakemodel.cpp b/shared/lib_windwakemodel.cpp
@@ -49,6 +49,23 @@ bool windTurbine::setPowerCurve(std::vector<double> windSpeeds, std::vector<doub
 	return 1;
 }
 
+bool windTurbine::setCtCurve(std::vector<double> thrustCoeffCurve)
+{
+    if (powerCurveWS.empty())
+    {
+        errDetails = "Power curve must be set before thrust curve.";
+        return 0;
+    }
+    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;
@@ -104,9 +121,9 @@ void windTurbine::turbinePower(double windVelocity, double airDensity, double *t
 
 	// Find power from turbine power curve
 	double out_pwr = 0.0;
+    int j = 1; // an index for the correct wind speed in the power curve for interpolations
 	if ((windVelocity > densityCorrectedWS[0]) && (windVelocity < densityCorrectedWS[powerCurveArrayLength - 1]))
 	{
-		int j = 1;
 		while (densityCorrectedWS[j] <= windVelocity)
 			j++; // find first m_adPowerCurveWS > windVelocity
 
@@ -132,6 +149,16 @@ 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)
+        {
+            //interpolate to right value of Ct based on wind speed, like for power curve
+            //can use the same index for the wind speed array that was determined above
+            *thrustCoefficient = util::interpolate(densityCorrectedWS[j - 1], ctCurve[j - 1], densityCorrectedWS[j], ctCurve[j], windVelocity);
+        }
+
 	} // out_pwr > (rated power * 0.001)
 
 	return;
@@ -214,7 +241,7 @@ 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);
 
-	double k = wakeDecayCoefficient;
+	double k = wakeDecayConstant;
 
 	double dRadiusOfWake = dRadiusUpstream + (k * distDownwind); // radius of circle formed by wake from upwind rotor
 	double dAreaOverlap = circle_overlap(distCrosswind, dRadiusDownstream, dRadiusOfWake);

diff --git a/shared/lib_windwakemodel.h b/shared/lib_windwakemodel.h
@@ -47,10 +47,11 @@ 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,
+                        ctCurve;                // vector that stores the optional coefficient of thrust curve input
 	double cutInSpeed;
 	double previousAirDensity;
 public:
@@ -68,8 +69,10 @@ class windTurbine
 		hubHeight = -999;
 		rotorDiameter = -999;
         previousAirDensity = physics::AIR_DENSITY_SEA_LEVEL;
+        ctCurve = { 0.0 }; // set to a default value of length 1 to mean that it's not assigned, and check for that length before using it
 	}
 	bool setPowerCurve(std::vector<double> windSpeeds, std::vector<double> powerOutput);
+    bool setCtCurve(std::vector<double> thrustCoeffCurve);
 
 	double tipSpeedRatio(double windSpeed);
 
@@ -151,17 +154,23 @@ class simpleWakeModel : public wakeModelBase{
 class parkWakeModel : public wakeModelBase{
 private:
 	double rotorDiameter;
-	double wakeDecayCoefficient = 0.07,
+	double wakeDecayConstant = 0.07,
 		   minThrustCoeff = 0.02;
 	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){ nTurbines = numberOfTurbinesInFarm; wTurbine = wt; }
+    parkWakeModel(size_t numberOfTurbinesInFarm, windTurbine* wt, double wdc)
+    {
+        nTurbines = numberOfTurbinesInFarm;
+        wTurbine = wt;
+        setWakeDecayConstant(wdc);
+    }
 	virtual ~parkWakeModel() {};
 	std::string getModelName() override { return "Park"; }
 	void setRotorDiameter(double d){ rotorDiameter = d; }
+    void setWakeDecayConstant(double w) { wakeDecayConstant = w; }
 	void wakeCalculations(
 		/*INPUTS*/
 		const double airDensity,					// not used in this model
@@ -302,4 +311,4 @@ class constantWakeModel : public wakeModelBase
             ) override;
 };
 
-#endif
+#endif