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
|
#
# This demo shows how to use bezier splines to define NACA four
# series airfoils and complex variables to define Joukowski
# Airfoils. It will be expanded after overplotting in implemented
# to plot Coefficient of Pressure as well.
# Alex Woo, Dec. 1992
#
# The definitions below follows: "Bezier presentation of airfoils",
# by Wolfgang Boehm, Computer Aided Geometric Design 4 (1987) pp 17-22.
#
# Gershon Elber, Nov. 1992
#
# m = percent camber
# p = percent chord with maximum camber
print "NACA four series airfoils by bezier splines"
print "Will add pressure distribution later with Overplotting"
mm = 0.6
# NACA6xxx
thick = 0.09
# nine percent NACAxx09
pp = 0.4
# NACAx4xx
# Combined this implies NACA6409 airfoil
#
# Airfoil thickness function.
#
set xlabel "NACA6409 -- 9% thick, 40% max camber, 6% camber"
x0 = 0.0
y0 = 0.0
x1 = 0.0
y1 = 0.18556
x2 = 0.03571
y2 = 0.34863
x3 = 0.10714
y3 = 0.48919
x4 = 0.21429
y4 = 0.58214
x5 = 0.35714
y5 = 0.55724
x6 = 0.53571
y6 = 0.44992
x7 = 0.75000
y7 = 0.30281
x8 = 1.00000
y8 = 0.01050
#
# Directly defining the order 8 Bezier basis function for a faster evaluation.
#
bez_d4_i0(x) = (1 - x)**4
bez_d4_i1(x) = 4 * (1 - x)**3 * x
bez_d4_i2(x) = 6 * (1 - x)**2 * x**2
bez_d4_i3(x) = 4 * (1 - x)**1 * x**3
bez_d4_i4(x) = x**4
bez_d8_i0(x) = (1 - x)**8
bez_d8_i1(x) = 8 * (1 - x)**7 * x
bez_d8_i2(x) = 28 * (1 - x)**6 * x**2
bez_d8_i3(x) = 56 * (1 - x)**5 * x**3
bez_d8_i4(x) = 70 * (1 - x)**4 * x**4
bez_d8_i5(x) = 56 * (1 - x)**3 * x**5
bez_d8_i6(x) = 28 * (1 - x)**2 * x**6
bez_d8_i7(x) = 8 * (1 - x) * x**7
bez_d8_i8(x) = x**8
m0 = 0.0
m1 = 0.1
m2 = 0.1
m3 = 0.1
m4 = 0.0
mean_y(t) = m0 * mm * bez_d4_i0(t) + \
m1 * mm * bez_d4_i1(t) + \
m2 * mm * bez_d4_i2(t) + \
m3 * mm * bez_d4_i3(t) + \
m4 * mm * bez_d4_i4(t)
p0 = 0.0
p1 = pp / 2
p2 = pp
p3 = (pp + 1) / 2
p4 = 1.0
mean_x(t) = p0 * bez_d4_i0(t) + \
p1 * bez_d4_i1(t) + \
p2 * bez_d4_i2(t) + \
p3 * bez_d4_i3(t) + \
p4 * bez_d4_i4(t)
z_x(x) = x0 * bez_d8_i0(x) + x1 * bez_d8_i1(x) + x2 * bez_d8_i2(x) + \
x3 * bez_d8_i3(x) + x4 * bez_d8_i4(x) + x5 * bez_d8_i5(x) + \
x6 * bez_d8_i6(x) + x7 * bez_d8_i7(x) + x8 * bez_d8_i8(x)
z_y(x, tk) = \
y0 * tk * bez_d8_i0(x) + y1 * tk * bez_d8_i1(x) + y2 * tk * bez_d8_i2(x) + \
y3 * tk * bez_d8_i3(x) + y4 * tk * bez_d8_i4(x) + y5 * tk * bez_d8_i5(x) + \
y6 * tk * bez_d8_i6(x) + y7 * tk * bez_d8_i7(x) + y8 * tk * bez_d8_i8(x)
#
# Given t value between zero and one, the airfoild curve is defined as
#
# c(t) = mean(t1(t)) +/- z(t2(t)) n(t1(t)),
#
# where n is the unit normal to the mean line. See the above paper for more.
#
# Unfortunately, the parametrization of c(t) is not the same for mean(t1)
# and z(t2). The mean line (and its normal) can assume linear function t1 = t,
# -1
# but the thickness z_y is, in fact, a function of z_x (t). Since it is
# not obvious how to compute this inverse function analytically, we instead
# replace t in c(t) equation above by z_x(t) to get:
#
# c(z_x(t)) = mean(z_x(t)) +/- z(t) n(z_x(t)),
#
# and compute and display this instead. Note we also ignore n(t) and assumes
# n(t) is constant in the y direction,
#
airfoil_y1(t, thick) = mean_y(z_x(t)) + z_y(t, thick)
airfoil_y2(t, thick) = mean_y(z_x(t)) - z_y(t, thick)
airfoil_y(t) = mean_y(z_x(t))
airfoil_x(t) = mean_x(z_x(t))
unset grid
unset zeroaxis
set parametric
set xrange [-0.1:1.1]
set yrange [-0.1:.7]
set trange [ 0.0:1.0]
set title "NACA6409 Airfoil"
plot airfoil_x(t), airfoil_y(t) title "mean line" w l lt 2, \
airfoil_x(t), airfoil_y1(t, thick) title "upper surface" w l lt 1, \
airfoil_x(t), airfoil_y2(t, thick) title "lower surface" w l lt 1
pause -1 "Press Return"
mm = 0.0
pp = .5
thick = .12
set title "NACA0012 Airfoil"
set xlabel "12% thick, no camber -- classical test case"
plot airfoil_x(t), airfoil_y(t) title "mean line" w l lt 2, \
airfoil_x(t), airfoil_y1(t, thick) title "upper surface" w l lt 1, \
airfoil_x(t), airfoil_y2(t, thick) title "lower surface" w l lt 1
pause -1 "Press Return"
set title ""
set xlab ""
set key box
set parametric
set samples 100
set isosamples 10
set style data lines
set style function lines
print "Joukowski Airfoil using Complex Variables"
set title "Joukowski Airfoil using Complex Variables" offset 0,0
#set timestamp
set yrange [-.2 : 1.8]
set trange [0: 2*pi]
set xrange [-.6:.6]
Zeta(t) = -eps + (a+eps)*exp(t*{0,1})
eta(t) = Zeta(t) + a*a/Zeta(t)
eps = 0.06
a =.250
set xlabel "eps = 0.06 real"
plot real(eta(t)),imag(eta(t))
pause -1 "Press Return"
eps = 0.06*{1,-1}
set xlabel "eps = 0.06 + i0.06"
plot real(eta(t)),imag(eta(t))
pause -1 "Press Return"
reset
|
