windSimSuite.bem

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 """A blade element momentum node. The main class 23 is described in L{BemNode}. 24 """ 25 import numpy as N 26 27 from mbdyn.nodes import StructuralNode 28 from mbdyn.elements import Body 29 from mbdyn.elts.force import BindingsForce 30 31 from windSimSuite.common import BasicObject 32

33 -def get_force_old(fkey, sec1, sec2):

34 """Integrate the force between two sections according 35 to the force key, fkey, 'normal' or 'tangential'""" 36 delta_r = sec2.radial_position - sec1.radial_position 37 acoeff = (sec2.press[fkey] - sec1.press[fkey]) / delta_r 38 bcoeff = (sec1.press[fkey] * sec2.radial_position - \ 39 sec2.press[fkey] * sec1.radial_position) / delta_r 40 square_delta_r = sec2.radial_position**2 - \ 41 sec1.radial_position**2 42 return 0.5 * acoeff * square_delta_r + bcoeff * delta_r

43 44

45 -def get_force(fkey, sec1, sec2):

46 """Return the torque contribution between 47 the two sections""" 48 delta_r = sec2.radial_position - sec1.radial_position 49 acoeff = (sec2.press[fkey] - sec1.press[fkey]) / delta_r 50 bcoeff = (sec1.press[fkey] * sec2.radial_position - \ 51 sec2.press[fkey] * sec1.radial_position) / delta_r 52 square_delta_r = sec2.radial_position**2 - \ 53 sec1.radial_position**2 54 threep_delta_r = sec2.radial_position**3 - \ 55 sec1.radial_position**3 56 return 1./3. * acoeff * threep_delta_r + \ 57 0.5 * bcoeff * square_delta_r

58 59

60 -class _ForceArea:

61 """The resulting force from 2 aerodynamics sections. 62 """ 63

64 - def __init__(self, lower_section, upper_section):

65 self.lower_section = lower_section 66 self.upper_section = upper_section

67

68 - def get_force(self, fkey):

69 """Return the force value according to the C{fkey}, 70 normal or tangential.""" 71 return get_force(fkey, self.lower_section, self.upper_section)

72 73

74 -class _CalculationSection:

75 """A section only used for calculating the load in one area 76 according to its radial position and two others 77 aerodynamic sections.""" 78

79 - def __init__(self, radial_position):

80 self.radial_position = radial_position 81 self.press = {"normal" : 0., 82 "tangential" : 0.}

83

84 - def fill(self, fkey, sec1, sec2):

85 """Fill the pressure from the two sections according 86 to 'fkey'""" 87 slope = (sec2.press[fkey] - sec1.press[fkey]) / \ 88 (sec2.radial_position - sec1.radial_position) 89 deltar = self.radial_position - sec1.radial_position 90 self.press[fkey] = slope * deltar + sec1.press[fkey]

91 92

93 -class _BoundaryForceArea:

94 """The top class for calculating the force area 95 at the boundaries of the L{BemNode} range of action 96 """ 97

98 - def __init__(self, radius, lower_section, upper_section):

99 self.section = _CalculationSection(radius) 100 self.lower_section = lower_section 101 self.upper_section = upper_section

102 103

104 -class _LowerForceArea(_BoundaryForceArea):

105 """The lower boundary of the L{BemNode} range of action, 106 closer to the root. 107 """ 108

109 - def get_force(self, fkey):

110 """Return the force contribution of the area""" 111 self.section.fill(fkey, self.lower_section, self.upper_section) 112 return get_force(fkey, self.section, self.upper_section)

113 114

115 -class _UpperForceArea(_BoundaryForceArea):

116 """The upper boundary of the L{BemNode} range of action, 117 closer to the tip. 118 """ 119

120 - def get_force(self, fkey):

121 """Return the force contribution of the area""" 122 self.section.fill(fkey, self.lower_section, self.upper_section) 123 return get_force(fkey, self.lower_section, self.section)

124 125

126 -class BemForce(BindingsForce):

127 """A custom bindings force that finds its value 128 from the aerodynamic sections, organised 129 into force areas. 130 """ 131

132 - def __init__(self, name):

133 BindingsForce.__init__(self, name) 134 self.areas = [] 135 self.forces = {"tangential" : 0., 136 "normal" : 0.} 137 self.torque = {} 138 self.tangential_force = self.forces["tangential"] 139 self.normal_force = self.forces["normal"] 140 141 self.xcontrib = 0. 142 self.ycontrib = 0. 143 self._load = N.zeros((3, 1)) 144 145 self.will_save("tangential_force", 146 "normal_force")

147 148 # To be removed as well

149 - def init_results(self):

150 """Init the results""" 151 self.common_init_results()

152

153 - def add_area(self, area):

154 """Add a force area to the list""" 155 self.areas.append(area)

156

157 - def add_lower_area(self, area):

158 """Add the lower force area to the list""" 159 self.areas.insert(0, area)

160

161 - def add_upper_area(self, area):

162 """Add the upper force area to the list""" 163 self.areas.append(area)

164

165 - def init_loading(self):

166 """Initialize the loading for the MBDyn input file""" 167 self._get_loads() 168 norm = N.sqrt(self.tangential_force ** 2 + \ 169 self.normal_force ** 2) 170 self.xcontrib = self.normal_force / norm 171 self.ycontrib = - self.tangential_force / norm 172 173 self.set_direction(self.xcontrib, 174 self.ycontrib, 175 0.) 176 self.set_amplitude(0.)

