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