Merge pull request #5238 from rmcdermo/master

FDS Source: implement WALE near wall turbulence model

Source/cons.f90

@@ -10,7 +10,7 @@ MODULE GLOBAL_CONSTANTS
 INTEGER, PARAMETER :: GAS_SPECIES=2,AEROSOL_SPECIES=3                                       ! For SPECIES%MODE
 INTEGER, PARAMETER :: EXPLICIT_EULER=1,RK2=2,RK3=3,RK2_RICHARDSON=4                         ! For COMBUSTION_ODE
 INTEGER, PARAMETER :: EXTINCTION_1=1,EXTINCTION_2=2,EXTINCTION_3=3                          ! For EXTINCT_MOD
-INTEGER, PARAMETER :: NO_TURB_MODEL=0,CONSMAG=1,DYNSMAG=2,DEARDORFF=3,VREMAN=4,RNG=5        ! Turbulence model
+INTEGER, PARAMETER :: NO_TURB_MODEL=0,CONSMAG=1,DYNSMAG=2,DEARDORFF=3,VREMAN=4,RNG=5,WALE=6 ! Turbulence model
 INTEGER, PARAMETER :: CONVECTIVE_FLUX_BC=-1,NET_FLUX_BC=0,SPECIFIED_TEMPERATURE=1,&
                       NO_CONVECTION=2,THERMALLY_THICK=3,INFLOW_OUTFLOW=4,&
                       INTERPOLATED_BC=6,THERMALLY_THICK_HT3D=7,&
@@ -214,7 +214,8 @@ MODULE GLOBAL_CONSTANTS
 ! Miscellaneous integer constants
 
 INTEGER :: ICYC,ICYC_RESTART=0,NFRAMES,PERIODIC_TEST=0,TURB_MODEL=0,SLIP_CONDITION=2,FISHPAK_BC(3)=-1,&
-           STOP_AT_ITER=0,HT3D_TEST=0,STOPFDS=-1,WALL_INCREMENT=2,WALL_INCREMENT_HT3D=1
+           STOP_AT_ITER=0,HT3D_TEST=0,STOPFDS=-1,WALL_INCREMENT=2,WALL_INCREMENT_HT3D=1,&
+           NEAR_WALL_TURB_MODEL=1
 
 ! Clocks for output file dumps
 

Source/read.f90

@@ -1604,7 +1604,8 @@ SUBROUTINE READ_MISC
 USE MATH_FUNCTIONS, ONLY: GET_RAMP_INDEX
 REAL(EB) :: CORIOLIS_VECTOR(3)=0._EB,FORCE_VECTOR(3)=0._EB,MAXIMUM_VISIBILITY
 CHARACTER(LABEL_LENGTH) :: ASSUMED_GAS_TEMPERATURE_RAMP,RAMP_FVX_T,RAMP_FVY_T,RAMP_FVZ_T,RAMP_GX,RAMP_GY,RAMP_GZ,&
-                           RAMP_U0_T,RAMP_V0_T,RAMP_W0_T,RAMP_U0_Z,RAMP_V0_Z,RAMP_W0_Z,RAMP_TMP0_Z,TURBULENCE_MODEL
+                           RAMP_U0_T,RAMP_V0_T,RAMP_W0_T,RAMP_U0_Z,RAMP_V0_Z,RAMP_W0_Z,RAMP_TMP0_Z,TURBULENCE_MODEL,&
+                           NEAR_WALL_TURBULENCE_MODEL
 
 NAMELIST /MISC/ AGGLOMERATION,AEROSOL_AL2O3,ALLOW_SURFACE_PARTICLES,ALLOW_UNDERSIDE_PARTICLES,ASSUMED_GAS_TEMPERATURE,&
                 ASSUMED_GAS_TEMPERATURE_RAMP,&
@@ -1622,7 +1623,7 @@ SUBROUTINE READ_MISC
                 IBLANK_SMV,IMMERSED_BOUNDARY_METHOD,INITIAL_UNMIXED_FRACTION,&
                 LAPSE_RATE,LES_FILTER_WIDTH,MAX_CHEMISTRY_ITERATIONS,&
                 MAX_LEAK_PATHS,MAXIMUM_VISIBILITY,MEAN_FORCING,MPI_TIMEOUT,&
