windSimSuite.main

1 #!/usr/bin/env python 2 # -*- coding: utf-8 -*- 3 # 4 # This file is part of MBDyn sim suite. 5 # Copyright (C) 2007 André ESPAZE, as part of a Master thesis supervised by 6 # Martin O.L.Hansen (DTU) and Nicolas Chauvat (Logilab) 7 8 # MBDyn sim suite is free software; you can redistribute it and/or modify 9 # it under the terms of the GNU General Public License as published by 10 # the Free Software Foundation; either version 2 of the License, or 11 # (at your option) any later version. 12 # 13 # This program is distributed in the hope that it will be useful, 14 # but WITHOUT ANY WARRANTY; without even the implied warranty of 15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 # GNU General Public License for more details. 17 # 18 # You should have received a copy of the GNU General Public License 19 # along with this program; if not, write to the Free Software 20 # Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 21 # 22 """The wind turbine and simulation definitions. This module 23 derived from the main module of the mbdyn package to provide 24 a simulation where a L{WindTurbine} can be added. 25 """ 26 import numpy as N 27 28 from windSimSuite.common import BasicObject 29 from windSimSuite.rotor import Rotor 30 from windSimSuite.nacelle import Nacelle 31 from windSimSuite.tower import Tower 32 33 from mbdyn.main import Simulation as SimulationBase 34 from mbdyn.main import PyMBDynFile as PyMBDynFileBase 35 36 37 CLASS = {"rotor" : Rotor, 38 "nacelle" : Nacelle, 39 "tower" : Tower} 40 41

42 -class WindTurbine(BasicObject):

43 """The wind turbine manipulating all the others components 44 """ 45

46 - def __init__(self, name="Wind turbine"):

47 BasicObject.__init__(self, name) 48 self.children_names = [] 49 50 self.tower = None 51 self.has_tower = False 52 self.nacelle = None 53 self.has_nacelle = False 54 self.rotor = None 55 self.has_rotor = False 56 self.has_update = False 57 58 self.mbdyn_coupling = True 59 self.own_para_names = ["children_names", "has_tower", 60 "has_nacelle", "has_rotor"]

61

62 - def set_tower(self, tower):

63 """Set the C{Tower} instance from L{windSimSuite.tower}""" 64 self.has_tower = True 65 self.tower = tower 66 self.children_names += ["tower"]

67

68 - def set_rotor(self, rotor):

69 """Set the C{Rotor} instance from L{windSimSuite.rotor}""" 70 self.has_rotor = True 71 self.rotor = rotor 72 self.children_names += ["rotor"]

73

74 - def set_nacelle(self, nacelle):

75 """Set the C{Nacelle} instance from L{windSimSuite.nacelle}""" 76 self.has_nacelle = True 77 self.nacelle = nacelle 78 self.children_names += ["nacelle"]

79

80 - def set_mbdyn_coupling(self, boolean):

81 """If False, only the unsteady computation will be run""" 82 self.mbdyn_coupling = boolean

83

84 - def _update_reference_frames(self):

85 """The reference frames for unsteady calculations 86 are updated according to the user values.""" 87 if self.has_nacelle: 88 self.rotor.prepend_references(self.nacelle.references.refs)

89

90 - def update_data(self):

91 """Update the attributes stored on the wind turbine on 92 all the objects that belong to it. 93 For example the wing angle information, set on the rotor, 94 is now transmitted to the blades. The reference frames 95 used for the unsteady calculation are also updated in 96 case MBDyn is not used.""" 97 if not self.has_update: 98 if not self.mbdyn_coupling: 99 self._update_reference_frames() 100 if self.has_rotor: 101 self.rotor.update_data() 102 self.has_update = True

103

104 - def create(self):

105 """Init the coupling between MBDyn and the unsteady BEM code. 106 MBDyn objects are created by being attributes of wind turbine 107 components. According to their type, they are stored in C{refs}, 108 C{nodes} or C{elts}. 109 Until now, this method is in development and is supposed to be 110 defined by a specific model. The reference frames, nodes and 111 elements defined will then be added on the simulation. 112 """ 113 pass

114

115 - def add_on_simulation(self, simu):

116 """Add the wind turbine objects add all its components, 117 reference frames, nodes and elements on the simulation. 118 Each method will scan the C{refs}, C{nodes} and C{elts} list 119 attributes.""" 120 for component_name in ["tower", "nacelle", "rotor"]: 121 if getattr(self, "has_" + component_name): 122 getattr(self, component_name).add_on_simulation(simu)

123

124 - def init_results(self):

