git.maemo.org Git - gnuplot/blob - demo/airfoil.dem

   1 #
   2 # $Id: airfoil.dem,v 1.10 2006/06/30 02:17:21 sfeam Exp $
   3 #
   4 # This demo shows how to use bezier splines to define NACA four
   5 # series airfoils and complex variables to define Joukowski
   6 # Airfoils.  It will be expanded after overplotting in implemented
   7 # to plot Coefficient of Pressure as well.
   8 #               Alex Woo, Dec. 1992
   9 #
  10 # The definitions below follows: "Bezier presentation of airfoils",
  11 # by Wolfgang Boehm, Computer Aided Geometric Design 4 (1987) pp 17-22.
  12 #
  13 #                               Gershon Elber, Nov. 1992
  14 #
  15 # m = percent camber
  16 # p = percent chord with maximum camber
  17 print  "NACA four series airfoils by bezier splines"
  18 print  "Will add pressure distribution later with Overplotting"
  19 mm = 0.6
  20 # NACA6xxx
  21 thick = 0.09
  22 # nine percent  NACAxx09
  23 pp = 0.4
  24 # NACAx4xx
  25 # Combined this implies NACA6409 airfoil
  26 #
  27 # Airfoil thickness function.
  28 #
  29 set xlabel "NACA6409 -- 9% thick, 40% max camber, 6% camber"
  30 x0 = 0.0
  31 y0 = 0.0
  32 x1 = 0.0
  33 y1 = 0.18556
  34 x2 = 0.03571
  35 y2 = 0.34863
  36 x3 = 0.10714
  37 y3 = 0.48919
  38 x4 = 0.21429
  39 y4 = 0.58214
  40 x5 = 0.35714
  41 y5 = 0.55724
  42 x6 = 0.53571
  43 y6 = 0.44992
  44 x7 = 0.75000
  45 y7 = 0.30281
  46 x8 = 1.00000
  47 y8 = 0.01050
  48 #
  49 # Directly defining the order 8 Bezier basis function for a faster evaluation.
  50 #
  51 bez_d4_i0(x) =     (1 - x)**4
  52 bez_d4_i1(x) = 4 * (1 - x)**3 * x
  53 bez_d4_i2(x) = 6 * (1 - x)**2 * x**2
  54 bez_d4_i3(x) = 4 * (1 - x)**1 * x**3
  55 bez_d4_i4(x) =                  x**4
  56
  57 bez_d8_i0(x) =      (1 - x)**8
  58 bez_d8_i1(x) =  8 * (1 - x)**7 * x
  59 bez_d8_i2(x) = 28 * (1 - x)**6 * x**2
  60 bez_d8_i3(x) = 56 * (1 - x)**5 * x**3
  61 bez_d8_i4(x) = 70 * (1 - x)**4 * x**4
  62 bez_d8_i5(x) = 56 * (1 - x)**3 * x**5
  63 bez_d8_i6(x) = 28 * (1 - x)**2 * x**6
  64 bez_d8_i7(x) =  8 * (1 - x)    * x**7
  65 bez_d8_i8(x) =                   x**8
  66
  67
  68 m0 = 0.0
  69 m1 = 0.1
  70 m2 = 0.1
  71 m3 = 0.1
  72 m4 = 0.0
  73 mean_y(t) = m0 * mm * bez_d4_i0(t) + \
  74             m1 * mm * bez_d4_i1(t) + \
  75             m2 * mm * bez_d4_i2(t) + \
  76             m3 * mm * bez_d4_i3(t) + \
  77             m4 * mm * bez_d4_i4(t)
  78
  79 p0 = 0.0
  80 p1 = pp / 2
  81 p2 = pp
  82 p3 = (pp + 1) / 2
  83 p4 = 1.0
  84 mean_x(t) = p0 * bez_d4_i0(t) + \
  85             p1 * bez_d4_i1(t) + \
  86             p2 * bez_d4_i2(t) + \
  87             p3 * bez_d4_i3(t) + \
  88             p4 * bez_d4_i4(t)
  89
  90 z_x(x) = x0 * bez_d8_i0(x) + x1 * bez_d8_i1(x) + x2 * bez_d8_i2(x) + \
  91          x3 * bez_d8_i3(x) + x4 * bez_d8_i4(x) + x5 * bez_d8_i5(x) + \
  92          x6 * bez_d8_i6(x) + x7 * bez_d8_i7(x) + x8 * bez_d8_i8(x)
  93
  94 z_y(x, tk) = \
  95    y0 * tk * bez_d8_i0(x) + y1 * tk * bez_d8_i1(x) + y2 * tk * bez_d8_i2(x) + \
  96    y3 * tk * bez_d8_i3(x) + y4 * tk * bez_d8_i4(x) + y5 * tk * bez_d8_i5(x) + \
  97    y6 * tk * bez_d8_i6(x) + y7 * tk * bez_d8_i7(x) + y8 * tk * bez_d8_i8(x)
  98
  99 #
 100 # Given t value between zero and one, the airfoild curve is defined as
 101 #
 102 #                       c(t) = mean(t1(t)) +/- z(t2(t)) n(t1(t)),
 103 #
 104 # where n is the unit normal to the mean line. See the above paper for more.
 105 #
 106 # Unfortunately, the parametrization of c(t) is not the same for mean(t1)
 107 # and z(t2). The mean line (and its normal) can assume linear function t1 = t,
 108 #                                                     -1
 109 # but the thickness z_y is, in fact, a function of z_x  (t). Since it is
 110 # not obvious how to compute this inverse function analytically, we instead
 111 # replace t in c(t) equation above by z_x(t) to get:
 112 #
 113 #                       c(z_x(t)) = mean(z_x(t)) +/- z(t) n(z_x(t)),
 114 #
 115 # and compute and display this instead. Note we also ignore n(t) and assumes
 116 # n(t) is constant in the y direction,
 117 #
 118
 119 airfoil_y1(t, thick) = mean_y(z_x(t)) + z_y(t, thick)
 120 airfoil_y2(t, thick) = mean_y(z_x(t)) - z_y(t, thick)
 121 airfoil_y(t) = mean_y(z_x(t))
 122 airfoil_x(t) = mean_x(z_x(t))
 123 unset grid
 124 unset zeroaxis
 125 set parametric
 126 set xrange [-0.1:1.1]
 127 set yrange [-0.1:.7]
 128 set trange [ 0.0:1.0]
 129 set title "NACA6409 Airfoil"
 130 plot airfoil_x(t), airfoil_y(t) title "mean line" w l lt 2, \
 131      airfoil_x(t), airfoil_y1(t, thick) title "upper surface" w l lt 1, \
 132      airfoil_x(t), airfoil_y2(t, thick) title "lower surface" w l lt 1
 133 pause -1 "Press Return"
 134 mm = 0.0
 135 pp = .5
 136 thick = .12
 137 set title "NACA0012 Airfoil"
 138 set xlabel "12% thick, no camber -- classical test case"
 139 plot airfoil_x(t), airfoil_y(t) title "mean line" w l lt 2, \
 140      airfoil_x(t), airfoil_y1(t, thick) title "upper surface" w l lt 1, \
 141      airfoil_x(t), airfoil_y2(t, thick) title "lower surface" w l lt 1
 142 pause -1 "Press Return"
 143 set title ""
 144 set xlab ""
 145 set key box
 146 set parametric
 147 set samples 100
 148 set isosamples 10
 149 set style data lines
 150 set style function lines
 151 print  "Joukowski Airfoil using Complex Variables"
 152 set title "Joukowski Airfoil using Complex Variables" offset 0,0
 153 set time
 154 set yrange [-.2 : 1.8]
 155 set trange [0: 2*pi]
 156 set xrange [-.6:.6]
 157 zeta(t) = -eps + (a+eps)*exp(t*{0,1})
 158 eta(t) = zeta(t) + a*a/zeta(t)
 159 eps = 0.06
 160 a =.250
 161 set xlabel "eps = 0.06 real"
 162 plot real(eta(t)),imag(eta(t))
 163 pause -1 "Press Return"
 164 eps = 0.06*{1,-1}
 165 set xlabel "eps = 0.06 + i0.06"
 166 plot real(eta(t)),imag(eta(t))
 167 pause -1 "Press Return"
 168 reset