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 16236 : 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 16236 : if (state.dataWindTurbine->GetInputFlag) {
129 1 : GetWindTurbineInput(state);
130 1 : state.dataWindTurbine->GetInputFlag = false;
131 : }
132 :
133 16236 : if (GeneratorIndex == 0) {
134 3 : WindTurbineNum = Util::FindItemInList(GeneratorName, state.dataWindTurbine->WindTurbineSys);
135 3 : if (WindTurbineNum == 0) {
136 0 : ShowFatalError(state, format("SimWindTurbine: Specified Generator not one of Valid Wind Turbine Generators {}", GeneratorName));
137 : }
138 3 : GeneratorIndex = WindTurbineNum;
139 : } else {
140 16233 : WindTurbineNum = GeneratorIndex;
141 16233 : int NumWindTurbines = (int)state.dataWindTurbine->WindTurbineSys.size();
142 16233 : 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 16233 : 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 16236 : InitWindTurbine(state, WindTurbineNum);
159 :
160 16236 : CalcWindTurbine(state, WindTurbineNum, RunFlag);
161 :
162 16236 : ReportWindTurbine(state, WindTurbineNum);
163 16236 : }
164 :
165 16236 : 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 16236 : GeneratorPower = state.dataWindTurbine->WindTurbineSys(GeneratorIndex).Power;
183 16236 : GeneratorEnergy = state.dataWindTurbine->WindTurbineSys(GeneratorIndex).Energy;
184 :
185 : // Thermal energy is ignored
186 16236 : ThermalPower = 0.0;
187 16236 : ThermalEnergy = 0.0;
188 16236 : }
189 :
190 1 : 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 3 : static std::string const CurrentModuleObject("Generator:WindTurbine");
205 1 : Real64 constexpr SysEffDefault(0.835); // Default value of overall system efficiency
206 1 : Real64 constexpr MaxTSR(12.0); // Maximum tip speed ratio
207 1 : Real64 constexpr DefaultPC(0.25); // Default power coefficient
208 1 : Real64 constexpr MaxPowerCoeff(0.59); // Maximum power coefficient
209 1 : Real64 constexpr DefaultH(50.0); // Default of height for local wind speed
210 :
211 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
212 1 : 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 1 : Array1D_string cAlphaArgs; // Alpha input items for object
219 1 : Array1D_string cAlphaFields; // Alpha field names
220 1 : Array1D_string cNumericFields; // Numeric field names
221 1 : Array1D<Real64> rNumericArgs; // Numeric input items for object
222 1 : Array1D_bool lAlphaBlanks; // Logical array, alpha field input BLANK = .TRUE.
223 1 : Array1D_bool lNumericBlanks; // Logical array, numeric field input BLANK = .TRUE.
224 :
225 : // Initializations and allocations
226 1 : state.dataInputProcessing->inputProcessor->getObjectDefMaxArgs(state, CurrentModuleObject, NumArgs, NumAlphas, NumNumbers);
227 1 : cAlphaArgs.allocate(NumAlphas);
228 1 : cAlphaFields.allocate(NumAlphas);
229 1 : cNumericFields.allocate(NumNumbers);
230 1 : rNumericArgs.dimension(NumNumbers, 0.0);
231 1 : lAlphaBlanks.dimension(NumAlphas, true);
232 1 : lNumericBlanks.dimension(NumNumbers, true);
233 :
234 1 : int NumWindTurbines = state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, CurrentModuleObject);
235 :
236 1 : state.dataWindTurbine->WindTurbineSys.allocate(NumWindTurbines);
237 :
238 4 : for (WindTurbineNum = 1; WindTurbineNum <= NumWindTurbines; ++WindTurbineNum) {
239 :
240 6 : state.dataInputProcessing->inputProcessor->getObjectItem(state,
241 : CurrentModuleObject,
242 : WindTurbineNum,
243 3 : state.dataIPShortCut->cAlphaArgs,
244 : NumAlphas,
245 3 : state.dataIPShortCut->rNumericArgs,
246 : NumNumbers,
247 : IOStat,
248 : lNumericBlanks,
249 : lAlphaBlanks,
250 : cAlphaFields,
251 : cNumericFields);
252 :
253 3 : ErrorObjectHeader eoh{routineName, CurrentModuleObject, state.dataIPShortCut->cAlphaArgs(1)};
254 :
255 3 : Util::IsNameEmpty(state, state.dataIPShortCut->cAlphaArgs(1), CurrentModuleObject, ErrorsFound);
256 :
257 3 : auto &windTurbine = state.dataWindTurbine->WindTurbineSys(WindTurbineNum);
258 :
259 3 : windTurbine.Name = state.dataIPShortCut->cAlphaArgs(1); // Name of wind turbine
260 :
261 3 : if (lAlphaBlanks(2)) {
262 0 : windTurbine.availSched = Sched::GetScheduleAlwaysOn(state);
263 3 : } 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 3 : windTurbine.rotorType =
269 3 : static_cast<RotorType>(getEnumValue(WindTurbine::RotorNamesUC, Util::makeUPPER(state.dataIPShortCut->cAlphaArgs(3))));
270 3 : 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 3 : windTurbine.controlType =
286 3 : static_cast<ControlType>(getEnumValue(WindTurbine::ControlNamesUC, Util::makeUPPER(state.dataIPShortCut->cAlphaArgs(4))));
287 3 : 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 3 : windTurbine.RatedRotorSpeed = state.dataIPShortCut->rNumericArgs(1); // Maximum rotor speed in rpm
302 3 : 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 3 : windTurbine.RotorDiameter = state.dataIPShortCut->rNumericArgs(2); // Rotor diameter in m
321 3 : 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 3 : windTurbine.RotorHeight = state.dataIPShortCut->rNumericArgs(3); // Overall height of the rotor
340 3 : 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 3 : windTurbine.NumOfBlade = state.dataIPShortCut->rNumericArgs(4); // Total number of blade
359 3 : 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 3 : windTurbine.RatedPower = state.dataIPShortCut->rNumericArgs(5); // Rated average power
370 3 : 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 3 : windTurbine.RatedWindSpeed = state.dataIPShortCut->rNumericArgs(6); // Rated wind speed
389 3 : 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 3 : windTurbine.CutInSpeed = state.dataIPShortCut->rNumericArgs(7); // Minimum wind speed for system operation
408 3 : 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 3 : windTurbine.CutOutSpeed = state.dataIPShortCut->rNumericArgs(8); // Minimum wind speed for system operation
427 3 : 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 3 : windTurbine.SysEfficiency = state.dataIPShortCut->rNumericArgs(9); // Overall wind turbine system efficiency
455 3 : 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 3 : windTurbine.MaxTipSpeedRatio = state.dataIPShortCut->rNumericArgs(10); // Maximum tip speed ratio
467 3 : 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 3 : 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 3 : windTurbine.MaxPowerCoeff = state.dataIPShortCut->rNumericArgs(11); // Maximum power coefficient
496 3 : 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 3 : 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 3 : windTurbine.LocalAnnualAvgWS = state.dataIPShortCut->rNumericArgs(12); // Local wind speed annually averaged
525 3 : 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 3 : windTurbine.HeightForLocalWS = state.dataIPShortCut->rNumericArgs(13); // Height of local meteorological station
544 3 : 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 3 : windTurbine.ChordArea = state.dataIPShortCut->rNumericArgs(14); // Chord area of a single blade for VAWTs
569 3 : 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 3 : windTurbine.DragCoeff = state.dataIPShortCut->rNumericArgs(15); // Blade drag coefficient
588 3 : 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 3 : windTurbine.LiftCoeff = state.dataIPShortCut->rNumericArgs(16); // Blade lift coefficient
607 3 : 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 3 : windTurbine.PowerCoeffs[0] = state.dataIPShortCut->rNumericArgs(17); // Empirical power coefficient C1
626 3 : if (lNumericBlanks(17)) {
627 2 : windTurbine.PowerCoeffs[0] = 0.0;
628 : }
629 3 : windTurbine.PowerCoeffs[1] = state.dataIPShortCut->rNumericArgs(18); // Empirical power coefficient C2
630 3 : if (lNumericBlanks(18)) {
631 2 : windTurbine.PowerCoeffs[1] = 0.0;
632 : }
633 3 : windTurbine.PowerCoeffs[2] = state.dataIPShortCut->rNumericArgs(19); // Empirical power coefficient C3
634 3 : if (lNumericBlanks(19)) {
635 2 : windTurbine.PowerCoeffs[2] = 0.0;
636 : }
637 3 : windTurbine.PowerCoeffs[3] = state.dataIPShortCut->rNumericArgs(20); // Empirical power coefficient C4
638 3 : if (lNumericBlanks(20)) {
639 2 : windTurbine.PowerCoeffs[3] = 0.0;
640 : }
641 3 : windTurbine.PowerCoeffs[4] = state.dataIPShortCut->rNumericArgs(21); // Empirical power coefficient C5
642 3 : if (lNumericBlanks(21)) {
643 2 : windTurbine.PowerCoeffs[4] = 0.0;
644 : }
645 3 : windTurbine.PowerCoeffs[5] = state.dataIPShortCut->rNumericArgs(22); // Empirical power coefficient C6
646 3 : if (lNumericBlanks(22)) {
647 2 : windTurbine.PowerCoeffs[5] = 0.0;
648 : }
649 : }
650 :
651 1 : cAlphaArgs.deallocate();
652 1 : cAlphaFields.deallocate();
653 1 : cNumericFields.deallocate();
654 1 : rNumericArgs.deallocate();
655 1 : lAlphaBlanks.deallocate();
656 1 : lNumericBlanks.deallocate();
657 :
658 1 : if (ErrorsFound) {
659 0 : ShowFatalError(state, format("{} errors occurred in input. Program terminates.", CurrentModuleObject));
660 : }
661 :
662 4 : for (WindTurbineNum = 1; WindTurbineNum <= NumWindTurbines; ++WindTurbineNum) {
663 3 : auto &windTurbine = state.dataWindTurbine->WindTurbineSys(WindTurbineNum);
664 6 : SetupOutputVariable(state,
665 : "Generator Produced AC Electricity Rate",
666 : Constant::Units::W,
667 3 : windTurbine.Power,
668 : OutputProcessor::TimeStepType::System,
669 : OutputProcessor::StoreType::Average,
670 3 : windTurbine.Name);
671 6 : SetupOutputVariable(state,
672 : "Generator Produced AC Electricity Energy",
673 : Constant::Units::J,
674 3 : windTurbine.Energy,
675 : OutputProcessor::TimeStepType::System,
676 : OutputProcessor::StoreType::Sum,
677 3 : windTurbine.Name,
678 : Constant::eResource::ElectricityProduced,
679 : OutputProcessor::Group::Plant,
680 : OutputProcessor::EndUseCat::WindTurbine);
681 6 : SetupOutputVariable(state,
682 : "Generator Turbine Local Wind Speed",
683 : Constant::Units::m_s,
684 3 : windTurbine.LocalWindSpeed,
685 : OutputProcessor::TimeStepType::System,
686 : OutputProcessor::StoreType::Average,
687 3 : windTurbine.Name);
688 6 : SetupOutputVariable(state,
689 : "Generator Turbine Local Air Density",
690 : Constant::Units::kg_m3,
691 3 : windTurbine.LocalAirDensity,
692 : OutputProcessor::TimeStepType::System,
693 : OutputProcessor::StoreType::Average,
694 3 : windTurbine.Name);
695 6 : SetupOutputVariable(state,
696 : "Generator Turbine Tip Speed Ratio",
697 : Constant::Units::None,
698 3 : windTurbine.TipSpeedRatio,
699 : OutputProcessor::TimeStepType::System,
700 : OutputProcessor::StoreType::Average,
701 3 : windTurbine.Name);
702 3 : switch (windTurbine.rotorType) {
703 2 : case RotorType::HorizontalAxis: {
704 4 : SetupOutputVariable(state,
705 : "Generator Turbine Power Coefficient",
706 : Constant::Units::None,
707 2 : windTurbine.PowerCoeff,
708 : OutputProcessor::TimeStepType::System,
709 : OutputProcessor::StoreType::Average,
710 2 : windTurbine.Name);
711 2 : } break;
712 1 : case RotorType::VerticalAxis: {
713 2 : SetupOutputVariable(state,
714 : "Generator Turbine Chordal Component Velocity",
715 : Constant::Units::m_s,
716 1 : windTurbine.ChordalVel,
717 : OutputProcessor::TimeStepType::System,
718 : OutputProcessor::StoreType::Average,
719 1 : windTurbine.Name);
720 2 : SetupOutputVariable(state,
721 : "Generator Turbine Normal Component Velocity",
722 : Constant::Units::m_s,
723 1 : windTurbine.NormalVel,
724 : OutputProcessor::TimeStepType::System,
725 : OutputProcessor::StoreType::Average,
726 1 : windTurbine.Name);
727 2 : SetupOutputVariable(state,
728 : "Generator Turbine Relative Flow Velocity",
729 : Constant::Units::m_s,
730 1 : windTurbine.RelFlowVel,
731 : OutputProcessor::TimeStepType::System,
732 : OutputProcessor::StoreType::Average,
733 1 : windTurbine.Name);
734 2 : SetupOutputVariable(state,
735 : "Generator Turbine Attack Angle",
736 : Constant::Units::deg,
737 1 : windTurbine.AngOfAttack,
738 : OutputProcessor::TimeStepType::System,
739 : OutputProcessor::StoreType::Average,
740 1 : windTurbine.Name);
741 1 : } break;
742 0 : default:
743 0 : break;
744 : }
745 : }
746 1 : }
747 :
748 16236 : void InitWindTurbine(EnergyPlusData &state, int const WindTurbineNum)
749 : {
750 :
751 : // SUBROUTINE INFORMATION:
752 : // AUTHOR Daeho Kang
753 : // DATE WRITTEN Oct 2009
754 : // MODIFIED Linda K. Lawrie, December 2009 for reading stat file
755 :
756 : // PURPOSE OF THIS SUBROUTINE:
757 : // This subroutine reads monthly average wind speed from stat file and then
758 : // determines annual average wind speed. Differences between this TMY wind speed
759 : // and local wind speed that the user inputs are then factored.
760 : // IF the user has no local wind data and does not enter the local wind speed to be factored,
761 : // then the factor of 1 is assigned, so that wind speed estimated
762 : // at the particular rotor height is used with no factorization.
763 : // It also initializes module variables at each time step.
764 :
765 : static char const TabChr('\t'); // Tab character
766 :
767 16236 : Array1D<Real64> MonthWS(12);
768 :
769 : // Estimate average annual wind speed once
770 16236 : if (state.dataWindTurbine->MyOneTimeFlag) {
771 1 : Real64 AnnualTMYWS = 0.0;
772 1 : if (FileSystem::fileExists(state.files.inStatFilePath.filePath)) {
773 2 : auto statFile = state.files.inStatFilePath.open(state, "InitWindTurbine");
774 1 : bool wsStatFound = false; // logical noting that wind stats were found
775 255 : while (statFile.good()) { // end of file
776 255 : auto lineIn = statFile.readLine();
777 : // reconcile line with different versions of stat file
778 255 : size_t lnPtr = index(lineIn.data, "Wind Speed");
779 255 : if (lnPtr == std::string::npos) {
780 254 : continue;
781 : }
782 : // have hit correct section.
783 8 : while (statFile.good()) { // find daily avg line
784 8 : lineIn = statFile.readLine();
785 8 : lnPtr = index(lineIn.data, "Daily Avg");
786 8 : if (lnPtr == std::string::npos) {
787 7 : continue;
788 : }
789 : // tab delimited file
790 1 : lineIn.data.erase(0, lnPtr + 10);
791 1 : MonthWS = 0.0;
792 1 : wsStatFound = true;
793 1 : bool warningShown = false; // true if the <365 warning has already been shown
794 13 : for (int mon = 1; mon <= 12; ++mon) {
795 12 : lnPtr = index(lineIn.data, TabChr);
796 12 : if (lnPtr != 1) {
797 12 : if ((lnPtr == std::string::npos) || (!stripped(lineIn.data.substr(0, lnPtr)).empty())) {
798 12 : if (lnPtr != std::string::npos) {
799 12 : bool error = false;
800 12 : MonthWS(mon) = Util::ProcessNumber(lineIn.data.substr(0, lnPtr), error);
801 :
802 12 : if (error) {
803 : // probably should throw some error here
804 : }
805 :
806 12 : lineIn.data.erase(0, lnPtr + 1);
807 : }
808 : } else { // blank field
809 0 : if (!warningShown) {
810 0 : ShowWarningError(state,
811 0 : format("InitWindTurbine: read from {} file shows <365 days in weather file. Annual average "
812 : "wind speed used will be inaccurate.",
813 0 : state.files.inStatFilePath.filePath));
814 0 : lineIn.data.erase(0, lnPtr + 1);
815 0 : warningShown = true;
816 : }
817 : }
818 : } else { // two tabs in succession
819 0 : if (!warningShown) {
820 0 : ShowWarningError(state,
821 0 : format("InitWindTurbine: read from {} file shows <365 days in weather file. Annual average wind "
822 : "speed used will be inaccurate.",
823 0 : state.files.inStatFilePath.filePath));
824 0 : lineIn.data.erase(0, lnPtr + 1);
825 0 : warningShown = true;
826 : }
827 : }
828 : }
829 1 : break;
830 : }
831 1 : if (wsStatFound) {
832 1 : break;
833 : }
834 255 : }
835 1 : if (wsStatFound) {
836 1 : AnnualTMYWS = sum(MonthWS) / 12.0;
837 : } else {
838 0 : ShowWarningError(
839 : state, "InitWindTurbine: stat file did not include Wind Speed statistics. TMY Wind Speed adjusted at the height is used.");
840 : }
841 1 : } else { // No stat file
842 0 : ShowWarningError(state, "InitWindTurbine: stat file missing. TMY Wind Speed adjusted at the height is used.");
843 : }
844 :
845 : // assign this value to all the wind turbines once here
846 4 : for (auto &wt : state.dataWindTurbine->WindTurbineSys) {
847 3 : wt.AnnualTMYWS = AnnualTMYWS;
848 1 : }
849 :
850 1 : state.dataWindTurbine->MyOneTimeFlag = false;
851 : }
852 :
853 16236 : auto &windTurbine = state.dataWindTurbine->WindTurbineSys(WindTurbineNum);
854 :
855 : // Factor differences between TMY wind data and local wind data once
856 16236 : if (windTurbine.AnnualTMYWS > 0.0 && windTurbine.WSFactor == 0.0 && windTurbine.LocalAnnualAvgWS > 0) {
857 : // Convert the annual wind speed to the local wind speed at the height of the local station, then factor
858 3 : Real64 LocalTMYWS = windTurbine.AnnualTMYWS * state.dataEnvrn->WeatherFileWindModCoeff *
859 3 : std::pow(windTurbine.HeightForLocalWS / state.dataEnvrn->SiteWindBLHeight, state.dataEnvrn->SiteWindExp);
860 3 : windTurbine.WSFactor = LocalTMYWS / windTurbine.LocalAnnualAvgWS;
861 : }
862 : // Assign factor of 1.0 if no stat file or no input of local average wind speed
863 16236 : if (windTurbine.WSFactor == 0.0) {
864 0 : windTurbine.WSFactor = 1.0;
865 : }
866 :
867 : // Do every time step initialization
868 16236 : windTurbine.Power = 0.0;
869 16236 : windTurbine.TotPower = 0.0;
870 16236 : windTurbine.PowerCoeff = 0.0;
871 16236 : windTurbine.TipSpeedRatio = 0.0;
872 16236 : windTurbine.ChordalVel = 0.0;
873 16236 : windTurbine.NormalVel = 0.0;
874 16236 : windTurbine.RelFlowVel = 0.0;
875 16236 : windTurbine.AngOfAttack = 0.0;
876 16236 : windTurbine.TanForce = 0.0;
877 16236 : windTurbine.NorForce = 0.0;
878 16236 : windTurbine.TotTorque = 0.0;
879 16236 : }
880 :
881 16236 : void CalcWindTurbine(EnergyPlusData &state,
882 : int const WindTurbineNum, // System is on
883 : [[maybe_unused]] bool const RunFlag // System is on
884 : )
885 : {
886 :
887 : // SUBROUTINE INFORMATION:
888 : // AUTHOR Daeho Kang
889 : // DATE WRITTEN October 2009
890 :
891 : // REFERENCES:
892 : // Sathyajith Mathew. 2006. Wind Energy: Fundamental, Resource Analysis and Economics. Springer,
893 : // Chap. 2, pp. 11-15
894 : // Mazharul Islam, David S.K. Ting, and Amir Fartaj. 2008. Aerodynamic Models for Darrieus-type sSraight-bladed
895 : // Vertical Axis Wind Turbines. Renewable & Sustainable Energy Reviews, Volume 12, pp.1087-1109
896 :
897 : using DataEnvironment::OutBaroPressAt;
898 : using DataEnvironment::OutDryBulbTempAt;
899 : using DataEnvironment::OutWetBulbTempAt;
900 : using Psychrometrics::PsyRhoAirFnPbTdbW;
901 : using Psychrometrics::PsyWFnTdbTwbPb;
902 :
903 16236 : Real64 constexpr MaxTheta(90.0); // Maximum of theta
904 16236 : Real64 constexpr MaxDegree(360.0); // Maximum limit of outdoor air wind speed in m/s
905 16236 : Real64 constexpr SecInMin(60.0);
906 :
907 : Real64 LocalWindSpeed; // Ambient wind speed at the specific height in m/s
908 : Real64 RotorH; // Height of the rotor in m
909 : Real64 RotorD; // Diameter of the rotor in m
910 : Real64 LocalHumRat; // Local humidity ratio of the air at the specific height
911 : Real64 LocalAirDensity; // Local density at the specific height in m
912 : Real64 PowerCoeff; // Power coefficient
913 : Real64 SweptArea; // Swept area of the rotor in m2
914 16236 : Real64 WTPower(0.0); // Total Power generated by the turbine in the quasi-steady state in Watts
915 : Real64 Power; // Actual power of the turbine in Watts
916 : Real64 TipSpeedRatio; // Tip speed ratio (TSR)
917 : Real64 TipSpeedRatioAtI; // Tip speed ratio at the ith time step
918 : Real64 AzimuthAng; // Azimuth angle of blades in degree
919 : Real64 ChordalVel; // Chordal velocity component in m/s
920 : Real64 NormalVel; // Normal velocity component in m/s
921 : Real64 AngOfAttack; // Angle of attack of a single blade in degree
922 : Real64 RelFlowVel; // Relative flow velocity component in m/s
923 : Real64 TanForce; // Tangential force in N.m
924 : Real64 NorForce; // Normal force in N.m
925 : Real64 RotorVel; // Rotor velocity in m/s
926 : Real64 AvgTanForce; // Average tangential force in N.m
927 : Real64 Constant; // Constants within integrand of tangential force
928 : Real64 IntRelFlowVel; // Integration of relative flow velocity
929 : Real64 TotTorque; // Total torque for the number of blades
930 : Real64 Omega; // Angular velocity of rotor in rad/s
931 : Real64 TanForceCoeff; // Tangential force coefficient
932 : Real64 NorForceCoeff; // Normal force coefficient
933 : Real64 Period; // Period of sine and cosine functions
934 : Real64 C1; // Empirical power coefficient C1
935 : Real64 C2; // Empirical power coefficient C2
936 : Real64 C3; // Empirical power coefficient C3
937 : Real64 C4; // Empirical power coefficient C4
938 : Real64 C5; // Empirical power coefficient C5
939 : Real64 C6; // Empirical power coefficient C6
940 : Real64 LocalTemp; // Ambient air temperature at the height in degree C
941 : Real64 LocalPress; // Local atmospheric pressure at the particular height in Pa
942 : Real64 InducedVel; // Induced velocity on the rotor in m/s
943 : // unused REAL(r64) :: SysEfficiency ! Overall wind turbine system efficiency including generator and inverter
944 : Real64 MaxPowerCoeff; // Maximum power coefficient
945 : Real64 RotorSpeed; // Speed of rotors
946 :
947 16236 : auto &windTurbine = state.dataWindTurbine->WindTurbineSys(WindTurbineNum);
948 : // Estimate local velocity and density
949 16236 : RotorH = windTurbine.RotorHeight;
950 16236 : RotorD = windTurbine.RotorDiameter;
951 16236 : RotorSpeed = windTurbine.RatedRotorSpeed;
952 16236 : LocalTemp = OutDryBulbTempAt(state, RotorH);
953 16236 : LocalPress = OutBaroPressAt(state, RotorH);
954 16236 : LocalHumRat = PsyWFnTdbTwbPb(state, LocalTemp, OutWetBulbTempAt(state, RotorH), LocalPress);
955 16236 : LocalAirDensity = PsyRhoAirFnPbTdbW(state, LocalPress, LocalTemp, LocalHumRat);
956 16236 : LocalWindSpeed = DataEnvironment::WindSpeedAt(state, RotorH);
957 16236 : LocalWindSpeed /= windTurbine.WSFactor;
958 :
959 : // Check wind conditions for system operation
960 16236 : if (windTurbine.availSched->getCurrentVal() > 0 && LocalWindSpeed > windTurbine.CutInSpeed && LocalWindSpeed < windTurbine.CutOutSpeed) {
961 :
962 : // System is on
963 16209 : Period = 2.0 * Constant::Pi;
964 16209 : Omega = (RotorSpeed * Period) / SecInMin;
965 16209 : SweptArea = (Constant::Pi * pow_2(RotorD)) / 4;
966 16209 : TipSpeedRatio = (Omega * (RotorD / 2.0)) / LocalWindSpeed;
967 :
968 : // Limit maximum tip speed ratio
969 16209 : if (TipSpeedRatio > windTurbine.MaxTipSpeedRatio) {
970 10815 : TipSpeedRatio = windTurbine.MaxTipSpeedRatio;
971 : }
972 :
973 16209 : switch (windTurbine.rotorType) {
974 10806 : case RotorType::HorizontalAxis: { // Horizontal axis wind turbine
975 10806 : MaxPowerCoeff = windTurbine.MaxPowerCoeff;
976 : // Check if empirical constants are available
977 10806 : C1 = windTurbine.PowerCoeffs[0];
978 10806 : C2 = windTurbine.PowerCoeffs[1];
979 10806 : C3 = windTurbine.PowerCoeffs[2];
980 10806 : C4 = windTurbine.PowerCoeffs[3];
981 10806 : C5 = windTurbine.PowerCoeffs[4];
982 10806 : C6 = windTurbine.PowerCoeffs[5];
983 :
984 10806 : Real64 const LocalWindSpeed_3(pow_3(LocalWindSpeed));
985 10806 : if (C1 > 0.0 && C2 > 0.0 && C3 > 0.0 && C4 >= 0.0 && C5 > 0.0 && C6 > 0.0) {
986 : // Analytical approximation
987 : // Maximum power, i.e., rotor speed is at maximum, and pitch angle is zero
988 : // TipSpeedRatioAtI = 1.0 / ( ( 1.0 / ( TipSpeedRatio + 0.08 * PitchAngle ) ) - ( 0.035 / ( pow_3( PitchAngle ) + 1.0 ) ) );
989 : // //Tuned PitchAngle is zero
990 5403 : TipSpeedRatioAtI = TipSpeedRatio / (1.0 - (TipSpeedRatio * 0.035));
991 : // PowerCoeff = C1 * ( ( C2 / TipSpeedRatioAtI ) - ( C3 * PitchAngle ) - ( C4 * std::pow( PitchAngle, 1.5 ) ) - C5 ) * (
992 : // std::exp( -( C6 / TipSpeedRatioAtI ) ) ); //Tuned PitchAngle is zero
993 5403 : PowerCoeff = C1 * ((C2 / TipSpeedRatioAtI) - C5) * std::exp(-(C6 / TipSpeedRatioAtI));
994 5403 : if (PowerCoeff > MaxPowerCoeff) {
995 0 : PowerCoeff = MaxPowerCoeff;
996 : }
997 5403 : WTPower = 0.5 * LocalAirDensity * PowerCoeff * SweptArea * LocalWindSpeed_3;
998 : } else { // Simple approximation
999 5403 : WTPower = 0.5 * LocalAirDensity * SweptArea * LocalWindSpeed_3 * MaxPowerCoeff;
1000 5403 : PowerCoeff = MaxPowerCoeff;
1001 : }
1002 : // Maximum of rated power
1003 10806 : if (LocalWindSpeed >= windTurbine.RatedWindSpeed || WTPower > windTurbine.RatedPower) {
1004 0 : WTPower = windTurbine.RatedPower;
1005 0 : PowerCoeff = WTPower / (0.5 * LocalAirDensity * SweptArea * LocalWindSpeed_3);
1006 : }
1007 : // Recalculated Cp at the rated power
1008 10806 : windTurbine.PowerCoeff = PowerCoeff;
1009 10806 : } break;
1010 5403 : case RotorType::VerticalAxis: { // Vertical axis wind turbine
1011 5403 : RotorVel = Omega * (RotorD / 2.0);
1012 : // Recalculated omega, if TSR is greater than the maximum
1013 5403 : if (TipSpeedRatio >= windTurbine.MaxTipSpeedRatio) {
1014 5403 : RotorVel = LocalWindSpeed * windTurbine.MaxTipSpeedRatio;
1015 5403 : Omega = RotorVel / (RotorD / 2.0);
1016 : }
1017 :
1018 5403 : AzimuthAng = MaxDegree / windTurbine.NumOfBlade;
1019 : // Azimuth angle between zero and 90 degree
1020 5403 : if (AzimuthAng > MaxTheta) { // Number of blades is 2 or 3
1021 5403 : AzimuthAng -= MaxTheta;
1022 5403 : if (AzimuthAng == MaxTheta) { // 2 blades
1023 0 : AzimuthAng = 0.0;
1024 : }
1025 0 : } else if (AzimuthAng == MaxTheta) { // 4 blades
1026 0 : AzimuthAng = 0.0;
1027 : }
1028 :
1029 5403 : InducedVel = LocalWindSpeed * 2.0 / 3.0;
1030 : // Velocity components
1031 5403 : Real64 const sin_AzimuthAng = std::sin(AzimuthAng * Constant::DegToRad);
1032 5403 : Real64 const cos_AzimuthAng = std::cos(AzimuthAng * Constant::DegToRad);
1033 5403 : ChordalVel = RotorVel + InducedVel * cos_AzimuthAng;
1034 5403 : NormalVel = InducedVel * sin_AzimuthAng;
1035 5403 : RelFlowVel = std::sqrt(pow_2(ChordalVel) + pow_2(NormalVel));
1036 :
1037 : // Angle of attack
1038 5403 : AngOfAttack = std::atan((sin_AzimuthAng / ((RotorVel / LocalWindSpeed) / (InducedVel / LocalWindSpeed) + cos_AzimuthAng)));
1039 :
1040 : // Force coefficients
1041 5403 : Real64 const sin_AngOfAttack = std::sin(AngOfAttack * Constant::DegToRad);
1042 5403 : Real64 const cos_AngOfAttack = std::cos(AngOfAttack * Constant::DegToRad);
1043 5403 : TanForceCoeff = std::abs(windTurbine.LiftCoeff * sin_AngOfAttack - windTurbine.DragCoeff * cos_AngOfAttack);
1044 5403 : NorForceCoeff = windTurbine.LiftCoeff * cos_AngOfAttack + windTurbine.DragCoeff * sin_AngOfAttack;
1045 :
1046 : // Net tangential and normal forces
1047 5403 : Real64 const RelFlowVel_2(pow_2(RelFlowVel));
1048 5403 : Real64 const density_fac(0.5 * LocalAirDensity * windTurbine.ChordArea * RelFlowVel_2);
1049 5403 : TanForce = TanForceCoeff * density_fac;
1050 5403 : NorForce = NorForceCoeff * density_fac;
1051 5403 : Constant = (1.0 / Period) * (TanForce / RelFlowVel_2);
1052 :
1053 : // Relative flow velocity is the only function of theta in net tangential force
1054 : // Integral of cos(theta) on zero to 2pi goes to zero
1055 : // Integrate constants only
1056 5403 : IntRelFlowVel = pow_2(RotorVel) * Period + pow_2(InducedVel) * Period;
1057 :
1058 : // Average tangential force on a single blade
1059 5403 : AvgTanForce = Constant * IntRelFlowVel;
1060 5403 : TotTorque = windTurbine.NumOfBlade * AvgTanForce * (RotorD / 2.0);
1061 5403 : WTPower = TotTorque * Omega;
1062 :
1063 : // Check if power produced is greater than maximum or rated power
1064 5403 : if (WTPower > windTurbine.RatedPower) {
1065 5403 : WTPower = windTurbine.RatedPower;
1066 : }
1067 :
1068 5403 : windTurbine.ChordalVel = ChordalVel;
1069 5403 : windTurbine.NormalVel = NormalVel;
1070 5403 : windTurbine.RelFlowVel = RelFlowVel;
1071 5403 : windTurbine.TanForce = TanForce;
1072 5403 : windTurbine.NorForce = NorForce;
1073 5403 : windTurbine.TotTorque = TotTorque;
1074 5403 : } break;
1075 0 : default: {
1076 0 : assert(false);
1077 : } break;
1078 : }
1079 :
1080 16209 : if (WTPower > windTurbine.RatedPower) {
1081 0 : WTPower = windTurbine.RatedPower;
1082 : }
1083 :
1084 : // Actual power generated by the wind turbine system
1085 16209 : Power = WTPower * windTurbine.SysEfficiency;
1086 :
1087 16209 : windTurbine.Power = Power;
1088 16209 : windTurbine.TotPower = WTPower;
1089 16209 : windTurbine.LocalWindSpeed = LocalWindSpeed;
1090 16209 : windTurbine.LocalAirDensity = LocalAirDensity;
1091 16209 : windTurbine.TipSpeedRatio = TipSpeedRatio;
1092 :
1093 : } else { // System is off
1094 27 : windTurbine.Power = 0.0;
1095 27 : windTurbine.TotPower = 0.0;
1096 27 : windTurbine.PowerCoeff = 0.0;
1097 27 : windTurbine.LocalWindSpeed = LocalWindSpeed;
1098 27 : windTurbine.LocalAirDensity = LocalAirDensity;
1099 27 : windTurbine.TipSpeedRatio = 0.0;
1100 27 : windTurbine.ChordalVel = 0.0;
1101 27 : windTurbine.NormalVel = 0.0;
1102 27 : windTurbine.RelFlowVel = 0.0;
1103 27 : windTurbine.AngOfAttack = 0.0;
1104 27 : windTurbine.TanForce = 0.0;
1105 27 : windTurbine.NorForce = 0.0;
1106 27 : windTurbine.TotTorque = 0.0;
1107 : }
1108 16236 : }
1109 :
1110 16236 : void ReportWindTurbine(EnergyPlusData &state, int const WindTurbineNum)
1111 : {
1112 : // SUBROUTINE INFORMATION:
1113 : // AUTHOR Daeho Kang
1114 : // DATE WRITTEN October 2009
1115 :
1116 : // PURPOSE OF THIS SUBROUTINE:
1117 : // This subroutine fills remaining report variables.
1118 :
1119 16236 : Real64 TimeStepSysSec = state.dataHVACGlobal->TimeStepSysSec;
1120 16236 : auto &windTurbine = state.dataWindTurbine->WindTurbineSys(WindTurbineNum);
1121 :
1122 16236 : windTurbine.Energy = windTurbine.Power * TimeStepSysSec;
1123 16236 : }
1124 :
1125 : //*****************************************************************************************
1126 :
1127 : } // namespace WindTurbine
1128 :
1129 : } // namespace EnergyPlus
|