125 """Suppose to initialize all the wind turbine object 126 results. Until now, results are only collected on 127 the rotor.""" 128 if self.has_rotor: 129 self.rotor.init_results()

130

131 - def update(self):

132 """Supposed to udpate some of the wind turbine 133 parameters according to the MBDyn new state. 134 For the same reason as L{create}, this method is in 135 development and this is the model that is supposed 136 to supply that method.""" 137 pass

138

139 - def collect_parameters(self):

140 """Store the wind turbine parameters in a dictionary""" 141 self.collect_own_parameters() 142 for children_name in self.children_names: 143 children = getattr(self, children_name) 144 children.collect_parameters() 145 self.para[children_name] = children.para

146

147 - def set_parameters(self, para):

148 """Set the wind turbine parameters from the C{para} dictionary""" 149 self.set_own_parameters(para) 150 for children_name in self.children_names: 151 children = CLASS[children_name]() 152 children.set_parameters(para[children_name]) 153 setattr(self, children_name, children)

154 155

156 -class AerodynamicsSimulation:

157 """The unsteady Blade Element Momentum method 158 applied on the wind turbine from the simulation. 159 """ 160

161 - def __init__(self):

162 # Overidden by the main Simulation 163 self.turbine = None 164 self.wind_speed = None 165 self.air_density = None 166 self.shear = 0.2 167 168 self.current_time = 0. 169 self.time_step = 0. 170 self.init_time = 1. 171 self.induction_fmax = 0.5 172 self.adjust_factor = 0.6 173 174 self.has_init_unsteady = False 175 self._try_enable_dynamic_wake = self._enable_dynamic_wake

176

177 - def set_init_time(self, init_time):

178 """Set the time to run before starting 179 to save the results and applying the dynamic wake""" 180 self.init_time = init_time

181

182 - def set_wind_speed(self, value):

183 """Set the wind spee value in M{m/s}""" 184 self.wind_speed = value

185

186 - def set_air_density(self, value):

187 """Set the air density in M{kg/m^3}.""" 188 self.air_density = value

189

190 - def set_shear_value(self, value):

191 """Set the shear value for the wind profile. 192 The default is 0.2, will depend of the site roughness""" 193 self.shear = value

194

195 - def get_velocity_at_height(self, height):

196 """Return the velocity at a particular height from 197 the wind profile""" 198 if height <= self.turbine.tower.aero_height: 199 return self.wind_speed * \ 200 (height/self.turbine.tower.aero_height)**self.shear 201 else: 202 return self.wind_speed

203

204 - def get_velocity_at(self, pt_array):

205 """Return the velocity for a point in the absolute 206 referential. The wind tower model is applied. 207 208 @type pt_array: a Numpy array 209 @param pt_array: a position vector of shape 3x1 210 211 In the referential used, C{pt_array[2]} 212 will be the height along the M{z} axis. 213 214 @rtype: a Numpy array 215 @return: the wind vector of size 3x1 216 """ 217 xval = pt_array[0][0] 218 yval = pt_array[1][0] 219 height = pt_array[2][0] 220 221 wind_speed = self.get_velocity_at_height(height) 222 #if pt_array[0] < self.turbine.tower.correction_height: 223 if height < self.turbine.tower.aero_height: 224 tower_radius = self.turbine.tower.get_radius_at(height) 225 pt_radius = N.sqrt(xval**2 + yval**2) 226 cos_theta = xval / pt_radius 227 sin_theta = yval / pt_radius 228 radius_velocity = wind_speed * \ 229 (1. - (tower_radius/pt_radius)**2)*cos_theta 230 theta_velocity = - wind_speed * \ 231 (1. + (tower_radius/pt_radius)**2)*sin_theta 232 x_velocity = radius_velocity * cos_theta - \ 233 theta_velocity * sin_theta 234 y_velocity = radius_velocity * sin_theta + \ 235 theta_velocity * cos_theta 236 # Correction 237 #y_velocity *= self.turbine.rotor.get_velocity_correction() 238 return N.array([ [x_velocity], 239 [y_velocity], 240 [0.] ]) 241 else: 242 return N.array([ [wind_speed], 243 [0.], 244 [0.] ])

245

246 - def enable_dynamic_wake(self, enable):

247 """If C{True}, will enable the dynamic wake""" 248 self._enable_dynamic_wake(enable) 249 if enable: 250 self._try_enable_dynamic_wake = self._enable_dynamic_wake 251 else: 252 self._try_enable_dynamic_wake = self._no_dynamic_wake

253

254 - def _no_dynamic_wake(self, enable):

255 """The method executed for an absent dynamic wake""" 256 pass

257

258 - def _enable_dynamic_wake(self, enable):

