c ************************************************************** c * * c * this program computes velocity and pressure coefficient * c * around an airfoil, whose contour is approximated by * c * vortex panels of linearly varying strength. * c * * c ************************************************************** c c parameter (npts = 13) dimension xb(npts),yb(npts),x(npts),y(npts-1),s(npts-1), . sine(npts-1),cosine(npts-1),theta(npts-1),v(npts-1), . cp(npts-1),gama(npts),rhs(npts),cn1(npts-1,npts-1), . cn2(npts-1,npts-1),ct1(npts-1,npts-1), . ct2(npts-1,npts-1),an(npts,npts),at(npts-1,npts), . indx(npts) c c Specify coordinates (xb,yb) of boundary points on airfoil c surface. The last point coincides with the first. c c c Nominal data c data xb 1 /1.0000,.933, .750, .500, .250, .067, 0. , .067, .250, .500, 2 .750, .933, 1.00/ data yb 1 /.0000,-.005,-.017,-.033,-.042,-.033,0.,.045,.076, 2 .072,.044,.013, 0./ c c 18 panel Wortmann airfoil data c c data xb/ c . 0.1000E+01, 0.8889E+00, 0.7778E+00, 0.6667E+00, 0.5556E+00, c . 0.4444E+00, 0.3333E+00, 0.2222E+00, 0.1111E+00, 0.0000E+00, c . 0.1111E+00, 0.2222E+00, 0.3333E+00, 0.4444E+00, 0.5556E+00, c . 0.6667E+00, 0.7778E+00, 0.8889E+00, 0.1000E+01/ c data yb/ c . 0.0000E+00, 0.2364E-01, 0.2745E-01, 0.2148E-01, 0.9967E-02, c . -0.3737E-02, -0.1508E-01, -0.2161E-01, -0.2209E-01, 0.0000E+00, c . 0.8570E-01, 0.1112E+00, 0.1211E+00, 0.1199E+00, 0.1095E+00, c . 0.9087E-01, 0.6605E-01, 0.3765E-01, 0.0000E+00/ c m = npts - 1 mp1 = npts pi = 4.*atan(1.) alpha = 8.*pi/180. c c Coordinates (x,y) of control point and panel length s are c computed for each of the vortex panels. rhs represents c the right-hand side of eq. (47). c do 1 i = 1,m ip1 = i+1 x(i) = .5*(xb(i)+xb(ip1)) y(i) = .5*(yb(i)+yb(ip1)) s(i) = sqrt((xb(ip1)-xb(i))**2+(yb(ip1)-yb(i))**2) theta(i) = atan2(yb(ip1)-yb(i),xb(ip1)-xb(i)) sine(i) = sin(theta(i)) cosine(i) = cos(theta(i)) 1 rhs(i) = sin(theta(i)-alpha) do 3 i = 1,m do 3 j = 1,m if (i.eq.j) go to 2 a = -(x(i)-xb(j))*cosine(j)-(y(i)-yb(j))*sine(j) b = (x(i)-xb(j))**2+(y(i)-yb(j))**2 c = sin(theta(i)-theta(j)) d = cos(theta(i)-theta(j)) e = (x(i)-xb(j))*sine(j)-(y(i)-yb(j))*cosine(j) f = alog(1.+s(j)*(s(j)+2.*a)/b) g = atan2(e*s(j),b+a*s(j)) p = (x(i)-xb(j))*sin(theta(i)-2.*theta(j)) . +(y(i)-yb(j))*cos(theta(i)-2.*theta(j)) q = (x(i)-xb(j))*cos(theta(i)-2.*theta(j)) . -(y(i)-yb(j))*sin(theta(i)-2.*theta(j)) cn2(i,j) = d+.5*q*f/s(j)-(a*c+d*e)*g/s(j) cn1(i,j) = .5*d*f+c*g-cn2(i,j) ct2(i,j) = c+.5*p*f/s(j)+(a*d-c*e)*g/s(j) ct1(i,j) = .5*c*f-d*g-ct2(i,j) go to 3 2 cn1(i,j) = -1.0 cn2(i,j) = 1.0 ct1(i,j) = .5*pi ct2(i,j) = .5*pi 3 continue c c compute influence coefficients in eqs. (47) and (49), c respectively. c do 4 i = 1,m an(i,1) = cn1(i,1) an(i,mp1) = cn2(i,m) at(i,1) = ct1(i,1) at(i,mp1) = ct2(i,m) do 4 j = 2,m an(i,j) = cn1(i,j)+cn2(i,j-1) 4 at(i,j) = ct1(i,j)+ct2(i,j-1) an(mp1,1) = 1.0 an(mp1,mp1) = 1.0 do 5 j = 2,m 5 an(mp1,j) = 0. rhs(mp1) = 0. c c Solve eq.(47) for dimensionless strengths gama using c LU decomposition. Then compute and print dimensionless c velocity and pressure coefficient at control points. c write(6,55)m,alpha*180./pi 55 format('1'/30x,i2,' panels, alpha =',f5.2,' deg'/) write(6,6) 6 format(11x,'i',4x,'x(i)',4x,'y(i)',4x,'theta(i)', 1 3x,'s(i)',3x,'gama(i)',3x,'v(i)',6x,'cp(i)'/10x, 2 '---',3x,'----',4x,'----',4x,'--------',3x,'----', 3 3x,'-------',3x,'----',6x,'-----') c c Solve the system of equations with an LU decomposition and c back substitution c call ludcmp(an,mp1,mp1,indx,d) call lubksb(an,mp1,mp1,indx,rhs) c do 105 i = 1 , mp1 105 gama(i) = rhs(i) c do 8 i = 1,m v(i) = cos(theta(i)-alpha) do 7 j = 1,mp1 v(i) = v(i)+at(i,j)*gama(j) 7 cp(i) = 1.-v(i)**2 8 write(6,9)i,x(i),y(i),theta(i),s(i),gama(i),v(i),cp(i) c 9 format(10x,i2,f8.4,f9.4,f10.4,f8.4,2f9.4,f10.4) write(6,10)mp1,gama(mp1) 10 format(10x,i2,35x,f9.4/) c c cfx=0.0 cfy=0.0 cm=0.0 c do 95 i=1,m dx=xb(i+1)-xb(i) dy=yb(i+1)-yb(i) cfx=cfx+cp(i)*dy cfy=cfy-cp(i)*dx cm=cm+cp(i)*(dx*x(i)+dy*y(i)) 95 continue c cd=cfx*cos(alpha)+cfy*sin(alpha) cl=cfy*cos(alpha)-cfx*sin(alpha) c write(6,110)cd,cl,cm 110 format(' cd =',f8.5,' cl =',f8.5,' cm =',f8.5) c c stop end c c c lu decomposition and backsubstitution c taken from "numerical recipes" c june 6, 1990 c c ...................................................................... c "ludcmp()"; given an "n" by "n" matrix "a", with physical dimension c "np", this routine replaces it by the lu decomposition of a rowwise c permutation of itself. "a" and "n" are input. "a" is output; "indx" c is an output vector which records the row permutation effected by the c partial pivoting; "d" is output as + or - 1 depending on whether the c number of row interchanges was even or odd, respectively. this c routine is used in combination with "lubksb()" to solve linear c equations or invert a matrix. c subroutine ludcmp(a, n, np, indx, d) parameter (nmax = 100, tiny = 1.0e-20) dimension a(np, np), indx(n), vv(nmax) d = 1.0 do 12 i = 1, n aamax = 0.0 do 11 j = 1, n if (abs(a(i,j)) .gt. aamax) aamax = abs(a(i,j)) 11 continue if (aamax .eq. 0.0) pause 'singular matrix' vv(i) = 1.0 / aamax 12 continue do 19 j = 1, n do 14 i = 1, j - 1 sum = a(i,j) do 13 k = 1, i - 1 sum = sum - a(i,k) * a(k,j) 13 continue a(i,j) = sum 14 continue aamax = 0.0 do 16 i = j, n sum = a(i,j) do 15 k = 1, j - 1 sum = sum - a(i,k) * a(k,j) 15 continue a(i,j) = sum dum = vv(i) * abs(sum) if (dum .ge. aamax) then imax = i aamax = dum endif 16 continue if (j .ne. imax) then do 17 k = 1, n dum = a(imax,k) a(imax,k) = a(j,k) a(j,k) = dum 17 continue d = -d vv(imax) = vv(j) endif indx(j) = imax if (a(j,j) .eq. 0.0) a(j,j) = tiny if (j .ne. n) then dum = 1.0 / a(j,j) do 18 i = j + 1, n a(i,j) = a(i,j) * dum 18 continue endif 19 continue return end c c ..................................................................... c "lubksb()" solves the set of "n" linear equations [a] * {x} = {b}. c here "a" is input, not as the matrix "a" but rather as its lu c decomposition, determined by the routine "ludcmp()". "indx" is c input as the permutation vector returned by "ludcmp()". "b" is c input as the right-hand side vector, and returns with the solution c vector {x}. "a", "n", "np" and "indx" are not modified by this c routine and can be left in place for successive calls with different c right-hand sides "b". this routine takes into account the possibility c that "b" will begin with many zero elements, so it is efficient for c use in matrix inversion. c subroutine lubksb(a, n, np, indx, b) dimension a(np, np), indx(n), b(n) ii = 0 do 12 i = 1, n ll = indx(i) sum = b(ll) b(ll) = b(i) if (ii .ne. 0) then do 11 j = ii, i - 1 sum = sum - a(i,j) * b(j) 11 continue else if (sum .ne. 0.0) then ii = i endif b(i) = sum 12 continue do 14 i = n, 1, -1 sum = b(i) if (i .lt. n) then do 13 j = i + 1, n sum = sum - a(i,j) * b(j) 13 continue endif b(i) = sum / a(i,i) 14 continue return end