-
Notifications
You must be signed in to change notification settings - Fork 42
/
Panel_Method_Geometry.py
109 lines (92 loc) · 7.91 KB
/
Panel_Method_Geometry.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
# PANEL METHOD GEOMETRY
# Written by: JoshTheEngineer
# Started: 01/24/19
# Updated: 01/24/19 - Started code in MATLAB
# - Works as expected
# 02/03/19 - Transferred code to Python
# - Works as expected
import numpy as np
import math as m
import matplotlib.pyplot as plt
from LOAD_AIRFOIL_SELIG import LOAD_AIRFOIL_SELIG
# %% CREATE/LOAD GEOMETRY
# Load circle or airfoil
#load = 'Airfoil' # Use the airfoil geometry
load = 'Circle' # Use the circle geometry
AoA = 0 # Angle of attack [deg]
if (load == 'Airfoil'): # If airfoil is to be loaded
fileName = 'goe623' # Specify the airfoil filename
XB, YB = LOAD_AIRFOIL_SELIG(fileName) # Load X and Y data points of airfoil
elif (load == 'Circle'): # If circle is to be used
numB = 9 # Number of boundary points
tO = 22.5 # Angle offset [deg]
# Angles used to compute boundary points
theta = np.linspace(0,360,numB) # Angles to compute boundary points [deg]
theta = theta + tO # Add angle offset [deg]
theta = theta*(np.pi/180) # Convert angle to radians [rad]
# Boundary points
XB = np.cos(theta) # X value of boundary points
YB = np.sin(theta) # Y value of boundary points
# Number of panels
numPan =len(XB)-1 # Number of panels
# Check for direction of points
edge = np.zeros(numPan) # Initialize edge check value
for i in range(numPan): # Loop over all panels
edge[i] = (XB[i+1]-XB[i])*(YB[i+1]+YB[i]) # Compute edge value for each panel
sumEdge = np.sum(edge) # Sum all panel edge values
# If panels are CCW, flip them (don't if CW)
if (sumEdge < 0): # If sum is negative
print('Points are counter-clockwise. Flipping.\n') # Display message in console
XB = np.flipud(XB) # Flip the X boundary points array
YB = np.flipud(YB) # Flip the Y boundary points array
elif (sumEdge > 0): # If sum is positive
print('Points are clockwise. Not flipping.\n') # Do nothing, display message in consolve
# %% COMPUTE GEOMETRIC VARIABLES
# Initialize variables
XC = np.zeros(numPan) # Initialize X control points
YC = np.zeros(numPan) # Initialize Y control points
S = np.zeros(numPan) # Initialize panel lengths
phi = np.zeros(numPan) # Initialize panel orientation angles
# 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 control point coordinate
YC[i] = 0.5*(YB[i]+YB[i+1]) # Y control point coordinate
dx = XB[i+1]-XB[i] # Panel X length
dy = YB[i+1]-YB[i] # Panel Y length
S[i] = (dx**2 + dy**2)**0.5 # Panel length
phi[i] = m.atan2(dy,dx) # Panel orientation angle [rad]
if (phi[i] < 0): # If panel orientation is negative
phi[i] = phi[i] + 2*np.pi # Add 2pi to the panel angle
# Compute angle of panel normal w.r.t. horizontal and include AoA
delta = phi + (np.pi/2) # Compute panel normal angle [rad]
beta = delta - (AoA*(np.pi/180)) # Angle between freestream and panel normal [rad]
# %% PLOTTING
# Dashed circle defined
T = np.linspace(0,2*np.pi,1000) # Angle array to compute circle
x = np.cos(T) # Circle X points
y = np.sin(T) # Circle Y points
# Plot the paneled geometry
fig = plt.figure(1) # Create figure
plt.cla() # Get ready for plotting
if (load == 'Circle'): # If circle is selected
plt.plot(x,y,'k--') # Plot actual circle outline
plt.fill(XB,YB,'k') # Plot polygon (circle or airfoil)
X = np.zeros(2) # Initialize panel X variable
Y = np.zeros(2) # Initialize panel Y variable
for i in range(numPan): # Loop over all panels
X[0] = XC[i] # Panel starting X point
X[1] = XC[i] + S[i]*np.cos(delta[i]) # Panel ending X point
Y[0] = YC[i] # Panel starting Y point
Y[1] = YC[i] + S[i]*np.sin(delta[i]) # Panel ending Y point
if (i == 0): # For first panel
plt.plot(X,Y,'b-',label='First Panel') # Plot the first panel normal vector
elif (i == 1): # For second panel
plt.plot(X,Y,'g-',label='Second Panel') # Plot the second panel normal vector
else: # For every other panel
plt.plot(X,Y,'r-') # Plot the panel normal vector
plt.xlabel('X-Axis') # Set X-label
plt.ylabel('Y-Axis') # Set Y-label
plt.title('Panel Geometry') # Set title
plt.axis('equal') # Set axes equal
plt.legend() # Plot legend
plt.show() # Display plot