259 """If C{True}, enable the dynamic wake on each aerodynamic 260 section. A private method, the user needs to use 261 L{enable_dynamic_wake}.""" 262 for blade in self.turbine.rotor.blades: 263 for section in blade.sections: 264 if enable: 265 section.calculate_induced_velocity = \ 266 section.apply_dynamic_wake 267 else: 268 section.calculate_induced_velocity = \ 269 section.disable_dynamic_wake

270

271 - def init_aero_unsteady(self):

272 """Initialize the aerodynamics calculation of the unsteady 273 BEM code""" 274 self.turbine.rotor.init_unsteady() 275 init_time_list = N.arange(0., 276 self.init_time + self.time_step, 277 self.time_step) 278 self._try_enable_dynamic_wake(False) 279 start_angle = self.turbine.rotor.wing_angle["rad"] 280 281 for time in init_time_list: 282 self.run_unsteady_at_time(time) 283 284 self.turbine.rotor.set_wing_angle(start_angle, "rad") 285 self._try_enable_dynamic_wake(True) 286 self.current_time = -self.time_step 287 self.has_init_unsteady = True

288

289 - def init_unsteady(self):

290 """Will be defined by the derived class""" 291 pass

292

293 - def save_unsteady_state_at(self, time):

294 """Will be defined by the derived class""" 295 pass

296

297 - def run_unsteady_for(self, time_length, nb_pts=100):

298 """Run the unsteady Blade Element Momentum code 299 for the C{time_lenght}. C{nb_pts} is an integer that 300 gives the number of calculation points for one revolution. 301 This method is used only when there is no coupling with MBDyn, 302 this is else the MBDyn solver that gives the current time.""" 303 revolution_period = 2*N.pi / \ 304 self.turbine.rotor.rotational_speed["rad/s"] 305 self.time_step = revolution_period / nb_pts 306 if not self.has_init_unsteady: 307 self.init_unsteady() 308 final_time = time_length + self.current_time 309 time_list = N.arange(self.current_time + self.time_step, 310 final_time + self.time_step, 311 self.time_step) 312 # Should be called from MBDyn but will not in that case 313 # Strange notation too hightlight that 'init_results' 314 # can not be defined in that class 315 getattr(self, "init_results")() 316 self.turbine.init_results() 317 318 for time in time_list: 319 self.run_unsteady_at_time(time) 320 self.save_unsteady_state_at(time) 321 self.current_time = time

322

323 - def run_unsteady_at_time(self, time):

324 """Run the unsteady BEM code calculations on the whole 325 turbine at a given time. At the end, nothing is saved""" 326 rotor = self.turbine.rotor 327 rotor.set_wing_angle_at_time(time) 328 vec_to_shaft = self.turbine.tower.get_abs_vector() + \ 329 self.turbine.nacelle.get_abs_vector() 330 for blade in rotor.blades: 331 for section in blade.sections: 332 matrix = blade.references.get_matrix_from("pitch") 333 abs_rvector = N.asarray(vec_to_shaft + \ 334 matrix * section.rvector) 335 abs_velocity = self.get_velocity_at(abs_rvector) 336 section.wind_speed = matrix.T * abs_velocity 337 section.calculate_relative_velocity( 338 rotor.rotational_speed, 339 blade.angle["cone"], 340 blade.references["pitch"].matrix 341 ) 342 section.calculate_pressures(self.air_density, 343 self.time_step, 344 blade.angle["pitch"]) 345 prandtlf = rotor.get_prandtl_factor(blade, section) 346 section.calculate_qs_induced_velocity(prandtlf, 347 self.air_density, 348 rotor.nb_blades) 349 section.calculate_induced_velocity(self.induction_fmax, 350 self.adjust_factor, 351 self.wind_speed, 352 self.time_step, 353 blade.radius) 354 blade.compute_forces_and_power(rotor.rotational_speed) 355 rotor.update_induced_velocity() 356 rotor.get_forces_and_power()

357

358 - def save_unsteady(self):

359 """Save the unsteady BEM code on the rotor.""" 360 self.turbine.rotor.save()

361 362

363 -class Simulation(AerodynamicsSimulation, SimulationBase):

364 """A specific simulation for the wind turbine. 365 The gathering between the aerodynamics calculation and 366 MBDyn. 367 """ 368

369 - def __init__(self, name="mbdyn simu", title="Mbdyn used from Python", 370 keep=False, overwrite=False, mbdyn_coupling=True):

