Line data Source code
1 : // EnergyPlus, Copyright (c) 1996-2023, The Board of Trustees of the University of Illinois,
2 : // The Regents of the University of California, through Lawrence Berkeley National Laboratory
3 : // (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge
4 : // National Laboratory, managed by UT-Battelle, Alliance for Sustainable Energy, LLC, and other
5 : // contributors. All rights reserved.
6 : //
7 : // NOTICE: This Software was developed under funding from the U.S. Department of Energy and the
8 : // U.S. Government consequently retains certain rights. As such, the U.S. Government has been
9 : // granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable,
10 : // worldwide license in the Software to reproduce, distribute copies to the public, prepare
11 : // derivative works, and perform publicly and display publicly, and to permit others to do so.
12 : //
13 : // Redistribution and use in source and binary forms, with or without modification, are permitted
14 : // provided that the following conditions are met:
15 : //
16 : // (1) Redistributions of source code must retain the above copyright notice, this list of
17 : // conditions and the following disclaimer.
18 : //
19 : // (2) Redistributions in binary form must reproduce the above copyright notice, this list of
20 : // conditions and the following disclaimer in the documentation and/or other materials
21 : // provided with the distribution.
22 : //
23 : // (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory,
24 : // the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be
25 : // used to endorse or promote products derived from this software without specific prior
26 : // written permission.
27 : //
28 : // (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form
29 : // without changes from the version obtained under this License, or (ii) Licensee makes a
30 : // reference solely to the software portion of its product, Licensee must refer to the
31 : // software as "EnergyPlus version X" software, where "X" is the version number Licensee
32 : // obtained under this License and may not use a different name for the software. Except as
33 : // specifically required in this Section (4), Licensee shall not use in a company name, a
34 : // product name, in advertising, publicity, or other promotional activities any name, trade
35 : // name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly
36 : // similar designation, without the U.S. Department of Energy's prior written consent.
37 : //
38 : // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
39 : // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
40 : // AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
41 : // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
42 : // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
43 : // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
44 : // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
45 : // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
46 : // POSSIBILITY OF SUCH DAMAGE.
47 :
48 : // C++ Headers
49 : #include <cassert>
50 : #include <cmath>
51 :
52 : // ObjexxFCL Headers
53 : #include <ObjexxFCL/Array.functions.hh>
54 : #include <ObjexxFCL/Fmath.hh>
55 : #include <ObjexxFCL/string.functions.hh>
56 :
57 : // EnergyPlus Headers
58 : #include <EnergyPlus/CommandLineInterface.hh>
59 : #include <EnergyPlus/Data/EnergyPlusData.hh>
60 : #include <EnergyPlus/DataEnvironment.hh>
61 : #include <EnergyPlus/DataHVACGlobals.hh>
62 : #include <EnergyPlus/DataIPShortCuts.hh>
63 : #include <EnergyPlus/FileSystem.hh>
64 : #include <EnergyPlus/General.hh>
65 : #include <EnergyPlus/InputProcessing/InputProcessor.hh>
66 : #include <EnergyPlus/OutputProcessor.hh>
67 : #include <EnergyPlus/Psychrometrics.hh>
68 : #include <EnergyPlus/ScheduleManager.hh>
69 : #include <EnergyPlus/StringUtilities.hh>
70 : #include <EnergyPlus/UtilityRoutines.hh>
71 : #include <EnergyPlus/WindTurbine.hh>
72 :
73 : namespace EnergyPlus {
74 :
75 : // (ref: Object: Generator:WindTurbine)
76 :
77 : namespace WindTurbine {
78 : // Module containing the data for wind turbine system
79 :
80 : // MODULE INFORMATION:
81 : // AUTHOR Daeho Kang
82 : // DATE WRITTEN October 2009
83 : // MODIFIED na
84 : // RE-ENGINEERED na
85 :
86 : // PURPOSE OF THIS MODULE:
87 : // This module is to calculate the electrical power output that wind turbine systems produce.
88 : // Both horizontal and vertical axis wind turbine systems are modeled.
89 :
90 : // REFERENCES:
91 : // Sathyajith Mathew. 2006. Wind Energy: Fundamental, Resource Analysis and Economics. Springer,
92 : // Chap. 2, pp. 11-15
93 : // Mazharul Islam, David S.K. Ting, and Amir Fartaj. 2008. Aerodynamic Models for Darrieus-type sSraight-bladed
94 : // Vertical Axis Wind Turbines. Renewable & Sustainable Energy Reviews, Volume 12, pp.1087-1109
95 :
96 : constexpr std::array<std::string_view, static_cast<int>(ControlType::Num)> ControlNamesUC{
97 : "FIXEDSPEEDFIXEDPITCH",
98 : "FIXEDSPEEDVARIABLEPITCH",
99 : "VARIABLESPEEDFIXEDPITCH",
100 : "VARIABLESPEEDVARIABLEPITCH",
101 : };
102 :
103 : constexpr std::array<std::string_view, static_cast<int>(RotorType::Num)> RotorNamesUC{
104 : "HORIZONTALAXISWINDTURBINE",
105 : "VERTICALAXISWINDTURBINE",
106 : };
107 :
108 16344 : void SimWindTurbine(EnergyPlusData &state,
109 : [[maybe_unused]] GeneratorType const GeneratorType, // Type of Generator
110 : std::string const &GeneratorName, // User specified name of Generator
111 : int &GeneratorIndex, // Generator index
112 : bool const RunFlag, // ON or OFF
113 : [[maybe_unused]] Real64 const WTLoad // Electrical load on WT (not used)
114 : )
115 : {
116 :
117 : // SUBROUTINE INFORMATION:
118 : // AUTHOR Daeho Kang
119 : // DATE WRITTEN October 2009
120 : // MODIFIED na
121 : // RE-ENGINEERED na
122 :
123 : // PURPOSE OF THIS SUBROUTINE:
124 : // This subroutine manages the simulation of wind turbine component.
125 : // This drivers manages the calls to all of the other drivers and simulation algorithms.
126 :
127 : // Using/Aliasing
128 :
129 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
130 : int WindTurbineNum;
131 : // Obtains and allocates heat balance related parameters from input
132 :
133 16344 : if (state.dataWindTurbine->GetInputFlag) {
134 1 : GetWindTurbineInput(state);
135 1 : state.dataWindTurbine->GetInputFlag = false;
136 : }
137 :
138 16344 : if (GeneratorIndex == 0) {
139 3 : WindTurbineNum = UtilityRoutines::FindItemInList(GeneratorName, state.dataWindTurbine->WindTurbineSys);
140 3 : if (WindTurbineNum == 0) {
141 0 : ShowFatalError(state, "SimWindTurbine: Specified Generator not one of Valid Wind Turbine Generators " + GeneratorName);
142 : }
143 3 : GeneratorIndex = WindTurbineNum;
144 : } else {
145 16341 : WindTurbineNum = GeneratorIndex;
146 16341 : int NumWindTurbines = (int)state.dataWindTurbine->WindTurbineSys.size();
147 16341 : if (WindTurbineNum > NumWindTurbines || WindTurbineNum < 1) {
148 0 : ShowFatalError(state,
149 0 : format("SimWindTurbine: Invalid GeneratorIndex passed={}, Number of Wind Turbine Generators={}, Generator name={}",
150 : WindTurbineNum,
151 : NumWindTurbines,
152 0 : GeneratorName));
153 : }
154 16341 : if (GeneratorName != state.dataWindTurbine->WindTurbineSys(WindTurbineNum).Name) {
155 0 : ShowFatalError(state,
156 0 : format("SimMWindTurbine: Invalid GeneratorIndex passed={}, Generator name={}, stored Generator Name for that index={}",
157 : WindTurbineNum,
158 : GeneratorName,
159 0 : state.dataWindTurbine->WindTurbineSys(WindTurbineNum).Name));
160 : }
161 : }
162 :
163 16344 : InitWindTurbine(state, WindTurbineNum);
164 :
165 16344 : CalcWindTurbine(state, WindTurbineNum, RunFlag);
166 :
167 16344 : ReportWindTurbine(state, WindTurbineNum);
168 16344 : }
169 :
170 16344 : void GetWTGeneratorResults(EnergyPlusData &state,
171 : [[maybe_unused]] GeneratorType const GeneratorType, // Type of Generator
172 : int const GeneratorIndex, // Generator number
173 : Real64 &GeneratorPower, // Electrical power
174 : Real64 &GeneratorEnergy, // Electrical energy
175 : Real64 &ThermalPower,
176 : Real64 &ThermalEnergy)
177 : {
178 :
179 : // SUBROUTINE INFORMATION:
180 : // AUTHOR B. Griffith
181 : // DATE WRITTEN Aug. 2008
182 : // MODIFIED D Kang, October 2009 for Wind Turbine
183 : // RE-ENGINEERED na
184 :
185 : // PURPOSE OF THIS SUBROUTINE:
186 : // This subroutine provides a "get" method to collect results for individual electic load centers.
187 :
188 16344 : GeneratorPower = state.dataWindTurbine->WindTurbineSys(GeneratorIndex).Power;
189 16344 : GeneratorEnergy = state.dataWindTurbine->WindTurbineSys(GeneratorIndex).Energy;
190 :
191 : // Thermal energy is ignored
192 16344 : ThermalPower = 0.0;
193 16344 : ThermalEnergy = 0.0;
194 16344 : }
195 :
196 1 : void GetWindTurbineInput(EnergyPlusData &state)
197 : {
198 :
199 : // SUBROUTINE INFORMATION:
200 : // AUTHOR Daeho Kang
201 : // DATE WRITTEN October 2009
202 : // MODIFIED na
203 : // RE-ENGINEERED na
204 :
205 : // PURPOSE OF THIS SUBROUTINE:
206 : // This subroutine gets input data for wind turbine components
207 : // and stores it in the wind turbine data structure.
208 :
209 : // Using/Aliasing
210 :
211 : using ScheduleManager::GetScheduleIndex;
212 :
213 : // SUBROUTINE PARAMETER DEFINITIONS:
214 1 : static std::string const CurrentModuleObject("Generator:WindTurbine");
215 1 : Real64 constexpr SysEffDefault(0.835); // Default value of overall system efficiency
216 1 : Real64 constexpr MaxTSR(12.0); // Maximum tip speed ratio
217 1 : Real64 constexpr DefaultPC(0.25); // Default power coefficient
218 1 : Real64 constexpr MaxPowerCoeff(0.59); // Maximum power coefficient
219 1 : Real64 constexpr DefaultH(50.0); // Default of height for local wind speed
220 :
221 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
222 1 : bool ErrorsFound(false); // If errors detected in input
223 : int WindTurbineNum; // Wind turbine number
224 : int NumAlphas; // Number of Alphas for each GetobjectItem call
225 : int NumNumbers; // Number of Numbers for each GetobjectItem call
226 : int NumArgs;
227 : int IOStat;
228 2 : Array1D_string cAlphaArgs; // Alpha input items for object
229 2 : Array1D_string cAlphaFields; // Alpha field names
230 2 : Array1D_string cNumericFields; // Numeric field names
231 2 : Array1D<Real64> rNumericArgs; // Numeric input items for object
232 2 : Array1D_bool lAlphaBlanks; // Logical array, alpha field input BLANK = .TRUE.
233 2 : Array1D_bool lNumericBlanks; // Logical array, numeric field input BLANK = .TRUE.
234 :
235 : // Initializations and allocations
236 1 : state.dataInputProcessing->inputProcessor->getObjectDefMaxArgs(state, CurrentModuleObject, NumArgs, NumAlphas, NumNumbers);
237 1 : cAlphaArgs.allocate(NumAlphas);
238 1 : cAlphaFields.allocate(NumAlphas);
239 1 : cNumericFields.allocate(NumNumbers);
240 1 : rNumericArgs.dimension(NumNumbers, 0.0);
241 1 : lAlphaBlanks.dimension(NumAlphas, true);
242 1 : lNumericBlanks.dimension(NumNumbers, true);
243 :
244 1 : int NumWindTurbines = state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, CurrentModuleObject);
245 :
246 1 : state.dataWindTurbine->WindTurbineSys.allocate(NumWindTurbines);
247 :
248 4 : for (WindTurbineNum = 1; WindTurbineNum <= NumWindTurbines; ++WindTurbineNum) {
249 :
250 9 : state.dataInputProcessing->inputProcessor->getObjectItem(state,
251 : CurrentModuleObject,
252 : WindTurbineNum,
253 3 : state.dataIPShortCut->cAlphaArgs,
254 : NumAlphas,
255 3 : state.dataIPShortCut->rNumericArgs,
256 : NumNumbers,
257 : IOStat,
258 : lNumericBlanks,
259 : lAlphaBlanks,
260 : cAlphaFields,
261 : cNumericFields);
262 3 : UtilityRoutines::IsNameEmpty(state, state.dataIPShortCut->cAlphaArgs(1), CurrentModuleObject, ErrorsFound);
263 :
264 3 : auto &windTurbine = state.dataWindTurbine->WindTurbineSys(WindTurbineNum);
265 :
266 3 : windTurbine.Name = state.dataIPShortCut->cAlphaArgs(1); // Name of wind turbine
267 :
268 3 : windTurbine.Schedule = state.dataIPShortCut->cAlphaArgs(2); // Get schedule
269 3 : if (lAlphaBlanks(2)) {
270 0 : windTurbine.SchedPtr = DataGlobalConstants::ScheduleAlwaysOn;
271 : } else {
272 3 : windTurbine.SchedPtr = GetScheduleIndex(state, state.dataIPShortCut->cAlphaArgs(2));
273 3 : if (windTurbine.SchedPtr == 0) {
274 0 : ShowSevereError(state,
275 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cAlphaFields(2) + "=\"" +
276 0 : state.dataIPShortCut->cAlphaArgs(2) + "\" not found.");
277 0 : ErrorsFound = true;
278 : }
279 : }
280 : // Select rotor type
281 3 : windTurbine.rotorType = static_cast<RotorType>(
282 6 : getEnumerationValue(WindTurbine::RotorNamesUC, UtilityRoutines::MakeUPPERCase(state.dataIPShortCut->cAlphaArgs(3))));
283 3 : if (windTurbine.rotorType == RotorType::Invalid) {
284 0 : if (state.dataIPShortCut->cAlphaArgs(3).empty()) {
285 0 : windTurbine.rotorType = RotorType::HorizontalAxis;
286 : } else {
287 0 : ShowSevereError(state,
288 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cAlphaFields(3) + "=\"" +
289 0 : state.dataIPShortCut->cAlphaArgs(3) + "\".");
290 0 : ErrorsFound = true;
291 : }
292 : }
293 :
294 : // Select control type
295 3 : windTurbine.controlType = static_cast<ControlType>(
296 6 : getEnumerationValue(WindTurbine::ControlNamesUC, UtilityRoutines::MakeUPPERCase(state.dataIPShortCut->cAlphaArgs(4))));
297 3 : if (windTurbine.controlType == ControlType::Invalid) {
298 0 : if (state.dataIPShortCut->cAlphaArgs(4).empty()) {
299 0 : windTurbine.controlType = ControlType::VariableSpeedVariablePitch;
300 : } else {
301 0 : ShowSevereError(state,
302 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cAlphaFields(4) + "=\"" +
303 0 : state.dataIPShortCut->cAlphaArgs(4) + "\".");
304 0 : ErrorsFound = true;
305 : }
306 : }
307 :
308 3 : windTurbine.RatedRotorSpeed = state.dataIPShortCut->rNumericArgs(1); // Maximum rotor speed in rpm
309 3 : if (windTurbine.RatedRotorSpeed <= 0.0) {
310 0 : if (lNumericBlanks(1)) {
311 0 : ShowSevereError(state,
312 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cNumericFields(1) +
313 : " is required but input is blank.");
314 : } else {
315 0 : ShowSevereError(state,
316 0 : format("{}=\"{}\" invalid {}=[{:.2R}] must be greater than zero.",
317 : CurrentModuleObject,
318 0 : state.dataIPShortCut->cAlphaArgs(1),
319 : cNumericFields(1),
320 0 : rNumericArgs(1)));
321 : }
322 0 : ErrorsFound = true;
323 : }
324 :
325 3 : windTurbine.RotorDiameter = state.dataIPShortCut->rNumericArgs(2); // Rotor diameter in m
326 3 : if (windTurbine.RotorDiameter <= 0.0) {
327 0 : if (lNumericBlanks(2)) {
328 0 : ShowSevereError(state,
329 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cNumericFields(2) +
330 : " is required but input is blank.");
331 : } else {
332 0 : ShowSevereError(state,
333 0 : format("{}=\"{}\" invalid {}=[{:.1R}] must be greater than zero.",
334 : CurrentModuleObject,
335 0 : state.dataIPShortCut->cAlphaArgs(1),
336 : cNumericFields(2),
337 0 : rNumericArgs(2)));
338 : }
339 0 : ErrorsFound = true;
340 : }
341 :
342 3 : windTurbine.RotorHeight = state.dataIPShortCut->rNumericArgs(3); // Overall height of the rotor
343 3 : if (windTurbine.RotorHeight <= 0.0) {
344 0 : if (lNumericBlanks(3)) {
345 0 : ShowSevereError(state,
346 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cNumericFields(3) +
347 : " is required but input is blank.");
348 : } else {
349 0 : ShowSevereError(state,
350 0 : format("{}=\"{}\" invalid {}=[{:.1R}] must be greater than zero.",
351 : CurrentModuleObject,
352 0 : state.dataIPShortCut->cAlphaArgs(1),
353 : cNumericFields(3),
354 0 : rNumericArgs(3)));
355 : }
356 0 : ErrorsFound = true;
357 : }
358 :
359 3 : windTurbine.NumOfBlade = state.dataIPShortCut->rNumericArgs(4); // Total number of blade
360 3 : if (windTurbine.NumOfBlade == 0) {
361 0 : ShowSevereError(state,
362 0 : format("{}=\"{}\" invalid {}=[{:.0R}] must be greater than zero.",
363 : CurrentModuleObject,
364 0 : state.dataIPShortCut->cAlphaArgs(1),
365 : cNumericFields(4),
366 0 : rNumericArgs(4)));
367 0 : ErrorsFound = true;
368 : }
369 :
370 3 : windTurbine.RatedPower = state.dataIPShortCut->rNumericArgs(5); // Rated average power
371 3 : if (windTurbine.RatedPower == 0.0) {
372 0 : if (lNumericBlanks(5)) {
373 0 : ShowSevereError(state,
374 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cNumericFields(5) +
375 : " is required but input is blank.");
376 : } else {
377 0 : ShowSevereError(state,
378 0 : format("{}=\"{}\" invalid {}=[{:.2R}] must be greater than zero.",
379 : CurrentModuleObject,
380 0 : state.dataIPShortCut->cAlphaArgs(1),
381 : cNumericFields(5),
382 0 : rNumericArgs(5)));
383 : }
384 0 : ErrorsFound = true;
385 : }
386 :
387 3 : windTurbine.RatedWindSpeed = state.dataIPShortCut->rNumericArgs(6); // Rated wind speed
388 3 : if (windTurbine.RatedWindSpeed == 0.0) {
389 0 : if (lNumericBlanks(6)) {
390 0 : ShowSevereError(state,
391 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cNumericFields(6) +
392 : " is required but input is blank.");
393 : } else {
394 0 : ShowSevereError(state,
395 0 : format("{}=\"{}\" invalid {}=[{:.2R}] must be greater than zero.",
396 : CurrentModuleObject,
397 0 : state.dataIPShortCut->cAlphaArgs(1),
398 : cNumericFields(6),
399 0 : rNumericArgs(6)));
400 : }
401 0 : ErrorsFound = true;
402 : }
403 :
404 3 : windTurbine.CutInSpeed = state.dataIPShortCut->rNumericArgs(7); // Minimum wind speed for system operation
405 3 : if (windTurbine.CutInSpeed == 0.0) {
406 0 : if (lNumericBlanks(7)) {
407 0 : ShowSevereError(state,
408 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cNumericFields(7) +
409 : " is required but input is blank.");
410 : } else {
411 0 : ShowSevereError(state,
412 0 : format("{}=\"{}\" invalid {}=[{:.2R}] must be greater than zero.",
413 : CurrentModuleObject,
414 0 : state.dataIPShortCut->cAlphaArgs(1),
415 : cNumericFields(7),
416 0 : rNumericArgs(7)));
417 : }
418 0 : ErrorsFound = true;
419 : }
420 :
421 3 : windTurbine.CutOutSpeed = state.dataIPShortCut->rNumericArgs(8); // Minimum wind speed for system operation
422 3 : if (windTurbine.CutOutSpeed == 0.0) {
423 0 : if (lNumericBlanks(8)) {
424 0 : ShowSevereError(state,
425 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cNumericFields(8) +
426 : " is required but input is blank.");
427 0 : } else if (windTurbine.CutOutSpeed <= windTurbine.RatedWindSpeed) {
428 0 : ShowSevereError(state,
429 0 : format("{}=\"{}\" invalid {}=[{:.2R}] must be greater than {}=[{:.2R}].",
430 : CurrentModuleObject,
431 0 : state.dataIPShortCut->cAlphaArgs(1),
432 : cNumericFields(8),
433 : rNumericArgs(8),
434 : cNumericFields(6),
435 0 : rNumericArgs(6)));
436 : } else {
437 0 : ShowSevereError(state,
438 0 : format("{}=\"{}\" invalid {}=[{:.2R}] must be greater than zero",
439 : CurrentModuleObject,
440 0 : state.dataIPShortCut->cAlphaArgs(1),
441 : cNumericFields(8),
442 0 : rNumericArgs(8)));
443 : }
444 0 : ErrorsFound = true;
445 : }
446 :
447 3 : windTurbine.SysEfficiency = state.dataIPShortCut->rNumericArgs(9); // Overall wind turbine system efficiency
448 3 : if (lNumericBlanks(9) || windTurbine.SysEfficiency == 0.0 || windTurbine.SysEfficiency > 1.0) {
449 0 : windTurbine.SysEfficiency = SysEffDefault;
450 0 : ShowWarningError(state,
451 0 : format("{}=\"{}\" invalid {}=[{:.2R}].",
452 : CurrentModuleObject,
453 0 : state.dataIPShortCut->cAlphaArgs(1),
454 : cNumericFields(9),
455 0 : state.dataIPShortCut->rNumericArgs(9)));
456 0 : ShowContinueError(state, format("...The default value of {:.3R} was assumed. for {}", SysEffDefault, cNumericFields(9)));
457 : }
458 :
459 3 : windTurbine.MaxTipSpeedRatio = state.dataIPShortCut->rNumericArgs(10); // Maximum tip speed ratio
460 3 : if (windTurbine.MaxTipSpeedRatio == 0.0) {
461 0 : if (lNumericBlanks(10)) {
462 0 : ShowSevereError(state,
463 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cNumericFields(10) +
464 : " is required but input is blank.");
465 : } else {
466 0 : ShowSevereError(state,
467 0 : format("{}=\"{}\" invalid {}=[{:.2R}] must be greater than zero.",
468 : CurrentModuleObject,
469 0 : state.dataIPShortCut->cAlphaArgs(1),
470 : cNumericFields(10),
471 0 : rNumericArgs(10)));
472 : }
473 0 : ErrorsFound = true;
474 : }
475 3 : if (windTurbine.SysEfficiency > MaxTSR) {
476 0 : windTurbine.SysEfficiency = MaxTSR;
477 0 : ShowWarningError(state,
478 0 : format("{}=\"{}\" invalid {}=[{:.2R}].",
479 : CurrentModuleObject,
480 0 : state.dataIPShortCut->cAlphaArgs(1),
481 : cNumericFields(10),
482 0 : state.dataIPShortCut->rNumericArgs(10)));
483 0 : ShowContinueError(state, format("...The default value of {:.1R} was assumed. for {}", MaxTSR, cNumericFields(10)));
484 : }
485 :
486 3 : windTurbine.MaxPowerCoeff = state.dataIPShortCut->rNumericArgs(11); // Maximum power coefficient
487 3 : if (windTurbine.rotorType == RotorType::HorizontalAxis && windTurbine.MaxPowerCoeff == 0.0) {
488 0 : if (lNumericBlanks(11)) {
489 0 : ShowSevereError(state,
490 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cNumericFields(11) +
491 : " is required but input is blank.");
492 : } else {
493 0 : ShowSevereError(state,
494 0 : format("{}=\"{}\" invalid {}=[{:.2R}] must be greater than zero.",
495 : CurrentModuleObject,
496 0 : state.dataIPShortCut->cAlphaArgs(1),
497 : cNumericFields(11),
498 0 : rNumericArgs(11)));
499 : }
500 0 : ErrorsFound = true;
501 : }
502 3 : if (windTurbine.MaxPowerCoeff > MaxPowerCoeff) {
503 0 : windTurbine.MaxPowerCoeff = DefaultPC;
504 0 : ShowWarningError(state,
505 0 : format("{}=\"{}\" invalid {}=[{:.2R}].",
506 : CurrentModuleObject,
507 0 : state.dataIPShortCut->cAlphaArgs(1),
508 : cNumericFields(11),
509 0 : state.dataIPShortCut->rNumericArgs(11)));
510 0 : ShowContinueError(state, format("...The default value of {:.2R} will be used. for {}", DefaultPC, cNumericFields(11)));
511 : }
512 :
513 3 : windTurbine.LocalAnnualAvgWS = state.dataIPShortCut->rNumericArgs(12); // Local wind speed annually averaged
514 3 : if (windTurbine.LocalAnnualAvgWS == 0.0) {
515 0 : if (lNumericBlanks(12)) {
516 0 : ShowWarningError(state,
517 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cNumericFields(12) +
518 : " is necessary for accurate prediction but input is blank.");
519 : } else {
520 0 : ShowSevereError(state,
521 0 : format("{}=\"{}\" invalid {}=[{:.2R}] must be greater than zero.",
522 : CurrentModuleObject,
523 0 : state.dataIPShortCut->cAlphaArgs(1),
524 : cNumericFields(12),
525 0 : rNumericArgs(12)));
526 0 : ErrorsFound = true;
527 : }
528 : }
529 :
530 3 : windTurbine.HeightForLocalWS = state.dataIPShortCut->rNumericArgs(13); // Height of local meteorological station
531 3 : if (windTurbine.HeightForLocalWS == 0.0) {
532 0 : if (windTurbine.LocalAnnualAvgWS == 0.0) {
533 0 : windTurbine.HeightForLocalWS = 0.0;
534 : } else {
535 0 : windTurbine.HeightForLocalWS = DefaultH;
536 0 : if (lNumericBlanks(13)) {
537 0 : ShowWarningError(state,
538 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cNumericFields(13) +
539 : " is necessary for accurate prediction but input is blank.");
540 0 : ShowContinueError(state, format("...The default value of {:.2R} will be used. for {}", DefaultH, cNumericFields(13)));
541 : } else {
542 0 : ShowSevereError(state,
543 0 : format("{}=\"{}\" invalid {}=[{:.2R}] must be greater than zero.",
544 : CurrentModuleObject,
545 0 : state.dataIPShortCut->cAlphaArgs(1),
546 : cNumericFields(13),
547 0 : rNumericArgs(13)));
548 0 : ErrorsFound = true;
549 : }
550 : }
551 : }
552 :
553 3 : windTurbine.ChordArea = state.dataIPShortCut->rNumericArgs(14); // Chord area of a single blade for VAWTs
554 3 : if (windTurbine.rotorType == RotorType::VerticalAxis && windTurbine.ChordArea == 0.0) {
555 0 : if (lNumericBlanks(14)) {
556 0 : ShowSevereError(state,
557 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cNumericFields(14) +
558 : " is required but input is blank.");
559 : } else {
560 0 : ShowSevereError(state,
561 0 : format("{}=\"{}\" invalid {}=[{:.2R}] must be greater than zero.",
562 : CurrentModuleObject,
563 0 : state.dataIPShortCut->cAlphaArgs(1),
564 : cNumericFields(14),
565 0 : rNumericArgs(14)));
566 : }
567 0 : ErrorsFound = true;
568 : }
569 :
570 3 : windTurbine.DragCoeff = state.dataIPShortCut->rNumericArgs(15); // Blade drag coefficient
571 3 : if (windTurbine.rotorType == RotorType::VerticalAxis && windTurbine.DragCoeff == 0.0) {
572 0 : if (lNumericBlanks(15)) {
573 0 : ShowSevereError(state,
574 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cNumericFields(15) +
575 : " is required but input is blank.");
576 : } else {
577 0 : ShowSevereError(state,
578 0 : format("{}=\"{}\" invalid {}=[{:.2R}] must be greater than zero.",
579 : CurrentModuleObject,
580 0 : state.dataIPShortCut->cAlphaArgs(1),
581 : cNumericFields(15),
582 0 : rNumericArgs(15)));
583 : }
584 0 : ErrorsFound = true;
585 : }
586 :
587 3 : windTurbine.LiftCoeff = state.dataIPShortCut->rNumericArgs(16); // Blade lift coefficient
588 3 : if (windTurbine.rotorType == RotorType::VerticalAxis && windTurbine.LiftCoeff == 0.0) {
589 0 : if (lNumericBlanks(16)) {
590 0 : ShowSevereError(state,
591 0 : CurrentModuleObject + "=\"" + state.dataIPShortCut->cAlphaArgs(1) + "\" invalid " + cNumericFields(16) +
592 : " is required but input is blank.");
593 : } else {
594 0 : ShowSevereError(state,
595 0 : format("{}=\"{}\" invalid {}=[{:.2R}] must be greater than zero.",
596 : CurrentModuleObject,
597 0 : state.dataIPShortCut->cAlphaArgs(1),
598 : cNumericFields(16),
599 0 : rNumericArgs(16)));
600 : }
601 0 : ErrorsFound = true;
602 : }
603 :
604 3 : windTurbine.PowerCoeffs[0] = state.dataIPShortCut->rNumericArgs(17); // Empirical power coefficient C1
605 3 : if (lNumericBlanks(17)) {
606 2 : windTurbine.PowerCoeffs[0] = 0.0;
607 : }
608 3 : windTurbine.PowerCoeffs[1] = state.dataIPShortCut->rNumericArgs(18); // Empirical power coefficient C2
609 3 : if (lNumericBlanks(18)) {
610 2 : windTurbine.PowerCoeffs[1] = 0.0;
611 : }
612 3 : windTurbine.PowerCoeffs[2] = state.dataIPShortCut->rNumericArgs(19); // Empirical power coefficient C3
613 3 : if (lNumericBlanks(19)) {
614 2 : windTurbine.PowerCoeffs[2] = 0.0;
615 : }
616 3 : windTurbine.PowerCoeffs[3] = state.dataIPShortCut->rNumericArgs(20); // Empirical power coefficient C4
617 3 : if (lNumericBlanks(20)) {
618 2 : windTurbine.PowerCoeffs[3] = 0.0;
619 : }
620 3 : windTurbine.PowerCoeffs[4] = state.dataIPShortCut->rNumericArgs(21); // Empirical power coefficient C5
621 3 : if (lNumericBlanks(21)) {
622 2 : windTurbine.PowerCoeffs[4] = 0.0;
623 : }
624 3 : windTurbine.PowerCoeffs[5] = state.dataIPShortCut->rNumericArgs(22); // Empirical power coefficient C6
625 3 : if (lNumericBlanks(22)) {
626 2 : windTurbine.PowerCoeffs[5] = 0.0;
627 : }
628 : }
629 :
630 1 : cAlphaArgs.deallocate();
631 1 : cAlphaFields.deallocate();
632 1 : cNumericFields.deallocate();
633 1 : rNumericArgs.deallocate();
634 1 : lAlphaBlanks.deallocate();
635 1 : lNumericBlanks.deallocate();
636 :
637 1 : if (ErrorsFound) ShowFatalError(state, CurrentModuleObject + " errors occurred in input. Program terminates.");
638 :
639 4 : for (WindTurbineNum = 1; WindTurbineNum <= NumWindTurbines; ++WindTurbineNum) {
640 3 : auto &windTurbine = state.dataWindTurbine->WindTurbineSys(WindTurbineNum);
641 6 : SetupOutputVariable(state,
642 : "Generator Produced AC Electricity Rate",
643 : OutputProcessor::Unit::W,
644 : windTurbine.Power,
645 : OutputProcessor::SOVTimeStepType::System,
646 : OutputProcessor::SOVStoreType::Average,
647 3 : windTurbine.Name);
648 6 : SetupOutputVariable(state,
649 : "Generator Produced AC Electricity Energy",
650 : OutputProcessor::Unit::J,
651 : windTurbine.Energy,
652 : OutputProcessor::SOVTimeStepType::System,
653 : OutputProcessor::SOVStoreType::Summed,
654 : windTurbine.Name,
655 : _,
656 : "ElectricityProduced",
657 : "WINDTURBINE",
658 : _,
659 3 : "Plant");
660 6 : SetupOutputVariable(state,
661 : "Generator Turbine Local Wind Speed",
662 : OutputProcessor::Unit::m_s,
663 : windTurbine.LocalWindSpeed,
664 : OutputProcessor::SOVTimeStepType::System,
665 : OutputProcessor::SOVStoreType::Average,
666 3 : windTurbine.Name);
667 6 : SetupOutputVariable(state,
668 : "Generator Turbine Local Air Density",
669 : OutputProcessor::Unit::kg_m3,
670 : windTurbine.LocalAirDensity,
671 : OutputProcessor::SOVTimeStepType::System,
672 : OutputProcessor::SOVStoreType::Average,
673 3 : windTurbine.Name);
674 6 : SetupOutputVariable(state,
675 : "Generator Turbine Tip Speed Ratio",
676 : OutputProcessor::Unit::None,
677 : windTurbine.TipSpeedRatio,
678 : OutputProcessor::SOVTimeStepType::System,
679 : OutputProcessor::SOVStoreType::Average,
680 3 : windTurbine.Name);
681 3 : switch (windTurbine.rotorType) {
682 2 : case RotorType::HorizontalAxis: {
683 4 : SetupOutputVariable(state,
684 : "Generator Turbine Power Coefficient",
685 : OutputProcessor::Unit::None,
686 : windTurbine.PowerCoeff,
687 : OutputProcessor::SOVTimeStepType::System,
688 : OutputProcessor::SOVStoreType::Average,
689 2 : windTurbine.Name);
690 2 : } break;
691 1 : case RotorType::VerticalAxis: {
692 2 : SetupOutputVariable(state,
693 : "Generator Turbine Chordal Component Velocity",
694 : OutputProcessor::Unit::m_s,
695 : windTurbine.ChordalVel,
696 : OutputProcessor::SOVTimeStepType::System,
697 : OutputProcessor::SOVStoreType::Average,
698 1 : windTurbine.Name);
699 2 : SetupOutputVariable(state,
700 : "Generator Turbine Normal Component Velocity",
701 : OutputProcessor::Unit::m_s,
702 : windTurbine.NormalVel,
703 : OutputProcessor::SOVTimeStepType::System,
704 : OutputProcessor::SOVStoreType::Average,
705 1 : windTurbine.Name);
706 2 : SetupOutputVariable(state,
707 : "Generator Turbine Relative Flow Velocity",
708 : OutputProcessor::Unit::m_s,
709 : windTurbine.RelFlowVel,
710 : OutputProcessor::SOVTimeStepType::System,
711 : OutputProcessor::SOVStoreType::Average,
712 1 : windTurbine.Name);
713 2 : SetupOutputVariable(state,
714 : "Generator Turbine Attack Angle",
715 : OutputProcessor::Unit::deg,
716 : windTurbine.AngOfAttack,
717 : OutputProcessor::SOVTimeStepType::System,
718 : OutputProcessor::SOVStoreType::Average,
719 1 : windTurbine.Name);
720 1 : } break;
721 0 : default:
722 0 : break;
723 : }
724 : }
725 1 : }
726 :
727 16344 : void InitWindTurbine(EnergyPlusData &state, int const WindTurbineNum)
728 : {
729 :
730 : // SUBROUTINE INFORMATION:
731 : // AUTHOR Daeho Kang
732 : // DATE WRITTEN Oct 2009
733 : // MODIFIED Linda K. Lawrie, December 2009 for reading stat file
734 : // RE-ENGINEERED na
735 :
736 : // PURPOSE OF THIS SUBROUTINE:
737 : // This subroutine reads monthly average wind speed from stat file and then
738 : // determines annual average wind speed. Differences between this TMY wind speed
739 : // and local wind speed that the user inputs are then factored.
740 : // IF the user has no local wind data and does not enter the local wind speed to be factored,
741 : // then the factor of 1 is assigned, so that wind speed estimated
742 : // at the particular rotor height is used with no factorization.
743 : // It also initializes module variables at each time step.
744 :
745 : static char const TabChr('\t'); // Tab character
746 :
747 : int mon; // loop counter
748 : bool wsStatFound; // logical noting that wind stats were found
749 : bool warningShown; // true if the <365 warning has already been shown
750 32688 : Array1D<Real64> MonthWS(12);
751 : Real64 LocalTMYWS; // Annual average wind speed at the rotor height
752 :
753 : // Estimate average annual wind speed once
754 16344 : if (state.dataWindTurbine->MyOneTimeFlag) {
755 1 : wsStatFound = false;
756 1 : Real64 AnnualTMYWS = 0.0;
757 1 : if (FileSystem::fileExists(state.files.inStatFilePath.filePath)) {
758 2 : auto statFile = state.files.inStatFilePath.open(state, "InitWindTurbine");
759 509 : while (statFile.good()) { // end of file
760 255 : auto lineIn = statFile.readLine();
761 : // reconcile line with different versions of stat file
762 255 : auto lnPtr = index(lineIn.data, "Wind Speed");
763 255 : if (lnPtr == std::string::npos) continue;
764 : // have hit correct section.
765 15 : while (statFile.good()) { // find daily avg line
766 8 : lineIn = statFile.readLine();
767 8 : lnPtr = index(lineIn.data, "Daily Avg");
768 8 : if (lnPtr == std::string::npos) continue;
769 : // tab delimited file
770 1 : lineIn.data.erase(0, lnPtr + 10);
771 1 : MonthWS = 0.0;
772 1 : wsStatFound = true;
773 1 : warningShown = false;
774 13 : for (mon = 1; mon <= 12; ++mon) {
775 12 : lnPtr = index(lineIn.data, TabChr);
776 12 : if (lnPtr != 1) {
777 12 : if ((lnPtr == std::string::npos) || (!stripped(lineIn.data.substr(0, lnPtr)).empty())) {
778 12 : if (lnPtr != std::string::npos) {
779 12 : bool error = false;
780 12 : MonthWS(mon) = UtilityRoutines::ProcessNumber(lineIn.data.substr(0, lnPtr), error);
781 :
782 12 : if (error) {
783 : // probably should throw some error here
784 : }
785 :
786 12 : lineIn.data.erase(0, lnPtr + 1);
787 : }
788 : } else { // blank field
789 0 : if (!warningShown) {
790 0 : ShowWarningError(
791 : state,
792 0 : "InitWindTurbine: read from " + state.files.inStatFilePath.filePath.string() +
793 : " file shows <365 days in weather file. Annual average wind speed used will be inaccurate.");
794 0 : lineIn.data.erase(0, lnPtr + 1);
795 0 : warningShown = true;
796 : }
797 : }
798 : } else { // two tabs in succession
799 0 : if (!warningShown) {
800 0 : ShowWarningError(state,
801 0 : "InitWindTurbine: read from " + state.files.inStatFilePath.filePath.string() +
802 : " file shows <365 days in weather file. Annual average wind speed used will be inaccurate.");
803 0 : lineIn.data.erase(0, lnPtr + 1);
804 0 : warningShown = true;
805 : }
806 : }
807 : }
808 1 : break;
809 : }
810 1 : if (wsStatFound) break;
811 : }
812 1 : if (wsStatFound) {
813 1 : AnnualTMYWS = sum(MonthWS) / 12.0;
814 : } else {
815 0 : ShowWarningError(
816 : state, "InitWindTurbine: stat file did not include Wind Speed statistics. TMY Wind Speed adjusted at the height is used.");
817 : }
818 : } else { // No stat file
819 0 : ShowWarningError(state, "InitWindTurbine: stat file missing. TMY Wind Speed adjusted at the height is used.");
820 : }
821 :
822 : // assign this value to all the wind turbines once here
823 4 : for (auto &wt : state.dataWindTurbine->WindTurbineSys) {
824 3 : wt.AnnualTMYWS = AnnualTMYWS;
825 : }
826 :
827 1 : state.dataWindTurbine->MyOneTimeFlag = false;
828 : }
829 :
830 16344 : auto &windTurbine = state.dataWindTurbine->WindTurbineSys(WindTurbineNum);
831 :
832 : // Factor differences between TMY wind data and local wind data once
833 16344 : if (windTurbine.AnnualTMYWS > 0.0 && windTurbine.WSFactor == 0.0 && windTurbine.LocalAnnualAvgWS > 0) {
834 : // Convert the annual wind speed to the local wind speed at the height of the local station, then factor
835 6 : LocalTMYWS = windTurbine.AnnualTMYWS * state.dataEnvrn->WeatherFileWindModCoeff *
836 3 : std::pow(windTurbine.HeightForLocalWS / state.dataEnvrn->SiteWindBLHeight, state.dataEnvrn->SiteWindExp);
837 3 : windTurbine.WSFactor = LocalTMYWS / windTurbine.LocalAnnualAvgWS;
838 : }
839 : // Assign factor of 1.0 if no stat file or no input of local average wind speed
840 16344 : if (windTurbine.WSFactor == 0.0) windTurbine.WSFactor = 1.0;
841 :
842 : // Do every time step initialization
843 16344 : windTurbine.Power = 0.0;
844 16344 : windTurbine.TotPower = 0.0;
845 16344 : windTurbine.PowerCoeff = 0.0;
846 16344 : windTurbine.TipSpeedRatio = 0.0;
847 16344 : windTurbine.ChordalVel = 0.0;
848 16344 : windTurbine.NormalVel = 0.0;
849 16344 : windTurbine.RelFlowVel = 0.0;
850 16344 : windTurbine.AngOfAttack = 0.0;
851 16344 : windTurbine.TanForce = 0.0;
852 16344 : windTurbine.NorForce = 0.0;
853 16344 : windTurbine.TotTorque = 0.0;
854 16344 : }
855 :
856 16344 : void CalcWindTurbine(EnergyPlusData &state,
857 : int const WindTurbineNum, // System is on
858 : [[maybe_unused]] bool const RunFlag // System is on
859 : )
860 : {
861 :
862 : // SUBROUTINE INFORMATION:
863 : // AUTHOR Daeho Kang
864 : // DATE WRITTEN Octorber 2009
865 : // MODIFIED na
866 : // RE-ENGINEERED na
867 :
868 : // REFERENCES:
869 : // Sathyajith Mathew. 2006. Wind Energy: Fundamental, Resource Analysis and Economics. Springer,
870 : // Chap. 2, pp. 11-15
871 : // Mazharul Islam, David S.K. Ting, and Amir Fartaj. 2008. Aerodynamic Models for Darrieus-type sSraight-bladed
872 : // Vertical Axis Wind Turbines. Renewable & Sustainable Energy Reviews, Volume 12, pp.1087-1109
873 :
874 : using DataEnvironment::OutBaroPressAt;
875 : using DataEnvironment::OutDryBulbTempAt;
876 : using DataEnvironment::OutWetBulbTempAt;
877 : using Psychrometrics::PsyRhoAirFnPbTdbW;
878 : using Psychrometrics::PsyWFnTdbTwbPb;
879 : using ScheduleManager::GetCurrentScheduleValue;
880 :
881 16344 : Real64 constexpr MaxTheta(90.0); // Maximum of theta
882 16344 : Real64 constexpr MaxDegree(360.0); // Maximum limit of outdoor air wind speed in m/s
883 16344 : Real64 constexpr SecInMin(60.0);
884 :
885 : Real64 LocalWindSpeed; // Ambient wind speed at the specific height in m/s
886 : Real64 RotorH; // Height of the rotor in m
887 : Real64 RotorD; // Diameter of the rotor in m
888 : Real64 LocalHumRat; // Local humidity ratio of the air at the specific height
889 : Real64 LocalAirDensity; // Local density at the specific height in m
890 : Real64 PowerCoeff; // Power coefficient
891 : Real64 SweptArea; // Swept area of the rotor in m2
892 16344 : Real64 WTPower(0.0); // Total Power generated by the turbine in the quasi-steady state in Watts
893 : Real64 Power; // Actual power of the turbine in Watts
894 : Real64 TipSpeedRatio; // Tip speed ratio (TSR)
895 : Real64 TipSpeedRatioAtI; // Tip speed ratio at the ith time step
896 : Real64 AzimuthAng; // Azimuth angle of blades in degree
897 : Real64 ChordalVel; // Chordal velocity component in m/s
898 : Real64 NormalVel; // Normal velocity component in m/s
899 : Real64 AngOfAttack; // Angle of attack of a single blade in degree
900 : Real64 RelFlowVel; // Relative flow velocity component in m/s
901 : Real64 TanForce; // Tangential force in N.m
902 : Real64 NorForce; // Normal force in N.m
903 : Real64 RotorVel; // Rotor velocity in m/s
904 : Real64 AvgTanForce; // Average tangential force in N.m
905 : Real64 Constant; // Constants within integrand of tangential force
906 : Real64 IntRelFlowVel; // Integration of relative flow velocity
907 : Real64 TotTorque; // Total torque for the number of blades
908 : Real64 Omega; // Angular velocity of rotor in rad/s
909 : Real64 TanForceCoeff; // Tnagential force coefficient
910 : Real64 NorForceCoeff; // Normal force coefficient
911 : Real64 Period; // Period of sine and cosine functions
912 : Real64 C1; // Empirical power coefficient C1
913 : Real64 C2; // Empirical power coefficient C2
914 : Real64 C3; // Empirical power coefficient C3
915 : Real64 C4; // Empirical power coefficient C4
916 : Real64 C5; // Empirical power coefficient C5
917 : Real64 C6; // Empirical power coefficient C6
918 : Real64 LocalTemp; // Ambient air temperature at the height in degree C
919 : Real64 LocalPress; // Local atmospheric pressure at the particular height in Pa
920 : Real64 InducedVel; // Induced velocity on the rotor in m/s
921 : // unused REAL(r64) :: SysEfficiency ! Overall wind turbine system efficiency including generator and inverter
922 : Real64 MaxPowerCoeff; // Maximum power coefficient
923 : Real64 RotorSpeed; // Speed of rotors
924 :
925 16344 : auto &windTurbine = state.dataWindTurbine->WindTurbineSys(WindTurbineNum);
926 : // Estimate local velocity and density
927 16344 : RotorH = windTurbine.RotorHeight;
928 16344 : RotorD = windTurbine.RotorDiameter;
929 16344 : RotorSpeed = windTurbine.RatedRotorSpeed;
930 16344 : LocalTemp = OutDryBulbTempAt(state, RotorH);
931 16344 : LocalPress = OutBaroPressAt(state, RotorH);
932 16344 : LocalHumRat = PsyWFnTdbTwbPb(state, LocalTemp, OutWetBulbTempAt(state, RotorH), LocalPress);
933 16344 : LocalAirDensity = PsyRhoAirFnPbTdbW(state, LocalPress, LocalTemp, LocalHumRat);
934 16344 : LocalWindSpeed = DataEnvironment::WindSpeedAt(state, RotorH);
935 16344 : LocalWindSpeed /= windTurbine.WSFactor;
936 :
937 : // Check wind conditions for system operation
938 32661 : if (GetCurrentScheduleValue(state, windTurbine.SchedPtr) > 0 && LocalWindSpeed > windTurbine.CutInSpeed &&
939 16317 : LocalWindSpeed < windTurbine.CutOutSpeed) {
940 :
941 : // System is on
942 16317 : Period = 2.0 * DataGlobalConstants::Pi;
943 16317 : Omega = (RotorSpeed * Period) / SecInMin;
944 16317 : SweptArea = (DataGlobalConstants::Pi * pow_2(RotorD)) / 4;
945 16317 : TipSpeedRatio = (Omega * (RotorD / 2.0)) / LocalWindSpeed;
946 :
947 : // Limit maximum tip speed ratio
948 16317 : if (TipSpeedRatio > windTurbine.MaxTipSpeedRatio) {
949 10887 : TipSpeedRatio = windTurbine.MaxTipSpeedRatio;
950 : }
951 :
952 16317 : switch (windTurbine.rotorType) {
953 10878 : case RotorType::HorizontalAxis: { // Horizontal axis wind turbine
954 10878 : MaxPowerCoeff = windTurbine.MaxPowerCoeff;
955 : // Check if empirical constants are available
956 10878 : C1 = windTurbine.PowerCoeffs[0];
957 10878 : C2 = windTurbine.PowerCoeffs[1];
958 10878 : C3 = windTurbine.PowerCoeffs[2];
959 10878 : C4 = windTurbine.PowerCoeffs[3];
960 10878 : C5 = windTurbine.PowerCoeffs[4];
961 10878 : C6 = windTurbine.PowerCoeffs[5];
962 :
963 10878 : Real64 const LocalWindSpeed_3(pow_3(LocalWindSpeed));
964 10878 : if (C1 > 0.0 && C2 > 0.0 && C3 > 0.0 && C4 >= 0.0 && C5 > 0.0 && C6 > 0.0) {
965 : // Analytical approximation
966 : // Maximum power, i.e., rotor speed is at maximum, and pitch angle is zero
967 : // TipSpeedRatioAtI = 1.0 / ( ( 1.0 / ( TipSpeedRatio + 0.08 * PitchAngle ) ) - ( 0.035 / ( pow_3( PitchAngle ) + 1.0 ) ) );
968 : // //Tuned PitchAngle is zero
969 5439 : TipSpeedRatioAtI = TipSpeedRatio / (1.0 - (TipSpeedRatio * 0.035));
970 : // PowerCoeff = C1 * ( ( C2 / TipSpeedRatioAtI ) - ( C3 * PitchAngle ) - ( C4 * std::pow( PitchAngle, 1.5 ) ) - C5 ) * (
971 : // std::exp( -( C6 / TipSpeedRatioAtI ) ) ); //Tuned PitchAngle is zero
972 5439 : PowerCoeff = C1 * ((C2 / TipSpeedRatioAtI) - C5) * std::exp(-(C6 / TipSpeedRatioAtI));
973 5439 : if (PowerCoeff > MaxPowerCoeff) {
974 0 : PowerCoeff = MaxPowerCoeff;
975 : }
976 5439 : WTPower = 0.5 * LocalAirDensity * PowerCoeff * SweptArea * LocalWindSpeed_3;
977 : } else { // Simple approximation
978 5439 : WTPower = 0.5 * LocalAirDensity * SweptArea * LocalWindSpeed_3 * MaxPowerCoeff;
979 5439 : PowerCoeff = MaxPowerCoeff;
980 : }
981 : // Maximum of rated power
982 10878 : if (LocalWindSpeed >= windTurbine.RatedWindSpeed || WTPower > windTurbine.RatedPower) {
983 0 : WTPower = windTurbine.RatedPower;
984 0 : PowerCoeff = WTPower / (0.5 * LocalAirDensity * SweptArea * LocalWindSpeed_3);
985 : }
986 : // Recalculated Cp at the rated power
987 10878 : windTurbine.PowerCoeff = PowerCoeff;
988 10878 : } break;
989 5439 : case RotorType::VerticalAxis: { // Vertical axis wind turbine
990 5439 : RotorVel = Omega * (RotorD / 2.0);
991 : // Recalculated omega, if TSR is greater than the maximum
992 5439 : if (TipSpeedRatio >= windTurbine.MaxTipSpeedRatio) {
993 5439 : RotorVel = LocalWindSpeed * windTurbine.MaxTipSpeedRatio;
994 5439 : Omega = RotorVel / (RotorD / 2.0);
995 : }
996 :
997 5439 : AzimuthAng = MaxDegree / windTurbine.NumOfBlade;
998 : // Azimuth angle between zero and 90 degree
999 5439 : if (AzimuthAng > MaxTheta) { // Number of blades is 2 or 3
1000 5439 : AzimuthAng -= MaxTheta;
1001 5439 : if (AzimuthAng == MaxTheta) { // 2 blades
1002 0 : AzimuthAng = 0.0;
1003 : }
1004 0 : } else if (AzimuthAng == MaxTheta) { // 4 blades
1005 0 : AzimuthAng = 0.0;
1006 : }
1007 :
1008 5439 : InducedVel = LocalWindSpeed * 2.0 / 3.0;
1009 : // Velocity components
1010 5439 : Real64 const sin_AzimuthAng(std::sin(AzimuthAng * DataGlobalConstants::DegToRadians));
1011 5439 : Real64 const cos_AzimuthAng(std::cos(AzimuthAng * DataGlobalConstants::DegToRadians));
1012 5439 : ChordalVel = RotorVel + InducedVel * cos_AzimuthAng;
1013 5439 : NormalVel = InducedVel * sin_AzimuthAng;
1014 5439 : RelFlowVel = std::sqrt(pow_2(ChordalVel) + pow_2(NormalVel));
1015 :
1016 : // Angle of attack
1017 5439 : AngOfAttack = std::atan((sin_AzimuthAng / ((RotorVel / LocalWindSpeed) / (InducedVel / LocalWindSpeed) + cos_AzimuthAng)));
1018 :
1019 : // Force coefficients
1020 5439 : Real64 const sin_AngOfAttack(std::sin(AngOfAttack * DataGlobalConstants::DegToRadians));
1021 5439 : Real64 const cos_AngOfAttack(std::cos(AngOfAttack * DataGlobalConstants::DegToRadians));
1022 5439 : TanForceCoeff = std::abs(windTurbine.LiftCoeff * sin_AngOfAttack - windTurbine.DragCoeff * cos_AngOfAttack);
1023 5439 : NorForceCoeff = windTurbine.LiftCoeff * cos_AngOfAttack + windTurbine.DragCoeff * sin_AngOfAttack;
1024 :
1025 : // Net tangential and normal forces
1026 5439 : Real64 const RelFlowVel_2(pow_2(RelFlowVel));
1027 5439 : Real64 const density_fac(0.5 * LocalAirDensity * windTurbine.ChordArea * RelFlowVel_2);
1028 5439 : TanForce = TanForceCoeff * density_fac;
1029 5439 : NorForce = NorForceCoeff * density_fac;
1030 5439 : Constant = (1.0 / Period) * (TanForce / RelFlowVel_2);
1031 :
1032 : // Relative flow velocity is the only function of theta in net tangential force
1033 : // Integral of cos(theta) on zero to 2pi goes to zero
1034 : // Integrate constants only
1035 5439 : IntRelFlowVel = pow_2(RotorVel) * Period + pow_2(InducedVel) * Period;
1036 :
1037 : // Average tangential force on a single blade
1038 5439 : AvgTanForce = Constant * IntRelFlowVel;
1039 5439 : TotTorque = windTurbine.NumOfBlade * AvgTanForce * (RotorD / 2.0);
1040 5439 : WTPower = TotTorque * Omega;
1041 :
1042 : // Check if power produced is greater than maximum or rated power
1043 5439 : if (WTPower > windTurbine.RatedPower) {
1044 5439 : WTPower = windTurbine.RatedPower;
1045 : }
1046 :
1047 5439 : windTurbine.ChordalVel = ChordalVel;
1048 5439 : windTurbine.NormalVel = NormalVel;
1049 5439 : windTurbine.RelFlowVel = RelFlowVel;
1050 5439 : windTurbine.TanForce = TanForce;
1051 5439 : windTurbine.NorForce = NorForce;
1052 5439 : windTurbine.TotTorque = TotTorque;
1053 5439 : } break;
1054 0 : default: {
1055 0 : assert(false);
1056 : } break;
1057 : }
1058 :
1059 16317 : if (WTPower > windTurbine.RatedPower) {
1060 0 : WTPower = windTurbine.RatedPower;
1061 : }
1062 :
1063 : // Actual power generated by the wind turbine system
1064 16317 : Power = WTPower * windTurbine.SysEfficiency;
1065 :
1066 16317 : windTurbine.Power = Power;
1067 16317 : windTurbine.TotPower = WTPower;
1068 16317 : windTurbine.LocalWindSpeed = LocalWindSpeed;
1069 16317 : windTurbine.LocalAirDensity = LocalAirDensity;
1070 16317 : windTurbine.TipSpeedRatio = TipSpeedRatio;
1071 :
1072 : } else { // System is off
1073 27 : windTurbine.Power = 0.0;
1074 27 : windTurbine.TotPower = 0.0;
1075 27 : windTurbine.PowerCoeff = 0.0;
1076 27 : windTurbine.LocalWindSpeed = LocalWindSpeed;
1077 27 : windTurbine.LocalAirDensity = LocalAirDensity;
1078 27 : windTurbine.TipSpeedRatio = 0.0;
1079 27 : windTurbine.ChordalVel = 0.0;
1080 27 : windTurbine.NormalVel = 0.0;
1081 27 : windTurbine.RelFlowVel = 0.0;
1082 27 : windTurbine.AngOfAttack = 0.0;
1083 27 : windTurbine.TanForce = 0.0;
1084 27 : windTurbine.NorForce = 0.0;
1085 27 : windTurbine.TotTorque = 0.0;
1086 : }
1087 16344 : }
1088 :
1089 16344 : void ReportWindTurbine(EnergyPlusData &state, int const WindTurbineNum)
1090 : {
1091 : // SUBROUTINE INFORMATION:
1092 : // AUTHOR Daeho Kang
1093 : // DATE WRITTEN October 2009
1094 : // MODIFIED na
1095 : // RE-ENGINEERED na
1096 :
1097 : // PURPOSE OF THIS SUBROUTINE:
1098 : // This subroutine fills remaining report variables.
1099 :
1100 16344 : auto &TimeStepSys = state.dataHVACGlobal->TimeStepSys;
1101 16344 : auto &windTurbine = state.dataWindTurbine->WindTurbineSys(WindTurbineNum);
1102 :
1103 16344 : windTurbine.Energy = windTurbine.Power * TimeStepSys * DataGlobalConstants::SecInHour;
1104 16344 : }
1105 :
1106 : //*****************************************************************************************
1107 :
1108 : } // namespace WindTurbine
1109 :
1110 2313 : } // namespace EnergyPlus
|
