-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathmain.m
226 lines (188 loc) · 7.29 KB
/
main.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
%
% ////// ///// //// ///// // // /////
% // // // // // // // // // //
% // // // // // // // //
% // //// // // // //// // // /////
% // // // // // // // // // //
% ///// // //// //// ////// // /////
%
%
% This program computes the incompressible flow around an airfoil
% with and without uniform blowing
%
% ------------------------------------------------------------------------------
% If you use this code and find it helpful please cite
%
% M.Reder, A. Stroh and D. Gatti, "Preliminary study of flow control via uniform
% blowing on airfoils with a boundary element method", Notes on Numerical Fluid
% Mechanics and Multidisciplinary Design, (submitted, 2018)
% ------------------------------------------------------------------------------
%
% Copyright (C) 2018 Dr. Davide Gatti & M. Reder
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <https://www.gnu.org/licenses/>.
%
% Contacts:
%
% Dr. Davide Gatti
% davide.gatti [at] kit.edu
% msc.davide.gatti [at] gmail.com
%
% Karlsruhe Institute of Technology
% Institute of Fluid Dynamics
% Kaiserstraße 10
% 76131 Karlsruhe
%
init % initialize: add paths etc.
parameters % load parameters for calculation
% inviscid solution
[InvRef.prf,InvRef.flo,InvRef.CoeffMatrix,InvRef.Uinv,InvRef.sges ] = InviscidSolution(prf,flo,eng);
% InvRef is a struct containing all information from inviscid solution
%-> can be committed by next use of airfoil function as input to get shorter calculation process
[sol,prf,flo,~,~,~,~]=airfoil(prf,flo,tri,blo,eng,InvRef);
PlotStuff(prf,flo.wake,sol, 'Cp');
PlotStuff(prf,flo.wake,sol, 'tau');
PlotProfile( prf, flo.wake,sol,flo,1)
%%
%---------------------------------------------------------
% set reference case without uniform blowing
CL_ref=sol.CL; CD_ref=sol.Cdrag; tr_ref=sol.tran.x;ratio_ref=CL_ref/CD_ref;
% use tripping at transition location of uncontrolled case
% ->neglect effect of blowing on transition, uncomment to consider transition changes
solRef=sol;prfRef=prf;
%tri.active=[true;true];
%tri.x=solRef.tran.x;
%flo.nkrit=9;
% turn on uniform blowing
blo.active=true;
[solB,prfB,flo,~,~,~,~]=airfoil(prf,flo,tri,blo,eng,InvRef);
% compare
BlowingComparison(prf,flo.wake,sol,flo,prfB,solB,1,'delta');
PlotProfile( prfB, flo.wake,solB,flo,1)
PlotStuff(prfB,flo.wake,solB, 'Cp');
PlotStuff(prfB,flo.wake,solB, 'tau');
% %% Plotting Functions avaivable
%
% % % plot quantities
% %------------------------------------------------------------------
%
% % quantities that can be plotted:
% %'tau','Cf','Cfint','Cp','delta' ,'U', 'CD','D','H12' ,'H32', 'Ret'
% % section: 1 - suction side, 2 - pressure side, 3 - wake, 4 - suction and pressure side , 5 - all sections
% % OverArclength=true -> plot over arclength s instead of x
% % section=4;
% % OverArclength=false;
% % Quantity='tau';
% % PlotStuff(prf,flo.wake,sol, Quantity,section,OverArclength);
% %------------------------------------------------------------------
%
% % default: section = 4, OverArclength=false;
% PlotStuff(prf,flo.wake,sol, 'tau');
%
%
% % Overview-plots of profiles
% %------------------------------------------------------------------
% % mode 1: plots Profil and wake with displacement thickness and shows CL, Cd usw.
% % mode 2: plots Profil and wake with the node positions used
% % mode 3: plots Profil with transition locations and blowing distribution in blowing case
% %------------------------------------------------------------------
% mode=1;
% PlotProfile(prf,flo.wake,sol, mode);
%
%
% % Comparison between blowing case and reference case without blowing
% %------------------------------------------------------------------
% % section 0: only plots overview with reduction in CL and Cd
% % section 1: plots overview + plots for suction side
% % section 2: plots overview + plots for pressure side
% % Quantities that can be plotted: 'Cf','Cp','delta','q','H12','H32','CD','D','Ret'
% %------------------------------------------------------------------
% section=1;Quantity='Cf';
% BlowingComparison(prf,flo.wake,sol,prfB,solB,section,Quantity);
%
% % Plot line forces
% %------------------------------------------------------------------
% % mode=0: sum of lower and upper airfoil part
% % mode=1: upper airfoil part
% % mode=2: lower airfoil part
% % if sol2 comitted -> plot both solutions
% %------------------------------------------------------------------
% PlotForces(prf,flo,sol,mode, solRef)
%
% % Plot stress/force vectors
% %------------------------------------------------------------------
% % each_nth: only plot each n-th node (default 1 -> every node)
% % total =true: plot the total stresses consisting of pressure and shear stresses
% % pressure=true: plot the pressure stresses
% % shear =true: plot the shear stresses
% % force = true : plots force instead of stress vectors (integrated over area)
% %------------------------------------------------------------------
% each_nth=3,total=true; pressure=false; shear=false; force=false;
% PlotStressVectors( prf,sol,flo,each_nth,total, pressure, shear, force )
return
%% Calculation of polar curves
%---------------------------------------------------------
% loop over alfas + write out for polar curves
step=1;start=0;
for k=start:step:12
flo.alfa= k *pi/180;
disp(['ALFA=',num2str(k)])
disp('--------------------------------------------')
% evaluation
[sol,prf,flo,~,~,~]=airfoil(prf,flo,tri,blo,eng);
if k==start
CP=sol.Cp(1:prf.N);
TAU=sol.tau(1:prf.N);
CL =sol.CL ;
Cnu =sol.Cnu;
Cdrag=sol.Cdrag;
Cd_p=sol.Cdrag-sol.Cnu;
alph=k;
xTT_up= [sol.tran.x(1), sol.xseparation(1), sol.xreattach(1) ];
xTT_lo= [sol.tran.x(2), sol.xseparation(2), sol.xreattach(2) ];
else
CP=[CP,sol.Cp(1:prf.N)];
TAU=[TAU,sol.tau(1:prf.N)];
CL=[CL;sol.CL];
Cnu=[Cnu;sol.Cnu];
Cdrag=[Cdrag;sol.Cdrag];
Cd_p=[Cd_p;sol.Cdrag-sol.Cnu];
alph=[alph;k];
xTT_up=[xTT_up;sol.tran.x(1), sol.xseparation(1), sol.xreattach(1) ];
xTT_lo=[xTT_lo;sol.tran.x(2), sol.xseparation(2), sol.xreattach(2) ];
end
end
%% Even distributed paremeters for blowing location
%
% N=20;
% xTop=linspace(-0.1,1,N);
% xBot=linspace(-0.1,1,N);
%
% tri.x=tr_ref;
% tri.active=[true; true];
% blo.active=true;
%
% for i=1:N
% for j=1:N
% blo.x= {[xTop(i)]*prf.c; % midpoint of blowing area
% [xBot(j)]*prf.c;};
% [sol,~,~,~,~,~]=airfoil(prf,flo,tri,blo,eng);
% CL(i,j)=sol.CL;
% CD(i,j)=sol.Cdrag;
% end
% end
%
% contourf(xTop,xBot,(real(CL(:,:))./CD)'/(CL_ref/CD_ref))
% xlabel('x^{TOP}')
% ylabel('x^{BOTTOM}')