Package meshpy :: Module naca
[hide private]
[frames] | no frames]

Source Code for Module meshpy.naca

  1  from __future__ import division 
  2   
  3  import numpy 
  4   
  5   
  6   
7 -class FourDigitsSymmetric:
8 - def __init__(self, thickness, edge_coeff):
9 self.thickness = thickness 10 self.edge_coeff = edge_coeff
11
12 - def __call__(self, x, side):
13 t = self.thickness 14 15 def y_upper(y): 16 return y
17 18 def y_lower(y): 19 return -y
20 21 def x_upper(x): 22 return x 23 24 def x_lower(x): 25 return x 26 27 y = t * 5 * (0.2969 * numpy.sqrt(x) + ((((-self.edge_coeff * x + 28 0.2843) * x - 0.3516) * x) - 0.126) * x) 29 30 if side == "upper": 31 return numpy.array([x_upper(x), y_upper(y)]) 32 elif side == "lower": 33 return numpy.array([x_lower(x), y_lower(y)]) 34 else: 35 raise ValueError("Neither upper nor lower side selected in the call.") 36 37 38 39 40
41 -class FourDigitsCambered:
42 - def __init__(self, thickness, max_camber, max_camber_pos, edge_coeff):
43 self.thickness = thickness 44 self.max_camber = max_camber 45 self.max_camber_pos = max_camber_pos 46 self.edge_coeff = edge_coeff
47
48 - def __call__(self, x, side):
49 t = self.thickness 50 m = self.max_camber 51 p = self.max_camber_pos 52 53 def y_upper(y_c, y, theta): 54 return y_c + y * numpy.cos(theta)
55 56 def y_lower(y_c, y, theta): 57 return y_c - y * numpy.cos(theta)
58 59 def x_upper(x, y, theta): 60 return x - y * numpy.sin(theta) 61 62 def x_lower(x, y, theta): 63 return x + y * numpy.sin(theta) 64 65 y = t * 5 * (0.2969 * numpy.sqrt(x) + ((((-self.edge_coeff * x + 66 0.2843) * x - 0.3516) * x) -0.126) * x) 67 68 if x <= p: 69 y_c = m * x / p ** 2 * (2 * p - x) 70 theta = numpy.arctan(2 * m / p ** 2 * (p - x)) 71 else: 72 y_c = m * (1 - x) / (1 - p) ** 2 * (1 + x - 2 * p) 73 theta = numpy.arctan(m / (1 - p) ** 2 * (2 * p - 2 * x)) 74 75 if side == "upper": 76 return numpy.array([x_upper(x, y, theta), y_upper(y_c, y, theta)]) 77 elif side == "lower": 78 return numpy.array([x_lower(x, y, theta), y_lower(y_c, y, theta)]) 79 else: 80 raise ValueError("Neither upper nor lower side selected in the call.") 81 82 83 84 85
86 -class FiveDigits:
87 - def __init__(self, thickness, m, k1, edge_coeff):
88 self.thickness = thickness 89 self.edge_coeff = edge_coeff 90 self.m = m 91 self.k1 = k1
92
93 - def __call__(self, x, side):
94 t = self.thickness 95 m = self.m 96 k1 = self.k1 97 98 def y_upper(y_c, y, theta): 99 return y_c + y * numpy.cos(theta)
100 101 def y_lower(y_c, y, theta): 102 return y_c - y * numpy.cos(theta)
103 104 def x_upper(x, y, theta): 105 return x - y * numpy.sin(theta) 106 107 def x_lower(x, y, theta): 108 return x + y * numpy.sin(theta) 109 110 y = t * 5 * (0.2969 * numpy.sqrt(x) + ((((-self.edge_coeff * x + 111 0.2843) * x - 0.3516) * x) -0.126) * x) 112 113 if x <= m: 114 y_c = k1 / 6 * x *((x - 3 * m) * x + m ** 2 * (3 - m)) 115 theta = numpy.arctan(k1 / 6 * ((3 * x - 6 * m) * x + 116 m ** 2 * (3 - m))) 117 else: 118 y_c = k1 * m ** 3 / 6 * (1 - x) 119 theta = numpy.arctan(-k1 * m ** 3 / 6) 120 121 if side == "upper": 122 return numpy.array([x_upper(x, y, theta), y_upper(y_c, y, theta)]) 123 elif side == "lower": 124 return numpy.array([x_lower(x, y, theta), y_lower(y_c, y, theta)]) 125 else: 126 raise ValueError("Neither upper nor lower side selected in the call.") 127 128 129 130 131
132 -def get_naca_points(naca_digits, number_of_points=100, 133 sharp_trailing_edge=True, 134 abscissa_map=lambda x: 0.03*x+0.97*x**2, 135 verbose=False):
136 """ 137 Return a list of coordinates of NACA 4-digit and 5-digit series 138 airfoils. 139 """ 140 141 if verbose: 142 def explain(*s): 143 print " ".join(str(s_i) for s_i in s)
144 else: 145 def explain(*s): 146 pass 147 148 explain("Airfoil: NACA-%s" %(naca_digits)) 149 150 if sharp_trailing_edge == True: 151 explain("Sharp trailing edge") 152 edge_coeff = 0.1036 153 else: 154 explain("Blunt trailing edge") 155 edge_coeff = 0.1015 156 157 158 raw_abscissae = numpy.linspace(0, 1, number_of_points, endpoint=True) 159 abscissae = numpy.empty_like(raw_abscissae) 160 for i in xrange(number_of_points): 161 abscissae[i] = abscissa_map(raw_abscissae[i]) 162 163 digits_int = int(naca_digits) 164 if len(naca_digits) == 4: 165 thickness = (digits_int % 100) 166 max_camber_pos = (digits_int % 1000) - thickness 167 max_camber = (digits_int % 10000) - max_camber_pos - thickness 168 169 thickness = thickness / 1e2 170 max_camber_pos = max_camber_pos / 1e3 171 max_camber = max_camber / 1e5 172 173 explain("Thickness:", thickness) 174 explain("Position of maximum camber:", max_camber_pos) 175 explain("Maximum camber:", max_camber) 176 177 if max_camber == 0 and max_camber_pos == 0: 178 explain("Symmetric 4-digit airfoil") 179 points = FourDigitsSymmetric(thickness, edge_coeff) 180 elif max_camber != 0 and max_camber_pos != 0: 181 explain("Cambered 4-digit airfoil") 182 points = FourDigitsCambered(thickness, max_camber, 183 max_camber_pos, edge_coeff) 184 else: 185 raise NotImplementedError("You must decide whether your airfoil shall be cambered or not!") 186 187 elif len(naca_digits) == 5: 188 thickness = (digits_int % 100) 189 max_camber_pos = (digits_int % 10000) - thickness 190 191 thickness = thickness / 1e2 192 max_camber_pos = max_camber_pos / 2e4 193 194 explain("Thickness:", thickness) 195 explain("Position of maximum camber:", max_camber_pos) 196 197 identifier = digits_int // 100 198 if identifier == 210: 199 m = 0.058 200 k1 = 361.4 201 elif identifier == 220: 202 m = 0.126 203 k1 = 51.64 204 elif identifier == 230: 205 m = 0.2025 206 k1 = 15.957 207 elif identifier == 240: 208 m = 0.29 209 k1 = 6.643 210 elif identifier == 250: 211 m = 0.391 212 k1 = 3.23 213 else: 214 raise NotImplementedError("5-digit series only implemented for " 215 "the first three digits in 210, 220, 230, 240, 250!") 216 217 explain("5-digit airfoil") 218 points = FiveDigits(thickness, m, k1, edge_coeff) 219 220 else: 221 raise NotImplementedError( 222 "Only the 4-digit and 5-digit series are implemented!") 223 224 points_upper = numpy.zeros((len(abscissae),2)) 225 points_lower = numpy.zeros((len(abscissae),2)) 226 227 for i in range(len(abscissae)): 228 points_upper[i] = points(abscissae[i], "upper") 229 points_lower[i] = points(abscissae[i], "lower") 230 231 if sharp_trailing_edge == True: 232 return list(points_upper)[1:-1] + list(points_lower[::-1]) 233 else: 234 return list(points_upper)[1:] + list(points_lower[::-1]) 235 236 237 238 239
240 -def write_points(points, filename):
241 file = open(filename, "w") 242 for pt in points: 243 print >>file, "\t".join(repr(p_comp) for p_comp in pt)
244 245 246 247
248 -def main():
249 from optparse import OptionParser 250 251 parser = OptionParser(usage="%prog AIRFOIL-ID") 252 parser.add_option("-o", "--output", 253 help="write ouput to FILE", metavar="FILE") 254 parser.add_option("-p", "--points", type="int", 255 help="generate N points", metavar="N") 256 parser.add_option("-s", "--sharp-trailing-edge", action="store_true") 257 parser.add_option("-u", "--uniform-distribution", action="store_true") 258 parser.add_option("-q", "--quiet", 259 action="store_false", dest="verbose", default=True, 260 help="Don't print status messages to stdout") 261 262 (options, args) = parser.parse_args() 263 264 if not args: 265 parser.print_help() 266 return 267 268 if options.points is None: 269 options.points = 100 270 271 digits = args[0] 272 points = generate_naca(digits, 273 number_of_points=options.points, 274 sharp_trailing_edge=options.sharp_trailing_edge, 275 uniform_distribution=options.uniform_distribution, 276 verbose=options.verbose) 277 278 if options.output is None: 279 options.output = "naca-%s.dat" % digits 280 281 print "Output file:", options.output 282 write_points(points, options.output)
283 284 285 286 287 if __name__ == "__main__": 288 main() 289