-                N_FIXED_CHEMISTRY_SUBSTEPS,N_INITIAL_PARTICLE_SUBSTEPS,NEW_MOMENTUM_NUDGING,&
+                N_FIXED_CHEMISTRY_SUBSTEPS,N_INITIAL_PARTICLE_SUBSTEPS,NEAR_WALL_TURBULENCE_MODEL,NEW_MOMENTUM_NUDGING,&
                 NOISE,NOISE_VELOCITY,NO_EVACUATION,NO_RAMPS,NORTHANGLE,&
                 OVERWRITE,PARTICLE_CFL_MAX,PARTICLE_CFL_MIN,PARTICLE_CFL,PERIODIC_TEST,PROFILING,POROUS_FLOOR,&
                 POTENTIAL_TEMPERATURE_CORRECTION,PR,PROJECTION,P_INF,PROCESS_ALL_MESHES,RAMP_FVX_T,RAMP_FVY_T,RAMP_FVZ_T,&
@@ -1722,6 +1723,7 @@ SUBROUTINE READ_MISC
 MAXIMUM_VISIBILITY   = 30._EB ! m
 VISIBILITY_FACTOR    = 3._EB
 TURBULENCE_MODEL     = 'null'
+NEAR_WALL_TURBULENCE_MODEL = 'null'
 MAX_LEAK_PATHS       = 200
 NORTHANGLE = -190.0_EB
 IF (N_MPI_PROCESSES<=50) VERBOSE = .TRUE.
@@ -1867,13 +1869,24 @@ SUBROUTINE READ_MISC
       TURB_MODEL=VREMAN
    CASE ('RNG')
       TURB_MODEL=RNG
+   CASE ('WALE')
+      TURB_MODEL=WALE
    CASE ('null')
       TURB_MODEL=NO_TURB_MODEL
    CASE DEFAULT
       WRITE(MESSAGE,'(A,A,A)')  'ERROR: TURBULENCE_MODEL, ',TRIM(TURBULENCE_MODEL),', is not recognized.'
       CALL SHUTDOWN(MESSAGE) ; RETURN
 END SELECT
 
+! Near wall eddy viscosity model
+
+SELECT CASE (TRIM(NEAR_WALL_TURBULENCE_MODEL))
+   CASE DEFAULT
+      NEAR_WALL_TURB_MODEL=CONSMAG
+   CASE ('WALE')
+      NEAR_WALL_TURB_MODEL=WALE
+END SELECT
+
 ! Extinction Model
 
 IF (TRIM(EXTINCTION_MODEL)/='null') THEN

Source/turb.f90

@@ -12,14 +12,13 @@ MODULE TURBULENCE
 IMPLICIT NONE
 PRIVATE
 
-REAL(EB) :: DUMMY(3)
-
 PUBLIC :: INIT_TURB_ARRAYS, VARDEN_DYNSMAG, WANNIER_FLOW, &
           WALL_MODEL, COMPRESSION_WAVE, VELTAN2D,VELTAN3D, &
           SYNTHETIC_TURBULENCE, SYNTHETIC_EDDY_SETUP, TEST_FILTER, EX2G3D, TENSOR_DIFFUSIVITY_MODEL, &
           TWOD_VORTEX_CERFACS, TWOD_VORTEX_UMD, LOGLAW_HEAT_FLUX_MODEL, ABL_HEAT_FLUX_MODEL, RNG_EDDY_VISCOSITY, &
           NS_ANALYTICAL_SOLUTION, NS_U_EXACT, NS_V_EXACT, NS_H_EXACT, SANDIA_DAT, SPECTRAL_OUTPUT, SANDIA_OUT, &
-          FILL_EDGES, NATURAL_CONVECTION_MODEL, FORCED_CONVECTION_MODEL, RAYLEIGH_HEAT_FLUX_MODEL, YUAN_HEAT_FLUX_MODEL
+          FILL_EDGES, NATURAL_CONVECTION_MODEL, FORCED_CONVECTION_MODEL, RAYLEIGH_HEAT_FLUX_MODEL, YUAN_HEAT_FLUX_MODEL, &
+          WALE_VISCOSITY
 
 CONTAINS
 
@@ -1091,27 +1090,48 @@ SUBROUTINE RNG_EDDY_VISCOSITY(MU_EFF,MU_DNS_RNG,RHO,STRAIN_RATE,DELTA)
 END SUBROUTINE RNG_EDDY_VISCOSITY
 
 
-! SUBROUTINE WALE_VISCOSITY(MU_T,G_IJ,S,DELTA,RHO_G)
+SUBROUTINE WALE_VISCOSITY(NU_T,G,S,DELTA)
+
+! Wall Adaptive Local Eddy (WALE) Viscosity
 
