1 from __future__ import division
2
3 import numpy
4
5
6
8 - def __init__(self, thickness, edge_coeff):
9 self.thickness = thickness
10 self.edge_coeff = edge_coeff
11
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
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
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
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
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
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
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