177 # Finally not working during the assembly 178 #self.set_amplitude(norm) 179

180 - def _get_loads(self):

181 """Fill the C{_load} attribute. 182 The normal force is along M{x}, the tangential force 183 along M{y}.""" 184 for key in ["normal", "tangential"]: 185 #self.forces[key] = 0. 186 self.torque[key] = 0. 187 for area in self.areas: 188 #self.forces[key] += area.get_force(key) 189 self.torque[key] += area.get_force(key) 190 self.forces[key] = self.torque[key] / self.node.radial_position 191 192 193 self.tangential_force = self.forces["tangential"] 194 self.normal_force = self.forces["normal"] 195 self._load[0][0] = self.normal_force 196 self._load[1][0] = - self.tangential_force

197

198 - def update(self, simu):

199 """Before a new BEM step, the aerodynamics forces are calculated 200 and are applied on the MBDyn force""" 201 self._get_loads() 202 rot_matrix = N.asmatrix(self.node.mbdyn_inst.get_rotation_matrix()) 203 self.mbdyn_inst.set_value(N.asarray(rot_matrix * self._load))

204 205

206 -class BemNode(BasicObject):

207 """A node used for applying the Blade Element Momentum method. 208 This node gathers a C{StructuralNode} and a C{Body} from 209 the MBdyn package. It could have also 210 a C{BemForce} if its radial position contains aerodynamic sections 211 on which pressures are calculated. The C{BemNode} are defined by the user 212 along the blade and are not twisted, it means that all the 213 C{StructuralNode} have initially the same rotation matrix. 214 However the C{BemNode} follows the blade pitching. 215 Moreover the rotation matrix of each node could change 216 if the blade starts to deform. 217 """ 218

219 - def __init__(self, name="bem node"):

220 BasicObject.__init__(self, name) 221 self.mass = 0. 222 self.radial_position = 0. 223 224 self.lower_radius = None 225 self.upper_radius = None 226 self.force = None 227 228 self.structural = StructuralNode(self.name) 229 self.structural.set_type("dynamic") 230 231 self.body = Body(self.name + "body") 232 self.body.attach_to_node(self.structural) 233 234 self.own_para_names += ["mass", 235 "radial_position"]

236

237 - def set_mass(self, mass):

238 """Applied the mass on the body""" 239 self.mass = mass 240 self.body.set_mass(mass)

241

242 - def set_radius(self, radius):

243 """Set the radial position of the BEM node""" 244 self.radial_position = radius 245 self.structural.radial_position = radius

246

247 - def create_force(self):

248 """Create the L{BemForce}. 249 Supposed to do more but there is still a problem 250 between a general interface and a specific model 251 with MBDyn objects""" 252 self.force = BemForce(self.name + " force")

253

254 - def create_force_areas_from(self, sections):

255 """Create the force areas from the sections. Only 256 used when a single bem node is used on the whole 257 blade""" 258 for idx in range(len(sections) - 1): 259 area = _ForceArea(sections[idx], 260 sections[idx + 1]) 261 self.force.add_area(area)

262

263 - def set_boundaries(self, lower_node_radius, upper_node_radius):

264 """Find the boundaries from the two other bem node limits""" 265 deltar = (self.radial_position - lower_node_radius) / 2. 266 self.lower_radius = lower_node_radius + deltar 267 268 deltar = (upper_node_radius - self.radial_position) / 2. 269 self.upper_radius = self.radial_position + deltar

270

271 - def find_sections(self, section_radius, sections):

272 """Find the relevant sections for the bem node""" 273 cond_lower = self.lower_radius <= section_radius 274 cond_upper = section_radius <= self.upper_radius 275 276 section_idx = list(N.where(cond_lower * cond_upper)[0]) 277 278 if len(section_idx) > 0: 279 self.create_force() 280 for idx in range(len(section_idx) - 1): 281 area = _ForceArea(sections[idx], 282 sections[idx + 1]) 283 self.force.add_area(area) 284 self._handle_boundary_areas(section_idx[0], 285 section_idx[-1], 286 sections)

287

288 - def _handle_boundary_areas(self, first_section_idx, 289 last_section_idx, sections):

290 """Create the boundary force areas for a C{BemNode}. 291 A boundary area exists only if there is not already 292 an aerodynamics section on the C{BemNode} border and if 293 an extra aerodynamics section exists above the border.""" 294 first_section = sections[first_section_idx] 295 if first_section.radial_position != self.lower_radius: 296 if (first_section_idx - 1) != 0: 297 area = _LowerForceArea(self.lower_radius, 298 sections[first_section_idx - 1], 299 first_section) 300 self.force.add_lower_area(area) 301 302 last_section = sections[last_section_idx] 303 if last_section.radial_position != self.upper_radius: 304 if (last_section_idx + 1) != (len(sections) - 1): 305 area = _UpperForceArea(self.upper_radius, 306 last_section, 307 sections[last_section_idx + 1]) 308 self.force.add_upper_area(area)

309

310 - def save(self):

311 """Save the bem node data""" 312 # To remove or not? 313 # could be interesting to see the amount of forces carried 314 # by each BEM node 315 if self.force: 316 self.force.save()

317

318 - def update(self):

319 """Update the forces before running MBDyn""" 320 if self.force: 321 self.force.update(None)

322

Source Code for Module windSimSuite.bem_node