-
Notifications
You must be signed in to change notification settings - Fork 1
/
airfoils.py
206 lines (170 loc) · 6.81 KB
/
airfoils.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
#===================================================#
# airfoils.py #
# #
# author: Guillaume Chauvat #
# email: [email protected] #
# created 2020-02-10, last modification: 2021-09-23 #
#===================================================#
import numpy as np
from scipy import interpolate
class Airfoil:
"""
Class representing an airfoil as closed 2D curve {t -> (x(t), y(t)) for t in [0, 1]} using cubic splines.
This allows for smooth interpolation of the coordinates as well as their derivatives.
"""
def __init__(self, x, y, n=10000, tol=1e-13, verbose=False):
self.n = n
self.tol = tol
self.verbose = verbose
self.x = x
self.y = y
(self.t, self.length) = self.compute_abscissas()
self.x_int = interpolate.CubicSpline(self.t, self.x)
self.y_int = interpolate.CubicSpline(self.t, self.y)
self.dx_int = self.x_int.derivative(1)
self.dy_int = self.y_int.derivative(1)
self.d2x_int = self.x_int.derivative(2)
self.d2y_int = self.y_int.derivative(2)
self.trigonometric_direction = self.is_trigonometric()
def compute_abscissas(self):
"""
Computes the normalized curvilinear abscissas t corresponding to the coordinates of the airfoil.
"""
# first, approximate the curve by linear interpolation to get an approximation of the curvilinear abscissa.
lengths = np.sqrt(np.diff(self.x)**2 + np.diff(self.y)**2)
s = np.concatenate(([0], lengths)).cumsum() # curvilinear abscissa
t = s/s[-1] # normalize
# compute the new curvilinear abscissa and iterate until convergence
tc1 = np.linspace(0, 1, self.n)
err = 1.
i = 0
if self.verbose:
print("Computing the curvilinear abscissa...")
while err > self.tol:
# get the high resolution points
x1 = interpolate.CubicSpline(t, self.x)(tc1)
y1 = interpolate.CubicSpline(t, self.y)(tc1)
# approximate the curvilinear abscissa based on a linear interpolation between high resolution points
lengths = np.sqrt(np.diff(x1)**2 + np.diff(y1)**2)
s1 = np.concatenate(([0], lengths)).cumsum()
t1 = s1/s1[-1]
# update the normalized curvilinear abscissa at coordinate points
t_new = interpolate.CubicSpline(tc1, t1)(t)
err = np.linalg.norm(t_new-t)
t = t_new
if self.verbose:
print(f"iteration {i}, err = {err}")
i = i + 1
return (t, s1[-1])
def pos(self, tc):
"""
Interpolates x: t->x(t) and y: t->y(t) at points in the array tc.
"""
return (self.x_int(tc), self.y_int(tc))
def curvilinear_abscissa(self, x0, t0, tol=1e-13):
"""
Returns the normalized curvilinear abscissa corresponding to a x/c location on the airfoil (where c is the chord).
The choice of t0 is important since there are always two possibilities for 0 < x < 1.
x0: location where the curvilinear abscissa must be computed
t0: first guess of the location
"""
# initial value
f0 = self.x_int(t0)-x0
# iterate with Newton-Raphson
i = 0
err = 1.0
if self.verbose:
print("Computing t...")
while np.any(err > tol):
df = self.dx_int(t0)
t1 = t0 - f0/df
f0 = self.x_int(t1)-x0
err = abs(t1-t0)
t0 = t1
if self.verbose:
print(f"iteration {i}, err = {err}")
i = i+1
return t1
def normal(self, t0):
"""
returns the normal vector to the surface at a given point
"""
dx = self.dx_int(t0)
dy = self.dy_int(t0)
l = np.sqrt(dx**2+dy**2)
nx = dy/l
ny = -dx/l
if not self.trigonometric_direction:
nx = -nx
ny = -ny
return (nx, ny)
def tangent(self, t0):
"""
returns the tangent vector to the surface at a given point
"""
dx = self.dx_int(t0)
dy = self.dy_int(t0)
l = np.sqrt(dx**2+dy**2)
nx = dx/l
ny = dy/l
return (nx, ny)
def dtangent(self, t0):
"""
returns the derivative of the tangent vector to the surface at a given point
"""
# ||dx**2 + dy**2|| is approximately equal to self.length, but we want to be very accurate and account for small deviations here
# so speed ~ self.length and fact1 << 1/self.length
dx = self.dx_int(t0)
dy = self.dy_int(t0)
d2x = self.d2x_int(t0)
d2y = self.d2y_int(t0)
speed = np.sqrt(dx**2+dy**2)
fact1 = (dx*d2x+dy*d2y)/speed**3
return (d2x/speed-fact1*dx, d2y/speed-fact1*dy)
def closest(self, x0, y0, t0=0.5, tol=1e-13, alpha=1):
"""
returns the closest point (x, y) on the surface to an arbitrary point (x0, y0)
"""
# t0 is the starting point
# the variable alpha can be chosen between 0 and 1 to sacrifice speed for stability
# find a zero of the function f: t -> (tangent(t)|x(t)-x0)
# f': t: -> 1 + (tangent'(t)|x(t)-x0)
t1 = t0
(x1, y1) = self.pos(t1)
(taux, tauy) = self.tangent(t1)
f1 = taux*(x1-x0)+tauy*(y1-y0)
err = 1
i = 0
if self.verbose:
print(f"Projecting point ({x0}, {y0}) on airfoil...")
while err > tol:
(dtaux, dtauy) = self.dtangent(t1)
dx = self.dx_int(t1)
dy = self.dy_int(t1)
df = dx*taux + dy*tauy + dtaux*(x1-x0) + dtauy*(y1-y0)
t2 = t1 - f1/df
err = abs(t2-t1)
t1 = (1-alpha)*t1+alpha*t2
(x1, y1) = self.pos(t1)
(taux, tauy) = self.tangent(t1)
f1 = taux*(x1-x0)+tauy*(y1-y0)
if self.verbose:
print(f"iteration {i}, t = {t1}, f = {f1}, err = {err}")
i = i+1
return (t1, x1, y1)
def is_trigonometric(self):
"""
checks whether the airfoil is described in the trigonometric direction
"""
dx1 = self.x[1]-self.x[0]
dx2 = self.x[-2]-self.x[0]
dy1 = self.y[1]-self.y[0]
dy2 = self.y[-2]-self.y[0]
return dx1*dy2-dx2*dy1 > 0
def rotate(x, y, angle, xr, yr):
"""
rotates a set of points around a point (xr, yr)
"""
x1 = xr + (x-xr)*np.cos(angle) + (y-yr)*np.sin(angle)
y1 = yr - (x-xr)*np.sin(angle) + (y-yr)*np.cos(angle)
return (x1, y1)