program qjairflow parameter(ix=41) real x(ix),axi(ix),zu(ix),zo(ix) real f(ix),fp(ix) real e(ix),ep(ix) real cpu(ix),cpo(ix) real vu(ix),vo(ix) real cpusd(ix),cposd(ix) real ve(ix),xx(ix),vgrad(ix),theta(ix),cfr(ix) complex un,im,z1,dzds,dfdzu,dfdzo open(unit=15,file='qjairflow.dat',form='formatted') open(unit=16,file='qjairflow.xzu',form='formatted') open(unit=17,file='qjairflow.xzo',form='formatted') open(unit=18,file='qjairflow.cpu',form='formatted') open(unit=19,file='qjairflow.cpo',form='formatted') open(unit=20,file='qjairflow.cfu',form='formatted') open(unit=21,file='qjairflow.cfo',form='formatted') open(unit=22,file='qjairflow.cpusd',form='formatted') open(unit=23,file='qjairflow.cposd',form='formatted') c*****constants eps=1.e-7 pi=2.*asin(1.) degrad=pi/180. dtet=pi/(ix-1) un=(1.,0.) im=(0.,1.) c*****geometric/flow coefficients Uref=1. read(15,*)dm read(15,*)em read(15,*)alphd delt=2.*dm epsi=4./(3.*sqrt(3.))*em+1.0e-3 alpha=alphd*degrad sina=sin(alpha) cosa=cos(alpha) teta=-dtet do 1 i=1,ix ic=ix+1-i teta=teta+dtet sint=sin(teta) cost=cos(teta) xi=.5*(1.+cost) x(ic)=xi zui=.5*(epsi*(1.-cost)*sint+delt*sint*sint) zoi=.5*(-epsi*(1.-cost)*sint+delt*sint*sint) zu(ic)=zui zo(ic)=zoi z1=(1.+epsi*(1.-cost)+delt*sint)*exp(im*teta) dfdzu=(exp(-im*alpha)+2.*(im*sina+un*(epsi*cosa-delt*sina)) & /(z1+1.))*(z1*z1/(z1+un*epsi-im*delt)**2) z1=(1.+epsi*(1.-cost)-delt*sint)*exp(-im*teta) dfdzo=(exp(-im*alpha)+2.*(im*sina+un*(epsi*cosa-delt*sina)) & /(z1+1.))*(z1*z1/(z1+un*epsi-im*delt)**2) c write(6,*)'dfdzu=',dfdzu,' dfdzo=',dfdzo cpui=-1.+abs(dfdzu)**2 cpoi=-1.+abs(dfdzo)**2 cpu(ic)=cpui cpo(ic)=cpoi dzds=un*(2.*sint)- & im*(2.*epsi*(cost-cost*cost+sint*sint)+4.*delt*sint*cost) anorm=abs(dzds) if(anorm.lt.eps)then dzds=(un-im*2.*delt)/sqrt(1.+4.*delt**2) else dzds=dzds/(anorm+eps) endif vu(ic)=real(dfdzu*dzds) dzds=un*(2.*sint)+ & im*(2.*epsi*(cost-cost*cost+sint*sint)-4.*delt*sint*cost) anorm=abs(dzds) if(anorm.lt.eps)then dzds=(un-im*2.*delt)/sqrt(1.+4.*delt**2) else dzds=dzds/(anorm+eps) endif vo(ic)=real(dfdzo*dzds) cpusd(ic)=8.*em*(1.-4.*xi/3.)/sqrt(3.) & +2.*alpha*sqrt((1.-xi)/(xi+eps))+16.*dm*sqrt(xi*(1.-xi)) cposd(ic)=8.*em*(1.-4.*xi/3.)/sqrt(3.) & -2.*alpha*sqrt((1.-xi)/(xi+eps))-16.*dm*sqrt(xi*(1.-xi)) 1 continue c*****output on listing write(6,*)'*************************' write(6,*)'dimensionless parameters:' write(6,*)' X=C*x' write(6,*)' Z=C*z' write(6,*)' P-Pinf=0.5*Rho*V**2*Cp' write(6,*)' L=0.5*Rho*V**2*C*Cl' write(6,*)' D=0.5*Rho*V**2*C*Cd' write(6,*)' CM,o=0.5*Rho*V**2*C**2*Cm,o' write(6,*) write(6,*)'*************' write(6,*)'profile data:' write(6,*)' relative camber=',dm write(6,*)' relative thickness=',em write(6,*)' angle of attack=',alpha,' (rd) =' & ,alphd,' (deg)' c*****initialization write(6,*) write(6,*)'do you want to print the mesh, geometry' & ,' and flow? y/n=1/0' read(5,*)ipr if(ipr.eq.1)write(6,*) if(ipr.eq.1)write(6,*)' i= ',' x(i)= ',' f(i)= ' & ,' fp(i)= ',' e(i)= ',' ep(i)= ',' vu(i)= ' & ,' vo(i)=' axii=.5*(1.-cos(.5*dtet)) fim=-4.*dm*axii*(1.-axii) eim=-.5*em*(1.+cos(.5*dtet))*sin(.5*dtet) axiim=-axii do 2 i=1,ix teti=(i-1)*dtet xi=.5*(1.-cos(teti)) axii=.5*(1.-cos(teti+.5*dtet)) x(i)=xi axi(i)=axii f(i)=4.*dm*axii*(1.-axii) fp(i)=(f(i)-fim)/(axii-axiim) e(i)=.5*em*(1.+cos(teti+.5*dtet))*sin(teti+.5*dtet) ep(i)=(e(i)-eim)/(axii-axiim) eim=e(i) e(i)=.5*em*(1.+cos(teti))*sin(teti) fim=f(i) f(i)=4.*dm*xi*(1.-xi) axiim=axii if(i.eq.ix-1)then fim=-fim eim=-3.*eim axiim=2.-axii endif if(ipr.eq.1)write(6,1000)i,x(i),f(i),fp(i),e(i),ep(i),vu(i) & ,vo(i) write(16,*)x(i),zu(i) write(17,*)x(i),zo(i) write(18,*)x(i),cpu(i) write(19,*)x(i),cpo(i) write(22,*)x(i),cpusd(i) write(23,*)x(i),cposd(i) 2 continue write(6,*) write(6,*)'****************************************' write(6,*)'exact inviscid aerodynamic coefficients:' Cl=2.*pi*((1.+epsi)*sina+delt*cosa) Cd=0. write(6,*)'Cl=',Cl,' Cd=',Cd write(6,*)' ' Cmo=-.5*pi*((1.+epsi)*cosa*sina+2.*delt) write(6,*)'Cm,o=',Cmo write(6,*) write(6,*)'*******************************************' write(6,*)'small disturbance aerodynamic coefficients:' Clsd=2.*pi*(alpha+delt) Cdsd=0. write(6,*)'Clsd=',Clsd,' Cdsd=',Cdsd write(6,*) Cmosd=-.5*pi*(alpha+2.*delt) write(6,*)'Cm,osd=',Cmosd c*****boundary layer correction write(6,*) write(6,*)' ************************' write(6,*)'**********BOUDARY LAYER CORRECTION************' write(6,*)' ************************' 100 continue write(6,*) write(6,*)'Reynolds Number Re=? (Re<1.E-7 or Re>1.E7 to exit)' read(5,*)Re if(Re.le.eps.or.Re.ge.1./eps)goto 200 write(6,*) write(6,*)'do you want to print the viscous flow? y/n=1/0' read(5,*)iprv if(iprv.eq.1)write(6,*) if(iprv.eq.1)write(6,*)'upper surface' do 3 i=1,ix ve(i)=vu(i) f(i)=f(i)+e(i) 3 continue call moran(iprv,ix,itrans,isep,Re,x,f,ve,xx,vgrad,theta,cfr) Cd=0. do 4 i=2,ix Cd=Cd+.5*(cfr(i-1)+cfr(i))*(x(i)-x(i-1)) if(i.gt.isep)then cfr(i)=0. endif write(20,*)x(i),cfr(i) 4 continue write(6,*)'upper surface itrans=',itrans,' isep=' & ,isep,' Cd=',Cd Cdt=Cd if(iprv.eq.1)write(6,*) if(iprv.eq.1)write(6,*) if(iprv.eq.1)write(6,*)'lower surface' do 5 i=1,ix ve(i)=vo(i) f(i)=f(i)-2.*e(i) 5 continue call moran(iprv,ix,itrans,isep,Re,x,f,ve,xx,vgrad,theta,cfr) Cd=0. do 6 i=2,ix Cd=Cd+.5*(cfr(i-1)+cfr(i))*(x(i)-x(i-1)) if(i.gt.isep)then cfr(i)=0. endif write(21,*)x(i),cfr(i) 6 continue write(6,*)'lower surface itrans=',itrans,' isep=' & ,isep,' Cd=',Cd write(6,*) Cdt=Cdt+Cd write(6,*)'Cd total=',Cdt goto 100 200 continue c*****files write(6,*) write(6,*)'******input files:' write(6,*)'qjairflow.dat :dm,em,alpha' write(6,*)'******output files:' write(6,*)'qjairflow.xzu :x, zu(x)' write(6,*)'qjairflow.xzo :x, zo(x)' write(6,*)'qjairflow.cpu :x, Cpu(x)' write(6,*)'qjairflow.cpo :x, Cpo(x)' write(6,*)'qjairflow.cfu :x, cfu(x)' write(6,*)'qjairflow.cfo :x, cfo(x)' write(6,*)'qjairflow.cpusd :x, Cpusd(x)' write(6,*)'qjairflow.cposd :x, Cposd(x)' 1000 format(1x,i3,2x,f8.4,2x,f8.4,2x,f8.4,2x,f8.4,2x,f8.4,2x,f8.4, & 2x,f8.4) end subroutine moran(iprv,ix,itrans,isep,Re,x,y,ve,xx,vgrad,theta & ,cfr) real x(ix),y(ix),ve(ix) real xx(ix),theta(ix),vgrad(ix),cfr(ix) real yy(2) itrans=ix isep=ix pi=3.1415926535 eps=1.e-12 xx(1)=0. do 100 i=2,ix dx=x(i)-x(i-1) dy=y(i)-y(i-1) xx(i)=xx(i-1)+sqrt(dx*dx+dy*dy) 100 continue v1=ve(3) x1=xx(3) v2=ve(1) x2=xx(1) do 110 i=1,ix v3=v1 x3=x1 v1=v2 x1=x2 if(i.lt.ix)then v2=ve(i+1) x2=xx(i+1) else v2=ve(ix-1) x2=xx(ix-1) v3=ve(ix-2) x3=xx(ix-2) endif fact=(x3-x1)/(x2-x1) vgrad(i)=((v2-v1)*fact-(v3-v1)/fact)/(x3-x2) 110 continue if(iprv.eq.1)write(6,1000) i=1 120 continue v0=ve(i) if(ve(i).lt.0.)then theta(i)=0. h=0. cf=0. cfr(i)=cf cfe=0. if(iprv.eq.1)write(6,1010)x(i),y(i),ve(i),vgrad(i),theta(i) & ,h,cf,cfe i=i+1 goto 120 endif dv=eps*vgrad(i) if(dv.lt.eps**2)dv=eps**2 theta(i)=sqrt(eps*(6.*v0**5+15.*v0**4*dv+20.*v0**3*dv**2 &+15.*v0**2*dv**3+6.*v0*dv**4+dv**5)/(6.*(v0+dv)**6)) 200 continue alam=theta(i)**2*vgrad(i)*Re if(alam.gt..1)alam=.1 if(alam.lt.-.0841)alam=-.0841 if(alam.lt.-.0842)goto 400 call thwats(alam,h,el) if(ve(i).lt.eps)then cf=2.*el/Re/theta(i)/eps else cf=2.*el/Re/theta(i)/ve(i) endif cfr(i)=cf cfe=.664/sqrt(Re*(x(i)+eps-x(1))/(x(ix)-x(1))) if(iprv.eq.1)write(6,1010)x(i),y(i),ve(i),vgrad(i),theta(i),h & ,cf,cfe i=i+1 if(i.gt.ix)return if(ve(i).lt.0.)goto 120 dth2ve6=.225*(ve(i)**5+ve(i-1)**5)*(xx(i)-xx(i-1))/Re theta(i)=sqrt(((theta(i-1)**2)*(ve(i-1)**6)+dth2ve6) & /(ve(i)**6)) Rex=Re*xx(i)*ve(i) Ret=Re*theta(i)*ve(i) Retmax=1.174*(1.+22400./Rex)*Rex**.46 if(Ret.lt.Retmax.or.i.le.1)goto 200 itrans=i h=1.47 i=itrans yy(2)=h1ofh(h) yy(1)=theta(i-1) 320 continue dx=xx(i)-xx(i-1) call runge2(ix,i,Re,dx,yy,ve,vgrad) theta(i)=yy(1) h=hofh1(yy(2)) rtheta=Re*ve(i)*theta(i) cf=cfturb(rtheta,h) cfr(i)=cf cfe=.0592/(Re*(x(i)+eps-x(1))/(x(ix)-x(1)))**.2 if(iprv.eq.1)write(6,1020)x(i),y(i),ve(i),vgrad(i),theta(i),h & ,cf,cfe if(h.gt.2.4)goto 410 i=i+1 if(i.le.ix)goto 320 return 400 continue isep=i if(iprv.eq.1)write(6,*)'laminar separation' return 410 continue isep=i if(iprv.eq.1)write(6,*)'turbulent separation' return 1000 format(//,9x,'x',8x,'y',7x,'ve',6x,'vdot',5x,'theta',8x,'h', & 8x,'cf',3x,'cf-flat plate',/) 1010 format('l',f10.5,f9.5,2e10.4,e10.4,e10.4,2e10.4) 1020 format('t',f10.5,f9.5,2e10.4,e10.4,e10.4,2e10.4) end subroutine thwats(alam,h,el) if(alam.lt.0.)then el=.22+1.402*alam+.018*alam/(.107+alam) h=2.088+.0731/(.14+alam) else el=.22+alam*(1.57-1.8*alam) h=2.61-alam*(3.75-5.24*alam) endif return end function h1ofh(h) if(h.gt.1.6)then h1ofh=3.3+1.5501*(h-.6778)**(-3.064) else h1ofh=3.3+.8234*(h-1.1)**(-1.287) endif return end function hofh1(h1) if(h1.lt.5.3)then if(h1.lt.3.3)then hofh1=3.0 else hofh1=.6778+1.1536*(h1-3.3)**(-.326) endif else hofh1=1.1+.86*(h1-3.3)**(-.777) endif return end function cfturb(rtheta,h) cfturb=.246*(10.**(-.678*h))*(rtheta)**(-.268) return end subroutine derivs(ix,i,Re,yp,yt,ve,vgrad) real ve(ix),vgrad(ix) real yp(2),yt(2) h1=h1ofh(hofh1(yt(2))) h=hofh1(h1) rtheta=Re*ve(i)*yt(1) yp(1)=-(h+2.)*yt(1)*vgrad(i)/ve(i) & +.5*cfturb(rtheta,h) yp(2)=-h1*(vgrad(i)/ve(i)+yp(1)/yt(1)) & +.0306*(h1-3.)**(-.6169)/yt(1) return end subroutine runge2(ix,i,Re,dx,yy,ve,vgrad) real ve(ix),vgrad(ix) real ys(2),yy(2) real yp(2),yt(2) do 1 j=1,2 yt(j)=yy(j) 1 continue call derivs(ix,i-1,Re,yp,yt,ve,vgrad) do 2 j=1,2 yt(j)=yy(j)+dx*yp(j) ys(j)=yy(j)+.5*dx*yp(j) 2 continue call derivs(ix,i,Re,yp,yt,ve,vgrad) do 3 j=1,2 yy(j)=ys(j)+.5*dx*yp(j) 3 continue return end