From c76b7bbc00de2cc4e334e6a1999829e310d2da27 Mon Sep 17 00:00:00 2001 From: mcgratta Date: Thu, 1 Aug 2024 17:18:57 -0400 Subject: [PATCH] FDS Source: Issue #13235. Move call to WALL_MODEL --- Source/main.f90 | 8 ++++---- Source/velo.f90 | 29 +++-------------------------- Source/wall.f90 | 12 +++++++++--- 3 files changed, 16 insertions(+), 33 deletions(-) diff --git a/Source/main.f90 b/Source/main.f90 index 9d5239d1eef..a7686978ee0 100644 --- a/Source/main.f90 +++ b/Source/main.f90 @@ -346,7 +346,7 @@ PROGRAM FDS DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX IF (TGA_SURF_INDEX>0) CYCLE - CALL COMPUTE_VISCOSITY(T_BEGIN,NM,APPLY_TO_ESTIMATED_VARIABLES=.FALSE.) ! needed here for KRES prior to mesh exchange + CALL COMPUTE_VISCOSITY(NM,APPLY_TO_ESTIMATED_VARIABLES=.FALSE.) ! needed here for KRES prior to mesh exchange ENDDO IF (MY_RANK==0 .AND. VERBOSE) CALL VERBOSE_PRINTOUT('Completed COMPUTE_VISCOSITY') @@ -366,7 +366,7 @@ PROGRAM FDS DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX IF (TGA_SURF_INDEX>0) CYCLE CALL MATCH_VELOCITY(NM) - CALL COMPUTE_VISCOSITY(T_BEGIN,NM,APPLY_TO_ESTIMATED_VARIABLES=.FALSE.) ! call again after mesh exchange + CALL COMPUTE_VISCOSITY(NM,APPLY_TO_ESTIMATED_VARIABLES=.FALSE.) ! call again after mesh exchange IF (SYNTHETIC_EDDY_METHOD) CALL SYNTHETIC_EDDY_SETUP(NM) CALL VISCOSITY_BC(NM,APPLY_TO_ESTIMATED_VARIABLES=.FALSE.) CALL VELOCITY_BC(T_BEGIN,NM,APPLY_TO_ESTIMATED_VARIABLES=.FALSE.) @@ -617,7 +617,7 @@ PROGRAM FDS COMPUTE_FINITE_DIFFERENCES_1: DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX CALL INSERT_ALL_PARTICLES(T,NM) - IF (.NOT.SOLID_PHASE_ONLY .AND. .NOT.FREEZE_VELOCITY) CALL COMPUTE_VISCOSITY(T,NM,APPLY_TO_ESTIMATED_VARIABLES=.FALSE.) + IF (.NOT.SOLID_PHASE_ONLY .AND. .NOT.FREEZE_VELOCITY) CALL COMPUTE_VISCOSITY(NM,APPLY_TO_ESTIMATED_VARIABLES=.FALSE.) CALL MASS_FINITE_DIFFERENCES(NM) ENDDO COMPUTE_FINITE_DIFFERENCES_1 @@ -797,7 +797,7 @@ PROGRAM FDS ! Finite differences for mass and momentum equations for the second half of the time step COMPUTE_FINITE_DIFFERENCES_2: DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX - IF (.NOT.SOLID_PHASE_ONLY .AND. .NOT.FREEZE_VELOCITY) CALL COMPUTE_VISCOSITY(T,NM,APPLY_TO_ESTIMATED_VARIABLES=.TRUE.) + IF (.NOT.SOLID_PHASE_ONLY .AND. .NOT.FREEZE_VELOCITY) CALL COMPUTE_VISCOSITY(NM,APPLY_TO_ESTIMATED_VARIABLES=.TRUE.) CALL MASS_FINITE_DIFFERENCES(NM) CALL DENSITY(T,DT,NM) IF (LEVEL_SET_MODE>0) CALL LEVEL_SET_FIRESPREAD(T,DT,NM) diff --git a/Source/velo.f90 b/Source/velo.f90 index 503d0d1961e..d74a13d3de0 100644 --- a/Source/velo.f90 +++ b/Source/velo.f90 @@ -21,23 +21,20 @@ MODULE VELO !> \brief Compute the viscosity of the gas. !> \callergraph !> \callgraph -!> \param T Current time (s). !> \param NM Mesh number. !> \param APPLY_TO_ESTIMATED_VARIABLES Flag indicating \f$\mu(T,\mathbf{u})\f$ or \f$\mu(T^*,\mathbf{u}^*)\f$ -SUBROUTINE COMPUTE_VISCOSITY(T,NM,APPLY_TO_ESTIMATED_VARIABLES) +SUBROUTINE COMPUTE_VISCOSITY(NM,APPLY_TO_ESTIMATED_VARIABLES) USE PHYSICAL_FUNCTIONS, ONLY: GET_VISCOSITY,GET_POTENTIAL_TEMPERATURE,GET_CONDUCTIVITY,GET_SPECIFIC_HEAT -USE TURBULENCE, ONLY: VARDEN_DYNSMAG,TEST_FILTER,FILL_EDGES,WALL_MODEL,WALE_VISCOSITY +USE TURBULENCE, ONLY: VARDEN_DYNSMAG,TEST_FILTER,FILL_EDGES,WALE_VISCOSITY USE MATH_FUNCTIONS, ONLY:EVALUATE_RAMP USE CC_SCALARS, ONLY : CC_COMPUTE_KRES,CC_COMPUTE_VISCOSITY,CUTFACE_VELOCITIES -REAL(EB), INTENT(IN) :: T INTEGER, INTENT(IN) :: NM LOGICAL, INTENT(IN) :: APPLY_TO_ESTIMATED_VARIABLES REAL(EB), ALLOCATABLE, DIMENSION(:) :: ZZ_GET REAL(EB) :: NU_EDDY,DELTA,KSGS,U2,V2,W2,AA,A_IJ(3,3),BB,B_IJ(3,3),& - DUDX,DUDY,DUDZ,DVDX,DVDY,DVDZ,DWDX,DWDY,DWDZ,SLIP_COEF,VEL_GAS,VEL_T,RAMP_T,TSI,& - VDF,WGT,T_NOW + DUDX,DUDY,DUDZ,DVDX,DVDY,DVDZ,DWDX,DWDY,DWDZ,VDF,WGT,T_NOW REAL(EB), PARAMETER :: RAPLUS=1._EB/26._EB INTEGER :: I,J,K,IIG,JJG,KKG,II,JJ,KK,IW,IOR,IC REAL(EB), POINTER, DIMENSION(:,:,:) :: RHOP=>NULL(),UP=>NULL(),VP=>NULL(),WP=>NULL(), & @@ -329,23 +326,6 @@ SUBROUTINE COMPUTE_VISCOSITY(T,NM,APPLY_TO_ESTIMATED_VARIABLES) CASE(SOLID_BOUNDARY) - IF (ABS(SF%T_IGN-T_BEGIN)<=SPACING(SF%T_IGN) .AND. SF%RAMP(TIME_VELO)%INDEX>=1) THEN - TSI = T - ELSE - TSI = T-SF%T_IGN - ENDIF - RAMP_T = EVALUATE_RAMP(TSI,SF%RAMP(TIME_VELO)%INDEX,TAU=SF%RAMP(TIME_VELO)%TAU) - VEL_T = RAMP_T*SQRT(SF%VEL_T(1)**2 + SF%VEL_T(2)**2) - - SELECT CASE(ABS(IOR)) - CASE(1) - VEL_GAS = SQRT( 0.25_EB*( (VV(IIG,JJG,KKG)+VV(IIG,JJG-1,KKG))**2 + (WW(IIG,JJG,KKG)+WW(IIG,JJG,KKG-1))**2 ) ) - CASE(2) - VEL_GAS = SQRT( 0.25_EB*( (UU(IIG,JJG,KKG)+UU(IIG-1,JJG,KKG))**2 + (WW(IIG,JJG,KKG)+WW(IIG,JJG,KKG-1))**2 ) ) - CASE(3) - VEL_GAS = SQRT( 0.25_EB*( (UU(IIG,JJG,KKG)+UU(IIG-1,JJG,KKG))**2 + (VV(IIG,JJG,KKG)+VV(IIG,JJG-1,KKG))**2 ) ) - END SELECT - IF (SIM_MODE/=DNS_MODE) THEN DELTA = LES_FILTER_WIDTH(IIG,JJG,KKG) SELECT CASE(SF%NEAR_WALL_TURB_MODEL) @@ -382,9 +362,6 @@ SUBROUTINE COMPUTE_VISCOSITY(T,NM,APPLY_TO_ESTIMATED_VARIABLES) IF (CELL(CELL_INDEX(II,JJ,KK))%SOLID) MU(II,JJ,KK) = MU(IIG,JJG,KKG) - CALL WALL_MODEL(SLIP_COEF,B2%U_TAU,B2%Y_PLUS,MU_DNS(IIG,JJG,KKG)/RHO(IIG,JJG,KKG),SF%ROUGHNESS,& - 0.5_EB/B1%RDN,VEL_GAS-VEL_T) - CASE(OPEN_BOUNDARY,MIRROR_BOUNDARY) MU(II,JJ,KK) = MU(IIG,JJG,KKG) diff --git a/Source/wall.f90 b/Source/wall.f90 index 8356c832851..4bdaa4b2132 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -27,12 +27,13 @@ MODULE WALL_ROUTINES SUBROUTINE WALL_BC(T,DT,NM) USE COMP_FUNCTIONS, ONLY: CURRENT_TIME -USE CC_SCALARS, ONLY : CFACE_THERMAL_GASVARS +USE CC_SCALARS, ONLY: CFACE_THERMAL_GASVARS +USE TURBULENCE, ONLY: WALL_MODEL REAL(EB) :: TNOW REAL(EB), INTENT(IN) :: T,DT INTEGER, INTENT(IN) :: NM LOGICAL :: CALL_HT_1D -REAL(EB) :: DT_BC +REAL(EB) :: DT_BC,SLIP_COEF INTEGER :: IW,IP,ICF,ITW TYPE(WALL_TYPE), POINTER :: WC TYPE(SURFACE_TYPE), POINTER :: SF @@ -83,7 +84,7 @@ SUBROUTINE WALL_BC(T,DT,NM) ENDIF ENDIF -!$OMP PARALLEL PRIVATE(IW,WC,SF,BC,B1,B2,LP,LPC,CFA,IP,ICF) +!$OMP PARALLEL PRIVATE(IW,WC,SF,BC,B1,B2,LP,LPC,CFA,IP,ICF,SLIP_COEF) ! Sweep through all WALL cells and apply boundary conditions @@ -109,6 +110,11 @@ SUBROUTINE WALL_BC(T,DT,NM) CALL CALCULATE_RHO_D_F(B1,BC,WALL_INDEX=IW) ENDIF + IF (WC%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. (ANY(SPECIES_MIXTURE%CONDENSATION_SMIX_INDEX>0).OR.DEPOSITION)) THEN + CALL WALL_MODEL(SLIP_COEF,B2%U_TAU,B2%Y_PLUS,MU_DNS(BC%IIG,BC%JJG,BC%KKG)/RHO(BC%IIG,BC%JJG,BC%KKG),SF%ROUGHNESS,& + 0.5_EB/B1%RDN,B1%U_TANG) + ENDIF + IF (DEPOSITION .AND. .NOT.INITIALIZATION_PHASE .AND. CORRECTOR .AND. .NOT.SOLID_PHASE_ONLY) THEN IF (WC%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. B1%NODE_INDEX==0 .AND. ABS(SF%VEL)