c...2-D Transonic Small Disturbance Equation Solver. c Flow over a circular arc airfoil parameter(idim=160,jdim=130,ndim=5000) dimension phi(idim,jdim),a(idim,jdim),mu(idim,jdim), & x(idim,jdim),y(idim,jdim),dx(idim,jdim),dy(idim,jdim), & amat(jdim),bmat(jdim),cmat(jdim),dmat(jdim), & dumm1(idim),dumm2(idim),rmax(ndim) real mu write(*,*) '2-D Transonic Small Disturbance Equation Solver' write(*,*) '(c) 2002, Stephen M. Ruffin' write(*,*) ' ' c...Input Parameters.................................................. write(*,*) 'Enter Mach Number' read(*,*) fsmach write(*,*) 'Enter circular-arc-airfoil thickness ratio' read(*,*) th nmax = 1000 imax = 74 jmax = 14 icst = 15 dxmin = 0.025 dymin = 0.03 rmin = 1.0e-3 gamma = 1.4 qinf = 1.0 rinf = 1.0 c..................................................................... open(unit=3,file='resid.dat',form='formatted') write(6,*) ' ',' N ',' RESIDmax ',' i at max ',' j at max ' write(3,*) ' ',' N ',' RESIDmax ',' i at max ',' j at max ' ainf = 1.0/fsmach pinf = 1.0/(gamma*fsmach**2) c...Set Initial Conditions do 10 i = 1,imax do 10 j = 1,jmax phi(i,j) = 1.0 10 continue c...Generate Grid call grid(dxmin,dymin,icst,icend,x,y,dx,dy,imax,jmax,idim,jdim) do 30 n = 1,nmax c...Compute A and MU call aandmu(phi,a,mu,dx,fsmach,gamma,qinf, & imax,jmax,idim,jdim) do 20 i = 3,imax-1 c...Fill interior of diagonal matricies call allmats(i,amat,bmat,cmat,dmat,phi,a,mu,dx,dy, & fsmach,gamma,qinf,imax,jmax,idim,jdim) c...Enforce B.C.'s through diagonal matricies call bcmat(i,amat,bmat,cmat,dmat,x,dy,th,qinf,icst,icend, & imax,jmax,idim,jdim) c...Solve tridiagonal equation set call scaltri(i,phi,amat,bmat,cmat,dmat,imax,jmax,idim,jdim) 20 continue c...Compute residual call resid(n,rmax,residm0,phi,a,mu,dx,dy,imax,jmax,idim,jdim,ndim) if(rmax(n).le.rmin) go to 31 30 continue 31 continue c...Write output call output(phi,x,y,dx,dy,fsmach,gamma,qinf,rinf,pinf, & dumm1,dumm2,n,rmax,imax,jmax,nmax,idim,jdim,ndim) close(3) stop end subroutine grid(dxmin,dymin,icst,icend, & x,y,dx,dy,imax,jmax,idim,jdim) dimension x(idim,jdim),y(idim,jdim), & dx(idim,jdim),dy(idim,jdim) y(1,1) = 0.0 y(1,2) = dymin do 10 j = 3,jmax y(1,j) = y(1,j-1) + 2.0*(j-2)*dymin 10 continue do 20 i=2,imax do 20 j=1,jmax y(i,j) = y(1,j) 20 continue icend = icst + 1.0/dxmin x(icst ,1) = 0.0 x(icst-1,1) = -dxmin x(icend ,1) = 1.0 x(icend+1,1) = 1.0 + dxmin do 30 i = icst-2,1,-1 x(i,1) = x(i+1,1) - 2.0*(icst-1-i)*dxmin 30 continue do 31 i = icst+1,icend-1 x(i,1) = x(i-1,1) + dxmin 31 continue do 32 i = icend+2,imax x(i,1) = x(i-1,1) + 2.0*(i-icend-1)*dxmin 32 continue do 33 i = 1,imax do 33 j = 2,jmax x(i,j) = x(i,1) 33 continue do 40 i = 1,imax do 40 j = 1,jmax-1 dy(i,j) = y(i,j+1)-y(i,j) 40 continue do 50 j=1,jmax do 50 i=1,imax-1 dx(i,j) = x(i+1,j)-x(i,j) 50 continue return end subroutine aandmu(phi,a,mu,dx,fsmach,gamma,qinf, & imax,jmax,idim,jdim) dimension phi(idim,jdim),a(idim,jdim),mu(idim,jdim), & dx(idim,jdim) real mu fsmach2 = fsmach**2 do 10 i = 1,imax do 10 j = 1,jmax if(i.eq.1) then dphidx = (phi(2,j)-phi(1,j))/dx(1,j) elseif(i.eq.imax) then dphidx = (phi(imax,j)-phi(imax-1,j))/dx(imax-1,j) else dphidx = (phi(i+1,j)-phi(i-1,j))/(dx(i,j)+dx(i-1,j)) endif a(i,j) = 1.0 - fsmach2 - (gamma+1.0)*fsmach2*dphidx/qinf if(a(i,j).gt.0.0) then mu(i,j) = 0.0 elseif(a(i,j).lt.0.0) then mu(i,j) = 1.0 else write(*,*) 'This cant be... A is exactly zero?' stop endif 10 continue return end subroutine allmats(i,amat,bmat,cmat,dmat,phi,a,mu,dx,dy, & fsmach,gamma,qinf,imax,jmax,idim,jdim) dimension amat(jdim),bmat(jdim),cmat(jdim),dmat(jdim), & phi(idim,jdim),a(idim,jdim),mu(idim,jdim), & dx(idim,jdim),dy(idim,jdim) real mu do 10 j = 2,jmax-1 dx1 = 2.0/(dx(i ,j)+dx(i-1,j )) dx2 = 2.0/(dx(i-1,j)+dx(i-2,j )) dy1 = 2.0/(dy(i ,j)+dy( i,j-1)) amat(j) = dy1/dy(i,j-1) bmat(j) = -(1.0-mu(i,j))*a(i,j)*dx1*(1.0/dx(i,j)+1.0/dx(i-1,j)) & +mu(i-1,j)*a(i-1,j)*dx2/dx(i-1,j) & -dy1*(1.0/dy(i,j)+1.0/dy(i,j-1)) cmat(j) = dy1/dy(i,j) dmat(j) = -(1.0-mu(i,j))*a(i,j)*dx1 & *(phi(i+1,j)/dx(i,j)+phi(i-1,j)/dx(i-1,j)) & -mu(i-1,j)*a(i-1,j)*dx2 &*(-phi(i-1,j)*(1.0/dx(i-1,j)+1.0/dx(i-2,j))+phi(i-2,j)/dx(i-2,j)) 10 continue return end subroutine output(phi,x,y,dx,dy,fsmach,gamma,qinf,rinf,pinf, & phis,up,n,rmax,imax,jmax,nmax,idim,jdim,ndim) dimension phi(idim,jdim),x(idim,jdim),y(idim,jdim), & dx(idim,jdim),dy(idim,jdim),phis(idim),up(idim), & rmax(ndim) character tab tab = char(9) gam1 = 0.5*(gamma-1.0) gam2 = gamma/(gamma-1.0) do 10 i=1,imax phis(i) = phi(i,2) 10 continue up(1) = (phis(2)-phis(1))/dx(1,2) up(imax) = (phis(imax)-phis(imax-1))/dx(imax-1,2) do 20 i=2,imax-1 up(i) = (phis(i+1)-phis(i-1))/(dx(i,2)+dx(i-1,2)) 20 continue open(unit=2,file='cp.dat',form='formatted') write(*,*) ' ' write(*,*) ' ',' i ',' X ',' Cp ' write(2,*) ' ',' i ',' X ',' Cp ' do 30 i=1,imax cps = -2.0*up(i)/qinf u = qinf + up(i) v = (phi(i,2)-phi(i,1))/(dy(i,1)) psisen = pinf*(1.0 & -gam1*fsmach**2*((u**2+v**2)/qinf**2-1.0))**gam2 cpsisen = (psisen-pinf)/(0.5*rinf*qinf**2) c write(*,*) 'x =',x(i,1),' cp =',cps,' cpisen =',cpsisen write(6,100) i,x(i,1),cps write(2,100) i,x(i,1),cps 30 continue 100 format(1x,i5,2(es11.3)) c open(unit=3,file='RMS.dat',form='formatted') c write(6,*) ' ',' N ',' RMSmax ' c write(3,*) ' ',' N ',' RMSmax ' c do 40 n = 1,nmax c write(6,200) n,rmax(n) c write(3,200) n,rmax(n) c 40 continue c 200 format(1x,i5,es11.3) write(1) imax,jmax write(1) ((x(i,j),i=1,imax),j=1,jmax), & ((y(i,j),i=1,imax),j=1,jmax) close(2) c close(3) return end subroutine scaltri(i,qmat,amat,bmat,cmat,dmat, & imax,jmax,idim,jdim) c...Scalar Tridiagonal Solver dimension qmat(idim,jdim), & amat(jdim),bmat(jdim),cmat(jdim),dmat(jdim) do 10 j = 2,jmax dumm = amat(j)/bmat(j-1) bmat(j) = bmat(j) - cmat(j-1)*dumm dmat(j) = dmat(j) - dmat(j-1)*dumm 10 continue qmat(i,jmax) = dmat(jmax)/bmat(jmax) do 20 j = jmax-1,1,-1 qmat(i,j) = (dmat(j)-cmat(j)*qmat(i,j+1))/bmat(j) 20 continue return end subroutine bcmat(i,amat,bmat,cmat,dmat,x,dy,th,qinf,icst,icend, & imax,jmax,idim,jdim) dimension x(idim,jdim),dy(idim,jdim), & amat(jdim),bmat(jdim),cmat(jdim),dmat(jdim) amat(jmax) = 0.0 bmat(jmax) = 1.0 cmat(jmax) = 0.0 dmat(jmax) = 1.0 if((i.ge.icst).and.(i.le.icend)) then tau = 0.5*th yc = (tau**2-0.25)/(2.0*tau) r2 = 0.25 + yc**2 dummx = x(i,1) - 0.5 dydx = -dummx/sqrt(r2-dummx**2) amat(1) = 0.0 bmat(1) = -1.0 cmat(1) = 1.0 dmat(1) = dy(i,1)*qinf*dydx else amat(1) = 0.0 bmat(1) = 1.0 cmat(1) = -1.0 dmat(1) = 0.0 endif return end subroutine resid(n,rmax,residm0,phi,a,mu,dx,dy,imax,jmax, &idim,jdim,ndim) dimension phi(idim,jdim),a(idim,jdim),mu(idim,jdim), & dx(idim,jdim),dy(idim,jdim),rmax(ndim) real mu residm = 0.0 if(n.eq.1) residm0 = 1.0 iresid = 0 jresid = 0 do 10 i=3,imax-1 do 10 j=2,jmax-1 dx1 = 2.0/(dx(i ,j)+dx(i-1,j )) dx2 = 2.0/(dx(i-1,j)+dx(i-2,j )) dy1 = 2.0/(dy(i ,j)+dy( i,j-1)) dumm1 = (1.0-mu(i,j))*a(i,j)*dx1 & *(phi(i+1,j)/dx(i,j)-phi(i,j)*(1.0/dx(i,j)+1.0/dx(i-1,j)) & +phi(i-1,j)/dx(i-1,j)) dumm2 = mu(i-1,j)*a(i-1,j)*dx2 & *(phi(i,j)/dx(i-1,j)-phi(i-1,j)*(1.0/dx(i-1,j)+1.0/dx(i-2,j)) & +phi(i-2,j)/dx(i-2,j)) dumm3 = dy1*(phi(i,j+1)/dy(i,j) & -phi(i,j)*(1.0/dy(i,j)+1.0/dy(i,j-1))+phi(i,j-1)/dy(i,j-1)) residt = abs(dumm1 + dumm2 + dumm3) residt = residt/residm0 if(residt.gt.residm) then residm = residt iresid = i jresid = j endif 10 continue if(n.eq.1) then residm0 = residm residm = 1.0 endif rmax(n) = residm write(6,200) n,residm,iresid,jresid write(3,200) n,residm,iresid,jresid 200 format(1x,i5,es11.3,i10,i10) return end