Source Code for Module windSimSuite.sections

  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 aerodynamics sections definition for the unsteady BEM code. 
 23  Each section will also need an L{Airfoil} definition that will 
 24  be made of L{Curve} instances. 
 25  """ 
 26  import math as M 
 27  import numpy as N 
 28   
 29  from windSimSuite.common import BasicObject, Angle, get_norm 
 30   
 31   
 32 -class Curve: 
 33      """Made of a list of angles and a list of coefficients.  
 34      The coefficient type will be either lift or drag. But 
 35      for each type  two different situations can exist: 
 36      fully attached or fully separated. 
 37      """ 
 38   
 39 -    def __init__(self, name="curve"): 
 40          self.name = name 
 41          self.angles = [] 
 42          self.coeffs = [] 
 43           
 44          self.array_angle = None 
 45     
 46 -    def _get_angle_inst(self, angle): 
 47          """Return the angle instance""" 
 48          if isinstance(angle, Angle): 
 49              return angle 
 50          elif type(angle) == tuple and len(angle) == 2: 
 51              return Angle(angle[0], angle[1]) 
 52          else: 
 53              return Angle(angle[0]) 
 54       
 55 -    def add_point(self, angle, coeff_value): 
 56          """Add a point to the curve. That is to say 
 57          an angle and its corresponding coefficient value""" 
 58          angle = self._get_angle_inst(angle) 
 59          self.angles.append(angle["rad"]) 
 60          self.coeffs.append(coeff_value) 
 61   
 62 -    def init(self): 
 63          """Initialize the process for getting 
 64          back a coefficient from an angle.""" 
 65          self.array_angle = N.asarray(self.angles) 
 66   
 67 -    def get_coefficient_from(self, attack_angle): 
 68          """Return a coefficient from a given angle of attack. 
 69           
 70          @type attack_angle: a float 
 71          @param attack_angle: the angle in radians. 
 72          """ 
 73          attack_indices = N.where(self.array_angle < attack_angle["rad"]) 
 74          aid = int(attack_indices[0][-1]) 
 75          slope = (self.coeffs[aid + 1] - self.coeffs[aid]) / \ 
 76                  (self.angles[aid + 1] - self.angles[aid]) 
 77          angle_diff = attack_angle["rad"] - self.angles[aid] 
 78          return slope * angle_diff + self.coeffs[aid] 
 79   
 80   
 81 -class Airfoil: 
 82      """An airfoil definition. It can be defined in two different manners. 
 83      In the first case, it can handle lift and drag coefficients according 
 84      to the angle of attack. In the second case, it can handle   
 85      lift fully attached, lift fully separated and drag coefficients 
 86      according to the angle of attack. The difference between the two 
 87      definitions will be set in the L{set_lift} method. 
 88      """ 
 89 -    def __init__(self, name="airfoil"): 
 90          self.name = name 
 91          self.curves = {} 
 92          self.coeffs = {} 
 93   
 94 -    def set_name(self, name): 
 95          """Set the airfoil name""" 
 96          self.name = name 
 97   
 98 -    def set_lift(self, lift_curve, separated_lift_curve=None): 
 99          """Set the lift L{Curve} instance(s).""" 
100          if separated_lift_curve == None: 
101              self.curves["lift"] = lift_curve 
102          else: 
103              self.curves["alift"] = lift_curve 
104              self.curves["slift"] = separated_lift_curve 
105   
106 -    def set_drag(self, drag_curve): 
107          """Set the drag L{Curve} instance.""" 
108          self.curves["drag"] = drag_curve 
109   
110 -    def init(self): 
111          """Initialize the curves""" 
112          for curve in self.curves.values(): 
113              curve.init() 
114   
115 -    def get_coefficient(self, key, attack_angle): 
116          """Return a coefficient from the curve key and an angle of attack""" 
117          return self.curves[key].get_coefficient_from(attack_angle) 
118   
119   
120 -class AerodynamicsSection: 
121      """The section for the aerodynamics calculations 
122      of the unsteady BEM code. 
123      """ 
124   
125 -    def __init__(self): 
126          self.airfoil = Airfoil() 
127          self.flow_angle = Angle(1., "rad") 
128          self.attack_angle = Angle(0., "deg") 
129           
130          self.radial_position = 0. 
131          self.twist = None 
132          self.chord = None 
133           
134          # 'cu' stands for current 
135          # 'int' for intermediate 
136          # 'qs' for quasi static 
137          self.induced_velocity = {} 
138          self.old_induced_velocity = {} 
139          for key in ["cu", "qs"]: 
140              self.induced_velocity[key] = N.zeros((3, 1)) 
141          for key in ["cu", "int", "qs"]: 
142              self.old_induced_velocity[key] = N.zeros((3, 1)) 
143           
144          self.rvector = N.zeros((3, 1)) 
145          self.wind_speed = N.zeros((3, 1)) 
146          self.velocity_prime = N.zeros((3, 1)) 
147          self.velocity_prime_norm = 0. 
148          self.relative_velocity = N.zeros((3, 1)) 
149   
150          self.prandtl_factor = 1. 
151          self.lift = 1. 
152           
153          # For the dynamic stall 
154          self.stall_coeff = {"lift" : 0.} 
155          self.old_attached_flow = 1. 
156          self.first_run = True 
157   
158          self.press = {"tangential" : 0., "normal" : 0.} 
159          self.tangential_press = 0. 
160          self.normal_press = 0. 
161   
162          self.rotor_plane_normal = N.array([ [-1.], 
163                                              [0.], 
164                                              [0.] ]) 
165          self._get_coefficients = self._get_direct_coefficients 
166       
167 -    def set_radius(self, radial_position): 
168          """Set the radial position of the section along 
169          the blade radius""" 
170          self.radial_position = radial_position 
171   
172 -    def set_chord(self, chord): 
173          """Set the chord of the section""" 
174          self.chord = chord 
175   
176 -    def set_twist(self, value, unit="deg"): 
177          """Set the twist angle for the section according 
178          to the blade plane. 
179           
180          @type value: a float 
181          @param value: the twist value 
182           
183          @type unit: a string 
184          @param unit: the 'rad' or 'deg' value""" 
185   
186          self.twist = Angle(value, unit) 
187   
188 -    def set_airfoil(self, airfoil): 
189          """Set the L{Airfoil} instance. In case the airfoil 
190          has a fully attached and fully separated lift curves, 
191          a dynamic stall model will be applied.""" 
192          if airfoil.curves.has_key("lift"): 
193              self._get_coefficients = self._get_direct_coefficients 
194          else: 
195              self._get_coefficients = self._get_dynamic_stall_coefficients 
196          self.airfoil = airfoil 
197   
198 -    def _get_direct_coefficients(self, time_step): 
199          """Return the lift and drag coefficient""" 
200          return self.airfoil.get_coefficient("lift", self.attack_angle), \ 
201                 self.airfoil.get_coefficient("drag", self.attack_angle) 
202   
203 -    def _apply_dynamic_stall(self, key, attached_coeff, 
204                                     separated_coeff, time_step): 
205          """Apply the dynamic stall model""" 
206          # TODO: Fixing the dynamic stall 
207          return attached_coeff 
208          time_constant = 4. * self.chord / \ 
209                          get_norm(self.relative_velocity) 
210          if attached_coeff == separated_coeff or self.first_run: 
211              attached_flow_st = 1. 
212              self.first_run = False 
213          else: 
214              attached_flow_st = (self.stall_coeff[key] - separated_coeff) / \ 
215                                 (attached_coeff - separated_coeff) 
216          attached_flow = attached_flow_st + \ 
217                          (self.old_attached_flow - attached_flow_st) * \ 
218                          M.exp(-time_step / time_constant) 
219          coeff = attached_flow * attached_coeff + \ 
220                      (1. - attached_flow) * separated_coeff 
221          # Update of the value for the next step 
222          self.old_attached_flow = attached_flow 
223          self.stall_coeff[key] = coeff 
224          return coeff 
225   
226 -    def _get_dynamic_stall_coefficients(self, time_step): 
227          """Return the lift and drag coefficients by applying 
228          a dynamic stall model""" 
229          acoeff = self.airfoil.get_coefficient("alift", self.attack_angle) 
230          scoeff = self.airfoil.get_coefficient("slift", self.attack_angle) 
231          lift_coeff = self._apply_dynamic_stall("lift", acoeff, 
232                                                 scoeff, time_step) 
233          drag_coeff = self.airfoil.get_coefficient("drag", self.attack_angle) 
234          return lift_coeff, drag_coeff 
235   
236 -    def init_unsteady(self): 
237          """Initialize the unsteady calculation""" 
238          self.calculate_rvector() 
239          self.airfoil.init() 
240   
241 -    def calculate_rvector(self): 
242          """Calculate the section vector""" 
243          self.rvector = N.array([ [0.], 
244                                   [0.], 
245                                   [self.radial_position] ]) 
246   
247 -    def calculate_relative_velocity(self, rotational_speed, cone_angle, 
248                                                              pitch_matrix): 
249          """Calculate the relative velocity from the velocity triangle""" 
250          rotational_vel = N.array([ [0.],  
251                                      [self.radial_position * \ 
252                                        rotational_speed["rad/s"] * \ 
253                                        N.cos(cone_angle["rad"])], 
254                                      [0.] ]) 
255          pitch_rotational_vel = N.asarray(pitch_matrix * rotational_vel) 
256          self.relative_velocity = self.wind_speed + pitch_rotational_vel + \ 
257                                   self.induced_velocity["cu"] 
258   
259 -    def calculate_pressures(self, air_density, time_step, pitch): 
260          """Calculate the aerodynamics pressures (in M{N/m}) at the  
261          section""" 
262          rvelocity_norm = get_norm(self.relative_velocity) 
263          self.flow_angle["rad"] = float(N.arctan(self.relative_velocity[0] / \ 
264                                                  self.relative_velocity[1])) 
265          self.attack_angle["rad"] = self.flow_angle["rad"] - self.twist["rad"] 
266                                             #(self.twist["rad"] + pitch["rad"]) 
267          lift_coeff, drag_coeff = self._get_coefficients(time_step) 
268          press = 0.5 * air_density * rvelocity_norm**2 * self.chord 
269          self.lift = press * lift_coeff 
270          drag = press * drag_coeff 
271          flow_sin = M.sin(self.flow_angle["rad"] + pitch["rad"]) 
272          flow_cos = M.cos(self.flow_angle["rad"] + pitch["rad"]) 
273          self.press["tangential"] = self.lift * flow_sin - drag * flow_cos 
274          self.press["normal"] = self.lift * flow_cos + drag * flow_sin 
275          self.tangential_press = self.press["tangential"] 
276          self.normal_press = self.press["normal"] 
277   
278 -    def calculate_qs_induced_velocity(self, prandlt_factor, air_density, 
279                                                              nb_blades): 
280          """Calculate the quasi static induced velocity, without taking 
281          into account the dynamic wake""" 
282          self.old_induced_velocity["qs"] = self.induced_velocity["qs"].copy() 
283          num = - nb_blades * self.lift 
284          dem = 4. * air_density * N.pi * \ 
285                self.radial_position * prandlt_factor 
286          scalar = N.dot(self.rotor_plane_normal.transpose(), 
287                         self.induced_velocity["cu"]) 
288          self.velocity_prime = self.wind_speed + \ 
289                                scalar * self.rotor_plane_normal 
290          self.velocity_prime_norm = get_norm(self.velocity_prime) 
291          self.induced_velocity["qs"][0] = num * N.cos(self.flow_angle["rad"]) \ 
292                                           / (dem * self.velocity_prime_norm) 
293          self.induced_velocity["qs"][1] = -num * N.sin(self.flow_angle["rad"]) \ 
294                                           / (dem * self.velocity_prime_norm) 
295    
296 -    def disable_dynamic_wake(self, induction_factor_max, adjust_factor, 
297                                     wind_speed, time_step, blade_radius): 
298          """Disable the dynamic wake by giving the reference of  
299          the quasi static induced velocity to the current induced  
300          velocity""" 
301          for key in ["cu", "int", "qs"]: 
302              self.old_induced_velocity[key] = self.induced_velocity["qs"].copy() 
303          self.induced_velocity["cu"] = self.induced_velocity["qs"].copy() 
304       
305 -    def apply_dynamic_wake(self, induction_factor_max, adjust_factor, 
306                                   wind_speed, time_step, blade_radius): 
307          """Calculate the current induced velocity by applying the two filters 
308          of the dynamic wake""" 
309          self.old_induced_velocity["cu"] = self.induced_velocity["cu"].copy() 
310          norm_wake_prime_velocity = get_norm(self.wind_speed + \ 
311                                              self.induced_velocity["cu"]) 
312          norm_wake_prime_velocity = self.velocity_prime_norm  
313          induction_factor = 1. - norm_wake_prime_velocity / wind_speed 
314          if induction_factor > induction_factor_max: 
315              induction_factor = induction_factor_max 
316          time_cst1 = 1.1 * blade_radius / \ 
317                      ((1. - 1.3 * induction_factor) * wind_speed) 
318          time_cst2 = (0.39 - 0.26 * (self.radial_position/blade_radius)**2) \ 
319                      * time_cst1 
320   
321          # First filter 
322          cst = self.induced_velocity["qs"] + \ 
323                adjust_factor * time_cst1 / time_step * \ 
324                (self.induced_velocity["qs"] - self.old_induced_velocity["qs"]) 
325          int_induced_velocity = cst + M.exp(-time_step / time_cst1) * \ 
326                                (self.old_induced_velocity["int"] - cst) 
327           
328          # Second filter and final value of the induced velocity 
329          self.induced_velocity["cu"] = int_induced_velocity + \ 
330                                        M.exp(-time_step / time_cst2) * \ 
331                                        (self.old_induced_velocity["cu"] - \ 
332                                         int_induced_velocity) 
333           
334          # Updating the value for the intermediate filter 
335          self.old_induced_velocity["int"] = int_induced_velocity.copy() 
336   
337   
338 -class Section(AerodynamicsSection, BasicObject): 
339      """The section definition. This class gathers 
340      the aerodynamics calculations and the saving tasks 
341      for retrieving the section status in another process, 
342      after calculations. 
343      """ 
344   
345 -    def __init__(self, name="section"): 
346          AerodynamicsSection.__init__(self) 
347          BasicObject.__init__(self, name) 
348           
349          self.own_para_names += ["radial_position", 
350                                  "chord", 
351                                  "twist", 
352                                  "airfoil"] 
353   
354 -    def collect_parameters(self): 
355          """Collect the section parameters into a dictionary""" 
356          self.collect_own_parameters() 
357   
358 -    def set_parameters(self, para): 
359          """Set the sections parameters from the C{para} 
360          dictionary""" 
361          self.set_own_parameters(para) 
362