windSimSuite.models.nrel

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 NREL turbine model. Still in heavy development that's why 23 it has an incomplete documentation.""" 24 import numpy as N 25 import mbdyn as M 26 import windSimSuite as W 27 28

29 -class Foundation:

30 """The wind turbine foundation. 31 """ 32

33 - def __init__(self, name="foundation"):

34 self.name = name 35 self.refs = W.List() 36 self.ref = None 37 38 self.phix = W.Angle(0., "deg") 39 self.phiy = W.Angle(0., "deg") 40 self.ground_level = 0.

41

42 - def set_lateral_misalignement(self, value, unit_key="deg"):

43 """The value is supposed to be entered in degree by default. 44 1. rotate laterally round global X-axis. 45 """ 46 self.phix[unit_key] = value

47

48 - def set_longitudinal_misalignement(self, value, unit_key="deg"):

49 """The value is supposed to be entered in degree by default. 50 2. rotate longitudinally round the (inclined) Y-axis. 51 """ 52 self.phiy[unit_key] = value

53

54 - def set_ground_level(self, value):

55 self.ground_level = value

56

57 - def create_references(self):

58 refx = M.ReferenceFrame("x " + self.name) 59 refx.set_orientation_matrix(one = (1., 0., 0.), 60 three = (0., 61 -N.sin(self.phix["rad"]), 62 N.cos(self.phix["rad"]))) 63 self.refs.add(refx) 64 65 refy = M.ReferenceFrame("y " + self.name) 66 refy.set_relative_from(refx) 67 refy.set_orientation_matrix(two = (0., 1., 0.), 68 three = (N.sin(self.phiy["rad"]), 69 0., 70 N.cos(self.phiy["rad"]))) 71 self.refs.add(refy) 72 73 ref = M.ReferenceFrame(self.name) 74 ref.set_relative_from(refy) 75 ref.set_position(0., 0., self.ground_level) 76 self.ref = ref 77 self.refs.add(ref)

78

79 - def create(self):

80 """Called by the tower""" 81 self.create_references()

82 83

84 -class Tower(W.Tower):

85 """Interface for building the 86 the tower of the wind turbine. 87 The needed parameters are: 88 89 - height 90 - mass 91 - stiffness as a dictionary 92 """ 93

94 - def __init__(self, name="tower"):

95 W.Tower.__init__(self, name) 96 self.foundation = Foundation() 97 98 self.stiffness = {} 99 self.height = 0. 100 self.mass = None 101 102 self.height_substructure = 0.

103

104 - def set_height(self, height):

105 self.height = height

106

107 - def create_references(self):

108 for ref in self.foundation.refs: 109 self.refs.add(ref) 110 111 base_ref = M.ReferenceFrame(self.name + " base") 112 base_ref.set_relative_from(self.foundation.ref) 113 base_ref.set_position(0., 0., self.height_substructure) 114 self.refs.add(base_ref, "base") 115 116 top_ref = M.ReferenceFrame(self.name + " mid") 117 top_ref.set_relative_from(base_ref) 118 top_ref.set_position(0., 0., 0.5 * self.height) 119 self.refs.add(top_ref, "mid") 120 121 top_ref = M.ReferenceFrame(self.name + " top") 122 top_ref.set_relative_from(base_ref) 123 top_ref.set_position(0., 0., self.height) 124 self.refs.add(top_ref, "top")

125

126 - def create_nodes(self):

127 nodes_desc = [ 128 ("base", "static", 0.), 129 ("mid", "dynamic", 0.5), 130 ("top", "dynamic", 1.) 131 ] 132 for name, ntype, factor in nodes_desc: 133 node = M.StructuralNode(self.name + " " + name) 134 node.set_relative_from(self.refs.get("base")) 135 com = "Node properties relative to %s" % \ 136 self.refs.get("base").name 137 node.set_comment(com) 138 node.set_type(ntype) 139 node.set_position(0., 0., factor * self.height) 140 self.nodes.add(node, name)

141

