Skip to content

Commit ac6f853

Browse files
EttoreSaettaEttoreSaetta
committed
OptimalPropeller script added.
1 parent c17aabd commit ac6f853

1 file changed

Lines changed: 371 additions & 0 deletions

File tree

‎SU2_PY/OptimalPropeller.py‎

Lines changed: 371 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,371 @@
1+
# ==============================================================================================
2+
# Name : OptimalPropeller
3+
# Author : Ettore Saetta, Lorenzo Russo, Renato Tognaccini
4+
# Theoretical and Applied Aerodynamic Research Group (TAARG),
5+
# University of Naples Federico II.
6+
# Version : 1.0.0 - Python
7+
# Date : 01/09/2020
8+
# Copyright :
9+
# Description : Compute the optimal load distribution along the propeller radius using
10+
# the inviscid theory of the optimal propeller.
11+
# Reference : Glauert H., Airplane Propellers, in Aerodynamic Theory, Ed. Durand W. F.,
12+
# Vol. IV, pp. 169 - 360, Springer, 1935.
13+
# Input : Interactive.
14+
# Output : ActuatorDisk.cfg, containing part of SU2 .cfg file.
15+
# ActuatorDisk.dat, containing propeller load distribution to be read by SU2_CFD.
16+
# Note : Python 3 or higher needed.
17+
# ==============================================================================================
18+
19+
import math
20+
import numpy as np
21+
import pylab as pl
22+
23+
##########################
24+
### Functions ###
25+
##########################
26+
def a_distribution (w0, Chi):
27+
28+
"""Function used to compute the value of the axial interference factor using the inviscid theory of the optimal propeller."""
29+
30+
a = (w0*pow(Chi,2))/(pow(Chi,2)+pow((1+(w0)),2))
31+
return a
32+
33+
def Print_external_file(CTrs, CPrs):
34+
35+
"""Function used to write the actuator disk input data file"""
36+
file = open('ActuatorDisk.dat', 'w')
37+
file.write('# Automatic generated actuator disk input data file using the Optimal Propeller code.\n')
38+
file.write('# Data file needed for the actuator disk VARIABLE_LOAD type.\n')
39+
file.write('# The load distribution is obtained using the inviscid theory of the optimal propeller\n')
40+
file.write('# using global data.\n')
41+
file.write('#\n')
42+
file.write('# The first three lines must be filled.\n')
43+
file.write('# An example of this file can be found in the TestCases directory.\n')
44+
file.write('#\n')
45+
file.write('# Author: Ettore Saetta, Lorenzo Russo, Renato Tognaccini.\n')
46+
file.write('# Theoretical and Applied Aerodynamic Research Group (TAARG),\n')
47+
file.write('# University of Naples Federico II\n')
48+
file.write('# -------------------------------------------------------------------------------------\n')
49+
file.write('#\n')
50+
file.write('MARKER_ACTDISK= \n')
51+
file.write('CENTER= \n')
52+
file.write('AXIS= \n')
53+
file.write('RADIUS= '+str(R)+'\n')
54+
file.write('ADV_RATIO= '+str(J)+'\n')
55+
file.write('NROW= '+str(Stations)+'\n')
56+
file.write('# rs=r/R dCT/drs dCP/drs dCR/drs\n')
57+
58+
for i in range(0, Stations):
59+
file.write(f' {r[i]:.7f} {CTrs[i]:.7f} {CPrs[i]:.7f} 0.0\n')
60+
61+
file.close()
62+
63+
##########################
64+
### Main ###
65+
##########################
66+
# Screen output
67+
print('------------------ Optimal Propeller vsn 1.0.0 ------------------')
68+
print('| Computation of the optimal dCT/dr and dCP/dr distributions. |')
69+
print('| Based on the inviscid theory of the optimal propeller. |')
70+
print('| |')
71+
print('| This code is used to generate the actuator disk input data |')
72+
print('| file needed for the VARIABLE_LOAD actuator disk type |')
73+
print('| implemented in SU2 7.0.7. |')
74+
print('| |')
75+
print('| Author: Ettore Saetta, Lorenzo Russo, Renato Tognaccini. |')
76+
print('| Theoretical and Applied Aerodynamic Research Group (TAARG), |')
77+
print('| University of Naples Federico II. |')
78+
print('-----------------------------------------------------------------')
79+
80+
print('')
81+
print('Warning: present version requires input in SI units.')
82+
print('')
83+
84+
# Number of radial stations in input.
85+
Stations = int(input('Number of radial stations: '))
86+
print('')
87+
88+
dStations = float(Stations)
89+
90+
# Resize the vectors using the number of radial stations.
91+
r = [None]*Stations
92+
chi = [None]*Stations
93+
dCp = [None]*Stations
94+
w = [None]*Stations
95+
a_new = [None]*Stations
96+
a_old = [None]*Stations
97+
a_0 = [None]*Stations
98+
ap_old = [None]*Stations
99+
dCt_new = [None]*Stations
100+
dCt_old = [None]*Stations
101+
dCt_0 = [None]*Stations
102+
DeltaP = [None]*Stations
103+
F = [None]*Stations
104+
a_optimal = [None]*Stations
105+
ap_optimal = [None]*Stations
106+
dCt_optimal = [None]*Stations
107+
108+
# Thrust coefficient in input.
109+
Ct = float(input('CT (Renard definition): '))
110+
print('')
111+
112+
# Propeller radius in input.
113+
R = float(input('R (propeller radius [m]): '))
114+
print('')
115+
116+
# Hub radius in input.
117+
rhub = float(input('r_hub (hub radius [m]): '))
118+
print('')
119+
120+
# Advance ratio in input.
121+
J = float(input('J (advance ratio): '))
122+
print('')
123+
124+
# Freestream velocity in input.
125+
Vinf = float(input('Vinf (m/s): '))
126+
print('')
127+
128+
# Asking if the tip loss Prandtl correction function needs to be used.
129+
Prandtl = input('Using tip loss Prandtl correction? (<y>/n): ')
130+
print('')
131+
132+
if Prandtl == 'y' or Prandtl == 'Y' or Prandtl == '':
133+
# Number of propeller blades in input.
134+
N = int(input('N (number of propeller blades): '))
135+
print('')
136+
137+
corr = True
138+
139+
else:
140+
corr = False
141+
142+
# Computation of the non-dimensional hub radius.
143+
rs_hub = rhub/R
144+
145+
# Computation of the non-dimensional radial stations.
146+
for i in range(1,Stations+1):
147+
r[i-1]=i/dStations
148+
if r[i-1] <= rs_hub:
149+
i_hub = i-1
150+
151+
# Computation of the propeller diameter.
152+
D = 2*R
153+
# Computation of the propeller angular velocity (Rounds/s).
154+
n = Vinf/(D*J)
155+
# Computation of the propeller angular velocity (Rad/s).
156+
Omega = n*2*math.pi
157+
158+
# Computation of the tip loss Prandtl correction function F.
159+
if corr == True:
160+
for i in range(0, Stations):
161+
F[i] = (2/math.pi)*math.acos(math.exp(-0.5*N*(1-r[i])*math.sqrt(1+pow(Omega*R/Vinf,2))))
162+
163+
else:
164+
for i in range(0, Stations):
165+
F[i] = 1.0
166+
167+
# Computation of the non-dimensional radius chi=Omega*r/Vinf.
168+
for i in range(0, Stations):
169+
chi[i] = Omega*r[i]*R/Vinf
170+
171+
eps = 5E-20
172+
# Computation of the propeller radial stations spacing.
173+
h = (1.0/Stations)
174+
175+
# Computation of the first try induced velocity distribution.
176+
for i in range(0, Stations):
177+
w[i] = (2/math.pow(Vinf,2))*((-1/Vinf)+math.sqrt(1+((math.pow(D,4)*(Ct)*math.pow(n,2))/(math.pow(Vinf,2)*math.pi*r[i]))))
178+
179+
# Computation of the first try Lagrange moltiplicator.
180+
w_0 = 0.0
181+
for i in range(0, Stations):
182+
w_0 += w[i]
183+
184+
w_0 = w_0/(Vinf*Stations)
185+
186+
# Computation of the first try axial interference factor distribution.
187+
for i in range(0, Stations):
188+
a_0[i] = a_distribution(w_0*F[i],chi[i])
189+
190+
# Computation of the thrust coefficient distribution
191+
for i in range(0, Stations):
192+
dCt_0[i]= math.pi*J*J*r[i]*(1+a_0[i])*a_0[i]
193+
194+
# Computation of the total thrust coefficient.
195+
Ct_0 = 0.0
196+
for i in range(i_hub, Stations):
197+
Ct_0 += h*dCt_0[i]
198+
199+
# Compute the error with respect to the thrust coefficient given in input.
200+
err_0 = Ct_0 - Ct
201+
print('CONVERGENCE HISTORY:')
202+
print(err_0)
203+
204+
# Computation of the second try Lagrange moltiplicator.
205+
w_old = w_0 + 0.1
206+
207+
# Computation of the second try axial interference factor distribution.
208+
for i in range(0, Stations):
209+
a_old[i] = a_distribution(w_old*F[i],chi[i])
210+
211+
# Computation of the thrust coefficient distribution
212+
for i in range(0, Stations):
213+
dCt_old[i]= math.pi*J*J*r[i]*(1+a_old[i])*a_old[i]
214+
215+
# Computation of the total thrust coefficient.
216+
Ct_old = 0.0
217+
for i in range(i_hub, Stations):
218+
Ct_old += h*dCt_old[i]
219+
220+
# Compute the error with respect to the thrust coefficient given in input.
221+
err_old = Ct_old - Ct
222+
print(err_old)
223+
224+
##########################
225+
### Iterations ###
226+
##########################
227+
# Iterate using the false position methods.
228+
# Based on the error from the thrust coefficient given in input.
229+
iteration = 2
230+
err_new = err_old
231+
while math.fabs(err_new) >= eps and err_0 != err_old:
232+
233+
iteration += 1
234+
235+
# Computation of the new Lagrange moltiplicator value based on the false position method.
236+
w_new = (w_old*err_0 - w_0*err_old)/(err_0 - err_old)
237+
238+
# Computation of the new axial interference factor distribution.
239+
for i in range(0, Stations):
240+
a_new[i] = a_distribution(w_new*F[i],chi[i])
241+
242+
# Computation of the new thrust coefficient distribution.
243+
for i in range(0, Stations):
244+
dCt_new[i]= math.pi*J*J*r[i]*(1+a_new[i])*a_new[i]
245+
246+
# Computation of the new total thrust coefficient.
247+
Ct_new = 0.0
248+
for i in range(i_hub, Stations):
249+
Ct_new += h*dCt_new[i]
250+
251+
# Computation of the total thrust coefficient error with respect to the input value.
252+
err_new = Ct_new - Ct
253+
254+
print(err_new)
255+
256+
# Updating the stored values for the next iteration.
257+
err_0 = err_old
258+
err_old = err_new
259+
260+
w_0 = w_old
261+
w_old = w_new
262+
263+
# Computation of the correct axial and rotational interference factors (a and ap).
264+
for i in range(0, Stations):
265+
a_optimal[i] = a_distribution(w_new*F[i],chi[i])
266+
ap_optimal[i] = (w_new*F[i])*((1+w_new*F[i])/(chi[i]*chi[i]+math.pow(1+w_new*F[i],2)))
267+
268+
# Computation of the correct thrust coefficient distribution.
269+
for i in range(0, Stations):
270+
dCt_optimal[i] = math.pi*J*J*r[i]*(1+a_optimal[i])*a_optimal[i]
271+
272+
# Computation of the correct power coefficient distribution.
273+
for i in range(0, Stations):
274+
dCp[i] = (R*4*math.pi/(math.pow(n,3)*math.pow(D,5)))*(math.pow(Vinf,3)*math.pow(1+a_optimal[i],2)*a_optimal[i]*r[i]*R+math.pow(Omega,2)*Vinf*(1+a_optimal[i])*math.pow(ap_optimal[i],2)*math.pow(r[i]*R,3))
275+
276+
##########################
277+
### Check Results ###
278+
##########################
279+
# Computation of the total power coefficient.
280+
Cp = 0.0
281+
for i in range(i_hub, Stations):
282+
Cp += h*dCp[i]
283+
284+
# Computation of the total thrust coefficient.
285+
Ct_optimal = 0.0
286+
for i in range(i_hub, Stations):
287+
Ct_optimal += h*dCt_optimal[i]
288+
289+
# Computation of the static pressure jump distribution.
290+
for i in range(0, Stations):
291+
DeltaP[i] = (dCt_optimal[i])*(2*Vinf*Vinf)/(J*J*math.pi*r[i])
292+
293+
# Computation of the thrust over density (T) using the static pressure jump distribution.
294+
T = 0.0
295+
for i in range(i_hub, Stations):
296+
T += 2*math.pi*r[i]*math.pow(R,2)*h*DeltaP[i]
297+
298+
# Computation of the thrust coefficient using T.
299+
Ct_Renard = (T)/(math.pow(n,2)*math.pow(D,4))
300+
301+
# Computation of the efficiency.
302+
eta = J*(Ct_optimal/Cp)
303+
304+
# Screen output used to check that everything worked correcty.
305+
print('%%%%%%%%%%%%%%%%%%%%%%%%% CHECK OUTPUT VALUES %%%%%%%%%%%%%%%%%%%%%%%%%')
306+
print(f' dCT distribution integral: {Ct_optimal:.4f}')
307+
print(f' dCT computed using the static pressure jump: {Ct_Renard:.4f}')
308+
print(f' dCP distribution integral: {Cp:.4f}')
309+
print(f' Thrust over Density (T/rho): {T:.4f} [N*m^3/kg]')
310+
print(f' Efficiency eta: {eta:.4f}')
311+
print(f' w0/Vinf: {w_new:.4f}')
312+
print('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
313+
314+
##########################
315+
### File Writing ###
316+
##########################
317+
# Write the actuator disk configuration file
318+
file = open('ActuatorDisk.cfg', 'w')
319+
320+
file.write('% Automatic generated actuator disk configuration file.\n')
321+
file.write('%\n')
322+
file.write('% The first two elements of MARKER_ACTDISK must be filled.\n')
323+
file.write('% An example of this file can be found in the TestCases directory.\n')
324+
file.write('%\n')
325+
file.write('% Author: Ettore Saetta, Lorenzo Russo, Renato Tognaccini.\n')
326+
file.write('% Theoretical and Applied Aerodynamic Research Group (TAARG),\n')
327+
file.write('% University of Naples Federico II\n')
328+
file.write('\n')
329+
file.write('ACTDISK_TYPE = VARIABLE_LOAD\n')
330+
file.write('ACTDISK_FILENAME = ActuatorDisk.dat\n')
331+
file.write('MARKER_ACTDISK = ( , , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)\n')
332+
file.close()
333+
334+
print('SU2 file generated!')
335+
336+
# Write the actuator disk data file.
337+
# This is the actuator disk input data file.
338+
Print_external_file(dCt_optimal, dCp)
339+
340+
##########################
341+
### Plots ###
342+
##########################
343+
# Automatically plot the computed propeller performance.
344+
345+
f1 = pl.figure(1)
346+
pl.plot(r, dCt_optimal, 'r', markersize=4, label='$\\frac{dCT}{d\overline{r}}$')
347+
pl.plot(r, dCp, 'k', markersize=4, label='$\\frac{dCP}{d\overline{r}}$')
348+
pl.grid(True)
349+
pl.legend(numpoints=3)
350+
pl.xlabel('$\overline{r}$')
351+
pl.ylabel('')
352+
pl.title("Load Distribution")
353+
354+
f1 = pl.figure(2)
355+
pl.plot(chi, a_optimal, 'r', markersize=4, label='$a$')
356+
pl.plot(chi, ap_optimal, 'k', markersize=4, label='$a^1$')
357+
pl.grid(True)
358+
pl.legend(numpoints=3)
359+
pl.xlabel('$\chi$')
360+
pl.ylabel('')
361+
pl.title("Interference Factors")
362+
363+
if corr == True:
364+
f1 = pl.figure(3)
365+
pl.plot(r, F, 'k', markersize=4)
366+
pl.grid(True)
367+
pl.xlabel('$\overline{r}$')
368+
pl.ylabel('$F(\overline{r})$')
369+
pl.title("Tip Loss Prandtl Correction Function")
370+
371+
pl.show()

0 commit comments

Comments
 (0)