forked from jte0419/Panel_Methods
-
Notifications
You must be signed in to change notification settings - Fork 0
/
VP_Airfoil.py
389 lines (337 loc) · 27.9 KB
/
VP_Airfoil.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
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
# VORTEX PANEL METHOD - SINGLE AIRFOIL
# Written by: JoshTheEngineer
# YouTube : www.youtube.com/joshtheengineer
# Website : www.joshtheengineer.com
# Started : 12/09/18 - In MATLAB
# Updated : 02/03/19 - Transferred from MATLAB to Python
# - Works as expected
# 02/09/20 - Added DAT airfoil loading option with XFOIL function
# Notes : This code is not optimized, but is instead written in such a way
# that it is easy to follow along with my YouTube video derivations
#
# Functions Needed:
# - XFOIL.py
# - COMPUTE_KL_VPM.py
# - STREAMLINE_VPM.py
# - COMPUTE_CIRCULATION.py
#
# Programs Needed:
# - xfoil.exe
#
# Folder Needed:
# - Airfoil_DAT_Selig: folder containing all Selig-format airfoils
#
# References
# - [1]: Panel Method Geometry
# Link: https://www.youtube.com/watch?v=kIqxbd937PI
# - [2]: Normal Geometric Integral VPM, K(ij)
# Link: https://www.youtube.com/watch?v=5lmIv2CUpoc
# - [3]: Tangential Geometric Integral VPM, L(ij)
# Link: https://www.youtube.com/watch?v=IxWJzwIG_gY
# - [4]: Streamline Geometric Integral VPM, Nx(pj) and Ny(pj)
# Link: https://www.youtube.com/watch?v=TBwBnW87hso
# - [5]: Solving the System of Equations (VPM)
# Link: https://www.youtube.com/watch?v=ep7vPzGYsbw
# - [6]: How To Compute Circulation
# Link: https://www.youtube.com/watch?v=b8EnhiSjL3o
# - [7]: UIUC Airfoil Database: Download All Files using Python
# Link: https://www.youtube.com/watch?v=nILo18DlqAo
# - [8]: Python code for downloading Selig airfoil DAT files
# Link: http://www.joshtheengineer.com/2019/01/30/uiuc-airfoil-database-file-download/
import numpy as np
import math as math
import matplotlib.pyplot as plt
from matplotlib import path
from XFOIL import XFOIL
from COMPUTE_KL_VPM import COMPUTE_KL_VPM
from STREAMLINE_VPM import STREAMLINE_VPM
from COMPUTE_CIRCULATION import COMPUTE_CIRCULATION
# %% KNOWNS
# Flag to specify creating or loading airfoil
flagAirfoil = [1, # Create specified NACA airfoil in XFOIL
0] # Load Selig-format airfoil from directory
# User-defined knowns
Vinf = 1 # Freestream velocity [] (just leave this at 1)
AoA = 0 # Angle of attack [deg]
NACA = '0012' # NACA airfoil to load [####]
# Convert angle of attack to radians
AoAR = AoA*(np.pi/180) # Angle of attack [rad]
# Plotting flags
flagPlot = [0, # Airfoil with panel normal vectors
0, # Geometry boundary pts, control pts, first panel, second panel
1, # Cp vectors at airfoil surface panels
1, # Pressure coefficient comparison (XFOIL vs. VPM)
0, # Airfoil streamlines
0] # Pressure coefficient contour
# %% XFOIL - CREATE/LOAD AIRFOIL
# PPAR menu options
PPAR = ['170', # "Number of panel nodes"
'4', # "Panel bunching paramter"
'1.5', # "TE/LE panel density ratios"
'1', # "Refined area/LE panel density ratio"
'1 1', # "Top side refined area x/c limits"
'1 1'] # "Bottom side refined area x/c limits"
# Call XFOIL function to obtain the following:
# - Airfoil coordinates
# - Pressure coefficient along airfoil surface
# - Lift, drag, and moment coefficients
xFoilResults = XFOIL(NACA, PPAR, AoA, flagAirfoil)
# Separate out results from XFOIL function results
afName = xFoilResults[0] # Airfoil name
xFoilX = xFoilResults[1] # X-coordinate for Cp result
xFoilY = xFoilResults[2] # Y-coordinate for Cp result
xFoilCP = xFoilResults[3] # Pressure coefficient
XB = xFoilResults[4] # Boundary point X-coordinate
YB = xFoilResults[5] # Boundary point Y-coordinate
xFoilCL = xFoilResults[6] # Lift coefficient
xFoilCD = xFoilResults[7] # Drag coefficient
xFoilCM = xFoilResults[8] # Moment coefficient
# Number of boundary points and panels
numPts = len(XB) # Number of boundary points
numPan = numPts - 1 # Number of panels (control points)
# %% CHECK PANEL DIRECTIONS - FLIP IF NECESSARY
# Check for direction of points
edge = np.zeros(numPan) # Initialize edge value array
for i in range(numPan): # Loop over all panels
edge[i] = (XB[i+1]-XB[i])*(YB[i+1]+YB[i]) # Compute edge values
sumEdge = np.sum(edge) # Sum all edge values
# If panels are CCW, flip them (don't if CW)
if (sumEdge < 0): # If panels are CCW
XB = np.flipud(XB) # Flip the X-data array
YB = np.flipud(YB) # Flip the Y-data array
# %% PANEL METHOD GEOMETRY - REF [1]
# Initialize variables
XC = np.zeros(numPan) # Initialize control point X-coordinate array
YC = np.zeros(numPan) # Initialize control point Y-coordinate array
S = np.zeros(numPan) # Initialize panel length array
phi = np.zeros(numPan) # Initialize panel orientation angle array [deg]
# Find geometric quantities of the airfoil
for i in range(numPan): # Loop over all panels
XC[i] = 0.5*(XB[i]+XB[i+1]) # X-value of control point
YC[i] = 0.5*(YB[i]+YB[i+1]) # Y-value of control point
dx = XB[i+1]-XB[i] # Change in X between boundary points
dy = YB[i+1]-YB[i] # Change in Y between boundary points
S[i] = (dx**2 + dy**2)**0.5 # Length of the panel
phi[i] = math.atan2(dy,dx) # Angle of panel (positive X-axis to inside face)
if (phi[i] < 0): # Make all panel angles positive [rad]
phi[i] = phi[i] + 2*np.pi
# Compute angle of panel normal w.r.t. horizontal and include AoA
delta = phi + (np.pi/2) # Angle from positive X-axis to outward normal vector [rad]
beta = delta - AoAR # Angle between freestream vector and outward normal vector [rad]
beta[beta > 2*np.pi] = beta[beta > 2*np.pi] - 2*np.pi # Make all panel angles between 0 and 2pi [rad]
# %% COMPUTE VORTEX PANEL STRENGTHS - REF [5]
# Geometric integral (normal [K] and tangential [L])
# - Refs [2] and [3]
K, L = COMPUTE_KL_VPM(XC,YC,XB,YB,phi,S) # Compute geometric integrals
# Populate A matrix
A = np.zeros([numPan,numPan]) # Initialize the A matrix
for i in range(numPan): # Loop over all i panels
for j in range(numPan): # Loop over all j panels
if (i == j): # If the panels are the same
A[i,j] = 0 # Set A equal to zero
else: # If panels are not the same
A[i,j] = -K[i,j] # Set A equal to negative geometric integral
# Populate b array
b = np.zeros(numPan) # Initialize the b array
for i in range(numPan): # Loop over all panels
b[i] = -Vinf*2*np.pi*np.cos(beta[i]) # Compute RHS array
# Satisfy the Kutta condition
pct = 100 # Panel replacement percentage
panRep = int((pct/100)*numPan)-1 # Replace this panel with Kutta condition equation
if (panRep >= numPan): # If we specify the last panel
panRep = numPan-1 # Set appropriate replacement panel index
A[panRep,:] = 0 # Set all colums of the replaced panel equation to zero
A[panRep,0] = 1 # Set first column of replaced panel equal to 1
A[panRep,numPan-1] = 1 # Set last column of replaced panel equal to 1
b[panRep] = 0 # Set replaced panel value in b array equal to zero
# Compute gamma values
gamma = np.linalg.solve(A,b) # Compute all vortex strength values
# %% COMPUTE PANEL VELOCITIES AND PRESSURE COEFFICIENTS
# Compute velocities
Vt = np.zeros(numPan) # Initialize tangential velocity array
Cp = np.zeros(numPan) # Initialize pressure coefficient array
for i in range(numPan): # Loop over all i panels
addVal = 0 # Reset summation value to zero
for j in range(numPan): # Loop over all j panels
addVal = addVal - (gamma[j]/(2*np.pi))*L[i,j] # Sum all tangential vortex panel terms
Vt[i] = Vinf*np.sin(beta[i]) + addVal + gamma[i]/2 # Compute tangential velocity by adding uniform flow and i=j terms
Cp[i] = 1 - (Vt[i]/Vinf)**2 # Compute pressure coefficient
# %% COMPUTE LIFT AND MOMENT COEFFICIENTS
# Compute normal and axial force coefficients
CN = -Cp*S*np.sin(beta) # Normal force coefficient []
CA = -Cp*S*np.cos(beta) # Axial force coefficient []
# Compute lift, drag, and moment coefficients
CL = sum(CN*np.cos(AoAR)) - sum(CA*np.sin(AoAR)) # Decompose axial and normal to lift coefficient []
CM = sum(Cp*(XC-0.25)*S*np.cos(phi)) # Moment coefficient []
# Print the results to the Console
print("======= RESULTS =======")
print("Lift Coefficient (CL)")
print(" K-J : %2.8f" % (2*sum(gamma*S))) # From Kutta-Joukowski lift equation
print(" VPM : %2.8f" % CL) # From this VPM code
print(" XFOIL: %2.8f" % xFoilCL) # From XFOIL program
print("Moment Coefficient (CM)")
print(" VPM : %2.8f" % CM) # From this VPM code
print(" XFOIL: %2.8f" % xFoilCM) # From XFOIL program
# %% COMPUTE STREAMLINES - REF [4]
if (flagPlot[4] == 1 or flagPlot[5] == 1): # If we are plotting streamlines or pressure coefficient contours
# Grid parameters
nGridX = 150 # X-grid for streamlines and contours
nGridY = 150 # Y-grid for streamlines and contours
xVals = [-0.5, 1.5] # X-grid extents [min, max]
yVals = [-0.3, 0.3] # Y-grid extents [min, max]
# Streamline parameters
slPct = 25 # Percentage of streamlines of the grid
Ysl = np.linspace(yVals[0],yVals[1],int((slPct/100)*nGridY)) # Create array of Y streamline starting points
Xsl = xVals[0]*np.ones(len(Ysl)) # Create array of X streamline starting points
XYsl = np.vstack((Xsl.T,Ysl.T)).T # Concatenate X and Y streamline starting points
# Generate the grid points
Xgrid = np.linspace(xVals[0],xVals[1],nGridX) # X-values in evenly spaced grid
Ygrid = np.linspace(yVals[0],yVals[1],nGridY) # Y-values in evenly spaced grid
XX, YY = np.meshgrid(Xgrid,Ygrid) # Create meshgrid from X and Y grid arrays
# Initialize velocities
Vx = np.zeros([nGridX,nGridY]) # Initialize X velocity matrix
Vy = np.zeros([nGridX,nGridY]) # Initialize Y velocity matrix
# Path to figure out if grid point is inside polygon or not
AF = np.vstack((XB.T,YB.T)).T # Concatenate XB and YB geometry points
afPath = path.Path(AF) # Create a path for the geometry
# Solve for grid point X and Y velocities
for m in range(nGridX): # Loop over X-grid points
for n in range(nGridY): # Loop over Y-grid points
XP = XX[m,n] # Current iteration's X grid point
YP = YY[m,n] # Current iteration's Y grid point
Nx, Ny = STREAMLINE_VPM(XP,YP,XB,YB,phi,S) # Compute Nx and Ny geometric integrals
# Check if grid points are in object
# - If they are, assign a velocity of zero
if afPath.contains_points([(XP,YP)]): # If (XP,YP) is in the body
Vx[m,n] = 0 # Set X-velocity equal to zero
Vy[m,n] = 0 # Set Y-velocity equal to zero
else:
Vx[m,n] = Vinf*np.cos(AoAR) + sum(-gamma*Nx/(2*np.pi)) # Compute X-velocity
Vy[m,n] = Vinf*np.sin(AoAR) + sum(-gamma*Ny/(2*np.pi)) # Compute Y-velocity
# Compute grid point velocity magnitude and pressure coefficient
Vxy = np.sqrt(Vx**2 + Vy**2) # Compute magnitude of velocity vector []
CpXY = 1 - (Vxy/Vinf)**2 # Pressure coefficient []
# %% CIRCULATION AND VORTEX STRENGTH CHECK
if (flagPlot[4] == 1 or flagPlot[5] == 1): # If we are plotting streamlines or Cp contours
# Compute circulation
aa = 0.75 # Ellipse horizontal half-length
bb = 0.25 # Ellipse vertical half-length
x0 = 0.5 # Ellipse center X-coordinate
y0 = 0 # Ellipse center Y-coordinate
numT = 5000 # Number of points on ellipse
Circulation, xC, yC, VxC, VyC = COMPUTE_CIRCULATION(aa,bb,x0,y0, # Compute circulation around ellipse
numT,Vx,Vy,Xgrid,Ygrid)
# Print values to Console
print("======= CIRCULATION RESULTS =======")
print("Sum of L : %2.8f" % sum(gamma*S)) # Print sum of vortex strengths
print("Circulation: %2.8f" % Circulation) # Print circulation
print("Lift Coef : %2.8f" % (2.0*Circulation)) # Lift coefficient from K-J equation
# %% PLOTTING
# FIGURE: Airfoil with panel normal vectors
if (flagPlot[0] == 1):
fig = plt.figure(1) # Create the figure
plt.cla() # Clear the axes
plt.fill(XB,YB,'k') # Plot the airfoil
X = np.zeros(2) # Initialize 'X'
Y = np.zeros(2) # Initialize 'Y'
for i in range(numPan): # Loop over all panels
X[0] = XC[i] # Set X start of panel orientation vector
X[1] = XC[i] + S[i]*np.cos(delta[i]) # Set X end of panel orientation vector
Y[0] = YC[i] # Set Y start of panel orientation vector
Y[1] = YC[i] + S[i]*np.sin(delta[i]) # Set Y end of panel orientation vector
if (i == 0): # If it's the first panel index
plt.plot(X,Y,'b-',label='First Panel') # Plot normal vector for first panel
elif (i == 1): # If it's the second panel index
plt.plot(X,Y,'g-',label='Second Panel') # Plot normal vector for second panel
else: # If it's neither the first nor second panel index
plt.plot(X,Y,'r-') # Plot normal vector for all other panels
plt.xlabel('X Units') # Set X-label
plt.ylabel('Y Units') # Set Y-label
plt.title('Panel Geometry') # Set title
plt.axis('equal') # Set axes equal
plt.legend() # Display legend
plt.show() # Display plot
# FIGURE: Geometry with the following indicated:
# - Boundary points, control points, first panel, second panel
if (flagPlot[1] == 1):
fig = plt.figure(2) # Create figure
plt.cla() # Get ready for plotting
plt.plot(XB,YB,'k-') # Plot airfoil panels
plt.plot([XB[0], XB[1]],[YB[0], YB[1]],'b-',label='First Panel') # Plot first panel
plt.plot([XB[1], XB[2]],[YB[1], YB[2]],'g-',label='Second Panel') # Plot second panel
plt.plot(XB,YB,'ko',markerfacecolor='k',label='Boundary Pts') # Plot boundary points (black circles)
plt.plot(XC,YC,'ko',markerfacecolor='r',label='Control Pts') # Plot control points (red circles)
plt.xlabel('X Units') # Set X-label
plt.ylabel('Y Units') # Set Y-label
plt.axis('equal') # Set axes equal
plt.legend() # Display legend
plt.show() # Display plot
# FIGURE: Cp vectors at airfoil control points
if (flagPlot[2] == 1):
fig = plt.figure(3) # Create figure
plt.cla() # Get ready for plotting
Cps = np.absolute(Cp*0.15) # Scale and make positive all Cp values
X = np.zeros(2) # Initialize X values
Y = np.zeros(2) # Initialize Y values
for i in range(len(Cps)): # Loop over all panels
X[0] = XC[i] # Control point X-coordinate
X[1] = XC[i] + Cps[i]*np.cos(delta[i]) # Ending X-value based on Cp magnitude
Y[0] = YC[i] # Control point Y-coordinate
Y[1] = YC[i] + Cps[i]*np.sin(delta[i]) # Ending Y-value based on Cp magnitude
if (Cp[i] < 0): # If pressure coefficient is negative
plt.plot(X,Y,'r-') # Plot as a red line
elif (Cp[i] >= 0): # If pressure coefficient is zero or positive
plt.plot(X,Y,'b-') # Plot as a blue line
plt.fill(XB,YB,'k') # Plot the airfoil as black polygon
plt.xlabel('X Units') # Set X-label
plt.ylabel('Y Units') # Set Y-label
plt.gca().set_aspect('equal') # Set aspect ratio equal
plt.show() # Show the plot
# FIGURE: Pressure coefficient comparison (XFOIL vs. VPM)
if (flagPlot[3] == 1):
fig = plt.figure(4) # Create figure
plt.cla() # Get ready for plotting
midIndX = int(np.floor(len(xFoilCP)/2)) # Airfoil middle index for XFOIL data
midIndS = int(np.floor(len(Cp)/2)) # Airfoil middle index for VPM data
plt.plot(xFoilX[0:midIndX],xFoilCP[0:midIndX], # Plot Cp for upper surface of airfoil from XFoil
'b-',label='XFOIL Upper')
plt.plot(xFoilX[midIndX+1:len(xFoilX)],xFoilCP[midIndX+1:len(xFoilX)], # Plot Cp for lower surface of airfoil from XFoil
'r-',label='XFOIL Lower')
plt.plot(XC[midIndS+1:len(XC)],Cp[midIndS+1:len(XC)], # Plot Cp for upper surface of airfoil from panel method
'ks',markerfacecolor='b',label='VPM Upper')
plt.plot(XC[0:midIndS],Cp[0:midIndS], # Plot Cp for lower surface of airfoil from panel method
'ks',markerfacecolor='r',label='VPM Lower')
plt.xlim(0,1) # Set X-limits
plt.xlabel('X Coordinate') # Set X-label
plt.ylabel('Cp') # Set Y-label
plt.title('Pressure Coefficient') # Set title
plt.show() # Display plot
plt.legend() # Display legend
plt.gca().invert_yaxis() # Invert Cp (Y) axis
# FIGURE: Airfoil streamlines
if (flagPlot[4] == 1):
fig = plt.figure(5) # Create figure
plt.cla() # Get ready for plotting
np.seterr(under="ignore") # Ignore underflow error message
plt.streamplot(XX,YY,Vx,Vy, linewidth=0.5, density=40, color='r', # Plot streamlines
arrowstyle='-', start_points=XYsl)
plt.clim(vmin=0, vmax=2)
plt.fill(XB,YB,'k') # Plot airfoil as black polygon
plt.xlabel('X Units') # Set X-label
plt.ylabel('Y Units') # Set Y-label
plt.gca().set_aspect('equal') # Set axes equal
plt.xlim(xVals) # Set X-limits
plt.ylim(yVals) # Set Y-limits
plt.show() # Display plot
# FIGURE: Pressure coefficient contour
if (flagPlot[5] == 1):
fig = plt.figure(6) # Create figure
plt.cla() # Get ready for plotting
plt.contourf(XX,YY,CpXY,500,cmap='jet') # Plot contour
plt.fill(XB,YB,'k') # Plot airfoil as black polygon
plt.xlabel('X Units') # Set X-label
plt.ylabel('Y Units') # Set Y-label
plt.gca().set_aspect('equal') # Set axes equal
plt.xlim(xVals) # Set X-limits
plt.ylim(yVals) # Set Y-limits
plt.show() # Display plot