371 AerodynamicsSimulation.__init__(self) 372 SimulationBase.__init__(self, name, title, keep, overwrite) 373 self.mbdyn_coupling = mbdyn_coupling 374 375 self.has_turbine = False 376 self.turbine_para = None 377 378 self.children_names = ["references"] 379 self.own_para_names += ["has_turbine", 380 "children_names", 381 "turbine_para", 382 "wind_speed", 383 "air_density"]

384

385 - def set_turbine(self, turbine):

386 """Set the L{WindTurbine} instance""" 387 turbine.set_mbdyn_coupling(self.mbdyn_coupling) 388 turbine.update_data() 389 390 self.turbine = turbine 391 self.turbine_para = turbine.para 392 self.children_names += ["turbine"] 393 self.has_turbine = True

394

395 - def init_unsteady(self):

396 """Init the unsteady BEM code by running it 397 for the 'init_time'""" 398 if self.mbdyn_coupling: 399 # FIX ME: Calculation of the time step for 400 # initializing the unsteady BEM 401 self.time_step = 0.027 402 # May be to small, not a good solution 403 #self.time_step = self.integrator["time step"] 404 self.init_aero_unsteady()

405

406 - def init(self):

407 """Initialize the simulation according to WraptMBDyn. 408 The Python object get references to the object from 409 the bindings module, interfacing the C++ 410 ones. This method is called by the C{run_full} method 411 of the C{SimulationBase}. It also overwrites the C{init} 412 of the C{SimulationBase}.""" 413 if not self.has_init: 414 self._init_coupling() 415 self.turbine.add_on_simulation(self) 416 self._init_wrapt_mbdyn() 417 self._init_wind_simulation() 418 self.has_init = True 419 else: 420 print "'init()' already called on the Simulation"

421

422 - def _init_coupling(self):

423 """Init the unsteady code and stabilize the aerodynamics loads. 424 Then create the wind turbine with MBDyn objects, that take into 425 account the ones from the unsteady BEM code. 426 Then run the MBDyn simulation, first L{init} will be called. 427 Then L{save} and L{update} are called at each MBDyn time step.""" 428 self.init_unsteady() 429 self.turbine.create()

430

431 - def _init_wind_simulation(self):

432 """The Python objects are getting a reference of the C++ ones. 433 The results are also initialized.""" 434 self.turbine.init_results() 435 # "self.init_results" is called by init_simulation 436 self.init_simulation()

437

438 - def run(self):

439 """Run the full simulation""" 440 self.run_full()

441

442 - def update(self):

443 """Update the wind turbine simulation. 444 Called by the C{run_full()} method at each MBDyn iteration. 445 This method is doing the coupling between MBDyn and 446 the unsteady BEM code. First the state of the wind turbine 447 is updated according to the MBDyn objects that are attached 448 to it. Then the unsteady BEM code is run from this new 449 configuration. Finally the MBDyn elements are updated.""" 450 self.turbine.update() 451 self.run_unsteady_at_time(self.wm.solver.current_time) 452 self.update_elts()

453

454 - def save_unsteady_state_at(self, time):

455 """Save the unsteady state at a particular time. This function 456 is used only when C{mbdyn_coupling} is set to C{False}. Else 457 the time is saved from the MBDyn simulation.""" 458 self.time = time 459 self.save_results() 460 self.save_unsteady()

461

462 - def save(self):

463 """The save method called by C{WraptMBDyn}. This method is called 464 before L{update}, because the values are saved just after a successful 465 time step. L{save_mbdyn} will also save the simulation 466 status (by calling the C{save_results} method.""" 467 self.save_mbdyn() 468 self.save_unsteady()

469

470 - def collect_parameters(self):

471 """Collect the simulation parameters into a dictionary""" 472 self.collect_own_parameters() 473 self.refs.collect_parameters() 474 self.para["refs"] = self.refs.para 475 if self.has_turbine: 476 self.turbine.collect_parameters()

477

478 - def set_parameters(self, para):

479 """Set the simulation parameters from the C{para} dictionary""" 480 self.set_own_parameters(para) 481 self.refs.set_parameters(para["refs"]) 482 if self.has_turbine: 483 self.turbine = WindTurbine() 484 self.turbine.set_parameters(self.turbine_para) 485 self.turbine.add_on_simulation(self)

486 487

488 -class PyMBDynFile(PyMBDynFileBase):

489 """The container to save a wind simulation. 490 The only difference with its base class is now to 491 save a L{Simulation} class suited to wind engineering. 492 """ 493

494 - def __init__(self, name, mode="r", path=".", autoload=True):

495 PyMBDynFileBase.__init__(self, name, mode, path, autoload=False) 496 self.CLASS = Simulation 497 if (mode=="r") and autoload: 498 self.load()

499

Source Code for Module windSimSuite.main