forked from jte0419/Panel_Methods
-
Notifications
You must be signed in to change notification settings - Fork 0
/
SP_Circle.m
320 lines (276 loc) · 20.5 KB
/
SP_Circle.m
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
% SOURCE PANEL METHOD - CIRCLE
% Written by: JoshTheEngineer
% YouTube : www.youtube.com/joshtheengineer
% Website : www.joshtheengineer.com
% Started: 01/01/19
% Updated: 01/01/19 - Started code
% - Works as expected
% 01/11/19 - Added streamline plotting
% 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:
% - COMPUTE_IJ_SPM.m
% - STREAMLINE_SPM.m
%
% References
% - [1]: Panel Method Geometry
% Link: https://www.youtube.com/watch?v=kIqxbd937PI
% - [2]: Normal Geometric Integral SPM, I(ij)
% Link: https://www.youtube.com/watch?v=76vPudNET6U
% - [3]: Tangential Geometric Integral SPM, J(ij)
% Link: https://www.youtube.com/watch?v=JRHnOsueic8
% - [4]: Streamline Geometric Integral SPM, Mx(ij) and My(ij)
% Link: https://www.youtube.com/watch?v=BnPZjGCatcg
% - [5]: Solving the System of Equations
% Link: https://www.youtube.com/watch?v=ep7vPzGYsbw
clear;
clc;
%% KNOWNS
% User-defined knowns
Vinf = 1; % Freestream velocity
AoA = 0; % Angle of attack [deg]
numB = 67; % Number of boundary points (including endpoint)
tO = (360/(numB-1))/2; % Boundary point angle offset [deg]
% Plotting flags
flagPlot = [1; % Shape polygon with panel normal vectors
1; % Geometry boundary pts, control pts, first panel, second panel
1; % Analytical and SPM pressure coefficient plot
1; % Streamlines
1]; % Pressure coefficient contours
%% CREATE CIRCLE BOUNDARY POINTS
% Angles used to compute boundary points
theta = linspace(0,360,numB)'; % Create angles for computing boundary point locations [deg]
theta = theta + tO; % Add panel angle offset [deg]
theta = theta*(pi/180); % Convert from degrees to radians [rad]
% Boundary points
XB = cos(theta); % Compute boundary point X-coordinate (radius of 1)
YB = sin(theta); % Compute boundary point Y-coordinate (radius of 1)
% Number of panels
numPan = length(XB)-1; % Number of panels (control points)
%% CHECK PANEL DIRECTIONS - FLIP IF NECESSARY
% Check for direction of points
edge = zeros(numPan,1); % Initialize edge value array
for i = 1:1:numPan % Loop over all panels
edge(i) = (XB(i+1)-XB(i))*(YB(i+1)+YB(i)); % Compute edge value
end
sumEdge = sum(edge); % Sum of all edge values
% If panels are CCW, flip them (don't if CW)
if (sumEdge < 0) % If panels are CCW
XB = flipud(XB); % Flip the X-data array
YB = flipud(YB); % Flip the Y-data array
end
%% PANEL METHOD GEOMETRY - REF [1]
% Initialize variables
XC = zeros(numPan,1); % Initialize control point X-coordinate array
YC = zeros(numPan,1); % Initialize control point Y-coordinate array
S = zeros(numPan,1); % Initialize panel length array
phiD = zeros(numPan,1); % Initialize panel orientation angle array [deg]
% Find geometric quantities of airfoil
for i = 1:1: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
phiD(i) = atan2d(dy,dx); % Angle of the panel (positive X-axis to inside face)
if (phiD(i) < 0) % Make all panel angles positive
phiD(i) = phiD(i) + 360;
end
end
% Compute angle of panel normal w.r.t horizontal and include AoA
deltaD = phiD + 90; % Angle from positive X-axis to outward normal vector [deg]
betaD = deltaD - AoA; % Angle between freestream vector and outward normal vector [deg]
betaD(betaD > 360) = betaD(betaD > 360) - 360; % Make sure angles aren't greater than 360 [deg]
% Convert angles from [deg] to [rad]
phi = phiD.*(pi/180); % Convert from [deg] to [rad]
beta = betaD.*(pi/180); % Convert from [deg] to [rad]
%% COMPUTE SOURCE PANEL STRENGTHS - REF [5]
% Geometric integral (normal [I] and tangential [J])
% - Refs [2] and [3]
[I,J] = COMPUTE_IJ_SPM(XC,YC,XB,YB,phi,S); % Compute geometric integrals
% Populate A matrix
% - Simpler option: A = I + pi*eye(numPan,numPan);
A = zeros(numPan,numPan); % Initialize the A matrix
for i = 1:1:numPan % Loop over all i panels
for j = 1:1:numPan % Loop over all j panels
if (i == j) % If the panels are the same
A(i,j) = pi; % Set A equal to pi
else % If panels are not the same
A(i,j) = I(i,j); % Set A equal to geometric integral
end
end
end
% Populate b array
% - Simpler option: b = -Vinf*2*pi*cos(beta);
b = zeros(numPan,1); % Initialize the b array
for i = 1:1:numPan % Loop over all panels
b(i) = -Vinf*2*pi*cos(beta(i)); % Compute RHS array
end
% Compute source panel strengths (lambda) from system of equations
lambda = A\b; % Compute all source strength values
% Check the sum of the source strenghts
% - This should be very close to zero for a closed polygon
fprintf('Sum of L: %g\n',sum(lambda.*S)); % Print sum of all source strengths
%% COMPUTE PANEL VELOCITIES AND PRESSURE COEFFICIENTS
% Compute velocities
% - Simpler method: Vt = Vinf*sin(beta) + J*lambda/(2*pi);
% Cp = 1-(Vt/Vinf).^2;
Vt = zeros(numPan,1); % Initialize tangential velocity array
Cp = zeros(numPan,1); % Initialize pressure coefficient array
for i = 1:1:numPan % Loop over all i panels
addVal = 0; % Reset the summation value to zero
for j = 1:1:numPan % Loop over all j panels
addVal = addVal + (lambda(j)/(2*pi))*(J(i,j)); % Sum all tangential source panel terms
end
Vt(i) = Vinf*sin(beta(i)) + addVal; % Compute tangential velocity by adding uniform flow term
Cp(i) = 1-(Vt(i)/Vinf)^2; % Compute pressure coefficient
end
% Analytical angles and pressure coefficients
analyticTheta = linspace(0,2*pi,200)'; % Analytical theta angles [rad]
analyticCP = 1-4*sin(analyticTheta).^2; % Analytical pressure coefficient []
%% COMPUTE LIFT AND DRAG
% Compute normal and axial force coefficients
CN = -Cp.*S.*sin(beta); % Normal force coefficient []
CA = -Cp.*S.*cos(beta); % Axial force coefficient []
% Compute lift and drag coefficients
CL = sum(CN.*cosd(AoA)) - sum(CA.*sind(AoA)); % Decompose axial and normal to lift coefficient []
CD = sum(CN.*sind(AoA)) + sum(CA.*cosd(AoA)); % Decompose axial and normal to drag coefficient []
% Display lift and drag coefficients in command window
fprintf('CL : %g\n',CL); % Display lift coefficient (should be zero)
fprintf('CD : %g\n',CD); % Display drag coefficient (should be zero)
%% COMPUTE STREAMLINES - REF [4]
if (flagPlot(4) == 1 || flagPlot(5) == 1)
% Grid parameters
nGridX = 100; % X-grid for streamlines and contours
nGridY = 100; % Y-grid for streamlines and contours
xVals = [-1.5; 1.5]; % X-grid extents [min, max]
yVals = [-1.5; 1.5]; % Y-grid extents [min, max]
% Streamline parameters
stepsize = 0.01; % Step size for streamline propagation
maxVert = nGridX*nGridY*10; % Maximum vertices
slPct = 30; % Percentage of streamlines of the grid
Ysl = linspace(yVals(1),yVals(2),floor((slPct/100)*nGridY))'; % Create array of Y streamline starting points
% Generate the grid points
Xgrid = linspace(xVals(1),xVals(2),nGridX)'; % X-values in evenly spaced grid
Ygrid = linspace(yVals(1),yVals(2),nGridY)'; % Y-values in evenly spaced grid
[XX,YY] = meshgrid(Xgrid,Ygrid); % Create meshgrid from X and Y grid arrays
% Initialize velocities
Vx = zeros(nGridX,nGridY); % Initialize X velocity matrix
Vy = zeros(nGridX,nGridY); % Initialize Y velocity matrix
% Solve for grid point X and Y velocities
for m = 1:1:nGridX % Loop over X-grid points
for n = 1:1:nGridY % Loop over Y-grid points
XP = XX(m,n); % Isolate X point
YP = YY(m,n); % Isolate Y point
[Mx,My] = STREAMLINE_SPM(XP,YP,XB,YB,phi,S); % Compute streamline Mx and My values (Ref [4])
% Check if grid points are in object
% - If they are, assign a velocity of zero
[in,on] = inpolygon(XP,YP,XB,YB); % Find whether (XP,YP) is in polygon body
if (in == 1 || on == 1) % If (XP, YP) is in the polygon body
Vx(m,n) = 0; % X-velocity is zero
Vy(m,n) = 0; % Y-velocity is zero
else % If (XP,YP) is not in the polygon body
Vx(m,n) = Vinf*cosd(AoA) + sum(lambda.*Mx./(2*pi)); % Compute X-velocity
Vy(m,n) = Vinf*sind(AoA) + sum(lambda.*My./(2*pi)); % Compute Y-velocity
end
end
end
% Compute grid point velocity magnitude and pressure coefficient
Vxy = sqrt(Vx.^2 + Vy.^2); % Compute magnitude of velocity vector
CpXY = 1-(Vxy./Vinf).^2; % Pressure coefficient []
end
%% PLOTTING
% FIGURE: Shape polygon with panel normal vectors
if (flagPlot(1) == 1)
figure(1); % Create figure
cla; hold on; grid off; % Get ready for plotting
set(gcf,'Color','White'); % Set color to white
set(gca,'FontSize',12); % Set font size
fill(XB,YB,'k'); % Plot polygon
for i = 1:1:numPan % Loop over all panels
X(1) = XC(i); % Set X start of panel orientation vector
X(2) = XC(i) + S(i)*cosd(betaD(i)+AoA); % Set X end of panel orientation vector
Y(1) = YC(i); % Set Y start of panel orientation vector
Y(2) = YC(i) + S(i)*sind(betaD(i)+AoA); % Set Y end of panel orientation vector
plot(X,Y,'r-','LineWidth',3); % Plot panel normal vector
end
xlabel('X Units'); % Set X-label
ylabel('Y Units'); % Set Y-label
title('Shape andd normal vectors'); % Set title
axis equal; % Set axes equal
zoom reset; % Reset zoom
end
% FIGURE: Geometry with the following indicated:
% - Boundary points, control points, first panel, second panel
if (flagPlot(2) == 1)
figure(2); % Create figure
cla; hold on; grid on; % Get ready for plotting
set(gcf,'Color','White'); % Set color to white
set(gca,'FontSize',12); % Set font size
plot(XB,YB,'k-','LineWidth',3); % Plot panels
p1 = plot([XB(1) XB(2)],[YB(1) YB(2)],'g-','LineWidth',3); % Plot first panel
p2 = plot([XB(2) XB(3)],[YB(2) YB(3)],'m-','LineWidth',3); % Plot second panel
pB = plot(XB,YB,'ko','MarkerFaceColor','k','MarkerSize',10); % Plot boundary points
pC = plot(XC,YC,'ko','MarkerFaceColor','r','MarkerSize',10); % Plot control points
legend([pB,pC,p1,p2],... % Show legend
{'Boundary','Control','First Panel','Second Panel'});
ylim([-1 1]); % Set Y-limits
axis equal; % Set axes equal
xlabel('X Units'); % Set X-label
ylabel('Y Units'); % Set Y-label
title('Geometry and panels'); % Set title
zoom reset; % Reset zoom
end
% FIGURE: Analytical and SPM pressure coefficient plot
if (flagPlot(3) == 1)
figure(3); % Create figure
cla; hold on; grid on; % Get ready for plotting
set(gcf,'Color','White'); % Set color to white
set(gca,'FontSize',12); % Set font size
pA = plot(analyticTheta,analyticCP,'k-','LineWidth',3); % Plot analytical pressure coefficient
pC = plot(beta,Cp,'ks','MarkerFaceColor','r','MarkerSize',10); % Plot compute pressure coefficient
xlabel('Angle [rad]'); % Set X-label
ylabel('Cp'); % Set Y-label
xlim([0 2*pi]); % Set X-limits
ylim([-3.5 1.5]); % Set Y-limits
title('Analytical and pressure coefficient plot'); % Set title
legend([pA,pC],{'Analytical - p. 230 Libro Fluido','SPM'},'Location','S'); % Add legend
end
% FIGURE: Streamlines (and quiver if commented in)
if (flagPlot(4) == 1)
figure(4); % Create figure
cla; hold on; grid off; % Get ready for plotting
set(gcf,'Color','White'); % Set color to white
set(gca,'FontSize',12); % Set font size
axis equal; % Set axes equal
xlabel('X Units'); % Set X-label
ylabel('Y Units'); % Set Y-label
% quiver(XX,YY,Vx,Vy,'r'); % Plot velocity vectors
for i = 1:1:length(Ysl) % Loop over all streamlines
sl = streamline(XX,YY,Vx,Vy,xVals(1),Ysl(i),[stepsize,maxVert]); % streamline(X,Y,U,V,startx,starty)
set(sl,'LineWidth',3); % Change streamline width
end
fill(XB,YB,'k'); % Plot polygon
xlim(xVals); % Set X-limits
ylim(yVals); % Set Y-limits
title('Streamlines'); % Set title
zoom reset; % Reset zoom
end
% FIGURE: Pressure coefficient contours
if (flagPlot(5) == 1)
figure(5); % Create figure
cla; hold on; grid on; % Get ready for plotting
set(gcf,'Color','White'); % Set color to white
set(gca,'FontSize',12); % Set font size
axis equal; % Set axes equal
xlabel('X Units'); % Set X-label
ylabel('Y Units'); % Set Y-label
contourf(XX,YY,CpXY,100,'EdgeColor','none'); % Plot contour
fill(XB,YB,'k'); % Plot polygon
xlim(xVals); % Set X-limits
ylim(yVals); % Set Y-limits
title('Pressure coefficient contours'); % Set title
colorbar
zoom reset; % Reset zoom
end