c
REAL FUNCTION KINVIST(I)
INCLUDE 'farray'
INCLUDE 'mfmdat'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'satear'
INCLUDE 'grdear'
INCLUDE 'prpcmn'
COMMON/GENI/IGF1(2),NXNYST,NDIR,KDUMM,IGF2(4),NFM,IGF3(39),
1 ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IGF4(4)
1 /F0/IF01(29),L0XYDX,L0XYDY,IF02(4),L0XYXG,IF03,L0XYYG,IF04,
1 L0XYDXG,L0XYDYG,IF05(124),L0XYDZ,L0XYDZG,IF06(137)
1 /LRNTM1/L0WDIS,L0FMU,L0FONE,L0FTWO,L0REYN,L0REYT,L0UD1,
1 L0UD2,L0UD3,L0UD4
1 /LRNTM2/LBWDIS,LBFMU,LBFONE,LBFTWO,LBREYN,LBREYT,LBEPKE
1 /LRNTM3/L0UTAU,NUMWAL,L0WALL,L0DSKN,IVPRST
1 /MMKKE/ LBFOMG,L0FOMG
1 /DOMSIZ/SIZMAX
1 /CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
1 /MCRMXI/LBRATE,L0APXY,L0SUXY,L0APCL,L0SUCL,L0RATE
C COMMON/REALKE1/LBCMU,LBFIL1(10)
COMMON/REALKE2/L0CMU,LBFIL2(10)
COMMON/KWMOD2/L0KWF1(7),L0BF1,L0BF2,L0BF3,L0KWF2(7)
C
C List of turbulence models (according to IENUTA)
C -----------------------------------------------
c< laminar flow
c = 0 !< standard k-e model
c = 1 !< RNG high-Re k-e model
c = 2 !< Chen-Kim high-Re k-e model
c = 3 !< Lam-Bremhorst low-Re k-e model
c = 4 !< Chen-Kim low-Re k-e model
c = 6 !< Simple eddy-viscosity multiplier for low Reynolds numbers
c = 7 !< Two-scale high-Re k-e model
c = 8 !< Two-layer low-Re k-e model
c = 9 !< Another simple eddy-viscosity multiplier for low Reynolds numbers
c = 10 !< Wilcox (1988) high-Re k-omega model
c = 11 !< Wilcox (1988) low-Re k-omega model
c = 12 !< MMK high-Re k-e model
c = 13 !< Kato-Launder high-Re k-e model
c = 14 !< Realisable high-Re k-e model
c = 15 !< Wilcox (2008) high-Re k-omega model
c = 16 !< Wilcox (2008) low-Re k-omega model
c = 17 !< Menter (1992) high-Re k-omega model
c = 18 !< Menter (1992) low-Re k-omega model
C = 19 !< high-Re k-omega SST model
C = 20 !< low-Re k-omega SST model
c other turbulence models where IENUTA is not used:
c !< LVEL model
c !< model of constant effective viscosity
C
C.... KINVIST is turbulent contribution to total kinematic viscosity
IF(IGRND==-1) THEN ! use the PIL variable ENUT
KINVIST=PRPRTY
ELSEIF(IGRND==1) THEN ! linear function of length
KINVIST= ENUTA + ENUTB*F(L0LEN+I)
ELSEIF(IGRND==2) THEN ! Prandtl mixing-length formula
KINVIST= F(L0LEN+I)**2 * SQRT(ABS(F(L0GEN+I)))
ELSEIF(IGRND==3) THEN ! Prandtl-Kolmogorov formula
KINVIST= CMU*F(L0LEN+I)*SQRT(ABS(F(L0KE+I)))
ELSEIF(IGRND==4) THEN ! two-fluid model of turbulence
ABSDV = VSLPCL(I,XCYCZ,SLUNX1)
KINVIST= ENUTA*F(L0LEN+I)*ABSDV*F(L0R1+I)*F(L0R2+I)
ELSEIF(IGRND==5) THEN ! Harlow-Nakayama formula
KINVIST= CMUCD*F(L0KE+I)**2/AMAX1(F(L0EP+I),TINY)
ELSEIF(IGRND==6) THEN ! Saffman-Spalding formula
KINVIST= CMUCD*F(L0KE+I)/SQRT(AMAX1(F(L0VOSQ+I),TINY))
ELSEIF(IGRND==7) THEN ! Kolmogorov formula
GOMEGA = AMAX1(F(L0OMEG+I),TINY)
GKE = F(L0KE+I)
IF(IENUTA==15.OR.IENUTA==16) THEN ! Wilcox (2008) k-w model
CLIM = 7./8.
GGEN1 = SQRT(F(L0GEN+I)) ! S=Sqrt(2.*SijSij) = sqrt(GEN1)
GOMLIM = AMAX1(GOMEGA,CLIM*GGEN1/0.3)
KINVIST = GKE/GOMLIM
ELSE IF(IENUTA==19.OR.IENUTA==20) THEN ! k-w SST model
GGEN1 = SQRT(F(L0GEN+I)) ! S=Sqrt(2.*SijSij) = sqrt(GEN1)
GBF2 = F(L0BF2+I)
C GBF3 = F(L0BF3+I)
C GSF2 = GGEN1*GBF2*GBF3
GSF2 = GGEN1*GBF2
GOMLIM = AMAX1(0.3*GOMEGA,GSF2)
KINVIST = 0.3*GKE/GOMLIM
ELSE ! other k-w models
KINVIST= GKE/GOMEGA
ENDIF
ELSEIF(IGRND==8) THEN ! LVEL model (SKNFG2 returns the
VISLAM= F(L0VISL+I) ! skin-friction factor at the wall)
REYN = SQRT(ABSVSQ(I,XCYCZ,SLUNX1))*F(L0LEN+I)/(VISLAM+TINY)
IF(LBREYN>0) THEN
F(L0REYN+I)=REYN
ENDIF
IF(REYN<1.0) THEN ! very close to the wall there
KINVIST= 0.0 ! is no turbulent contribution
ELSE ! formula from Spalding's wall law
c TAUP= SKNFG2(REYN,5+ISKINA,ISKINB,.FALSE.,REYN,EEFF,SKINT)
c AKUP= AK/SQRT(TAUP) ! K times uplus
c SUM = ( (0.1666667*AKUP + 0.5) * AKUP + 1.0) * AKUP + 1.0
c KINVIST = VISLAM * ( EXP(AKUP) - SUM ) * AK / EWAL
c dbs 17.06.12 The above recourse to sknfg2 is not necessary; and indeed
c it is incorrect; for there is no clearly-connected
c near-surface cell where the skin friction can be calculated
c for circumstances sketched here.
c
c x <-- point for which effective viscosity
c is to be calculated
c ////////////////
c ///// solid ////
c ////////////////
c
c The following call to enutpl suffice:
c
KINVIST=VISLAM*ENUTPL(REYN)
ENDIF
ELSEIF(IGRND==9) THEN
KINVIST= ENUTA ! uniform value
ELSEIF(IGRND==10) THEN
KINVIST= VISCON*F(L0MIXL+I)**2*F(L0RATE+I) ! Formula used in some MFM
ENDIF ! cases
C.... Modifications for low Reynolds numbers:
IF(IENUTA==3 .OR.IENUTA==4) THEN ! Lam-Bremhorst K-E model
FMU = 1.0
REYT= AMAX1( F(L0KE+I)**2/(F(L0EP+I)*F(L0VISL+I)+TINY) , 1.E-9)
REYN= SQRT(ABS(F(L0KE+I)))*F(L0WDIS+I)/(F(L0VISL+I)+TINY)
IF(ABS(REYN)<800.) FMU= 1.0 - EXP(-0.0165*REYN)
FMU= AMAX1(AMIN1(FMU**2*(1.+20.5/REYT), 1.0), 1.E-9)
KINVIST= KINVIST*FMU
F(L0REYT+I)= REYT
F(L0REYN+I)= REYN
F(L0FMU +I)= FMU
ELSEIF(IENUTA==6) THEN ! a Simple low-Re modification of ENUT
C VisT = VisT_nom * amin1(1.0 , (A * VisT_nom / VisL ) ** B)
C where: VisT = turbulent kinematic viscosity to be used,
C VisT_nom = nominal (high-Reynolds-number) value of VisT
C VisL = laminar viscosity
C Suitable values for A and B may be 0.2 and 4.0
VISLAM= F(L0VISL+I)
FACTOR= ENUTA*KINVIST
IF(FACTOR0.0) THEN
VSLP = AMAX1(CFIPA, VSLPCL(I,XCYCZ,SLUNX1))
VISADD = ENUTC * F(L0R2+I) * (1.0 - F(L0R2+I)) * VSLP
IF(L0LEN/=0) THEN
VISADD = VISADD * F(L0LEN+I)
ENDIF
KINVIST = KINVIST + VISADD
ENDIF
ENDIF
END
C***********************************************************************
SUBROUTINE SLBVST(IPILOPT,dbgloc)
INCLUDE 'farray'
INCLUDE 'mfmdat'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'satear'
INCLUDE 'grdear'
INCLUDE 'prpcmn'
COMMON/VMSCMN/FL1CON
1 /LRNTM1/L0WDIS,L0FMU,L0FONE,L0FTWO,L0REYN,L0REYT,L0UD1,
1 L0UD2,L0UD3,L0UD4
1 /LRNTM2/LBWDIS,LBFMU,LBFONE,LBFTWO,LBREYN,LBREYT,LBEPKE
1 /MMKKE/ LBFOMG,L0FOMG
1 /VELCMN/L0UVW(6),L0UVW2(6) /FLPCMN/IFILP(30)
1 /CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
COMMON/TSKEMI/ LBKP,LBKT,LBET,LBVOSQ,LBOMEG
1 /TSKEMR/ GKTDKP
COMMON/REALKE1/LBCMU,LBFIL1(10)
COMMON/REALKE2/L0CMU,LBFIL2(10)
COMMON/KWMOD1/LBKWF1(7),LBBF1,LBBF2,LBBF3,LBSIGK,LBSIGW,LBKWF2(5)
COMMON/KWMOD2/LBKWF3(7),L0BF1,L0BF2,L0BF3,L0SIGK,L0SIGW,LBKWF4(5)
COMMON/GENI/IGF1(2),NXNYST,NDIR,KDUMM,IGF2(4),NFM,NWHOLE, IGSH,
1 IGF3(19),IPRL,
1 IBTAU,IGF4(16),ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,
1 IPRPS,IRADX,IRADY,IRADZ,IVFOL
1 /F0/ IF01(29),L0XYDX,L0XYDY,IF02(3),L0XYRV,L0XYXG,IF03,
1 L0XYYG,IF04,L0XYDXG,L0XYDYG,IF05(68),KZXCY,IF06(37),
1 L0AHZ,IF07(17),L0XYDZ,L0XYDZG,IF08(137)
COMMON /MCRMXI/LBRATE,L0APXY,L0SUXY,L0APCL,L0SUCL,L0RATE,L0FRCS,
1 L0ABSC,L0AV,L0BV,L0CV
COMMON/NAMFN/NAMFUN,NAMSUB
REAL KINVIST
LOGICAL DBGLOC,SLD,FL1CON,GRN,NEZ
CHARACTER*6 NAMFUN,NAMSUB
C
NAMSUB= 'SLBVST'
if(flag.or.dbgloc) call banner(1,'SLBVST',070519 )
IGR= 9
ISC= IPROP
C.... Call GROUND for user-set property:
IF(IPILOPT==0) THEN
IF(USEGRD) THEN
CALL GROUND
ENDIF
GO TO 800
ENDIF
IF(IPILOPT==-1) GO TO 700
C.... Set constants and other auxiliary variables:
C-----------------------------------------------------------------------
IF(IGRND==4.OR.IGRND==8.OR.(.NOT.ONEPHS.AND.ENUTC>0.0))THEN
IF(XCYCLE) XCYCZ=NEZ(F(KZXCY+IZ))
IF(STORE(3)) L0UVW(3)=L0F(3)
IF(STORE(5)) L0UVW(1)=L0F(5)
IF(STORE(7)) L0UVW(5)=L0F(7)
IF(.NOT.ONEPHS) THEN
IF(STORE(4)) L0UVW2(3)=L0F(4)
IF(STORE(6)) L0UVW2(1)=L0F(6)
IF(STORE(8)) L0UVW2(5)=L0F(8)
ENDIF
ENDIF
L0VISL = L0F(VISL)
IF(IGRND==3) THEN
IF(STORE(12)) L0KE=L0F(12)
IF(ENUTC>0.0.AND.STORE(10)) L0R2=L0F(10)
ELSEIF(IGRND==4) THEN
IF(STORE(9)) L0R1=L0F(9)
IF(SOLVE(10)) L0R2=L0F(10)
ELSEIF(IGRND==5) THEN
IF(STORE(12)) L0KE=L0F(12)
IF(STORE(13)) L0EP=L0F(13)
ELSEIF(IGRND==6)THEN
IF(STORE(12)) L0KE=L0F(12)
IF(LBVOSQ>0) L0VOSQ=L0F(LBVOSQ)
ELSEIF(IGRND==7) THEN
IF(STORE(12)) L0KE=L0F(12)
IF(LBOMEG>0) L0OMEG=L0F(LBOMEG)
ELSEIF(IGRND==8) THEN
L0VISL=L0F(VISL)
IF(XCYCLE) XCYCZ=NEZ(F(KZXCY+IZ))
IF(LBREYN/=0) L0REYN= L0F(LBREYN)
ELSEIF(IGRND==10) THEN
IF(LBRATE>0) L0RATE=L0F(LBRATE)
IF(LBMIXL>0) L0MIXL=L0F(LBMIXL)
ENDIF
IF(USEWDS) L0WDIS=L0F(LBWDIS)
IF(GENK) L0GEN=L0F(GEN1)
IF(LEN1>0 .OR. GRN(EL1)) L0LEN=L0F(LEN1)
IF(IENUTA==3 .OR.IENUTA==4) THEN
IF(LBREYT/=0) L0REYT= L0F(LBREYT)
IF(LBREYN/=0) L0REYN= L0F(LBREYN)
IF(LBFMU/=0) L0FMU = L0F(LBFMU)
ELSEIF(IENUTA==8) THEN
IF(LBREYN/=0) L0REYN= L0F(LBREYN)
IF(LBFMU/=0) L0FMU = L0F(LBFMU)
ELSEIF(IENUTA==11) THEN
IF(LBREYT/=0) L0REYT= L0F(LBREYT)
IF(LBFMU/=0) L0FMU = L0F(LBFMU)
ELSEIF(IENUTA==12) THEN
IF(LBFOMG/=0) L0FOMG = L0F(LBFOMG)
ELSEIF(IENUTA==14) THEN
IF(LBCMU/=0) L0CMU = L0F(LBCMU)
ELSEIF(IENUTA==19.OR.IENUTA==20) THEN
IF(LBBF2/=0) L0BF2 = L0F(LBBF2)
C IF(LBBF3/=0) L0BF3 = L0F(LBBF3)
ENDIF
C-----------------------------------------------------------------------
C.... Loop over slab to get and set cell properties:
700 IGRND=IPILOPT
IF(IPRPS==0) THEN
C.... One material only
DO 60 I= 1,NXNYST
60 F(KPROP+I)= KINVIST(I)
ELSE
C.... exclude solids
DO 70 I= 1,NXNYST
IF(SLD(I)) THEN
F(KPROP+I)= TINY
ELSE
GVIST = KINVIST(I)
F(KPROP+I)= KINVIST(I)
ENDIF
70 CONTINUE
ENDIF
C----------------------------------------------------------------------
C.... Corrections, debug print-out, and other property adjustments:
IF(IPILOPT==8 .AND. LBEPKE/=0) THEN
L0EPKE= L0F(LBEPKE)
DO 80 I= 1,NXNYST
IF(SLD(I)) THEN
F(L0EPKE+I)= 0.0
ELSE
RWDIS= 1./(F(L0WDIS+I)+TINY)
F(L0EPKE+I)= F(KPROP+I)*RWDIS*RWDIS
ENDIF
80 CONTINUE
ENDIF
C.... Call GREX to correct a property set above:
800 IF(USEGRX) CALL GREX3
C.... Call ALTPRP for an alternative property setting
IF(USEALT) CALL ALTPRP
C.... Call GROUND for the user to correct a property set above:
IF(USEGRD) THEN
IF(IPILOPT>0) CALL GROUND
ENDIF
if(flag.or.dbgloc) call banner(2,'SLBVST',070519 )
NAMSUB= 'SLBVST'
END
c-------------------------
SUBROUTINE INIVST
INCLUDE 'farray'
INCLUDE 'qdata'
INCLUDE 'satear'
INCLUDE 'satgrd'
INCLUDE 'prpcmn'
COMMON/TURB/CMU,CD,CMUCD,C1E,C2E,AK,EWAL,YPLUSC,REWC
1 /SPEPRC/MFM
1 /LRNTM1/L0WDIS,L0FMU,L0FONE,L0FTWO,L0REYN,L0REYT,L0UD1,
1 L0UD2,L0UD3,L0UD4/GENI/IGN1(49),ITEM1,ITEM2,IGN2(9)
1 /LRNTM2/LBWDIS,LBFMU,LBFONE,LBFTWO,LBREYN,LBREYT,LBEPKE
1 /LRNTM3/L0UTAU,NUMWAL,L0WALL,L0DSKN,IVPRST
1 /LRNTM4/RSCMCD,CDDAK2,TAUDKE,RTTDKE,AKC,EWC
COMMON/TSKEMI/ LBKP,LBKT,LBET,LBVOSQ,LBOMEG
1 /TSKEMR/ GKTDKP
1 /VELCMN/L0UVW(6),L0UVW2(6) /FLPCMN/IFILP(30)
COMMON/KWMOD3/CWWALL,SIGK1,SIGK2,SIGW1,SIGW2,CWA1,CWA2,CWB1,CWB2,
1 BETAST,CDSIG
LOGICAL GETREALS,KEMOD,KFMOD,SOLVE,KFHIRE
DIMENSION NAMPRP(30)
EQUIVALENCE (NAMPRP(1),DENST1)
C.... Preliminaries:
if(flag) call banner(1,'INIVST',120717 )
C.... Turbulence-model constants:
IF(.NOT.RSTM) THEN
IF(CMU<=0.0) CMU=0.5478
IF(CD<=0.0) CD= 0.1643
IF(C1E<=0.0) C1E=1.44
IF(C2E<=0.0) C2E=1.92
IF(AK<=0.0) AK= 0.41
IF(EWAL<=0.0) EWAL=8.6
IF(YPLUSC<=0.0) YPLUSC=11.126
KEMOD=SOLVE(12).AND.SOLVE(13)
KFMOD=SOLVE(12).AND.SOLVE(LBOMEG)
c
IF(GETREALS('TURCONST')) THEN
CMU=REALS(1) ; CD=REALS(2) ; C1E=REALS(3)
C2E=REALS(4) ; AK=REALS(5) ; EWAL=REALS(6)
ELSE
CALL GETSPD('KECONST','CMU',1,CMU,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','CD',1,CD,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','AK' ,1,AK,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','EWAL' ,1,EWAL,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','YPLUSC',1,YPLUSC,IDUM,.FALSE.,' ',IERR)
ENDIF
CMUCD=CMU*CD
IF(SCALWF.AND.ENUT/=GRND9) THEN
IF(.NOT.KEMOD) THEN ! only for KE models
SCALWF=.FALSE.
KFHIRE=IENUTA==10.OR.IENUTA==15.OR.IENUTA==17.OR.IENUTA==19
IF(KFMOD.AND.KFHIRE) SCALWF=.TRUE.! only high-RE
IF(ISKINB>0) THEN
SCALWF=.TRUE. ! allow for LVEL
CALL WRIT40('Scalable Wall Function allowed for LVEL')
ENDIF
ELSE
IF((IENUTA>2.AND.IENUTA<7).OR.(IENUTA>7.AND.IENUTA<10).OR.
1 IENUTA==11) THEN ! only high-RE
SCALWF=.FALSE.
ENDIF
ENDIF
IF(.NOT.SCALWF) THEN
CALL WRIT40('Scalable wall function switched off because')
CALL WRIT40('not compatible with current turbulence model')
ENDIF
ENDIF
REWC=YPLUSC**2. ! critical wall Reynolds number
RSCMCD=CMUCD**(-0.5) ; CDDAK2=2.0*CD/AK ; TAUDKE=SQRT(CMUCD)
RTTDKE=SQRT(TAUDKE) ; AKC=AK*RTTDKE ; EWC=EWAL*RTTDKE
C
C.... Standard K-E model:
IF(IENUTA==0) THEN
IF(.NOT.NULLPR) CALL PR_TURCON(IENUTA,
1 ' Standard k-e model constants ')
C.... Lam-Bremhorst low-Re K-E model:
ELSEIF(IENUTA==3) THEN
IF(.NOT.NULLPR) CALL PR_TURCON(IENUTA,
1 ' Lam-Bremhorst low-Re model constants ')
C.... Two-layer low-Re K-E model:
ELSEIF(IENUTA==8) THEN
IF(.NOT.NULLPR) CALL PR_TURCON(IENUTA,
1 ' Two-layer low-Re k-e model constants ')
ELSEIF(IENUTA==1) THEN
C.... RNG-derived K-E model:
CMUCD=0.0845 ; C1E=1.42 ; C2E=1.68
CD = CMUCD**0.75 ; CMU= CMUCD/CD
CALL GETSPD('KECONST','CMU',1,CMU,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','CD',1,CD,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,' ',IERR)
CMUCD=CMU*CD
CALL GXRNGM(0,0,FIXFLU)
IF(.NOT.NULLPR) CALL PR_TURCON(IENUTA,
1 ' RNG k-e turbulence model constants ')
ELSEIF(IENUTA==2 .OR. IENUTA==4) THEN
C.... Chen-Kim high & low-Re k-e models
C1E=1.15 ; C2E=1.9
CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,' ',IERR)
IF(.NOT.NULLPR) THEN
IF(IENUTA==2) CALL PR_TURCON(IENUTA,
1 ' Chen-Kim k-e turbulence model constants')
IF(IENUTA==4) CALL PR_TURCON(IENUTA,
1 ' Chen-Kim low-Re k-e model constants ')
ENDIF
ELSEIF(IENUTA==7) THEN
C.... Two-scale k-e model:
C1E=1.24 ; C2E=1.84 ; GC3E=0.21
CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','GC3E',1,GC3E,IDUM,.FALSE.,' ',IERR)
GKTDKP= AK*AK/(PRT(13)*SQRT(CMUCD)*(C2E-GC3E-C1E)) - 1.0
CALL GXTSKE(0,0,FIXFLU)
IF(.NOT.NULLPR) CALL PR_TURCON(IENUTA,
1 'Two-scale k-e model constants ')
ELSEIF(IENUTA==10 .OR. IENUTA==11) THEN
C.... Wilcox (1988) high & low-Re k-w models
C1E=5./9. ; C2E=0.075
CWWALL=6.0 ! low-Re wall constant for w
CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','CWWALL',1,CWWALL,IDUM,.FALSE.,
1 ' ',IERR)
KELIN= 0
IF(.NOT.NULLPR) THEN
IF(IENUTA==10) CALL PR_TURCON(IENUTA,
1 ' Wilcox (1988) k-w model constants ')
IF(IENUTA==11) CALL PR_TURCON(IENUTA,
1 'Wilcox (1988) low-Re k-w model constants')
ENDIF
ELSEIF(IENUTA==14) THEN
C.... Realisable high-Re K-E model:
C2E=1.9
CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,' ',IERR)
C... KELIN=3 must be used because other KELIN options use a constant
C CMU in the linearisation of the KE source/sink terms
KELIN=3
IF(.NOT.NULLPR) CALL PR_TURCON(IENUTA,
1 ' Realisable k-e model constants ')
ELSEIF(IENUTA==15 .OR. IENUTA==16) THEN
C.... Wilcox (2008) high & low-Re revised k-w models
CDSIG=1./8. ! cross-diffusion source-term coefficient
C1E=13./25. ; C2E=0.0708
CWWALL=6.0 ! low-Re wall constant for w
CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','CDSIG',1,CDSIG,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','CWWALL',1,CWWALL,IDUM,.FALSE.,
1 ' ',IERR)
KELIN= 0
IF(.NOT.NULLPR) THEN
IF(IENUTA==15) CALL PR_TURCON(IENUTA,
1 ' Wilcox (2008) k-w model constants ')
IF(IENUTA==16) CALL PR_TURCON(IENUTA,
1 'Wilcox (2008) low-Re k-w model constants')
ENDIF
ELSEIF(IENUTA==17.OR.IENUTA==18) THEN
C.... Menter (1992) high & low-Re Baseline k-w models
SIGK1 = 2.0 ; SIGW1 = 2.0 ! wall region limit
SIGK2 = 1.0 ; SIGW2 = 1./0.856 ! outer region limit
C1E=5./9. ; C2E=3./40.
CWWALL=6.0 ! low-Re wall constant for w
CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','CWWALL',1,CWWALL,IDUM,.FALSE.,
1 ' ',IERR)
BETAST=CMUCD
CWB1=C2E ! = 0.075
CWB2 = 0.0828 ; AK2 = AK*AK
CWA1 = CWB1/BETAST - AK2/(SIGW1*SQRT(BETAST)) ! CWA1 = 0.5532=C1E
CWA2 = CWB2/BETAST - AK2/(SIGW2*SQRT(BETAST)) ! CWA2 = 0.44
KELIN= 0
IF(.NOT.NULLPR) THEN
IF(IENUTA==17) CALL PR_TURCON(IENUTA,
1 ' Menter Baseline k-w model constants ')
IF(IENUTA==18) CALL PR_TURCON(IENUTA,
1 'Menter Baseline low-Re k-w model ')
ENDIF
ELSEIF(IENUTA==19.OR.IENUTA==20) THEN
C.... Menter (1992) high & low-Re k-w SST models
C Original Menter paper sets SIGK1=0.85 & SIGW1 = 0.65 ! check
SIGK1 = 2.0 ; SIGW1 = 2.0
SIGK2 = 1.0 ; SIGW2 = 1./0.856 ! outer region limit
C1E=5./9. ; C2E=3./40.
CWWALL=6.0 ! low-Re wall constant for w
CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','CWWALL',1,CWWALL,IDUM,.FALSE.,
1 ' ',IERR)
BETAST=CMUCD ! CMUCD =0.09
CWB1=C2E
CWB2 = 0.0828 ; AK2 = AK*AK
CWA1 = CWB1/BETAST - AK2/(SIGW1*SQRT(BETAST)) ! CWA1 = 0.5532
CWA2 = CWB2/BETAST - AK2/(SIGW2*SQRT(BETAST)) ! CWA2 = 0.44
KELIN= 0
IF(.NOT.NULLPR) THEN
IF(IENUTA==19) CALL PR_TURCON(IENUTA,
1 ' k-w SST turbulence model constants ')
IF(IENUTA==20) CALL PR_TURCON(IENUTA,
1 ' k-w SST low-Re model constants ')
ENDIF
ENDIF
C IF(.NOT.NULLPR) THEN
ELSE
IF(SCALWF) THEN
SCALWF=.FALSE.
CALL WRIT40('Scalable wall function switched off because')
CALL WRIT40('not compatible with RSTM turbulence model')
ENDIF
ENDIF
RSCMCD=CMUCD**(-0.5) ; CDDAK2=2.0*CD/AK ; TAUDKE=SQRT(CMUCD)
RTTDKE=SQRT(TAUDKE) ; AKC=AK*RTTDKE ; EWC=EWAL*RTTDKE
C.... Memory allocation for use in turbulence modelling:
IF(ENUT<0.0) THEN
CALL ADDL1C(L0UTAU,IENUTA==5,'l0utau')
IF(IENUTA==3 .OR. IENUTA==4 .OR. IENUTA==8 .OR.
1 IENUTA==11) THEN
NXNY=NX*NY
IF(LBFONE==0) CALL GXMAKE(L0FONE,NXNY,'FONE')
IF(LBFTWO==0) CALL GXMAKE(L0FTWO,NXNY,'FTWO')
IF(LBFMU==0) CALL GXMAKE(L0FMU,NXNY,'FMU ')
IF(LBREYN==0.AND.IENUTA/=11) CALL GXMAKE(L0REYN,NXNY,'RE')
IF(LBREYT==0.AND.IENUTA/=8) CALL GXMAKE(L0REYT,NXNY,'RET')
ENDIF
ENDIF
if(flag) call banner(2,'INIVST',070519 )
END
************************************************************************
REAL FUNCTION ENUTPL(REYN)
COMMON/TURB/CMU,CD,CMUCD,C1E,C2E,AK,EWAL,YPLUSC,REWC
COMMON/PLUSES/UPLUS,YPLUS ! this common appears also in sknfg2
c
RECEW=1.0/EWAL
REYNHI=REYN*1.001 ! To stop search when within 0.1% of
REYNLO=REYN*0.999 ! of specified Reynolds number
IF(REYN<=1.0) THEN
ENUTPL=0.0
YPLUS=SQRT(REYN)
UPLUS=YPLUS
RETURN
ELSEIF(REYN<=1.E1) THEN ! this block of conditions is designed to
UPSMALL=1.0 ! reduce the number of iterations.
UPBIG=3.2 ! The values set for upbig and upsmall
ELSEIF(REYN<=1.E2) THEN ! are valid only for ak=0.41. ewal=8.6
UPSMALL=3.1
UPBIG=10.0
ELSEIF(REYN<=1.E3) THEN
UPSMALL=9.0
UPBIG=16.0
ELSEIF(REYN<=1.E4) THEN
UPSMALL=15.0
UPBIG=21.0
ELSEIF(REYN<=1.E5) THEN
UPSMALL=20.0
UPBIG=26.0
ELSEIF(REYN<=1.E6) THEN
UPSMALL=25.
UPBIG=31.
ELSEIF(REYN<=1.E7) THEN
UPSMALL=30.0
UPBIG=36.0
ELSEIF(REYN<=1.E8) THEN
UPSMALL=35.0
UPBIG=42.0
ELSE
UPSMALL=41.0
UPBIG=100
ENDIF
DO I=1,100
NUMIT=I
UP=0.5*(UPBIG+UPSMALL)
AKUP=AK*UP
EXPKUP=EXP(AKUP)
YP=UP+RECEW*(EXPKUP -( 1 + AKUP*(1 + AKUP*(0.5 +AKUP*(0.1666667
1 + 0.0444444*AKUP)))))
RE=UP*YP
IF(REREYNLO) THEN ! but if the difference is small, don't bother
GOTO 10
ENDIF
UPSMALL=UP
ELSE ! up will need to decrease
IF(RE