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