From 0d4f43d355de79178b1142e9735902cf641670b6 Mon Sep 17 00:00:00 2001 From: Dimitri Sokolyuk Date: Mon, 11 May 2009 00:27:49 +0000 Subject: Xfoil 6.97 --- src/profil.f | 1034 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1034 insertions(+) create mode 100644 src/profil.f (limited to 'src/profil.f') diff --git a/src/profil.f b/src/profil.f new file mode 100644 index 0000000..91cb5eb --- /dev/null +++ b/src/profil.f @@ -0,0 +1,1034 @@ +C*********************************************************************** +C Module: profil.f +C +C Copyright (C) 2000 Mark Drela +C +C This program is free software; you can redistribute it and/or modify +C it under the terms of the GNU General Public License as published by +C the Free Software Foundation; either version 2 of the License, or +C (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program; if not, write to the Free Software +C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +C*********************************************************************** + + SUBROUTINE PRWALL(DSTAR,THETA,UO,RT,MS,CT, BB, + & DO, DO_DS, DO_TH, DO_UO, DO_RT, DO_MS, + & UI, UI_DS, UI_TH, UI_UO, UI_RT, UI_MS, + & HS, HS_DS, HS_TH, HS_UO, HS_RT, HS_MS, + & CF, CF_DS, CF_TH, CF_UO, CF_RT, CF_MS, + & CD, CD_DS, CD_TH, CD_UO, CD_RT, CD_MS, CD_CT ) + IMPLICIT REAL (A-H,M,O-Z) +C================================================================ +C Returns wall slip velocity and thickness of wall BL profile +C +C Input: +C DSTAR kinematic displacement thickness +C THETA kinematic momentum thickness +C RT momentum thickness based on ue and THETA +C MS Mach^2 based on ue +C +C UO uo/ue outer velocity; assumed = 1 in this version +C +C Output: +C BB outer profile exponent +C DO thickness of profile deck +C UI inner "slip" velocity +C CF wall skin friction +C================================================================ +C + PARAMETER (N=65) + DIMENSION ETA(N), UIP(N), UIP_DP(N), G(N), G_BB(N) +C +C---- pi/2 , 2/pi + DATA HPI, TOPI / 1.570796327 , 0.6366197723 / +C + DATA T, SQT / 0.28 , 0.5291502622 / +C +C---- TCON = ( atan(T^1/2) / T^1/2 - 1/(1+T)) / (2T) + 0.5/(1.0 + T) +C - atan(1/T^1/2) / 2T^1/2 + DATA TCON / -0.3864027035 / +C +C---- slip velocity coefficient + DATA AK / 0.09 / +C +C---- log-law constants + DATA VKAP, VB / 0.40 , 5.0 / +C + HK = DSTAR/THETA +C + UO = 1.0 + BB = 1.0 +C +C---- initialize variables + CALL CFT(HK,RT,MS,CF,CF_HK,CF_RT,CF_MS) + SGN = SIGN( 1.0 , CF ) + UT = SGN * SQRT(0.5*ABS(CF)) +C + UI = MIN( UT/AK * HPI , 0.90 ) + DO = HK*THETA / (1.0 - 0.5*(UO+UI)) +C + EBK = EXP(-VB*VKAP) +C + DO 1000 ITER=1, 12 +C + SGN = SIGN( 1.0 , UT ) +C +C------ set d+ = DP(UT DO ; RT TH) + DP = SGN * UT*RT*DO/THETA + DP_DO = SGN * UT*RT /THETA + DP_UT = SGN * RT*DO/THETA +C + DP_TH = -DP/THETA + DP_RT = SGN * UT *DO/THETA + DP_MS = 0.0 +C +C------ estimate inner profile edge velocity Ui+ using log-law + UPE = LOG(DP)/VKAP + VB +C +C------ converge exact Ui+ using Spalding formula + DO 10 ITUP=1, 5 + UK = UPE*VKAP + ARG = UK - VB*VKAP + EXU = EXP(ARG) + REZ = UPE + EXU - EBK*(1.0 + UK + UK*UK/2.0 + UK*UK*UK/6.0) + & - DP + DP_U = 1.0 + (EXU - EBK*(1.0 + UK + UK*UK/2.0))*VKAP +C + IF(ABS(REZ/DP) .LT. 1.0E-5) GO TO 11 +C + DUPE = -REZ/DP_U + UPE = UPE + DUPE + 10 CONTINUE + WRITE(*,*) 'PRWALL: Ue+ convergence failed, Res =', REZ/DP + 11 CONTINUE +C + UPE_DP = 1.0/DP_U +C +C 2 2 3 3 +C------ set d y+/du+ and d y+/du+ at BL edge + DP_UU = (EXU - EBK*(1.0 + UK))*VKAP**2 + DP_UUU = (EXU - EBK )*VKAP**3 +C +C------ set du+/dy+ at BL edge + UPD = 1.0/DP_U + UPD_DP = (-1.0/DP_U**3) * DP_UU +C +C 2 2 +C------ set d u+/dy+ at BL edge +CCC UPD_DP = (-1.0/DP_U**3) * DP_UU + UPDD = UPD_DP + UPDD_DP = (-1.0/DP_U**4) * DP_UUU + & + ( 3.0/DP_U**5) * DP_UU**2 +C +C------ set coefficients for Spalding profile correction polynomial + DC2 = 0.5*DP*DP*UPDD - DP*UPD + DC2_DP = DP*UPDD - UPD + & + 0.5*DP*DP*UPDD_DP - DP*UPD_DP +C + DC3 = -( DP*DP*UPDD - DP*UPD ) / 3.0 + DC3_DP = -(2.0 *DP*UPDD - UPD ) / 3.0 + & -( DP*DP*UPDD_DP - DP*UPD_DP) / 3.0 +C +C------ set outer profile amplitude DUO + DUO = UO - UT*(UPE + DC2 + DC3 ) + DUO_DP = - UT*(UPE_DP + DC2_DP + DC3_DP) +C + DUO_UT = - (UPE + DC2 + DC3 ) + & + DUO_DP*DP_UT + DUO_DO = DUO_DP*DP_DO +C + DUO_TH = DUO_DP*DP_TH + DUO_RT = DUO_DP*DP_RT + DUO_MS = DUO_DP*DP_MS +c +c write(*,*) 'dUo', duo, duo_dp, duo_ut, duo_do +c read(*,*) ddo, dut +c if(ddo.ne.0.0 .or. dut.ne.0.0) then +c do = do+ddo +c ut = ut+dut +c write(*,*) 'new', duo + duo_do*ddo + duo_ut*dut +c go to 666 +c endif +C +C------ set wake profile coefficients + BB1 = 3.0*(BB +2.0)*(BB+3.0)/(BB+7.0) + BB1_BB = 3.0*(BB*2.0+5.0 )/(BB+7.0) - BB1/(BB+7.0) + BB2 = -5.0*(BB +1.0)*(BB+3.0)/(BB+7.0) + BB2_BB = -5.0*(BB*2.0+4.0 )/(BB+7.0) - BB2/(BB+7.0) + BB3 = 2.0*(BB +1.0)*(BB+2.0)/(BB+7.0) + BB3_BB = 2.0*(BB*2.0+3.0 )/(BB+7.0) - BB3/(BB+7.0) +C +C------ fill eta coordinate and inner profile arrays +CCC EXUPE = EXP(UPE*VKAP - VB*VKAP) + EXUPE = EXU +C + DEXU = (EXUPE - EBK)/FLOAT(N-1) +C + I = 1 + UIP(I) = 0.0 + UIP_DP(I) = 0.0 + G(I) = 0.0 + G_BB(I) = 0.0 +C + DO 20 I=2, N +ccc EXU = EBK + DEXU*FLOAT(I-1) + EXU = EBK + (DEXU - 0.75*DEXU*FLOAT(N-I)/FLOAT(N-1)) + & *FLOAT(I-1) +C +CCC UK = UP*VKAP + UK = LOG(EXU) + VB*VKAP +C + UP = UK/VKAP +C +C-------- set "inverse" Spalding profile y+(u+) and derivatives + YP = UP + EXU - EBK*(1.0 + UK + UK*UK/2.0 + UK*UK*UK/6.0) + YP_U = 1.0 + (EXU - EBK*(1.0 + UK + UK*UK/2.0))*VKAP + YP_UU = (EXU - EBK*(1.0 + UK ))*VKAP**2 +C + ET = YP/DP +C +C-------- set final inner profile (fudged Spalding) + UIP(I) = UP + DC2 *ET**2 + DC3 *ET**3 + UIP_DP(I) = DC2_DP*ET**2 + DC3_DP*ET**3 +C +ccc UIPD(I) = 1.0/YP_U + 2.0*DC2 *ET + 3.0*DC3 *ET**2 +ccc UIPD_DP(I) = (-1.0/YP_U**3)*YPUU +ccc & + 2.0*DC2_DP*ET + 3.0*DC3_DP*ET**2 +C +C-------- set outer profile + ETB = ET**BB + ALE = LOG(ET) +C +ccc G(I) = 2.0*ETB - ETB**2 +ccc G_BB(I) = (2.0*ETB - 2.0*ETB**2)*ALE +C + G(I) = (BB1 *ET + BB2 *ET**2 + BB3 *ET**3)*ETB + G_BB(I) = (BB1_BB*ET + BB2_BB*ET**2 + BB3_BB*ET**3)*ETB + & + G(I)*ALE +C + ETA(I) = ET +C + 20 CONTINUE +C +C + DSN = 0.0 + DSN_DO = 0.0 + DSN_UT = 0.0 + DSN_BB = 0.0 +C + DSN_TH = 0.0 + DSN_RT = 0.0 + DSN_MS = 0.0 +C +C + THN = 0.0 + THN_DO = 0.0 + THN_UT = 0.0 + THN_BB = 0.0 +C + THN_TH = 0.0 + THN_RT = 0.0 + THN_MS = 0.0 +C +c TSN = 0.0 +c TSN_DO = 0.0 +c TSN_UT = 0.0 +c TSN_BB = 0.0 +c +c TSN_TH = 0.0 +c TSN_RT = 0.0 +c TSN_MS = 0.0 +C +C------ perform integration + DO 100 I=1, N-1 + DETA = ETA(I+1) - ETA(I) + GA = 0.5*(G(I+1) + G(I) ) + GA_BB = 0.5*(G_BB(I+1) + G_BB(I)) +C + UIPA = 0.5*(UIP(I+1) + UIP(I) ) + UIPA_DP = 0.5*(UIP_DP(I+1) + UIP_DP(I)) +C + U = UT*UIPA + DUO *GA + U_DP = UT*UIPA_DP +C + U_DO = DUO_DO*GA + U_DP*DP_DO + U_UT = UIPA + DUO_UT*GA + U_DP*DP_UT + U_BB = DUO *GA_BB +C + U_TH = DUO_TH*GA + U_DP*DP_TH + U_RT = DUO_RT*GA + U_DP*DP_RT + U_MS = DUO_MS*GA + U_DP*DP_MS +C +C + DSN = DSN + (1.0 - U )*DETA + DSN_DO = DSN_DO - U_DO *DETA + DSN_UT = DSN_UT - U_UT *DETA + DSN_BB = DSN_BB - U_BB *DETA +C + DSN_TH = DSN_TH - U_TH *DETA + DSN_RT = DSN_RT - U_RT *DETA + DSN_MS = DSN_MS - U_MS *DETA +C +C + THN = THN + (U - U*U) *DETA + THN_DO = THN_DO + (1.0 - 2.0*U)*U_DO*DETA + THN_UT = THN_UT + (1.0 - 2.0*U)*U_UT*DETA + THN_BB = THN_BB + (1.0 - 2.0*U)*U_BB*DETA +C + THN_TH = THN_TH + (1.0 - 2.0*U)*U_TH*DETA + THN_RT = THN_RT + (1.0 - 2.0*U)*U_RT*DETA + THN_MS = THN_MS + (1.0 - 2.0*U)*U_MS*DETA +C +c TSN = TSN + (U - U*U*U) *DETA +c TSN_DO = TSN_DO + (1.0 - 3.0*U*U)*U_DO*DETA +c TSN_UT = TSN_UT + (1.0 - 3.0*U*U)*U_UT*DETA +c TSN_BB = TSN_BB + (1.0 - 3.0*U*U)*U_BB*DETA +C +c TSN_TH = TSN_TH + (1.0 - 3.0*U*U)*U_TH*DETA +c TSN_RT = TSN_RT + (1.0 - 3.0*U*U)*U_RT*DETA +c TSN_MS = TSN_MS + (1.0 - 3.0*U*U)*U_MS*DETA +C + 100 CONTINUE +C +C------ set up 2x2 system for DO UT + REZ1 = DO*DSN - THETA*HK + A11 = DO*DSN_DO + DSN + A12 = DO*DSN_UT +cc A12 = DO*DSN_BB +C + REZ2 = DO*THN - THETA + A21 = DO*THN_DO + THN + A22 = DO*THN_UT +cc A22 = DO*THN_BB +C +cc IF(ABS(REZ1/THETA) .LT. 2.0E-5 .AND. +cc & ABS(REZ2/THETA) .LT. 2.0E-5 ) GO TO 1010 + IF(ABS(REZ1/THETA) .LT. 1.0E-3 .AND. + & ABS(REZ2/THETA) .LT. 1.0E-3 ) GO TO 1010 +C + DET = A11*A22 - A12*A21 + B11 = A22/DET + B12 = -A12/DET + B21 = -A21/DET + B22 = A11/DET +C + DDO = -(B11*REZ1 + B12*REZ2) + DUT = -(B21*REZ1 + B22*REZ2) +cc DBB = -(B21*REZ1 + B22*REZ2) +C + DMAX = MAX( ABS(DDO/DO) , ABS(DUT/0.05) ) +cc DMAX = MAX( ABS(DDO/DO) , ABS(DBB/BB ) ) + RLX = 1.0 + IF(DMAX.GT.0.5) RLX = 0.5/DMAX +C + DO = DO + RLX*DDO + UT = UT + RLX*DUT +cc BB = BB + RLX*DBB +c +cc write(*,*) iter, do, ut, rez1, rez2 +cc write(*,*) iter, do, bb, rez1, rez2 +C + 1000 CONTINUE +C + WRITE(*,*) 'PRWALL: Convergence failed. Res =', REZ1, REZ2 +C + 1010 CONTINUE +C +C +C +CCC REZ1 = DO*DSN - THETA*HK + Z1_HK = - THETA + Z1_TH = DO*DSN_TH - HK + Z1_RT = DO*DSN_RT + Z1_MS = DO*DSN_MS +C +CCC REZ2 = DO*THN - THETA + Z2_HK = 0.0 + Z2_TH = DO*THN_TH - 1.0 + Z2_RT = DO*THN_RT + Z2_MS = DO*THN_MS +C + DO_HK = -(B11*Z1_HK + B12*Z2_HK) + DO_TH = -(B11*Z1_TH + B12*Z2_TH) + DO_RT = -(B11*Z1_RT + B12*Z2_RT) + DO_MS = -(B11*Z1_MS + B12*Z2_MS) +C + UT_HK = 0.0 + UT_TH = 0.0 + UT_RT = 0.0 + UT_MS = 0.0 +C + BB_HK = 0.0 + BB_TH = 0.0 + BB_RT = 0.0 + BB_MS = 0.0 +C + UT_HK = -(B21*Z1_HK + B22*Z2_HK) + UT_TH = -(B21*Z1_TH + B22*Z2_TH) + UT_RT = -(B21*Z1_RT + B22*Z2_RT) + UT_MS = -(B21*Z1_MS + B22*Z2_MS) +C +cc BB_HK = -(B21*Z1_HK + B22*Z2_HK) +cc BB_TH = -(B21*Z1_TH + B22*Z2_TH) +cc BB_RT = -(B21*Z1_RT + B22*Z2_RT) +cc BB_MS = -(B21*Z1_MS + B22*Z2_MS) +C +C +C---- set and linearize Cf + CF = SGN*2.0*UT**2 + CF_UT = SGN*4.0*UT + CF_DO = 0.0 +C + CF_HK = CF_UT*UT_HK + CF_DO*DO_HK + CF_TH = CF_UT*UT_TH + CF_DO*DO_TH + CF_RT = CF_UT*UT_RT + CF_DO*DO_RT + CF_MS = CF_UT*UT_MS + CF_DO*DO_MS +C +C +C---- set and linearize "slip" velocity UI = UI( DUO(DO UT TH RT MS) ) + UI = UO - DUO + UI_UT = - DUO_UT + UI_DO = - DUO_DO +C + UI_HK = UI_UT*UT_HK + UI_DO*DO_HK + UI_TH = UI_UT*UT_TH + UI_DO*DO_TH - DUO_TH + UI_RT = UI_UT*UT_RT + UI_DO*DO_RT - DUO_RT + UI_MS = UI_UT*UT_MS + UI_DO*DO_MS - DUO_MS +C + RETURN + END ! PRWALL + + + + SUBROUTINE UWALL(TH,UO,DO,UI,RT,CF,BB, Y,U,N) +C------------------------------------------ +C Returns wall BL profile U(Y). +C +C Input: +C TH kinematic momentum thickness +C UO uo/ue outer velocity (= 1 for normal BL) +C DO BL thickness +C UI inner "slip" velocity +C RT momentum thickness based on ue and THETA +C CF wall skin friction +C BB outer profile exponent +C N number of profile array points +C +C Output: +C Y(i) normal coordinate array +C U(i) u/ue velocity profile array +C------------------------------------------- +C + IMPLICIT REAL (A-H,M,O-Z) + DIMENSION Y(N), U(N) + DATA HPI / 1.570796327 / + DATA AK / 0.09 / +C +C---- log-law constants + DATA VKAP, VB / 0.40 , 5.0 / +C + EBK = EXP(-VB*VKAP) +C + SGN = SIGN( 1.0 , CF ) + UT = SGN * SQRT(0.5*ABS(CF)) +C +C +C---- set d+ = DP(UT DO ; RT TH) + DP = SGN * UT*RT*DO/TH +C +C---- estimate inner profile edge velocity Ui+ using log-law + UPE = LOG(DP)/VKAP + VB +C +C---- converge exact Ui+ using Spalding formula + DO 10 ITUP=1, 5 + UK = UPE*VKAP + ARG = UK - VB*VKAP + EXU = EXP(ARG) + REZ = UPE + EXU - EBK*(1.0 + UK + UK*UK/2.0 + UK*UK*UK/6.0) + & - DP + DP_U = 1.0 + (EXU - EBK*(1.0 + UK + UK*UK/2.0))*VKAP +C + IF(ABS(REZ/DP) .LT. 1.0E-5) GO TO 11 +C + DUPE = -REZ/DP_U + UPE = UPE + DUPE + 10 CONTINUE + WRITE(*,*) 'UWALL: Ue+ convergence failed, Res =', REZ/DP + 11 CONTINUE +C +C 2 2 3 3 +C---- set d y+/du+ and d y+/du+ at BL edge + DP_UU = (EXU - EBK*(1.0 + UK))*VKAP**2 + DP_UUU = (EXU - EBK )*VKAP**3 +C +C 2 2 +C---- set du+/dy+ and d u+/dy+ at BL edge + UPD = 1.0/DP_U + UPDD = (-1.0/DP_U**3) * DP_UU +C +C---- set coefficients for Spalding profile correction polynomial + DC2 = 0.5*DP*DP*UPDD - DP*UPD + DC3 = -( DP*DP*UPDD - DP*UPD ) / 3.0 +C +C---- set outer profile amplitude DUO + DUO = UO - UT*(UPE + DC2 + DC3) +C +C + BB1 = 3.0*(BB +2.0)*(BB+3.0)/(BB+7.0) + BB2 = -5.0*(BB +1.0)*(BB+3.0)/(BB+7.0) + BB3 = 2.0*(BB +1.0)*(BB+2.0)/(BB+7.0) +C +c NE = (N*9)/10 + NE = N +C +C---- fill Y coordinate and U profile arrays +CCC EXUPE = EXP(UPE*VKAP - VB*VKAP) + EXUPE = EXU +C + DEXU = (EXUPE - EBK)/FLOAT(NE-1) +C + I = 1 + Y(I) = 0.0 + U(I) = 0.0 + DO 20 I=2, NE +ccc EXU = EBK + DEXU*FLOAT(I-1) + EXU = EBK + (DEXU - 0.75*DEXU*FLOAT(NE-I)/FLOAT(NE-1)) + & *FLOAT(I-1) +C +CCC UK = UP*VKAP + UK = LOG(EXU) + VB*VKAP +C + UP = UK/VKAP +C +C------ set "inverse" Spalding profile y+(u+) + YP = UP + EXU - EBK*(1.0 + UK + UK*UK/2.0 + UK*UK*UK/6.0) + YP_UP = 1.0 + (EXU - EBK*(1.0 + UK + UK*UK/2.0))*VKAP +C + ET = YP/DP +C +C------ set final inner profile (fudged Spalding) + UIP = UP + DC2*ET**2 + DC3*ET**3 +C +C------ set outer profile + SQE = SQRT(ET) + ETB = ET**BB +C +ccc G = 2.0*ETB - ETB**2 +C + G = (BB1 *ET + BB2 *ET**2 + BB3 *ET**3)*ETB +C + Y(I) = ET*DO + U(I) = UT*UIP + DUO*G +c + 20 CONTINUE +C +c DETA = 0.1 / FLOAT(N - NE - 1) +c DO 300 I=NE+1, N +c ETA = 1.0 + DETA*FLOAT(I-NE) +c Y(I) = DO*ETA +c U(I) = 1.0 +c 300 CONTINUE +C + RETURN + END ! UWALL + + + + SUBROUTINE FS(INORM,ISPEC,BSPEC,HSPEC,N,ETAE,GEO,ETA,F,U,S,DELTA) + IMPLICIT REAL (A-H,M,O-Z) + DIMENSION ETA(N), F(N), U(N), S(N) +C----------------------------------------------------- +C Routine for solving the Falkner-Skan equation. +C +C Input: +C ------ +C INORM 1: eta = y / sqrt(vx/Ue) "standard" Falkner-Skan coordinate +C 2: eta = y / sqrt(2vx/(m+1)Ue) Hartree's coordinate +C 3: eta = y / Theta momentum thickness normalized coordinate +C ISPEC 1: BU = x/Ue dUe/dx ( = "m") specified +C 2: H12 = Dstar/Theta specified +C BSPEC specified pressure gradient parameter (if ISPEC = 1) +C HSPEC specified shape parameter of U profile (if ISPEC = 2) +C N total number of points in profiles +C ETAE edge value of normal coordinate +C GEO exponential stretching factor for ETA: +C +C Output: +C ------- +C BSPEC calculated pressure gradient parameter (if ISPEC = 2) +C HSPEC calculated shape parameter of U profile (if ISPEC = 1) +C ETA normal BL coordinate +C F,U,S Falkner Skan profiles +C DELTA normal coordinate scale y = eta * Delta +C----------------------------------------------------- +C + PARAMETER (NMAX=257,NRMAX=3) + DIMENSION A(3,3,NMAX),B(3,3,NMAX),C(3,3,NMAX), + & R(3,NRMAX,NMAX) +C +C---- set number of righthand sides. + DATA NRHS / 3 / +C + ITMAX = 20 +C + IF(N.GT.NMAX) STOP 'FS: Array overflow.' +C + PI = 4.0*ATAN(1.0) +C +CCC if(u(n) .ne. 0.0) go to 9991 + +c +C---- initialize H or BetaU with empirical curve fits + IF(ISPEC.EQ.1) THEN + H = 2.6 + BU = BSPEC + ELSE + H = HSPEC + IF(H .LE. 2.17) THEN + WRITE(*,*) 'FS: Specified H too low' + H = 2.17 + ENDIF + BU = (0.058*(H-4.0)**2/(H-1.0) - 0.068) / (6.54*H - 14.07) * H**2 + IF(H .GT. 4.0) BU = MIN( BU , 0.0 ) + ENDIF +C +C---- initialize TN = Delta^2 Ue / vx + IF(INORM.EQ.3) THEN + TN = (6.54*H - 14.07) / H**2 + ELSE + TN = 1.0 + ENDIF +C +C---- set eta array + DETA = 1.0 + ETA(1) = 0.0 + DO 5 I=2, N + ETA(I) = ETA(I-1) + DETA + DETA = GEO*DETA + 5 CONTINUE +C + DO 6 I=1, N + ETA(I) = ETA(I) * ETAE/ETA(N) + 6 CONTINUE +C +C +C---- initial guess for profiles using a sine loop for U for half near wall + IF(H .LE. 3.0) THEN +C + IF(INORM.EQ.3) THEN + ETJOIN = 7.3 + ELSE + ETJOIN = 5.0 + ENDIF +C + EFAC = 0.5*PI/ETJOIN + DO 10 I=1, N + U(I) = SIN(EFAC*ETA(I)) + F(I) = 1.0/EFAC * (1.0 - COS(EFAC*ETA(I))) + S(I) = EFAC*COS(EFAC*ETA(I)) + IF(ETA(I) .GT. ETJOIN) GO TO 11 + 10 CONTINUE + 11 CONTINUE + IJOIN = I +C +C----- constant U for outer half + DO 12 I=IJOIN+1, N + U(I) = 1.0 + F(I) = F(IJOIN) + ETA(I) - ETA(IJOIN) + S(I) = 0. + 12 CONTINUE +C + ELSE +C + IF(INORM.EQ.3) THEN + ETJOIN1 = 0.0 + ETJOIN2 = 8.0 + IF(H .GT. 4.0) THEN + ETJOIN1 = H - 4.0 + ETJOIN2 = ETJOIN1 + 8.0 + ENDIF + ELSE + ETJOIN1 = 0.0 + ETJOIN2 = 8.0 + ENDIF +C + DO 13 I=1, N + U(I) = 0.0 + S(I) = 0.0 + F(I) = 0.0 + IF(ETA(I) .GE. ETJOIN1) GO TO 14 + 13 CONTINUE + 14 CONTINUE + IJOIN = I +C + EFAC = 0.5*PI/(ETJOIN2-ETJOIN1) + DO 15 I=IJOIN+1, N + EBAR = ETA(I) - ETJOIN1 + U(I) = 0.5 - 0.5*COS(2.0*EFAC*EBAR) + F(I) = 0.5*EBAR - 0.25/EFAC * SIN(2.0*EFAC*EBAR) + S(I) = EFAC*SIN(2.0*EFAC*EBAR) + IF(ETA(I) .GE. ETJOIN2) GO TO 16 + 15 CONTINUE + 16 CONTINUE + IJOIN = I +C +C----- constant U for outer half + DO 17 I=IJOIN+1, N + U(I) = 1.0 + F(I) = F(IJOIN) + ETA(I) - ETA(IJOIN) + S(I) = 0. + 17 CONTINUE +C + ENDIF +c + 9991 continue +C +C +C---- Newton iteration loop + DO 100 ITER=1, ITMAX +C +C------ zero out A,B,C blocks and righthand sides R + DO 20 I=1, N + DO 201 II=1,3 + DO 2001 III=1,3 + A(II,III,I) = 0. + B(II,III,I) = 0. + C(II,III,I) = 0. + 2001 CONTINUE + R(II,1,I) = 0. + R(II,2,I) = 0. + R(II,3,I) = 0. + 201 CONTINUE + 20 CONTINUE +C +C................................................... +C + A(1,1,1) = 1.0 + A(2,2,1) = 1.0 + A(3,2,N) = 1.0 + R(1,1,1) = F(1) + R(2,1,1) = U(1) + R(3,1,N) = U(N) - 1.0 +C + IF(INORM.EQ.2) THEN + BETU = 2.0*BU/(BU+1.0) + BETU_BU = (2.0 - BETU/(BU+1.0))/(BU+1.0) + BETN = 1.0 + BETN_BU = 0.0 + ELSE + BETU = BU + BETU_BU = 1.0 + BETN = 0.5*(1.0 + BU) + BETN_BU = 0.5 + ENDIF +C + DO 30 I=1,N-1 +C + DETA = ETA(I+1) - ETA(I) + R(1,1,I+1) = F(I+1) - F(I) - 0.5*DETA*(U(I+1)+U(I)) + R(2,1,I+1) = U(I+1) - U(I) - 0.5*DETA*(S(I+1)+S(I)) + R(3,1,I) = S(I+1) - S(I) + & + TN * ( BETN*DETA*0.5*(F(I+1)*S(I+1) + F(I)*S(I)) + & + BETU*DETA*(1.0 - 0.5*(U(I+1)**2 + U(I)**2)) ) +C + A(3,1,I) = TN * BETN*0.5*DETA*S(I) + C(3,1,I) = TN * BETN*0.5*DETA*S(I+1) + A(3,2,I) = -TN * BETU *DETA*U(I) + C(3,2,I) = -TN * BETU *DETA*U(I+1) + A(3,3,I) = TN * BETN*0.5*DETA*F(I) - 1.0 + C(3,3,I) = TN * BETN*0.5*DETA*F(I+1) + 1.0 +C + B(1,1,I+1) = -1.0 + A(1,1,I+1) = 1.0 + B(1,2,I+1) = -0.5*DETA + A(1,2,I+1) = -0.5*DETA +C + B(2,2,I+1) = -1.0 + A(2,2,I+1) = 1.0 + B(2,3,I+1) = -0.5*DETA + A(2,3,I+1) = -0.5*DETA +C + R(3,2,I) = TN + & * ( BETN_BU*DETA*0.5*(F(I+1)*S(I+1) + F(I)*S(I)) + & + BETU_BU*DETA*(1.0 - 0.5*(U(I+1)**2 + U(I)**2))) + R(3,3,I) = ( BETN*DETA*0.5*(F(I+1)*S(I+1) + F(I)*S(I)) + & + BETU*DETA*(1.0 - 0.5*(U(I+1)**2 + U(I)**2)) ) +C + 30 CONTINUE +C........................................................... +C +C +C---- solve Newton system for the three solution vectors + CALL B3SOLV(A,B,C,R,N,NRHS,NRMAX) +C +C +C---- calculate and linearize Dstar, Theta, in computational space + DSI = 0. + DSI1 = 0. + DSI2 = 0. + DSI3 = 0. +C + THI = 0. + THI1 = 0. + THI2 = 0. + THI3 = 0. +C + DO 40 I=1,N-1 + US = U(I) + U(I+1) + DETA = ETA(I+1) - ETA(I) +C + DSI = DSI + (1.0 - 0.5*US)*DETA + DSI_US = -0.5*DETA +C + THI = THI + (1.0 - 0.5*US)*0.5*US*DETA + THI_US = (0.5 - 0.5*US)*DETA +C + DSI1 = DSI1 + DSI_US*(R(2,1,I) + R(2,1,I+1)) + DSI2 = DSI2 + DSI_US*(R(2,2,I) + R(2,2,I+1)) + DSI3 = DSI3 + DSI_US*(R(2,3,I) + R(2,3,I+1)) +C + THI1 = THI1 + THI_US*(R(2,1,I) + R(2,1,I+1)) + THI2 = THI2 + THI_US*(R(2,2,I) + R(2,2,I+1)) + THI3 = THI3 + THI_US*(R(2,3,I) + R(2,3,I+1)) + 40 CONTINUE +C +C + IF(ISPEC.EQ.1) THEN +C +C----- set and linearize Bu = Bspec residual + R1 = BSPEC - BU + Q11 = 1.0 + Q12 = 0.0 +C + ELSE +C +C----- set and linearize H = Hspec residual + R1 = DSI - HSPEC*THI + & -DSI1 + HSPEC*THI1 + Q11 = -DSI2 + HSPEC*THI2 + Q12 = -DSI3 + HSPEC*THI3 +C + ENDIF +C +C + IF(INORM.EQ.3) THEN +C +C----- set and linearize normalized Theta = 1 residual + R2 = THI - 1.0 + & -THI1 + Q21 = -THI2 + Q22 = -THI3 +C + ELSE +C +C----- set eta scaling coefficient to unity + R2 = 1.0 - TN + Q21 = 0.0 + Q22 = 1.0 +C + ENDIF +C +C + DET = Q11*Q22 - Q12*Q21 + DBU = -(R1 *Q22 - Q12*R2 ) / DET + DTN = -(Q11*R2 - R1 *Q21) / DET +C +C +C---- calculate changes in F,U,S, and the max and rms change + RMAX = 0. + RMS = 0. + DO 50 I=1,N + DF = -R(1,1,I) - DBU*R(1,2,I) - DTN*R(1,3,I) + DU = -R(2,1,I) - DBU*R(2,2,I) - DTN*R(2,3,I) + DS = -R(3,1,I) - DBU*R(3,2,I) - DTN*R(3,3,I) +C + RMAX = MAX(RMAX,ABS(DF),ABS(DU),ABS(DS)) + RMS = DF**2 + DU**2 + DS**2 + RMS + 50 CONTINUE + RMS = SQRT(RMS/(3.0*FLOAT(N) + 3.0)) +C + RMAX = MAX(RMAX,ABS(DBU/0.5),ABS(DTN/TN)) +C +C---- set underrelaxation factor if necessary by limiting max change to 0.5 + RLX = 1.0 + IF(RMAX.GT.0.5) RLX = 0.5/RMAX +C +C---- update F,U,S + DO 60 I=1,N + DF = -R(1,1,I) - DBU*R(1,2,I) - DTN*R(1,3,I) + DU = -R(2,1,I) - DBU*R(2,2,I) - DTN*R(2,3,I) + DS = -R(3,1,I) - DBU*R(3,2,I) - DTN*R(3,3,I) +C + F(I) = F(I) + RLX*DF + U(I) = U(I) + RLX*DU + S(I) = S(I) + RLX*DS + 60 CONTINUE +C +C---- update BetaU and Theta + BU = BU + RLX*DBU + TN = TN + RLX*DTN +C +C---- check for convergence + IF(ITER.GT.3 .AND. RMS.LT.1.E-5) GO TO 105 +C + 100 CONTINUE + WRITE(*,*) 'FS: Convergence failed' +C + 105 CONTINUE +C + HSPEC = DSI/THI + BSPEC = BU +C + DELTA = SQRT(TN) +C + RETURN +C +C The + END + + + SUBROUTINE B3SOLV(A,B,C,R,N,NRHS,NRMAX) + IMPLICIT REAL (A-H,M,O-Z) + DIMENSION A(3,3,N), B(3,3,N), C(3,3,N), R(3,NRMAX,N) +C ********************************************************************** +C This routine solves a 3x3 block-tridiagonal system with an arbitrary +C number of righthand sides by a standard block elimination scheme. +C The solutions are returned in the Rj vectors. +C +C |A C ||d| |R..| +C |B A C ||d| |R..| +C | B . . ||.| = |R..| +C | . . C||.| |R..| +C | B A||d| |R..| +C Mark Drela 10 March 86 +C ********************************************************************** +C +CCC** Forward sweep: Elimination of lower block diagonal (B's). + DO 1 I=1, N +C + IM = I-1 +C +C------ don't eliminate B1 block because it doesn't exist + IF(I.EQ.1) GO TO 12 +C +C------ eliminate Bi block, thus modifying Ai and Ci blocks + DO 11 K=1, 3 + DO 111 L=1, 3 + A(K,L,I) = A(K,L,I) + & - (B(K,1,I)*C(1,L,IM) + B(K,2,I)*C(2,L,IM) + B(K,3,I)*C(3,L,IM)) + 111 CONTINUE + DO 112 L=1, NRHS + R(K,L,I) = R(K,L,I) + & - (B(K,1,I)*R(1,L,IM) + B(K,2,I)*R(2,L,IM) + B(K,3,I)*R(3,L,IM)) + 112 CONTINUE + 11 CONTINUE +C +C -1 +CCC---- multiply Ci block and righthand side Ri vectors by (Ai) +C using Gaussian elimination. +C + 12 DO 13 KPIV=1, 2 + KP1 = KPIV+1 +C +C-------- find max pivot index KX + KX = KPIV + DO 131 K=KP1, 3 + IF(ABS(A(K,KPIV,I)) .LE. ABS(A(KX,KPIV,I))) THEN + GO TO 131 + ELSE + GO TO 1311 + ENDIF + 1311 KX = K + 131 CONTINUE +C + IF(A(KX,KPIV,I).EQ.0.0) THEN + WRITE(*,*) 'Singular A block, i = ',I + STOP + ENDIF +C + PIVOT = 1.0/A(KX,KPIV,I) +C +C-------- switch pivots + A(KX,KPIV,I) = A(KPIV,KPIV,I) +C +C-------- switch rows & normalize pivot row + DO 132 L=KP1, 3 + TEMP = A(KX,L,I)*PIVOT + A(KX,L,I) = A(KPIV,L,I) + A(KPIV,L,I) = TEMP + 132 CONTINUE +C + DO 133 L=1, 3 + TEMP = C(KX,L,I)*PIVOT + C(KX,L,I) = C(KPIV,L,I) + C(KPIV,L,I) = TEMP + 133 CONTINUE +C + DO 134 L=1, NRHS + TEMP = R(KX,L,I)*PIVOT + R(KX,L,I) = R(KPIV,L,I) + R(KPIV,L,I) = TEMP + 134 CONTINUE +CB +C-------- forward eliminate everything + DO 135 K=KP1, 3 + DO 1351 L=KP1, 3 + A(K,L,I) = A(K,L,I) - A(K,KPIV,I)*A(KPIV,L,I) + 1351 CONTINUE + C(K,1,I) = C(K,1,I) - A(K,KPIV,I)*C(KPIV,1,I) + C(K,2,I) = C(K,2,I) - A(K,KPIV,I)*C(KPIV,2,I) + C(K,3,I) = C(K,3,I) - A(K,KPIV,I)*C(KPIV,3,I) + DO 1352 L=1, NRHS + R(K,L,I) = R(K,L,I) - A(K,KPIV,I)*R(KPIV,L,I) + 1352 CONTINUE + 135 CONTINUE +C + 13 CONTINUE +C +C------ solve for last row + IF(A(3,3,I).EQ.0.0) THEN + WRITE(*,*) 'Singular A block, i = ',I + STOP + ENDIF + PIVOT = 1.0/A(3,3,I) + C(3,1,I) = C(3,1,I)*PIVOT + C(3,2,I) = C(3,2,I)*PIVOT + C(3,3,I) = C(3,3,I)*PIVOT + DO 14 L=1, NRHS + R(3,L,I) = R(3,L,I)*PIVOT + 14 CONTINUE +C +C------ back substitute everything + DO 15 KPIV=2, 1, -1 + KP1 = KPIV+1 + DO 151 K=KP1, 3 + C(KPIV,1,I) = C(KPIV,1,I) - A(KPIV,K,I)*C(K,1,I) + C(KPIV,2,I) = C(KPIV,2,I) - A(KPIV,K,I)*C(K,2,I) + C(KPIV,3,I) = C(KPIV,3,I) - A(KPIV,K,I)*C(K,3,I) + DO 1511 L=1, NRHS + R(KPIV,L,I) = R(KPIV,L,I) - A(KPIV,K,I)*R(K,L,I) + 1511 CONTINUE + 151 CONTINUE + 15 CONTINUE + 1 CONTINUE +C +CCC** Backward sweep: Back substitution using upper block diagonal (Ci's). + DO 2 I=N-1, 1, -1 + IP = I+1 + DO 21 L=1, NRHS + DO 211 K=1, 3 + R(K,L,I) = R(K,L,I) + & - (R(1,L,IP)*C(K,1,I) + R(2,L,IP)*C(K,2,I) + R(3,L,IP)*C(K,3,I)) + 211 CONTINUE + 21 CONTINUE + 2 CONTINUE +C + RETURN + END ! B3SOLV -- cgit v1.2.3