program panel c c smith-hess panel method for single c element lifting airfoil in 2-d c incompressible flow c dimension z(100) pi=3.1415926585 c call indata call setup c 100 write(6,1000) 1000 format(///,' input alpha in degrees') read(5,*)alpha if (alpha.gt.90) go to 200 c cosalf=cos(alpha*pi/180.) sinalf=sin(alpha*pi/180.) call cofish(sinalf,cosalf) call gauss(1) call veldis(sinalf,cosalf) call fandm(sinalf,cosalf) go to 100 c 200 stop end c c********************************************************** c subroutine setup c common /bod/ nlower,nupper,nodtot,x(100),y(100), 2 costhe(100),sinthe(100) common /num/ pi,pi2inv c pi=3.1415926585 pi2inv=.5/pi c c set coordinates of nodes on body surface c write(6,1000) 1000 format(///' body shape'//4x,'x',9x,'y'/) npoints=nlower sign=-1.0 nstart=0 do 110 nsurf=1,2 do 100 n=1,npoints fract=float(n-1)/float(npoints) z=.5*(1.-cos(pi*fract)) i=nstart+n call body(z,sign,x(i),y(i)) write(6,1010)x(i),y(i) 1010 format(f8.4,f10.4) 100 continue npoints=nupper sign=1.0 nstart=nlower 110 continue nodtot=nupper+nlower x(nodtot+1)=x(1) y(nodtot+1)=y(1) c c set slopes of panels c do 200 i=1,nodtot dx=x(i+1)-x(i) dy=y(i+1)-y(i) dist=sqrt(dx*dx+dy*dy) sinthe(i)=dy/dist costhe(i)=dx/dist 200 continue c return end c c********************************************************** c subroutine body(z,sign,x,y) c c return coordinates of point on body surface c c z = node spacing parameter c x,y = cartesian coordinates c sign = +1. for upper surface, -1. for lower surface c common /par/ naca,tau,epsmax,ptmax c if (sign.lt.0.) z=1.-z call naca45(z,thick,camber,beta) x=z-sign*thick*sin(beta) y=camber+sign*thick*cos(beta) c return end c c********************************************************** c subroutine cofish(sinalf,cosalf) c c set coefficients of linear system c common /bod/ nlower,nupper,nodtot,x(100),y(100), 2 costhe(100),sinthe(100) common /cof/ a(101,111),kutta common /num/ pi,pi2inv c kutta=nodtot+1 c c initialize coefficients c do 90 j=1,kutta 90 a(kutta,j)=0.0 c c set vn=0. at midpoint of ith panel c do 120 i=1,nodtot xmid=.5*(x(i)+x(i+1)) ymid=.5*(y(i)+y(i+1)) a(i,kutta)=0.0 c c find contribution of jth panel c do 110 j=1,nodtot flog=0.0 ftan=pi if (j.eq.i) go to 100 dxj=xmid-x(j) dxjp=xmid-x(j+1) dyj=ymid-y(j) dyjp=ymid-y(j+1) flog=.5*alog((dxjp*dxjp+dyjp*dyjp)/(dxj*dxj+dyj*dyj)) ftan=atan2(dyjp*dxj-dxjp*dyj,dxjp*dxj+dyjp*dyj) 100 ctimtj=costhe(i)*costhe(j)+sinthe(i)*sinthe(j) stimtj=sinthe(i)*costhe(j)-sinthe(j)*costhe(i) a(i,j)=pi2inv*(ftan*ctimtj+flog*stimtj) b=pi2inv*(flog*ctimtj-ftan*stimtj) a(i,kutta)=a(i,kutta)+b if ((i.gt.1).and.(i.lt.nodtot)) go to 110 c c if ith panel touches trailing edge, add contribution c to kutta condition c a(kutta,j)=a(kutta,j)-b a(kutta,kutta)=a(kutta,kutta)+a(i,j) 110 continue c c fill in known sides c a(i,kutta+1)=sinthe(i)*cosalf-costhe(i)*sinalf 120 continue a(kutta,kutta+1)=-(costhe(1)+costhe(nodtot))*cosalf 2 -(sinthe(1)+sinthe(nodtot))*sinalf c return end c c********************************************************** c subroutine veldis(sinalf,cosalf) c c compute and print out pressure distribution c common /bod/ nlower,nupper,nodtot,x(100),y(100), 2 costhe(100),sinthe(100) common /cof/ a(101,111),kutta common /cpd/ cp(100) common /num/ pi,pi2inv dimension q(150) c write(6,1000) 1000 format(///' pressure distribution'//4x,'x',8x,'cp'/) c c retrieve solution from a-matrix c do 50 i=1,nodtot 50 q(i)=a(i,kutta+1) gamma=a(kutta,kutta+1) c c find vt and cp at midpoint of ith panel c do 130 i=1,nodtot xmid=.5*(x(i)+x(i+1)) ymid=.5*(y(i)+y(i+1)) vtang=cosalf*costhe(i)+sinalf*sinthe(i) c c add contributions of jth panel c do 120 j=1,nodtot flog=0.0 ftan=pi if (j.eq.i) go to 100 dxj=xmid-x(j) dxjp=xmid-x(j+1) dyj=ymid-y(j) dyjp=ymid-y(j+1) flog=.5*alog((dxjp*dxjp+dyjp*dyjp)/(dxj*dxj+dyj*dyj)) ftan=atan2(dyjp*dxj-dxjp*dyj,dxjp*dxj+dyjp*dyj) 100 ctimtj=costhe(i)*costhe(j)+sinthe(i)*sinthe(j) stimtj=sinthe(i)*costhe(j)-sinthe(j)*costhe(i) aa=pi2inv*(ftan*ctimtj+flog*stimtj) b=pi2inv*(flog*ctimtj-ftan*stimtj) vtang=vtang-b*q(j)+gamma*aa 120 continue cp(i)=1.-vtang*vtang write(6,1010)xmid,cp(i) 1010 format(f8.4,f10.4) 130 continue c return end c c********************************************************** c subroutine fandm(sinalf,cosalf) c c compute and print out cd,cl,cmle c common /bod/ nlower,nupper,nodtot,x(100),y(100), 2 costhe(100),sinthe(100) common /cpd/ cp(100) c cfx=0.0 cfy=0.0 cm=0.0 c do 100 i=1,nodtot xmid=.5*(x(i)+x(i+1)) ymid=.5*(y(i)+y(i+1)) dx=x(i+1)-x(i) dy=y(i+1)-y(i) cfx=cfx+cp(i)*dy cfy=cfy-cp(i)*dx cm=cm+cp(i)*(dx*xmid+dy*ymid) 100 continue cd=cfx*cosalf+cfy*sinalf cl=cfy*cosalf-cfx*sinalf write(6,1000)cd,cl,cm 1000 format(////' cd =',f8.5,' cl =',f8.5,' cm =',f8.5) c return end c c********************************************************** c subroutine gauss(nrhs) c c solution of linear algebraic system by c gaussian elimination with partial pivoting c c [a] = coefficient matrix c neqns = number of equations c nrhs = number of right-hand sides c c right-hand sides and solutions stored in c columns neqns+1 thru neqns+nrhs of a c common /cof/ a(101,111),neqns c np=neqns+1 ntot=neqns+nrhs c c gauss reduction c do 150 i=2,neqns c c search for largest entry in (i-1)th column c on or below main diagonal c im=i-1 imax=im amax=abs(a(im,im)) do 110 j=i,neqns if (amax.ge.abs(a(j,im))) go to 110 imax=j amax=abs(a(j,im)) 110 continue c c switch (i-1)th and imaxth equations c if (imax.ne.im) go to 140 do 130 j=im,ntot temp=a(im,j) a(im,j)=a(imax,j) a(imax,j)=temp 130 continue c c eliminate (i-1)th unknown from c ith thru neqnsth equations c 140 do 150 j=i,neqns r=a(j,im)/a(im,im) do 150 k=i,ntot 150 a(j,k)=a(j,k)-r*a(im,k) c c back substitution c do 220 k=np,ntot a(neqns,k)=a(neqns,k)/a(neqns,neqns) do 210 l=2,neqns i=neqns+1-l ip=i+1 do 200 j=ip,neqns 200 a(i,k)=a(i,k)-a(i,j)*a(j,k) 210 a(i,k)=a(i,k)/a(i,i) 220 continue c return end c c********************************************************** c subroutine indata c c set parameters of body shape, flow c situation, and node distribution c c user must input: c nlower = number of nodes on lower surface c nupper = number of nodes on upper surface c plus data on body c common /bod/ nlower,nupper,nodtot,x(100),y(100), 2 costhe(100),sinthe(100) common /par/ naca,tau,epsmax,ptmax c write(6,*)'input nlower,nupper' read(5,*)nlower,nupper write(6,*)'input naca number' read(5,*)naca c ieps=naca/1000 iptmax=naca/100-10*ieps itau=naca-1000*ieps-100*iptmax epsmax=ieps*0.01 ptmax=iptmax*0.1 tau=itau*0.01 if (ieps.lt.10) return ptmax=0.2025 epsmax=2.6595*ptmax**3 c return end c c********************************************************** c subroutine naca45(z,thick,camber,beta) c common /par/ naca,tau,epsmax,ptmax c c evaluate thickness and camber c for naca 4- or 5-digit airfoil c thick=0.0 if (z.lt.1.e-10) go to 100 thick=5.*tau*(.2969*sqrt(z)-z*(.126+z*(.3537 2 -z*(.2843-z*.1015)))) 100 if (epsmax.eq.0.) go to 130 if (naca.gt.9999) go to 140 if (z.gt.ptmax) go to 110 camber=epsmax/ptmax/ptmax*(2.*ptmax-z)*z dcamdx=2.*epsmax/ptmax/ptmax*(ptmax-z) go to 120 110 camber=epsmax/(1.-ptmax)**2*(1.+z-2.*ptmax)*(1.-z) dcamdx=2.*epsmax/(1.-ptmax)**2*(ptmax-z) 120 beta=atan(dcamdx) c return c 130 camber=0.0 beta=0.0 c return c 140 if (z.gt.ptmax) go to 150 w=z/ptmax camber=epsmax*w*((w-3.)*w+3.-ptmax) dcamdx=epsmax*3.*w*(1.-w)/ptmax go to 120 150 camber=epsmax*(1.-z) dcamdx=-epsmax go to 120 c end