From 3c65bf18c059760afa967cfe059ab5f4b5101876 Mon Sep 17 00:00:00 2001 From: sdhzhs Date: Sat, 31 Aug 2024 09:37:30 +0800 Subject: [PATCH] Add comment codes which implements low Re version of SA model --- lib/Condiff.f03 | 45 +++++++++++++++++++++++++-------------------- lib/Turvis.f03 | 9 ++++++++- lib/Wallfunc.f03 | 2 +- 3 files changed, 34 insertions(+), 22 deletions(-) diff --git a/lib/Condiff.f03 b/lib/Condiff.f03 index 4ec0ccb..4ed7c0c 100644 --- a/lib/Condiff.f03 +++ b/lib/Condiff.f03 @@ -3,7 +3,7 @@ Subroutine Condiff(scalar) implicit none integer i,j real(8) Dnow,Dnoe,Dnos,Dnon,Faw,Fae,Fks,Fkn,Fwallw,Fwalle -real(8) Xi,fnu1,fnu2,Sv,rm,gm,Ret,Dwplus,phi1,F1,betaistar,Fmt,alphaf,betai,Tmin,Dampk,Ymax +real(8) Xi,fnu1,fnu2,Sv,rm,gm,Ret,Dwplus,phi1,F1,betaistar,Fmt,alphaf,betai,Tmin,Dampk,Ymax,Ym real(8) aP,aW,aE,aS,aN,DF real(8),external:: interpl real(8) Fwall(Ib1:Ib2) @@ -12,14 +12,15 @@ Subroutine Condiff(scalar) real(8) St(Ic,Jc),Sm(Ic,Jc),fw1(Ic,Jc),alphastar(Ic,Jc),betastar(Ic,Jc),alpha(Ic,Jc),beta(Ic,Jc),Dwt(Ic,Jc),C3e(Ic,Jc) character(*) scalar character(6) wallfunutype,wallfunktype -logical(1) productlimit,sstlowre,sstcom,savortex +logical(1) productlimit,sstlowre,sstcom,saprodlimit wallfunutype='parvel' wallfunktype='loglaw' productlimit=.false. sstlowre=.false. sstcom=.false. -savortex=.false. +saprodlimit=.false. +Ym=11.225 !$OMP PARALLEL if(scalar=='U'.or.scalar=='V') then @@ -93,6 +94,15 @@ Subroutine Condiff(scalar) Ga=mu+mut/sigmatw !$OMP END WORKSHARE end if +if(Turmod=='ke') then +!$OMP WORKSHARE + Ymax=maxval(Ystar) +!$OMP END WORKSHARE +else +!$OMP WORKSHARE + Ymax=maxval(Yplus) +!$OMP END WORKSHARE +end if if(scalar=='T'.or.scalar=='Tk'.or.scalar=='Tw'.or.scalar=='Te') then !$OMP WORKSHARE St=sqrt(2*(Ux**2+Vy**2)+(Uy+Vx)**2) @@ -147,9 +157,11 @@ Subroutine Condiff(scalar) Xi=rho(i,j)*Tn(i,j)/mu(i,j) end if fnu1=Xi**3/(Xi**3+Cnu1**3) + !if(Ymax>Ym) fnu1=1.0 fnu2=1-Tn(i,j)/(mu(i,j)/rho(i,j)+Tn(i,j)*fnu1) + !if(Ymax>Ym) fnu2=0.0 Sv=abs(Uy(i,j)-Vx(i,j)) - if(savortex) Sv=Sv+2.0*min(0.0,St(i,j)-Sv) + if(saprodlimit) Sv=Sv+2.0*min(0.0,St(i,j)-Sv) Sm(i,j)=Sv+Tn(i,j)*fnu2/(kapa*d(i,j))**2 rm=Tn(i,j)/(Sm(i,j)*(kapa*d(i,j))**2) gm=rm+Cw2*(rm**6-rm) @@ -210,15 +222,6 @@ Subroutine Condiff(scalar) !$OMP WORKSHARE Tmin=minval(Tplus) !$OMP END WORKSHARE -if(Turmod=='ke') then -!$OMP WORKSHARE - Ymax=maxval(Ystar) -!$OMP END WORKSHARE -else -!$OMP WORKSHARE - Ymax=maxval(Yplus) -!$OMP END WORKSHARE -end if !$OMP WORKSHARE aM(1,:,:)=1 aM(2,:,:)=0 @@ -273,8 +276,9 @@ Subroutine Condiff(scalar) if(Tmptype=='fixed'.and.Tmin>=0) aP=aP+rho(i,j)*ustar(i)*DR(i)/Tplus(i) else if(scalar=='Tn') then aP=aP+2*Ds(i,j) - if(Ymax>10.and.fw1(i,j)>=0) aP=aP+rho(i,j)*Cw1*fw1(i,j)*Tn(i,j)/(kapa*d(i,j))**2*Jg(i,j)*dx*dy - if(Ymax<=10.and.fw1(i,j)>=0) aP=aP+rho(i,j)*Cw1*fw1(i,j)*Tn(i,j)/d(i,j)**2*Jg(i,j)*dx*dy + if(Ymax>Ym.and.fw1(i,j)>=0) aP=aP+rho(i,j)*Cw1*fw1(i,j)*Tn(i,j)/(kapa*d(i,j))**2*Jg(i,j)*dx*dy + if(Ymax<=Ym.and.fw1(i,j)>=0) aP=aP+rho(i,j)*Cw1*fw1(i,j)*Tn(i,j)/d(i,j)**2*Jg(i,j)*dx*dy + !if(fw1(i,j)>=0) aP=aP+rho(i,j)*Cw1*fw1(i,j)*Tn(i,j)/d(i,j)**2*Jg(i,j)*dx*dy else if(Turmod=='sst'.and.scalar=='Tk') then if(wallfunktype=='genlaw') then aP=aP+rho(i,j)*ustar(i)**3*Uplus(i)*Jg(i,j)*dx*dy/(Tk(i,j)*Yp(i)) @@ -303,8 +307,9 @@ Subroutine Condiff(scalar) if(scalar=='Te'.and.Turmod=='ke') aP=aP+C2e*rho(i,j)*Te(i,j)*Jg(i,j)*dx*dy/Tk(i,j) if(scalar=='Tn'.and.Walltreat=='lr'.and.fw1(i,j)>=0) aP=aP+rho(i,j)*Cw1*fw1(i,j)*Tn(i,j)/d(i,j)**2*Jg(i,j)*dx*dy if(scalar=='Tn'.and.Walltreat=='wf'.and.fw1(i,j)>=0) then - if(Ymax>10) aP=aP+rho(i,j)*Cw1*fw1(i,j)*Tn(i,j)/(kapa*d(i,j))**2*Jg(i,j)*dx*dy - if(Ymax<=10) aP=aP+rho(i,j)*Cw1*fw1(i,j)*Tn(i,j)/d(i,j)**2*Jg(i,j)*dx*dy + if(Ymax>Ym) aP=aP+rho(i,j)*Cw1*fw1(i,j)*Tn(i,j)/(kapa*d(i,j))**2*Jg(i,j)*dx*dy + if(Ymax<=Ym) aP=aP+rho(i,j)*Cw1*fw1(i,j)*Tn(i,j)/d(i,j)**2*Jg(i,j)*dx*dy + !aP=aP+rho(i,j)*Cw1*fw1(i,j)*Tn(i,j)/d(i,j)**2*Jg(i,j)*dx*dy end if if(scalar=='Tk'.and.Turmod=='sst') aP=aP+rho(i,j)*betastar(i,j)*Tw(i,j)*Jg(i,j)*dx*dy if(scalar=='Tw'.and.Turmod=='sst') aP=aP+rho(i,j)*beta(i,j)*Tw(i,j)*Jg(i,j)*dx*dy @@ -366,10 +371,10 @@ Subroutine Condiff(scalar) if(Proctrl=='com') then b(i,j)=Jg(i,j)*(U(i,j)*Px(i,j)+V(i,j)*Py(i,j))*dx*dy/ca+rho(i,j)*ustar(i)*Tf*DR(i)/Tplus(i) !b(i,j)=rho(i,j)*ustar(i)*Tf*DR(i)/Tplus(i) - if(visheat=='Y'.and.(Turmod=='sst'.or.Ymax<=10)) b(i,j)=b(i,j)+Jg(i,j)*(-2*(mu(i,j)+mut(i,j))*(Ux(i,j)+Vy(i,j))**2/3+(mu(i,j)+mut(i,j))*St(i,j)**2)*dx*dy/ca + if(visheat=='Y'.and.(Turmod=='sst'.or.Ymax<=Ym)) b(i,j)=b(i,j)+Jg(i,j)*(-2*(mu(i,j)+mut(i,j))*(Ux(i,j)+Vy(i,j))**2/3+(mu(i,j)+mut(i,j))*St(i,j)**2)*dx*dy/ca else if(Proctrl=='incom') then b(i,j)=rho(i,j)*ustar(i)*Tf*DR(i)/Tplus(i) - if(visheat=='Y'.and.(Turmod=='sst'.or.Ymax<=10)) b(i,j)=b(i,j)+Jg(i,j)*(mu(i,j)+mut(i,j))*St(i,j)**2*dx*dy/ca + if(visheat=='Y'.and.(Turmod=='sst'.or.Ymax<=Ym)) b(i,j)=b(i,j)+Jg(i,j)*(mu(i,j)+mut(i,j))*St(i,j)**2*dx*dy/ca end if if(Tmin<0) b(i,j)=b(i,j)-rho(i,j)*ustar(i)*T(i,j)*DR(i)/Tplus(i) end if @@ -416,7 +421,7 @@ Subroutine Condiff(scalar) else if(scalar=='Tn') then b(i,j)=Cb2*rho(i,j)*(Tnx(i,j)**2+Tny(i,j)**2)*Jg(i,j)*dx*dy/sigman+rho(i,j)*Cb1*Sm(i,j)*Tn(i,j)*Jg(i,j)*dx*dy if(fw1(i,j)<0) then - if(Walltreat=='wf'.and.Ymax>10) then + if(Walltreat=='wf'.and.Ymax>Ym) then b(i,j)=b(i,j)-rho(i,j)*Cw1*fw1(i,j)*(Tn(i,j)/(kapa*d(i,j)))**2*Jg(i,j)*dx*dy else b(i,j)=b(i,j)-rho(i,j)*Cw1*fw1(i,j)*(Tn(i,j)/d(i,j))**2*Jg(i,j)*dx*dy diff --git a/lib/Turvis.f03 b/lib/Turvis.f03 index 9769420..c25d22a 100644 --- a/lib/Turvis.f03 +++ b/lib/Turvis.f03 @@ -2,9 +2,15 @@ Subroutine Turvis use Aero2DCOM implicit none integer i,j -real(8) Xi,fnu1,Dwplus,phi1,F1,betai,St,phi2,F2,Ret,alphastar +real(8) Xi,fnu1,Dwplus,phi1,F1,betai,St,phi2,F2,Ret,alphastar,Ymax,Ym logical(1) sstlowre sstlowre=.false. + Ym=11.225 + if(Turmod=='ke') then + Ymax=maxval(Ystar) + else + Ymax=maxval(Yplus) + end if DO j=1,Jc DO i=1,Ic if(Turmod=='sa') then @@ -14,6 +20,7 @@ Subroutine Turvis Xi=rho(i,j)*Tn(i,j)/mu(i,j) end if fnu1=Xi**3/(Xi**3+Cnu1**3) + !if(Ymax>Ym) fnu1=1.0 mut(i,j)=rho(i,j)*Tn(i,j)*fnu1 else if(Turmod=='ke') then mut(i,j)=(rho(i,j)*Cu*Tk(i,j)**2)/Te(i,j) diff --git a/lib/Wallfunc.f03 b/lib/Wallfunc.f03 index 27ddc67..b53ce4a 100644 --- a/lib/Wallfunc.f03 +++ b/lib/Wallfunc.f03 @@ -60,7 +60,7 @@ Subroutine Wallfunc lamda=-0.01*(Pr(Ib1:Ib2,1)*Yplus)**4/(1+5.*Pr(Ib1:Ib2,1)**3*Yplus) Ymax=maxval(Yplus) DO i=Ib1,Ib2 - if(visheat=='Y'.and.Ymax>10) then + if(visheat=='Y'.and.Ymax>Ym) then Tplusl=Pr(i,1)*Uplusl(i) Tplust=Prt*(Uplust(i)+Prough(i)) !Tplust=Prt*Uplust(i)