aboutsummaryrefslogtreecommitdiff
path: root/orrs/src/osseq.f
diff options
context:
space:
mode:
authorDimitri Sokolyuk <demon@dim13.org>2009-05-11 00:27:49 +0000
committerDimitri Sokolyuk <demon@dim13.org>2009-05-11 00:27:49 +0000
commit0d4f43d355de79178b1142e9735902cf641670b6 (patch)
tree2ced2323f6351db2a51090b3fd13eb11f69ff53f /orrs/src/osseq.f
Xfoil 6.97
Diffstat (limited to 'orrs/src/osseq.f')
-rw-r--r--orrs/src/osseq.f407
1 files changed, 407 insertions, 0 deletions
diff --git a/orrs/src/osseq.f b/orrs/src/osseq.f
new file mode 100644
index 0000000..4633667
--- /dev/null
+++ b/orrs/src/osseq.f
@@ -0,0 +1,407 @@
+
+ PROGRAM OSSEQ
+C--------------------------------------------------------------------
+C Program for executing and displaying Orr-Sommerfeld solution
+C--------------------------------------------------------------------
+C
+ PARAMETER (NMAX=2001)
+ DIMENSION ETA(NMAX), F(NMAX), U(NMAX), S(NMAX)
+ DIMENSION UTR(NMAX), UTI(NMAX), UT(NMAX),
+ & VTR(NMAX), VTI(NMAX), VT(NMAX),
+ & WTR(NMAX), WTI(NMAX), WT(NMAX),
+ & CTR(NMAX), CTI(NMAX), CT(NMAX),
+ & PTR(NMAX), PTI(NMAX), PT(NMAX)
+ DIMENSION UU(NMAX), VV(NMAX), UV(NMAX), QQ(NMAX)
+ CHARACTER*1 ANS
+ CHARACTER*80 FNAME, ARGP1
+ DIMENSION XLIN(2), YLIN(2)
+C
+ DIMENSION AINPUT(10)
+ LOGICAL ERROR
+C
+ LST = 1
+ LRE = 1
+C
+ IDEV = 1
+ IDEVRP = 2
+ IPSLU = 0
+C
+ SIZE = 6.0
+ PAR = 0.75
+C
+ CALL PLINITIALIZE
+C
+C
+ N = 2001
+ ETAE = 16.0
+ GEO = 1.01
+C
+ CH = 0.021
+C
+ IF(N.GT.NMAX) STOP 'TEST: Array overflow.'
+C
+ CALL GETARG(1,ARGP1)
+ IF(ARGP1(1:1).EQ.' ') GO TO 50
+C
+ FNAME = ARGP1
+C
+C---- try formatted read first
+ OPEN(1,FILE=FNAME,STATUS='OLD',ERR=50)
+ READ(1,*,ERR=30) N, H
+ DO I=1, N
+ READ(1,*) ETA(I), U(I), S(I)
+ ENDDO
+ CLOSE(1)
+ ETAE = ETA(N)
+ GO TO 90
+C
+C---- now try unformatted read
+ 30 CONTINUE
+ OPEN(19,FILE=FNAME,STATUS='OLD',FORM='UNFORMATTED',ERR=50)
+ READ(19,ERR=50) N, H
+ READ(19) (ETA(I),I=1, N)
+ READ(19) (U(I) ,I=1, N)
+ READ(19) (S(I) ,I=1, N)
+ ETAE = ETA(N)
+ CLOSE(19)
+ GO TO 90
+C
+C---- no argument specified or read error... get Falkner-Skan parameter
+ 50 CONTINUE
+ WRITE(*,*) 'Enter Falkner-Skan parameter Beta (or H)'
+ READ (*,*) BETA
+C
+ IF(BETA .GT. 1.0) THEN
+ write(*,*) 'Enter ETAE, GEO'
+ read (*,*) etae, geo
+ H = BETA
+ CALL FS(3,2,BU,H,N,ETAE,GEO,ETA,F,U,S,DELTA)
+C
+ ELSE
+ write(*,*) 'Enter ETAE, GEO'
+ read (*,*) etae, geo
+ BU = BETA / (2.0 - BETA)
+ CALL FS(3,1,BU,H,N,ETAE,GEO,ETA,F,U,S,DELTA)
+C
+ ENDIF
+C---------------------
+C
+ 90 CONTINUE
+ EWT = 1.0/ETAE
+C
+ DETA = 1.0
+ IF(ETAE .GT. 16.01) DETA = 2.0
+ IF(ETAE .GT. 30.01) DETA = 5.0
+ ETAE = DETA * AINT( ETAE/DETA + 0.51 )
+C
+ CALL PLOPEN(0.7,IPSLU,IDEV)
+ CALL PLOTABS(0.5,1.0,-3)
+C
+ CALL NEWFACTOR(SIZE)
+ CALL GETCOLOR(ICOL0)
+C
+ CALL PLOT(0.5*PAR,0.0,-3)
+C
+ CALL NEWPEN(1)
+ CALL XAXIS(0.0,0.0,-PAR, 0.2*PAR,0.0, 0.2,0.7*CH,1)
+ CALL YAXIS(0.0,0.0,-1.0,DETA*EWT,0.0,DETA,0.7*CH,1)
+C
+ CALL NEWPEN(2)
+ XL = -4.0*CH
+ YL = (ETAE-1.5*DETA)*EWT - 0.5*CH
+ CALL PLCHAR(XL,YL,CH,'y/ ',0.0,3)
+ CALL PLMATH(XL,YL,CH,' q',0.0,3)
+C
+C
+ UWT = PAR
+ CALL NEWPEN(4)
+ CALL XYLINE(N,U,ETA,0.0,UWT,0.0,EWT,1)
+ CALL NEWPEN(3)
+ CALL XYSYMB(N,U,ETA,0.0,UWT,0.0,EWT,0.125*CH,1)
+C
+ CALL PLFLUSH
+C
+ RE = 100.0
+ OMEGAR = 0.1
+ ALPHAR = 2.0*OMEGAR
+ ALPHAI = 0.
+C
+ WRITE(*,*) 'Enter Rtheta1, Rtheta2:'
+ READ(*,*) RTH1, RTH2
+C
+ WRITE(*,*) 'Enter wr1, wr2:'
+ READ(*,*) OMG1, OMG2
+C
+ WRITE(*,*) 'Enter number of steps:'
+ READ(*,*) NPASS
+C
+ WRITE(*,*) 'Enter initial ar, ai:'
+ READ(*,*) ALPHAR, ALPHAI
+C
+ OPEN(19,FILE='a.dat',STATUS='unknown')
+ REWIND(19)
+C
+ DO 100 IPASS=1, NPASS
+ ITMAX = 20
+C
+ FRAC = FLOAT(IPASS-1) / FLOAT(NPASS-1)
+C
+ RE = RTH1 * EXP( LOG(RTH2/RTH1) * FRAC )
+ OMEGAR = OMG1 * EXP( LOG(OMG2/OMG1) * FRAC )
+C
+ OMEGAI = 0.0
+C
+C
+ ITLIM = ITMAX
+ CALL ORRS(LST,LRE,N,ETA,U,S, RE, ITLIM,
+ & ALPHAR,ALPHAI, OMEGAR,OMEGAI,
+ & UTR,UTI,VTR,VTI,WTR,WTI,CTR,CTI, DELMAX)
+C
+ CALL OSPRES(N,ETA,U, ALPHAR,ALPHAI, VTR,VTI, PTR,PTI )
+C
+ DO I=1, N
+ UT(I) = SQRT(UTR(I)**2 + UTI(I)**2)
+ VT(I) = SQRT(VTR(I)**2 + VTI(I)**2)
+ PT(I) = SQRT(PTR(I)**2 + PTI(I)**2)
+ UU(I) = 0.5*(UTR(I)*UTR(I) + UTI(I)*UTI(I))
+ VV(I) = 0.5*(VTR(I)*VTR(I) + VTI(I)*VTI(I))
+ UV(I) = 0.5*(UTR(I)*VTR(I) + UTI(I)*VTI(I))
+ QQ(I) = UU(I) + VV(I)
+ ENDDO
+C
+ QTHIK = 0.
+ DQTDX = 0.
+ UPRES = 0.
+ TWORK = 0.
+ DISS1 = 0.
+ DISS2 = 0.
+ DISS3 = 0.
+ PQINT = 0.
+ DO I = 2, N
+ UA = (U(I) + U(I-1))*0.5
+ DU = U(I) - U(I-1)
+ DY = ETA(I) - ETA(I-1)
+C
+ URA = (UTR(I) + UTR(I-1))*0.5
+ UIA = (UTI(I) + UTI(I-1))*0.5
+ VRA = (VTR(I) + VTR(I-1))*0.5
+ VIA = (VTI(I) + VTI(I-1))*0.5
+ WRA = (WTR(I) + WTR(I-1))*0.5
+ WIA = (WTI(I) + WTI(I-1))*0.5
+ PRA = (PTR(I) + PTR(I-1))*0.5
+ PIA = (PTI(I) + PTI(I-1))*0.5
+C
+ DUR = UTR(I) - UTR(I-1)
+ DUI = UTI(I) - UTI(I-1)
+ DVR = VTR(I) - VTR(I-1)
+ DVI = VTI(I) - VTI(I-1)
+ DWR = WTR(I) - WTR(I-1)
+ DWI = WTI(I) - WTI(I-1)
+C
+ QTHIK = QTHIK + 0.25*(UU(I)+UU(I-1)
+ & +VV(I)+VV(I-1))*UA*DY
+C
+ UDUDX = - (ALPHAI*URA + ALPHAR*UIA)*URA
+ & + (ALPHAR*URA - ALPHAI*UIA)*UIA
+ VDVDX = - (ALPHAI*VRA + ALPHAR*VIA)*VRA
+ & + (ALPHAR*VRA - ALPHAI*VIA)*VIA
+C
+ PDUDX = - (ALPHAI*URA + ALPHAR*UIA)*PRA
+ & + (ALPHAR*URA - ALPHAI*UIA)*PIA
+ UDPDX = - (ALPHAI*PRA + ALPHAR*PIA)*URA
+ & + (ALPHAR*PRA - ALPHAI*PIA)*UIA
+C
+ DQTDX = DQTDX + 0.5*(UDUDX + VDVDX)*DY * UA
+C
+ UPRES = UPRES - 0.5*(UDPDX + PDUDX)*DY
+C
+ TWORK = TWORK - 0.50*(UV(I)+UV(I-1))*DU
+C
+ DISS1 = DISS1 + ( ALPHAI*URA + ALPHAR*UIA )**2 * DY
+ & + ( ALPHAR*URA - ALPHAI*UIA )**2 * DY
+C
+ DISS2 = DISS2 + ( DVR**2 + DVI**2 ) / DY
+C
+ DISS3 = DISS3
+ & + 0.5 * ( DUR/DY - ALPHAI*VRA - ALPHAR*VIA )**2 * DY
+ & + 0.5 * ( DUI/DY + ALPHAR*VRA - ALPHAI*VIA )**2 * DY
+C
+ PQINT = PQINT
+ & - 0.5*URA*DWR - 0.5*(ALPHAI*WRA + ALPHAR*WIA)*VRA * DY
+ & - 0.5*UIA*DWI + 0.5*(ALPHAR*WRA - ALPHAI*WIA)*VIA * DY
+ ENDDO
+C
+ DISS1 = DISS1 / RE
+ DISS2 = DISS2 / RE
+ DISS3 = DISS3 / RE
+ PQINT = PQINT / RE
+C
+C
+ DQTDX = DQTDX / QTHIK
+ UPRES = UPRES / QTHIK
+ TWORK = TWORK / QTHIK
+ DISS1 = DISS1 / QTHIK
+ DISS2 = DISS2 / QTHIK
+ DISS3 = DISS3 / QTHIK
+ PQINT = PQINT / QTHIK
+C
+ DISS = DISS1 + DISS2 + DISS3
+C
+ WRITE(*,*)
+ WRITE(*,*) 'dEdx, P+Dx+D :',DQTDX,TWORK+UPRES+PQINT
+ WRITE(*,*) 'P Dx D e:', TWORK, UPRES, PQINT, -DISS
+ WRITE(*,*)
+C
+ WRITE(19,9944) RE, OMEGAR, DQTDX, TWORK, UPRES, PQINT
+ 9944 FORMAT(1X,8E14.5)
+C
+C
+ IF(IPASS.EQ.1) THEN
+ CALL SCALIT(N,UT,0.0,USF,ANN,NANN)
+ UWT = PAR*USF
+C
+ CALL SCALIT(N,VT,0.0,VSF,ANN,NANN)
+ VWT = PAR*VSF
+C
+ CALL SCALIT(N,PT,0.0,PSF,ANN,NANN)
+ PWT = PAR*PSF
+C
+ CALL SCALIT(N,QQ,0.0,TSF,ANN,NANN)
+ TWT = PAR*TSF
+C
+ EOFF = 0.
+ UOFF = 0.
+ POFF = 0.
+ TOFF = 0.
+C
+ PWT = UWT
+ ENDIF
+C
+C
+ CALL NEWPEN(3)
+C
+ XL = PAR + 5.0*CH
+C
+ YL = ETAE*EWT
+ CALL PLCHAR(XL ,YL , CH,'H = ',0.0, 8)
+ CALL PLNUMB(XL+8.0*CH,YL , CH, H ,0.0, 3)
+C
+ YL = YL - 3.5*CH
+ CALL PLCHAR(XL ,YL , CH,'Re = ',0.0, 8)
+ CALL PLMATH(XL+1.9*CH,YL-0.4*CH,0.8*CH, 'q' ,0.0, 1)
+ CALL PLNUMB(XL+8.0*CH,YL , CH, RE ,0.0,-1)
+C
+ YL = YL - 2.5*CH
+ CALL PLMATH(XL ,YL , CH,'w q/ = ',0.0, 8)
+ CALL PLCHAR(XL ,YL , CH,' U ',0.0, 8)
+ CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 'r' ,0.0, 1)
+ CALL PLNUMB(XL+8.0*CH,YL , CH, OMEGAR ,0.0, 5)
+C
+ YL = YL - 3.5*CH
+ CALL PLMATH(XL ,YL , CH,'a q = ',0.0, 8)
+ CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 'r' ,0.0, 1)
+ CALL PLNUMB(XL+8.0*CH,YL , CH, ALPHAR ,0.0, 5)
+C
+ YL = YL - 2.5*CH
+ CALL PLMATH(XL ,YL , CH,'a q = ',0.0, 8)
+ CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 'i' ,0.0, 1)
+ CALL PLNUMB(XL+8.0*CH,YL , CH, ALPHAI ,0.0, 5)
+C
+C
+C
+ XLIN(1) = -7.5*CH
+ XLIN(2) = -1.5*CH
+ YLIN(1) = 0.5*CH
+ YLIN(2) = 0.5*CH
+C
+ CALL NEWPEN(2)
+C
+ XL = PAR + 12.0*CH
+ YL = 0.50*ETAE*EWT
+C
+ CALL NEWCOLORNAME('red')
+ CALL XYLINE(N,UTR,ETA,UOFF,UWT,EOFF,EWT,2)
+ CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,2)
+ CALL PLCHAR(XL ,YL , CH,'u /U',0.0,4)
+ CALL PLMATH(XL ,YL , CH,' ` ',0.0,4)
+ CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 'r' ,0.0,1)
+C
+ YL = YL - 2.5*CH
+C
+ CALL NEWCOLORNAME('orange')
+ CALL XYLINE(N,UTI,ETA,UOFF,UWT,EOFF,EWT,3)
+ CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,3)
+ CALL PLCHAR(XL ,YL , CH,'u /U',0.0,4)
+ CALL PLMATH(XL ,YL , CH,' ` ',0.0,4)
+ CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 'i' ,0.0,1)
+C
+ YL = YL - 2.5*CH
+cC
+c CALL NEWCOLORNAME('yellow')
+c CALL XYLINE(N,UT ,ETA,UOFF,UWT,EOFF,EWT,1)
+c CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,1)
+c CALL PLCHAR(XL ,YL , CH,'u /U',0.0,4)
+c CALL PLMATH(XL ,YL , CH,' ` ',0.0,4)
+c CALL PLCHAR(XL-0.6*CH,YL ,0.9*CH,'| |' ,0.0,3)
+cC
+cC
+c YL = YL - 3.5*CH
+cC
+c CALL NEWCOLORNAME('violet')
+c CALL XYLINE(N,PTR,ETA,POFF,PWT,EOFF,EWT,5)
+c CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,5)
+c CALL PLCHAR(XL ,YL , CH,'p / U ',0.0,6)
+c CALL PLMATH(XL ,YL , CH,' ` r 2',0.0,6)
+c CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 'r' ,0.0,1)
+cC
+c YL = YL - 2.5*CH
+cC
+c CALL NEWCOLORNAME('blue')
+c CALL XYLINE(N,PTI,ETA,POFF,PWT,EOFF,EWT,6)
+c CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,6)
+c CALL PLCHAR(XL ,YL , CH,'p / U ',0.0,6)
+c CALL PLMATH(XL ,YL , CH,' ` r 2',0.0,6)
+c CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 'i' ,0.0,1)
+cC
+c YL = YL - 2.5*CH
+cC
+ CALL NEWCOLORNAME('cyan')
+ CALL XYLINE(N,PT ,ETA,POFF,PWT,EOFF,EWT,1)
+ CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,1)
+ CALL PLCHAR(XL ,YL , CH,'p / U ',0.0,6)
+ CALL PLMATH(XL ,YL , CH,' ` r 2',0.0,6)
+ CALL PLCHAR(XL-0.6*CH,YL ,0.9*CH,'| |' ,0.0,3)
+C
+ YL = YL - 3.5*CH
+C
+c CALL NEWCOLORNAME('green')
+c CALL XYLINE(N,QQ ,ETA,TOFF,TWT,EOFF,EWT,2)
+c CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,2)
+c CALL PLMATH(XL ,YL+0.2*CH, CH,' ___ ' ,0.0,5)
+c CALL PLMATH(XL ,YL , CH,' ` ` 2',0.0,8)
+c CALL PLCHAR(XL ,YL , CH,' q q /U ',0.0,8)
+cC
+c YL = YL - 2.5*CH
+cC
+c CALL NEWCOLORNAME('green')
+c CALL XYLINE(N,UV ,ETA,TOFF,-10.0*TWT,EOFF,EWT,1)
+c CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,1)
+c CALL PLMATH(XL ,YL+0.2*CH, CH,' ___ ',0.0,12)
+c CALL PLMATH(XL ,YL , CH,' ` ` 2 # ',0.0,12)
+c CALL PLCHAR(XL ,YL , CH,'-u v /U 10',0.0,12)
+cC
+ CALL NEWCOLOR(ICOL0)
+C
+ CALL PLFLUSH
+C
+ 100 CONTINUE
+ 101 CONTINUE
+C
+ CLOSE(19)
+C
+ CALL PLCLOSE
+ STOP
+ END
+
+
+