From 22b02b50f7fd5fb331fb692a99f7c4635f2b95f5 Mon Sep 17 00:00:00 2001 From: rmcdermo Date: Fri, 14 Jul 2017 17:25:50 -0400 Subject: [PATCH] FDS Source: implement WALE near wall turbulence model --- Source/cons.f90 | 5 ++-- Source/read.f90 | 17 ++++++++++++-- Source/turb.f90 | 62 ++++++++++++++++++++++++++++++++----------------- Source/velo.f90 | 52 ++++++++++++++++++++++++++++++++++++++--- 4 files changed, 108 insertions(+), 28 deletions(-) diff --git a/Source/cons.f90 b/Source/cons.f90 index 0cc8d3980d6..f7995cf7f5e 100644 --- a/Source/cons.f90 +++ b/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 diff --git a/Source/read.f90 b/Source/read.f90 index 658d119bdf9..920570e96c2 100644 --- a/Source/read.f90 +++ b/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,6 +1869,8 @@ 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 @@ -1874,6 +1878,15 @@ SUBROUTINE READ_MISC 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 diff --git a/Source/turb.f90 b/Source/turb.f90 index 0be72e58f42..bfae37fce6e 100644 --- a/Source/turb.f90 +++ b/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 diff --git a/Source/velo.f90 b/Source/velo.f90 index e8049c8930e..7ed60e0dac6 100644 --- a/Source/velo.f90 +++ b/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)