-! ! F. Nicoud, F. Ducros. Subgrid-scale stress modelling based on the square of the velocity gradient tensor.
-! ! Flow, Turbulence, and Combustion, Vol. 62, pp. 183-200, 1999.
+! F. Nicoud, F. Ducros. Subgrid-scale stress modelling based on the square of the velocity gradient tensor.
+! Flow, Turbulence, and Combustion, Vol. 62, pp. 183-200, 1999.
 
-! REAL(EB), INTENT(OUT) :: MU_T
-! REAL(EB), INTENT(IN) :: S,G_IJ(3,3),DELTA,RHO_G
-! REAL(EB) :: G2_IJ(3,3)
-! REAL(EB), PARAMETER :: C_M=0.65_EB ! C_M**2 = 10.6 * C_S**2
+REAL(EB), INTENT(OUT) :: NU_T
+REAL(EB), INTENT(IN) :: S,G(3,3),DELTA
+REAL(EB) :: G2(3,3),SD(3,3),SD2,ONTH_G2_KK
+INTEGER :: I,J
+REAL(EB), PARAMETER :: C_M=0.65_EB ! C_M**2 = 10.6 * C_S**2
 
-! ! compute square of velocity gradient tensor
+! compute inner product of velocity gradient tensor
 
-! DO J=1,3
-!    DO I=1,3
-!       G2_IJ(I,J) = DOT_PRODUCT(G_IJ(I,:),G_IJ(:,J))
-!    ENDDO
-! ENDDO
+DO J=1,3
+   DO I=1,3
+      G2(I,J) = DOT_PRODUCT(G(I,:),G(:,J))
+   ENDDO
+ENDDO
+ONTH_G2_KK = ONTH * ( G2(1,1) + G2(2,2) + G2(3,3) )
+
+SD(1,1) = G2(1,1) - ONTH_G2_KK
+SD(2,2) = G2(2,2) - ONTH_G2_KK
+SD(3,3) = G2(3,3) - ONTH_G2_KK
+SD(1,2) = 0.5_EB * ( G2(1,2) + G2(2,1) )
+SD(1,3) = 0.5_EB * ( G2(1,3) + G2(3,1) )
+SD(2,3) = 0.5_EB * ( G2(2,3) + G2(3,2) )
+SD(2,1) = SD(1,2)
+SD(3,1) = SD(1,3)
+SD(3,2) = SD(2,3)
+
+SD2 = 0._EB
+DO J=1,3
+   DO I=1,3
+      SD2 = SD2 + SD(I,J)*SD(I,J)
+   ENDDO
+ENDDO
 
-! MU_T = RHO_G * (C_M*DELTA)**2 * SD2**1.5_EB / ( S**5 + SD2**1.25_EB )
+NU_T = (C_M*DELTA)**2 * SD2**1.5_EB / ( S**5 + SD2**1.25_EB )
 
-! END SUBROUTINE WALE_VISCOSITY
+END SUBROUTINE WALE_VISCOSITY
 
 
 SUBROUTINE WALL_MODEL(SLIP_FACTOR,U_TAU,Y_PLUS,U,NU,DY,S)
