windSimSuite.models.tjaerborg

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 Tjaereborg turbine model. Still in very heavy development 23 but in progress. The documentation is incomplete as the design 24 is being done. This module manipulates however the objects from the 25 L{windSimSuite} and C{mbdyn} packages. 26 27 For more information on the simulated turbine, go to the 28 U{wind data site<http://www.winddata.com/>} done by 29 U{Kurt Schaldemose Hansen<http://www.dtu.dk/Service/Telefonbog.aspx?id=789&type=publications&lg=showcommon>}. 30 Navigate to 'Ressources Data -> 31 tjare -> description'. You can also 31 have a look to the link number 32.""" 32 import numpy as N 33 import mbdyn as M 34 import windSimSuite as W 35 36

37 -class Tower(W.Tower):

38 - def __init__(self):

39 W.Tower.__init__(self) 40 # The default values 41 self.top_radius = 2.125 42 self.middle_radius = 2.375 43 self.middle_height = 28. 44 self.calculate_radius_coeff() 45 46 self.height = 0.

47

48 - def calculate_radius_coeff(self):

49 self._radius_coeff = (self.middle_radius - self.top_radius) / \ 50 (self.aero_height - self.middle_height)

51

52 - def set_height(self, height):

53 self.set_aero_height(height) 54 self.set_correction_height(height + 2.) 55 self.height = height

56

57 - def get_radius_at(self, height):

58 """The diameter is going from middle radius at the middle height 59 to the top radius. For a given height, M{x}, the radius, M{y}, will be:: 60 61 y = (middle-radius - top_radius)/(height - middle_height)*(height - x) 62 + top_radius""" 63 if height < self.middle_height: 64 return self.middle_radius 65 else: 66 return self._radius_coeff * (self.aero_height - height) + \ 67 self.top_radius

68

69 - def create_references(self):

70 base_ref = M.ReferenceFrame(self.name + " base") 71 self.refs.add(base_ref, "base") 72 73 top_ref = M.ReferenceFrame(self.name + " mid") 74 top_ref.set_relative_from(base_ref) 75 top_ref.set_position(0., 0., 0.5 * self.height) 76 self.refs.add(top_ref, "mid") 77 78 top_ref = M.ReferenceFrame(self.name + " top") 79 top_ref.set_relative_from(base_ref) 80 top_ref.set_position(0., 0., self.height) 81 self.refs.add(top_ref, "top")

82 83

84 -class Nacelle(W.Nacelle):

85

86 - def __init__(self, name="Nacelle"):

87 W.Nacelle.__init__(self, name) 88 self._yaw = W.Angle(0., "deg") 89 self._tilt = W.Angle(0., "deg") 90 91 self.tower = None

92

93 - def set_yaw(self, value, unit_key="deg"):

94 self._yaw[unit_key] = value

95

96 - def set_tilt(self, value, unit_key="deg"):

97 self._tilt[unit_key] = value

98

99 - def set_tower(self, tower):

100 self.tower = tower

101

102 - def create_references_from(self, wing_angle):

103 yaw_ref = M.ReferenceFrame(self.name + " yaw") 104 yaw_ref.set_relative_from(self.tower.refs.get("top")) 105 yaw_ref.set_orientation_matrix(three=(0., 0., 1.), 106 one=(N.cos(self._yaw["rad"]), 107 N.sin(self._yaw["rad"]), 108 0.)) 109 self.refs.add(yaw_ref, "yaw") 110 111 tilt_ref = M.ReferenceFrame(self.name + " tilt") 112 tilt_ref.set_relative_from(yaw_ref) 113 tilt_ref.set_orientation_matrix(two=(0., 1., 0.), 114 one=(N.cos(self._tilt["rad"]), 115 0., 116 -N.sin(self._tilt["rad"]))) 117 self.refs.add(tilt_ref, "tilt") 118 119 ref = M.ReferenceFrame(self.name) 120 ref.set_relative_from(tilt_ref) 121 ref.set_position(-self.length, 0., 0.) 122 self.refs.add(ref, "fix")

123 124

125 -class HubNode(M.StructuralNode):

126 """The node of the hub used by the rotor to 127 find the rotational speed""" 128

129 - def get_matrix(self):

130 """Get the rotation matrix of the node, 131 used as a transformation matrix between the absolute reference frame 132 and the local reference frame, the one of the rotor""" 133 self.array = self.mbdyn_inst.get_rotation_matrix() 134 self.matrix = N.asmatrix(self.array)

135 136

137 -class GeneratorCouple(M.BindingsCouple):

138

139 - def do_as_init(self):

140 self._direction = N.array([ [-1.], 141 [0.], 142 [0.] ]) 143 self.amplitude = self.get_amplitude(self.simulation)

144

145 - def get_amplitude(self, simu):

146 """Return the amplitude of the torque. 147 The generator model goes at that place. 148 Supposed to be overwritten by the input file 149 and the user model. 150 """ 151 pass

152

153 - def update(self, simu):

154 self.amplitude = self.get_amplitude(simu) 155 value = self._direction * self.amplitude 156 rot_matrix = N.asmatrix(self.node.mbdyn_inst.get_rotation_matrix()) 157 self.mbdyn_inst.set_value(N.asarray(rot_matrix * value))

158 159

160 -class Hub(W.Hub):

161

162 - def __init__(self, name="nrel hub"):

163 W.Hub.__init__(self, name) 164 self.nacelle = None 165 self.node = HubNode() 166 self.torque = None

167

168 - def set_torque(self, torque):

169 self.torque = torque

170

171 - def set_nacelle(self, nacelle):

172 self.nacelle = nacelle

173

174 - def create_refs_from(self, wing_angle, rotational_speed):

175 """Create a first referential on the top of the fix node but 176 rotates its rotation axis and then place the referential 177 following the first blade""" 178 ref_tmp = M.ReferenceFrame("Hub tmp") 179 ref_tmp.set_relative_from(self.nacelle.refs.get("fix")) 180 ref_tmp.set_angular_velocity(rotational_speed["rad/s"], 0., 0.) 181 self.refs.add(ref_tmp, "hub tmp") 182 183 ref = M.ReferenceFrame("Hub") 184 ref.set_relative_from(ref_tmp) 185 ref.set_orientation_matrix(one=(1., 0., 0.), 186 two=(0., 187 N.cos(wing_angle["rad"]), 188 N.sin(wing_angle["rad"]))) 189 self.refs.add(ref, "hub")

190

191 - def create_nodes(self):

192 self.node.set_relative_from(self.refs.get("hub")) 193 self.nodes.add(self.node, "hub")

194

195 - def create_joints(self):

196 joint = M.RevoluteHinge("Hub hinge") 197 joint.set_relative_from(self.refs.get("hub")) 198 joint.set_node1(self.nacelle.nodes.get("fix")) 199 joint.set_node2(self.nodes.get("hub")) 200 joint.set_rotation_axis(two=(0., 1., 0.), 201 three=(1., 0., 0.)) 202 self.elts.add(joint)

203

204 - def create_load(self):

205 self.torque.set_relative_from(self.refs.get("hub")) 206 self.torque.set_direction(-1., 0., 0.) 207 self.torque.set_amplitude(0.) 208 self.torque.attach_to_node(self.nodes.get("hub")) 209 self.elts.add(self.torque)

210

211 -class PitchDriver(M.BindingsHinge):

212

213 - def do_as_init(self):

214 self._direction = N.array([ [0.], 215 [0.], 216 [1.] ])

217

218 - def set_pitch(self, blade_pitch):

219 self.pitch = blade_pitch

220

221 - def broken_update(self, simu):

222 self.mbdyn_inst.set_value(self.pitch["rad"] * self._direction)

223

224 - def update(self, simu):

225 pass

226 227

228 -class Blade(W.Blade):

229 - def __init__(self, idx=0):

230 W.Blade.__init__(self, idx) 231 self.stiffness = {} 232 self.flange_distance = 5. 233 self.pitch_driver = PitchDriver() 234 self.pitch_driver.set_pitch(self.angle["pitch"])

235

236 - def set_flange_distance(self, value):

237 self.flange_distance = value

238

239 - def create_base_ref_from(self, hub):

240 ref = M.ReferenceFrame(self.name + " base") 241 ref.set_relative_from(hub.refs.get("hub")) 242 offset_angle = self.offset_wing_angle["rad"] 243 ref.set_orientation_matrix(one=(1., 0., 0.), 244 two=(0., 245 N.cos(offset_angle), 246 N.sin(offset_angle))) 247 self.refs.add(ref, "base") 248 249 ref_flange = M.ReferenceFrame(self.name + " flange fix") 250 ref_flange.set_relative_from(ref) 251 ref_flange.set_position(0., 0., self.flange_distance) 252 self.refs.add(ref_flange, "flange fix") 253 254 #self.pitch_debug = N.pi/4. 255 self.pitch_debug = 0. 256 ref = M.ReferenceFrame(self.name + "flange") 257 ref.set_relative_from(ref_flange) 258 ref.set_orientation_matrix(one=(N.cos(self.pitch_debug), 259 N.sin(self.pitch_debug), 260 0.), 261 three=(0., 0., 1.)) 262 self.refs.add(ref, "flange")

263

264 - def create_flange_node(self):

265 node = M.StructuralNode(self.name + "flange") 266 node.set_relative_from(self.refs.get("flange")) 267 self.nodes.add(node, "flange")

268

269 - def add_bem_nodes(self):

270 """Add the StructuralNode of the bem nodes that exit 271 to the list of nodes""" 272 for bem_node in self.bem_nodes: 273 ref = M.ReferenceFrame(self.name + " " + bem_node.name) 274 ref.set_relative_from(self.refs.get("flange")) 275 ref.set_position(0., 276 0., 277 bem_node.radial_position) 278 self.refs.add(ref) 279 280 bem_node.structural.set_relative_from(ref) 281 bem_node.structural.set_comment("Node properties relative to %s" \ 282 % ref.name) 283 self.nodes.add(bem_node.structural)

284

285 - def add_bem_bodies(self):

286 for bem_node in self.bem_nodes: 287 self.elts.add(bem_node.body)

288

289 - def create_drive_hinge_from(self, hub):

290 joint = M.SphericalHinge(self.name + " flange spherical") 291 joint.set_relative_from(self.refs.get("flange")) 292 joint.set_node1(hub.nodes.get("hub")) 293 joint.set_node2(self.nodes.get("flange")) 294 self.elts.add(joint) 295 296 297 #joint = M.Prismatic() 298 #joint.set_relative_from(self.refs.get("flange")) 299 #joint.set_node1(hub.nodes.get("hub")) 300 #joint.set_node2(self.nodes.get("flange")) 301 302 joint = self.pitch_driver 303 joint.set_relative_from(self.refs.get("flange")) 304 joint.set_node1(hub.nodes.get("hub")) 305 joint.set_node2(self.nodes.get("flange")) 306 joint.set_orientation_axis(0., 0., 1.) 307 joint.set_value(self.pitch_debug) 308 309 self.elts.add(joint)

310

311 - def create_coincidence_joint(self):

312 joint = M.Coincidence(self.name) 313 joint.set_relative_from(self.bem_nodes[0].structural.ref) 314 joint.set_node1(self.nodes.get("flange")) 315 joint.set_node2(self.bem_nodes[0].structural) 316 self.elts.add(joint)

317

318 - def create_rigid_structure(self):

319 for idx in range(self.nb_bem_nodes - 1): 320 node1 = self.bem_nodes[idx] 321 node2 = self.bem_nodes[idx + 1] 322 323 joint = M.Coincidence(self.name + " " + str(idx)) 324 joint.set_relative_from(node1.structural.ref) 325 joint.set_node1(node1.structural) 326 joint.set_node2(node2.structural) 327 self.elts.add(joint)

328

329 - def create_beams(self):

330 beam = M.Beam3(self.name) 331 beam.set_relative_from(self.refs.get("root")) 332 beam.set_node1(self.nodes.get("root")) 333 beam.set_node2(self.nodes.get("mid")) 334 beam.set_node3(self.nodes.get("tip")) 335 336 beam.section1.set_orientation_matrix(M.EYE) 337 beam.section2.set_orientation_matrix(M.EYE) 338 self._set_stiffness(beam) 339 340 self.elts.add(beam, "beam")

341

342 - def _set_stiffness(self, beam):

343 """By default, the user will give the same constitutive law 344 for the upper and lower part: values are contained in the 345 'stifness' dictionary for a 6D diagonal matrix. As a an alternative, 346 he can manipulate the 'blades[i].beam.section1' instance and set what he 347 wants a law.""" 348 law = M.law.LinearElasticGeneric() 349 stiff = self.stiffness 350 law.set_stiffness(diag=(stiff["GA"], 351 stiff["GA"], 352 stiff["EA"], 353 stiff["EIx"], 354 stiff["EIy"], 355 stiff["GIz"])) 356 beam.section1.set_constitutive_law(law) 357 beam.section2.set_constitutive_law(law)

358

359 - def add_bem_forces(self):

360 """Init the coupling with the bem nodes and 361 set the bem force in the MBDyn elements.""" 362 # Set the sections for the BemNode 363 # Create the force 364 self.init_coupling() 365 for bem_node in self.bem_nodes: 366 if bem_node.force: 367 bem_node.force.set_relative_from(bem_node.structural.ref) 368 bem_node.force.attach_to_node(bem_node.structural) 369 bem_node.force.init_loading() 370 self.elts.add(bem_node.force)

371

372 - def update(self):

373 """Update""" 374 pass

375 376

377 -class Rotor(W.Rotor):

378

379 - def __init__(self, name="rotor"):

380 W.Rotor.__init__(self, name) 381 self.nacelle = None 382 383 self.rotational_speed_value = 0. 384 self.electrical_torque = 0. 385 self.electrical_power = 0. 386 self.will_save("rotational_speed_value", 387 "electrical_torque", 388 "electrical_power")

389

390 - def set_nacelle(self, nacelle):

391 self.hub.set_nacelle(nacelle)

392

393 - def create_references(self):

394 self.hub.create_refs_from(self.wing_angle, 395 self.rotational_speed) 396 for blade in self.blades: 397 blade.create_base_ref_from(self.hub)

398

399 - def create_nodes(self):

400 self.hub.create_nodes() 401 for blade in self.blades: 402 blade.create_flange_node() 403 blade.add_bem_nodes()

404

405 - def create_elements(self):

406 self.hub.create_joints() 407 self.hub.create_load() 408 for blade in self.blades: 409 blade.add_bem_bodies() 410 blade.create_drive_hinge_from(self.hub) 411 blade.create_coincidence_joint() 412 blade.add_bem_forces() 413 blade.create_rigid_structure()

414

415 - def create(self, use_beams):

416 self.hub.create(self.rotational_speed) 417 for blade in self.blades: 418 blade.create_base_ref_from(self.hub) 419 #blade.create_coincidence_joint_from(self.hub) 420 #blade.create_bodies() 421 if use_beams: 422 blade.create_beams()

423 #else: 424 # blade.create_rigid_structure() 425 #blade.create_forces() 426

427 - def save(self):

428 """Save the current status.""" 429 self.hub.node.get_matrix() 430 angv = self.hub.node.mbdyn_inst.get_angular_velocity() 431 local_angv = N.asarray(self.hub.node.matrix * angv) 432 # to activate 433 self.rotational_speed["rad/s"] = float(local_angv[0][0]) 434 435 self.rotational_speed_value = self.rotational_speed["rad/s"] 436 self.electrical_torque = self.hub.torque.amplitude 437 self.electrical_power = self.electrical_torque * \ 438 self.rotational_speed["rad/s"] 439 440 self.save_results() 441 442 for blade in self.blades: 443 blade.save()

444

445 - def update(self):

446 """Update the value at each time step, before running 447 the unsteady BEM code.""" 448 pass

449 450

451 -class WindTurbine(W.WindTurbine):

452 """The wind turbine of a model is mainly supposed to 453 describe two methods: create and update. 454 'create()' is a method called at the beginning to fill 455 the 'refs', 'nodes' and 'elts' lists with MBDyn objects. 456 'update()' is called at each time step and is supposed 457 to update all the wind turbine parameters, according to 458 the MBDyn objects status. Then the unsteady BEM code 459 will be applied. This method must not be 460 confused with the update method called by WraptMBDyn, 461 that is called to update the MBDyn objects before a new 462 time step. This step actually occurs after the unsteady BEM. 463 """ 464

465 - def __init__(self, name="tjaerborg", use_beams=True, 466 only_mbdyn_rotor=False):

467 W.WindTurbine.__init__(self, name) 468 self.use_beams = use_beams 469 if only_mbdyn_rotor: 470 self.create = self._create_only_rotor

471

472 - def _create_only_rotor(self):

473 """Create only the rotor with MBDyn components. 474 The main bearing fix node is clamped, only the tower 475 and nacelle referentials are created. A very useful 476 function to test only the rotor or develop 477 the nacelle and tower deformations step by step.""" 478 self.tower.create_references() 479 self.nacelle.set_tower(self.tower) 480 self.nacelle.create_references_from(self.rotor.wing_angle) 481 482 node = M.StructuralNode("Nacelle clamp") 483 node.set_relative_from(self.nacelle.refs.get("fix")) 484 self.nacelle.nodes.add(node, "fix") 485 486 joint = M.Clamp("Nacelle") 487 joint.set_relative_from(self.nacelle.refs.get("fix")) 488 joint.attach_to_node(node) 489 self.nacelle.elts.add(joint) 490 491 self.rotor.set_nacelle(self.nacelle) 492 self.rotor.create_references() 493 self.rotor.create_nodes() 494 self.rotor.create_elements()

495

496 - def create(self):

497 """Suppose to create the complete turbine. 498 Still in development""" 499 if self.has_tower and self.has_nacelle: 500 self.nacelle.set_tower(self.tower) 501 if self.has_nacelle and self.has_rotor: 502 self.rotor.set_nacelle(self.nacelle)

503

504 - def update(self):

505 """Called by the simulation at each MBDyn 506 time step""" 507 self.rotor.update()

508

Source Code for Module windSimSuite.models.tjaerborg