1 /*
2 * $Id: demo4.i,v 1.1 2005-09-18 22:05:55 dhmunro Exp $
3 * Airfoil demo
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 func demo4(mono=)
12 /* DOCUMENT demo4
13 or demo4, mono=1
14 solves for the flow past a 2D airfoil using Kutta-Jakowski theory.
15 The colors represent static pressure (darker is lower pressure),
16 red lines are streamlines.
17 Solutions for various angles of attack are shown by animation.
18 With the mono=1 keyword, only the streamlines are shown. (On a
19 monochrome terminal, the pressure makes the streamlines invisible.)
20 */
21 {
22 require, "movie.i";
23 local nx, ny;
24 imax= 50;
25 movie, display;
26 display, imax;
27 }
28
display(i)29 func display(i)
30 {
31 attack= 28.*double(imax-i)/(imax-1);
32 solve, attack, 2.0, 0.2;
33 if (i==1) {
34 extern cn, cx;
35 cn= 0.8*min(pressure);
36 cx= max(pressure);
37 }
38 limits, -2.8, 2.4, -2.4, 2.8;
39 plmesh, z.im,z.re;
40 if (!mono) plf, zncen(pressure), cmin=cn,cmax=cx;
41 plc, potential.im, marks=0,color="red", levs=span(-2.5,3.5,33);
42 plg,z.im(,1),z.re(,1), marks=0;
43 return i<imax;
44 }
45
solve(attack,chord,thick)46 func solve(attack, chord, thick)
47 {
48 extern w, z, potential, velocity, pressure; /* outputs */
49 attack*= pi/180.;
50 a= 0.5*chord;
51 w= get_mesh(a, attack);
52 /* the log(w) term adds just enough circulation to move the rear
53 stagnation point to the trailing edge -- the Kutta condition */
54 potential= jakowski(w, a) + 2i*a*sin(attack)*log(w);
55 dpdw= 1.-(a/w)^2 + 2i*a*sin(attack)/w;
56 emith= exp(-1i*attack);
57 a= a*emith;
58 b= thick*emith; /* could add camber here too someday? */
59 w-= b;
60 a-= b;
61 z= jakowski(w, a);
62 dzdw= 1.-(a/w)^2;
63 velocity= conj(dpdw/dzdw);
64 pressure= 0.5*(1.0-abs(velocity)^2);
65 }
66
get_mesh(a,attack)67 func get_mesh(a, attack)
68 {
69 /* get a mesh in the w-plane
70 -- the plane in which the airfoil is a circle
71 the mesh splits at the point which will become the trailing edge */
72 if (is_void(nx)) nx= 120;
73 if (is_void(ny)) ny= 30;
74 a= abs(a);
75 theta= span(1.e-9, 2*pi-1.e-9, nx) - attack;
76 r= a/span(1.0+1.e-9, 0.25, ny)(-,);
77 return r*exp(1i*theta);
78 }
79
jakowski(z,a)80 func jakowski(z, a)
81 {
82 /* Jakowski transform - circle of radius a into slot of length 4a */
83 return z + a*a/z;
84 }
85
ijakowski(z,a)86 func ijakowski(z, a)
87 {
88 /* Inverse Jakowski - slot back into circle (unused here) */
89 z*= 0.5;
90 a= complex(a);
91 sgn= sign((z*conj(a)).re);
92 return z + sgn(z);
93 }
94