@@ -1483,7 +1503,7 @@ REAL(EB) FUNCTION VELTAN2D(U_VELO,U_SURF,NN,DN,DIVU,GRADU,GRADP,TAU_IJ,DT,RRHO,M
 
 REAL(EB), INTENT(IN) :: U_VELO(2),U_SURF(2),NN(2),DN,DIVU,GRADU(2,2),GRADP(2),TAU_IJ(2,2),DT,RRHO,MU
 INTEGER, INTENT(IN) :: I_VEL
-REAL(EB) :: C(2,2),SS(2),SLIP_COEF,ETA,AA,BB,U_STRM_0, &
+REAL(EB) :: C(2,2),SS(2),SLIP_COEF,ETA,AA,BB,U_STRM_0,DUMMY(3),&
             U_STRM,U_NORM,U_STRM_WALL,U_NORM_WALL,DPDS,DUSDS,DUSDN,TSN,RDN
 INTEGER :: SUBIT
 
@@ -1562,7 +1582,7 @@ REAL(EB) FUNCTION VELTAN3D(U_VELO,U_SURF,NN,DN,DIVU,GRADU,GRADP,TAU_IJ,DT,RRHO,M
 
 REAL(EB), INTENT(IN) :: U_VELO(3),U_SURF(3),NN(3),DN,DIVU,GRADU(3,3),GRADP(3),TAU_IJ(3,3),DT,RRHO,MU,ROUGHNESS,U_INT
 INTEGER, INTENT(IN) :: I_VEL
-REAL(EB) :: C(3,3),SS(3),PP(3),SLIP_COEF,ETA,AA,BB,U_STRM_0,U_RELA(3), &
+REAL(EB) :: C(3,3),SS(3),PP(3),SLIP_COEF,ETA,AA,BB,U_STRM_0,U_RELA(3),DUMMY(3),&
             U_STRM,U_ORTH,U_NORM,DPDS,DUSDS,DUSDN,TSN,DUPDP,DUNDN,RDN
 INTEGER :: SUBIT,I,J
 
@@ -2094,7 +2114,7 @@ SUBROUTINE SANDIA_DAT(NM,FN_ISO)
 ! This exe generates a random velocity field with a spectrum that matches the
 ! Comte-Bellot/Corrsin 1971 experimental data.
 
-REAL(EB) :: MEANU,MEANV,MEANW
+REAL(EB) :: MEANU,MEANV,MEANW,DUMMY(3)
 INTEGER :: I,J,K,II,JJ,KK,FILE_NUM,IM,IW,IOR,IIO,JJO,KKO,NX,NX_BLOCK,I_MESH,J_MESH,K_MESH,MX
 INTEGER, INTENT(IN) :: NM
 CHARACTER(80), INTENT(IN) :: FN_ISO

Source/velo.f90

@@ -45,7 +45,7 @@ END SUBROUTINE COMPUTE_VELOCITY_FLUX
 SUBROUTINE COMPUTE_VISCOSITY(T,NM)
 
 USE PHYSICAL_FUNCTIONS, ONLY: GET_VISCOSITY,LES_FILTER_WIDTH_FUNCTION,GET_POTENTIAL_TEMPERATURE
-USE TURBULENCE, ONLY: VARDEN_DYNSMAG,TEST_FILTER,FILL_EDGES,WALL_MODEL,RNG_EDDY_VISCOSITY
+USE TURBULENCE, ONLY: VARDEN_DYNSMAG,TEST_FILTER,FILL_EDGES,WALL_MODEL,RNG_EDDY_VISCOSITY,WALE_VISCOSITY
 USE MATH_FUNCTIONS, ONLY:EVALUATE_RAMP
 REAL(EB), INTENT(IN) :: T
 INTEGER, INTENT(IN) :: NM
@@ -296,6 +296,34 @@ SUBROUTINE COMPUTE_VISCOSITY(T,NM)
          ENDDO
       ENDDO
 
+   CASE (WALE) SELECT_TURB
+
+      DO K=1,KBAR
+         DO J=1,JBAR
+            DO I=1,IBAR
+               IF (SOLID(CELL_INDEX(I,J,K))) CYCLE
+               DELTA = LES_FILTER_WIDTH_FUNCTION(DX(I),DY(J),DZ(K))
+               ! compute velocity gradient tensor
+               DUDX = RDX(I)*(UU(I,J,K)-UU(I-1,J,K))
+               DVDY = RDY(J)*(VV(I,J,K)-VV(I,J-1,K))
+               DWDZ = RDZ(K)*(WW(I,J,K)-WW(I,J,K-1))
+               DUDY = 0.25_EB*RDY(J)*(UU(I,J+1,K)-UU(I,J-1,K)+UU(I-1,J+1,K)-UU(I-1,J-1,K))
+               DUDZ = 0.25_EB*RDZ(K)*(UU(I,J,K+1)-UU(I,J,K-1)+UU(I-1,J,K+1)-UU(I-1,J,K-1))
+               DVDX = 0.25_EB*RDX(I)*(VV(I+1,J,K)-VV(I-1,J,K)+VV(I+1,J-1,K)-VV(I-1,J-1,K))
+               DVDZ = 0.25_EB*RDZ(K)*(VV(I,J,K+1)-VV(I,J,K-1)+VV(I,J-1,K+1)-VV(I,J-1,K-1))
+               DWDX = 0.25_EB*RDX(I)*(WW(I+1,J,K)-WW(I-1,J,K)+WW(I+1,J,K-1)-WW(I-1,J,K-1))
+               DWDY = 0.25_EB*RDY(J)*(WW(I,J+1,K)-WW(I,J-1,K)+WW(I,J+1,K-1)-WW(I,J-1,K-1))
+               A_IJ(1,1)=DUDX; A_IJ(1,2)=DUDY; A_IJ(1,3)=DUDZ
+               A_IJ(2,1)=DVDX; A_IJ(2,2)=DVDY; A_IJ(2,3)=DVDZ
+               A_IJ(3,1)=DWDX; A_IJ(3,2)=DWDY; A_IJ(3,3)=DWDZ
+
+               CALL WALE_VISCOSITY(NU_EDDY,A_IJ,STRAIN_RATE(I,J,K),DELTA)
+
+               MU(I,J,K) = MU_DNS(I,J,K) + RHOP(I,J,K)*NU_EDDY
+            ENDDO
+         ENDDO
+      ENDDO
+
 END SELECT SELECT_TURB
 
 ! Compute resolved kinetic energy per unit mass
@@ -355,8 +383,26 @@ SUBROUTINE COMPUTE_VISCOSITY(T,NM)
 
          IF (LES) THEN
             DELTA = LES_FILTER_WIDTH_FUNCTION(DX(IIG),DY(JJG),DZ(KKG))
-            VDF = 1._EB-EXP(-WC%Y_PLUS*RAPLUS)
-            NU_EDDY = (VDF*C_SMAGORINSKY*DELTA)**2*STRAIN_RATE(IIG,JJG,KKG)
+            SELECT CASE(NEAR_WALL_TURB_MODEL)
+               CASE DEFAULT ! Constant Smagorinsky with Van Driest damping
+                  VDF = 1._EB-EXP(-WC%Y_PLUS*RAPLUS)
+                  NU_EDDY = (VDF*C_SMAGORINSKY*DELTA)**2*STRAIN_RATE(IIG,JJG,KKG)
+               CASE(WALE)
+                  ! compute velocity gradient tensor
+                  DUDX = RDX(IIG)*(UU(IIG,JJG,KKG)-UU(IIG-1,JJG,KKG))
+                  DVDY = RDY(JJG)*(VV(IIG,JJG,KKG)-VV(IIG,JJG-1,KKG))
+                  DWDZ = RDZ(KKG)*(WW(IIG,JJG,KKG)-WW(IIG,JJG,KKG-1))
+                  DUDY = 0.25_EB*RDY(JJG)*(UU(IIG,JJG+1,KKG)-UU(IIG,JJG-1,KKG)+UU(IIG-1,JJG+1,KKG)-UU(IIG-1,JJG-1,KKG))
+                  DUDZ = 0.25_EB*RDZ(KKG)*(UU(IIG,JJG,KKG+1)-UU(IIG,JJG,KKG-1)+UU(IIG-1,JJG,KKG+1)-UU(IIG-1,JJG,KKG-1))
+                  DVDX = 0.25_EB*RDX(IIG)*(VV(IIG+1,JJG,KKG)-VV(IIG-1,JJG,KKG)+VV(IIG+1,JJG-1,KKG)-VV(IIG-1,JJG-1,KKG))
+                  DVDZ = 0.25_EB*RDZ(KKG)*(VV(IIG,JJG,KKG+1)-VV(IIG,JJG,KKG-1)+VV(IIG,JJG-1,KKG+1)-VV(IIG,JJG-1,KKG-1))
+                  DWDX = 0.25_EB*RDX(IIG)*(WW(IIG+1,JJG,KKG)-WW(IIG-1,JJG,KKG)+WW(IIG+1,JJG,KKG-1)-WW(IIG-1,JJG,KKG-1))
+                  DWDY = 0.25_EB*RDY(JJG)*(WW(IIG,JJG+1,KKG)-WW(IIG,JJG-1,KKG)+WW(IIG,JJG+1,KKG-1)-WW(IIG,JJG-1,KKG-1))
+                  A_IJ(1,1)=DUDX; A_IJ(1,2)=DUDY; A_IJ(1,3)=DUDZ
+                  A_IJ(2,1)=DVDX; A_IJ(2,2)=DVDY; A_IJ(2,3)=DVDZ
+                  A_IJ(3,1)=DWDX; A_IJ(3,2)=DWDY; A_IJ(3,3)=DWDZ
+                  CALL WALE_VISCOSITY(NU_EDDY,A_IJ,STRAIN_RATE(IIG,JJG,KKG),DELTA)
+            END SELECT
             IF (CELL_COUNTER(IIG,JJG,KKG)==0) MU(IIG,JJG,KKG) = 0._EB
             CELL_COUNTER(IIG,JJG,KKG) = CELL_COUNTER(IIG,JJG,KKG) + 1
             WGT = 1._EB/REAL(CELL_COUNTER(IIG,JJG,KKG),EB)