142 - def create_clamped_joint(self):

143 joint = M.Clamp(self.name + " " + "clamp") 144 joint.attach_to_node(self.nodes.get("base")) 145 joint.set_position("node") 146 joint.set_orientation_matrix("node") 147 self.elts.add(joint)

148

149 - def create_bodies(self):

150 Jmid = 1./12. * (0.5*self.mass) * (0.5*self.height)**2 151 Jtop = 1./3. * (0.25*self.mass) * (0.25*self.height)**2 152 bodies_desc = [ \ 153 # key, mass ratio, offset ratio, inertia 154 ("mid", 0.5, 0., Jmid), 155 ("top", 0.25, -0.125, Jtop) 156 ] 157 for key, mass_ratio, offset_ratio, inertia in bodies_desc: 158 body = M.Body(self.name + " body " + key) 159 body.attach_to_node(self.nodes.get(key)) 160 body.set_mass(mass_ratio * self.mass) 161 body.set_center_of_mass(0., 162 0., 163 offset_ratio * self.height) 164 body.set_inertia_matrix(diag=(inertia, 165 inertia, 166 1.)) 167 self.elts.add(body, "body " + key)

168

169 - def create_rigid_structure(self):

170 for key1, key2 in [("base", "mid"), ("mid", "top")]: 171 joint = M.Coincidence(self.name + key2) 172 joint.set_relative_from(self.refs.get(key2)) 173 joint.set_node1(self.nodes.get(key1)) 174 joint.set_node2(self.nodes.get(key2)) 175 self.elts.add(joint, key2)

176

177 - def create_beam(self):

178 beam = M.Beam3(self.name + " beam") 179 180 beam.set_node1(self.nodes.get("base")) 181 beam.set_node2(self.nodes.get("mid")) 182 beam.set_node3(self.nodes.get("top")) 183 184 beam.section1.set_orientation_matrix(M.EYE) 185 beam.section2.set_orientation_matrix(M.EYE) 186 self._set_stiffness(beam) 187 188 self.elts.add(beam, "beam")

189

190 - def _set_stiffness(self, beam):

191 """By default, the user will give the same constitutive law 192 for the upper and lower part: values are contained in the 193 'stiffness' dictionary for a 6D diagonal matrix. As a an alternative, 194 he can manipulate the 'tower.beam.section1' instance and set what he 195 wants as a law.""" 196 stif = self.stiffness 197 law = M.law.LinearElasticGeneric() 198 law.set_stiffness(diag=(stif["GA"], stif["GA"], stif["EA"], 199 stif["EIx"], stif["EIy"], stif["GIz"])) 200 beam.section1.set_constitutive_law(law) 201 beam.section2.set_constitutive_law(law)

202

203 - def create(self, use_beams):

204 """Call by the wind turbine tower when the simulation 205 is initialized""" 206 self.foundation.create() 207 self.create_references() 208 self.create_nodes() 209 self.create_clamped_joint() 210 #self.create_bodies() 211 if use_beams: 212 self.create_beam() 213 else: 214 self.create_rigid_structure()

215 216

217 -class MainBearing:

218 """The main bearing of the wind turbine, 219 part of the nacelle. 220 """ 221

222 - def __init__(self):

223 self.name = "main bearing" 224 self.refs = W.List() 225 self.nodes = W.List() 226 self.elts = W.List() 227 228 self.distance = 0. 229 self.nacelle = None

230

231 - def set_distance(self, value):

232 """The distance between the nacelle referential 233 and the main bearing center""" 234 self.distance = value

235

236 - def set_nacelle(self, nacelle):

237 self.nacelle = nacelle

238

239 - def create_references(self):

240 ref_fix = M.Reference(self.name + " fix") 241 ref_fix.set_relative_from(self.nacelle.refs.get("tilt")) 242 ref_fix.set_position(self.distance, 0., 0.) 243 self.refs.add(ref_fix, "fix")

244

245 - def create_nodes(self):

246 #for name, ntype in [("fix", "dynamic"), 247 # ("fix rotor", "static")]: 248 for name, ntype in [("fix", "dynamic")]: 249 node = M.StructuralNode(self.name + " " + name) 250 node.set_type(ntype) 251 com = "Node properties relative to %s" % \ 252 self.refs.get("fix").name 253 node.set_comment(com) 254 node.set_relative_from(self.refs.get("fix")) 255 self.nodes.add(node, name)

256 257 #self.nodes.get("fix").set_rotation_matrix(M.EYE) 258 # A rotation from X to Y, to make the rotor turns clock wise 259 #self.nodes.get("fix rotor").set_rotation_matrix(two=(0., 1., 0.), 260 # three=(-1., 0., 0.)) 261

262 - def create_joints(self):

263 """Clamp the nacelle node with the fix node 264 of the bearing a Coincidence joint. 265 """ 266 267 joint = M.Coincidence("Main bearing coincidence") 268 joint.set_relative_from(self.refs.get("fix")) 269 joint.set_node1(self.nacelle.nodes.get("nacelle")) 270 joint.set_node2(self.nodes.get("fix")) 271 self.elts.add(joint)

272 273

274 -class Nacelle(W.Nacelle):

275 - def __init__(self, name="nacelle"):

276 W.Nacelle.__init__(self, name) 277 self.main_bearing = MainBearing() 278 279 # Element, like joint are shared between the two 280 self.main_bearing.set_nacelle(self) 281 282 self.tower = None 283 self.ref = None 284 285 self.distance_to_rotor_axis = 0. 286 self.mass = None 287 self.gravity_center = (0., 0., 0.) 288 self.inertia = {"roll" : None, 289 "node" : None, 290 "yaw" : None} 291 self._yaw = W.Angle(0., "deg") 292 self._tilt = W.Angle(0., "deg")

293

294 - def set_distance_to_rotor_axis(self, value):

295 self.distance_to_rotor_axis = value

296

297 - def set_mass(self, value):

298 """Set the mass of the nacelle in kg""" 299 self.mass = value

300

301 - def set_gravity_center(self, x, y, z):

302 """Set the position of the center of gravity 303 relative to the node position, in [m]""" 304 self.gravity_center = (x, y, z)

305

306 - def set_inertia(self, Jroll, Jnode, Jyaw):

307 """Set the inertia round the axis of nacelle referential 308 referred to center of gravity, in [kg m²]""" 309 self.inertia = {"roll" : Jroll, 310 "node" : Jnode, 311 "yaw" : Jyaw}

312

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

314 self._yaw[unit_key] = value

315

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

317 self._tilt[unit_key] = value

318

319 - def set_tower(self, tower):

320 self.tower = tower

321

322 - def create_references(self):

323 yaw_ref = M.Reference(self.name + " yaw") 324 yaw_ref.set_relative_from(self.tower.refs.get("top")) 325 # yaw + pi for this angle because the x axe will now 326 # point upstream 327 yaw_ref.set_orientation_matrix(three=(0., 0., 1.), 328 one=(N.cos(self._yaw["rad"] + N.pi), 329 N.sin(self._yaw["rad"] + N.pi), 330 0.)) 331 self.refs.add(yaw_ref, "yaw") 332 333 ref = M.Reference(self.name) 334 ref.set_relative_from(yaw_ref) 335 ref.set_position(0., 0., self.distance_to_rotor_axis) 336 self.ref = ref 337 self.refs.add(ref) 338 339 tilt_ref = M.Reference(self.name + " tilt") 340 tilt_ref.set_relative_from(ref) 341 # Minus the tilt angle because a positive tilt is supposed 342 # to point the node upstream 343 tilt_ref.set_orientation_matrix(two=(0., 1., 0.), 344 one=(N.cos(-self._tilt["rad"]), 345 0., 346 -N.sin(-self._tilt["rad"]))) 347 self.refs.add(tilt_ref, "tilt")

348 349

350 - def create_node(self):

351 node = M.StructuralNode(self.name) 352 node.set_type("dynamic") 353 node.set_comment("Node properties relative to %s" % self.ref.name) 354 node.set_relative_from(self.ref) 355 self.nodes.add(node, "nacelle") 356 for node in self.main_bearing.nodes: 357 self.nodes.add(node)

358

359 - def create_elements(self):

360 """Create the nacelle elements: a body and 2 joints""" 361 self._create_body() 362 self._create_joints()

363

364 - def _create_body(self):

365 body = M.Body(self.name + " body") 366 body.set_mass(self.mass) 367 body.set_inertia_matrix(diag=(self.inertia["roll"], 368 self.inertia["node"], 369 self.inertia["yaw"])) 370 body.set_center_of_mass(*self.gravity_center) 371 body.attach_to_node(self.nodes.get("nacelle")) 372 self.elts.add(body, "body")

373

374 - def _create_joints(self):

375 """Clamp the tower top with the nacelle node""" 376 joint = M.Coincidence("nacelle") 377 joint.set_relative_from(self.ref) 378 joint.set_node1(self.tower.nodes.get("top")) 379 joint.set_node2(self.nodes.get("nacelle")) 380 self.elts.add(joint)

381

382 - def create(self, use_beams):

383 self.create_references() 384 self.main_bearing.create_references() 385 for ref in self.main_bearing.refs: 386 self.refs.add(ref) 387 self.main_bearing.create_nodes() 388 self.create_node() 389 #self.create_elements() 390 self._create_joints() 391 self.main_bearing.create_joints()

392 393

394 -class FakeTorque(M.BindingsCouple):

395

396 - def do_as_init(self):

397 self._value = N.array([ [self.norm, 398 0., 399 0.] ])

400

401 - def set_norm(self, norm):

402 self.norm = norm

403

404 - def set_node(self, node):

405 self.node = node

406

407 - def update(self, simu):

408 self.mbdyn_inst.set_value(self._value)

409 410

411 -class Hub(W.Hub):

412

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

414 W.Hub.__init__(self, name) 415 self.nacelle = None 416 self.use_fake_load = False

417

418 - def set_nacelle(self, nacelle):

419 self.nacelle = nacelle

420

421 - def create_refs(self):

422 ref_tmp = M.ReferenceFrame("Hub tmp") 423 ref_tmp.set_relative_from(self.nacelle.main_bearing.refs.get("fix")) 424 self.refs.add(ref_tmp, "hub tmp") 425 426 ref = M.ReferenceFrame("Hub") 427 ref.set_relative_from(ref_tmp) 428 ref.set_orientation_matrix(one=(-1, 0., 0.), 429 three=(0., 0., 1.)) 430 self.refs.add(ref, "hub")

431

432 - def create_nodes(self):

433 node = M.StructuralNode("Hub") 434 node.set_relative_from(self.refs.get("hub")) 435 self.nodes.add(node, "hub")

436

437 - def create_joints(self):

438 joint = M.RevoluteHinge("Hub hinge") 439 joint.set_relative_from(self.refs.get("hub")) 440 joint.set_node1(self.nacelle.main_bearing.nodes.get("fix")) 441 joint.set_node2(self.nodes.get("hub")) 442 joint.set_rotation_axis(two=(0., 1., 0.), 443 three=(1., 0., 0.)) 444 self.elts.add(joint)

445

446 - def create(self):

447 self.create_refs() 448 self.create_nodes() 449 self.create_joints()

450

451 - def create_fake_load(self):

452 torque = FakeTorque() 453 torque.set_relative_from(self.refs.get("tip")) 454 torque.set_direction(0., -1., 0.) 455 torque.set_node(self.nodes.get("tip")) 456 self.elts.add(torque)

457

458 - def create_forces(self):

459 if self.use_fake_load: 460 self.create_fake_load()

461 462 463

464 -class FakeLoad(M.BindingsForce):

465

466 - def do_as_init(self):

467 self.force = N.array([ [0.], 468 [-10.], 469 [0.] ])

470

471 - def set_node(self, node):

472 self.node = node

473

474 - def update(self, simu):

475 rot_matrix = N.asmatrix(self.node.mbdyn_inst.get_rotation_matrix()) 476 self.mbdyn_inst.set_value(N.asarray(rot_matrix * self.force))

477 478

479 -class Blade(W.Blade):

480 - def __init__(self, idx=0, use_fake_load=False):

481 W.Blade.__init__(self, idx) 482 self.mass = 1. 483 self.stiffness = {} 484 self.flange_distance = 5. 485 self.use_fake_load = use_fake_load

486

487 - def set_mass(self, mass):

488 self.mass = mass

489

490 - def will_use_fake_load(self, boolean):

491 self.use_fake_load = boolean

492

493 - def create_base_ref_from(self, hub):

494 ref = M.ReferenceFrame(self.name + "base") 495 ref.set_relative_from(hub.refs.get("hub")) 496 # To rotate clock wise by facing the rotor plane, 497 # the rotation angle needs to be negative in the hub referential 498 offset_angle = -self.angle["wing"]["rad"] 499 ref.set_position(0., 500 N.sin(offset_angle) * self.flange_distance, 501 N.cos(offset_angle) * self.flange_distance) 502 ref.set_orientation_matrix(one=(1., 0., 0.), 503 three=(0., 504 N.sin(offset_angle), 505 N.cos(offset_angle))) 506 self.refs.add(ref, "base")

507

508 - def create_nodes(self):

509 nodes_desc = [ 510 ("root", "static", 0.), 511 ("mid", "dynamic", 0.5), 512 ("tip", "dynamic", 1.) 513 ] 514 beam_length = self.radius - self.flange_distance 515 for key, ntype, rfactor in nodes_desc: 516 ref = M.ReferenceFrame(self.name + " " + key) 517 ref.set_relative_from(self.refs.get("base")) 518 ref.set_position(0., 519 0., 520 rfactor * beam_length) 521 self.refs.add(ref, key) 522 node = M.StructuralNode(self.name + " " + key) 523 node.set_relative_from(ref) 524 node.set_type(ntype) 525 node.set_comment("Node properties relative to %s" % ref.name) 526 self.nodes.add(node, key)

527

528 - def create_coincidence_joint_from(self, hub):

529 joint = M.Coincidence(self.name) 530 joint.set_relative_from(self.refs.get("root")) 531 joint.set_node1(hub.nodes.get("hub")) 532 joint.set_node2(self.nodes.get("root")) 533 self.elts.add(joint)

534

535 - def create_bodies(self):

536 for key in ["mid", "tip"]: 537 body = M.Body(self.name + " " + key) 538 body.attach_to_node(self.nodes.get(key)) 539 body.set_mass(self.mass) 540 #body.set_inertia_matrix(diag=(1., 1., 1.)) 541 self.elts.add(body, key)

542

543 - def create_rigid_structure(self):

544 for key1, key2 in [("root", "mid"), ("mid", "tip")]: 545 joint = M.Coincidence(self.name + key2) 546 joint.set_relative_from(self.refs.get(key2)) 547 joint.set_node1(self.nodes.get(key1)) 548 joint.set_node2(self.nodes.get(key2)) 549 self.elts.add(joint, key2)

550

551 - def create_beams(self):

552 beam = M.Beam3(self.name) 553 beam.set_relative_from(self.refs.get("root")) 554 beam.set_node1(self.nodes.get("root")) 555 beam.set_node2(self.nodes.get("mid")) 556 beam.set_node3(self.nodes.get("tip")) 557 558 beam.section1.set_orientation_matrix(M.EYE) 559 beam.section2.set_orientation_matrix(M.EYE) 560 self._set_stiffness(beam) 561 562 self.elts.add(beam, "beam")

563

564 - def _set_stiffness(self, beam):

565 """By default, the user will give the same constitutive law 566 for the upper and lower part: values are contained in the 567 'stifness' dictionary for a 6D diagonal matrix. As a an alternative, 568 he can manipulate the 'blades[i].beam.section1' instance and set what he 569 wants a law.""" 570 law = M.law.LinearElasticGeneric() 571 stiff = self.stiffness 572 law.set_stiffness(diag=(stiff["GA"], 573 stiff["GA"], 574 stiff["EA"], 575 stiff["EIx"], 576 stiff["EIy"], 577 stiff["GIz"])) 578 beam.section1.set_constitutive_law(law) 579 beam.section2.set_constitutive_law(law)

580

581 - def create_fake_load(self):

582 force = FakeLoad() 583 force.set_relative_from(self.refs.get("tip")) 584 force.set_direction(0., -1., 0.) 585 force.set_node(self.nodes.get("tip")) 586 self.elts.add(force)

587

588 - def create_forces(self):

589 if self.use_fake_load: 590 self.create_fake_load()

591 592 593

594 -def create_blades():

595 blades = [] 596 for idx in range(3): 597 blade = Blade(idx + 1) 598 blade.set_offset_angle(idx * 2.*N.pi/3., "rad") 599 #blade.set_offset_angle(idx * N.pi, "rad") 600 blades.append(blade) 601 return blades

602

603 -class Rotor(W.Rotor):

604 """The NREL turbine is made of 3 blades 605 """ 606

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

608 W.Rotor.__init__(self, name) 609 hub = Hub(name + " hub") 610 self.set_hub(hub) 611 self.nacelle = None

612

613 - def set_nacelle(self, nacelle):

614 self.hub.set_nacelle(nacelle)

615

616 - def create(self, use_beams):

617 self.hub.create() 618 for blade in self.blades: 619 blade.create_base_ref_from(self.hub) 620 blade.create_nodes() 621 blade.create_bodies() 622 if use_beams: 623 blade.create_beams() 624 else: 625 blade.create_rigid_structure() 626 blade.create_coincidence_joint_from(self.hub) 627 blade.create_forces()

628 #self.create_refs() 629 #self.create_nodes() 630 #self.create_bodies() 631 #self.create_joints() 632 633

634 -class WindTurbine(W.WindTurbine):

635 - def __init__(self, name="nrel turbine", use_beams=True, 636 only_mbdyn_rotor=False):

637 W.WindTurbine.__init__(self) 638 self.use_beams = use_beams 639 if only_mbdyn_rotor: 640 self.create = self._create_only_rotor

641

642 - def _create_only_rotor(self):

643 """Create only the rotor with MBDyn components. 644 The main bearing fix node is clamped, only the tower 645 and nacelle referentials are created. A very useful 646 function to test only the rotor or develop 647 the nacelle and tower deformations step by step.""" 648 self.tower.foundation.create() 649 self.tower.create_references() 650 self.nacelle.set_tower(self.tower) 651 self.nacelle.create_references() 652 self.nacelle.main_bearing.create_references() 653 for ref in self.nacelle.main_bearing.refs: 654 self.nacelle.refs.add(ref) 655 656 node = M.StructuralNode("Main bearing clamp") 657 node.set_relative_from(self.nacelle.main_bearing.refs.get("fix")) 658 self.nacelle.main_bearing.nodes.add(node, "fix") 659 self.nacelle.nodes.add(node, "fix") 660 661 joint = M.Clamp("Main bearing") 662 joint.attach_to_node(node) 663 self.nacelle.elts.add(joint) 664 665 self.rotor.set_nacelle(self.nacelle) 666 self.rotor.create(self.use_beams)

667

668 - def create(self):

669 if self.has_tower and self.has_nacelle: 670 self.nacelle.set_tower(self.tower) 671 if self.has_nacelle and self.has_rotor: 672 self.rotor.set_nacelle(self.nacelle) 673 674 for compo_name in ["tower", "nacelle", "rotor"]: 675 if getattr(self, "has_" + compo_name): 676 compo = getattr(self, compo_name) 677 compo.create(self.use_beams)

678

Source Code for Module windSimSuite.models.nrel