diff --git a/Source/func.f90 b/Source/func.f90 index d29d714ecc..6c95c51be5 100644 --- a/Source/func.f90 +++ b/Source/func.f90 @@ -1451,159 +1451,6 @@ END FUNCTION MP5 END SUBROUTINE GET_SCALAR_FACE_VALUE -SUBROUTINE GET_SCALAR_FACE_COEF(A,U,BF,I1,I2,J1,J2,K1,K2,IOR,LIMITER) - -USE GLOBAL_CONSTANTS, ONLY: CENTRAL_LIMITER,GODUNOV_LIMITER,SUPERBEE_LIMITER,CHARM_LIMITER - -REAL(EB), INTENT(IN) :: A(0:,0:,0:),U(-1:,-1:,-1:) -REAL(EB), INTENT(INOUT) :: BF(0:,0:,0:) -INTEGER, INTENT(IN) :: LIMITER,I1,I2,J1,J2,K1,K2,IOR -REAL(EB) :: R,B,DU_UP,DU_LOC -INTEGER :: I,J,K,IM1,JM1,KM1,IP1,JP1,KP1,IP2,JP2,KP2 - -SELECT CASE(LIMITER) - CASE(GODUNOV_LIMITER,CENTRAL_LIMITER) - BF=0._EB - RETURN -END SELECT - -SELECT CASE(IOR) - CASE(1) ; IM1=-1 ; JM1= 0 ; KM1= 0 ; IP1=1 ; JP1=0 ; KP1=0 ; IP2=2 ; JP2=0 ; KP2=0 - CASE(2) ; IM1= 0 ; JM1=-1 ; KM1= 0 ; IP1=0 ; JP1=1 ; KP1=0 ; IP2=0 ; JP2=2 ; KP2=0 - CASE(3) ; IM1= 0 ; JM1= 0 ; KM1=-1 ; IP1=0 ; JP1=0 ; KP1=1 ; IP2=0 ; JP2=0 ; KP2=2 -END SELECT - -!$OMP PARALLEL DEFAULT(NONE) IF(I2>1) & -!$OMP& SHARED(U,A,BF,I1,I2,J1,J2,K1,K2,IP1,JP1,KP1,IM1,JM1,KM1,IP2,JP2,KP2,LIMITER) & -!$OMP& PRIVATE(I,J,K,DU_UP,DU_LOC,B,R) - -SELECT CASE (LIMITER) - - CASE (SUPERBEE_LIMITER) - !$OMP DO COLLAPSE(3) SCHEDULE(STATIC) - DO K=K1,K2 - DO J=J1,J2 - DO I=I1,I2 - DU_LOC = U(I+IP1,J+JP1,K+KP1) - U(I,J,K) - IF (A(I,J,K) > 0._EB) THEN - DU_UP = U(I,J,K) - U(I+IM1,J+JM1,K+KM1) - ELSE - DU_UP = U(I+IP2,J+JP2,K+KP2) - U(I+IP1,J+JP1,K+KP1) - ENDIF - R = 0._EB ; B = 0._EB - IF (ABS(DU_LOC) > TWO_EPSILON_EB) R = DU_UP/DU_LOC - IF (R > TWO_EPSILON_EB) B = MAX(0._EB, MIN(2._EB*R,1._EB), MIN(R,2._EB)) - BF(I,J,K) = MAX(0._EB, MIN(B, BF(I,J,K))) - ENDDO - ENDDO - ENDDO - !$OMP END DO - - CASE (CHARM_LIMITER) - !$OMP DO COLLAPSE(3) SCHEDULE(STATIC) - DO K=K1,K2 - DO J=J1,J2 - DO I=I1,I2 - DU_LOC = U(I+IP1,J+JP1,K+KP1) - U(I,J,K) - IF (A(I,J,K) > 0._EB) THEN - DU_UP = U(I,J,K) - U(I+IM1,J+JM1,K+KM1) - ELSE - DU_UP = U(I+IP2,J+JP2,K+KP2) - U(I+IP1,J+JP1,K+KP1) - ENDIF - R = 0._EB ; B = 0._EB - IF (ABS(DU_UP) > TWO_EPSILON_EB) R = DU_LOC/DU_UP - IF (R > TWO_EPSILON_EB) B = R*(3._EB*R+1._EB)/((R+1._EB)**2) - BF(I,J,K) = MAX(0._EB, MIN(B, BF(I,J,K))) - ENDDO - ENDDO - ENDDO - !$OMP END DO - -END SELECT -!$OMP END PARALLEL - -END SUBROUTINE GET_SCALAR_FACE_COEF - - -SUBROUTINE GET_SCALAR_FACE_VALUE_NEW(A,U,F,BF,I1,I2,J1,J2,K1,K2,IOR,LIMITER) - -USE GLOBAL_CONSTANTS, ONLY: CENTRAL_LIMITER,GODUNOV_LIMITER,SUPERBEE_LIMITER,CHARM_LIMITER - -REAL(EB), INTENT(IN) :: A(0:,0:,0:),U(-1:,-1:,-1:),BF(0:,0:,0:) -REAL(EB), INTENT(OUT) :: F(0:,0:,0:) -INTEGER, INTENT(IN) :: LIMITER,I1,I2,J1,J2,K1,K2,IOR -REAL(EB) :: DU_UP,DU_LOC -INTEGER :: I,J,K,IM1,JM1,KM1,IP1,JP1,KP1,IP2,JP2,KP2 - -SELECT CASE(IOR) - CASE(1) ; IM1=-1 ; JM1= 0 ; KM1= 0 ; IP1=1 ; JP1=0 ; KP1=0 ; IP2=2 ; JP2=0 ; KP2=0 - CASE(2) ; IM1= 0 ; JM1=-1 ; KM1= 0 ; IP1=0 ; JP1=1 ; KP1=0 ; IP2=0 ; JP2=2 ; KP2=0 - CASE(3) ; IM1= 0 ; JM1= 0 ; KM1=-1 ; IP1=0 ; JP1=0 ; KP1=1 ; IP2=0 ; JP2=0 ; KP2=2 -END SELECT - -!$OMP PARALLEL IF(I2>1) -SELECT CASE(LIMITER) - CASE(CENTRAL_LIMITER) ! central differencing - !$OMP DO SCHEDULE(STATIC) - DO K=K1,K2 - DO J=J1,J2 - DO I=I1,I2 - F(I,J,K) = 0.5_EB*(U(I,J,K) + U(I+IP1,J+JP1,K+KP1)) - ENDDO - ENDDO - ENDDO - !$OMP END DO - CASE(GODUNOV_LIMITER) ! first-order upwinding - !$OMP DO SCHEDULE(STATIC) - DO K=K1,K2 - DO J=J1,J2 - DO I=I1,I2 - IF (A(I,J,K)>0._EB) THEN - F(I,J,K) = U(I,J,K) - ELSE - F(I,J,K) = U(I+IP1,J+JP1,K+KP1) - ENDIF - ENDDO - ENDDO - ENDDO - !$OMP END DO - CASE(SUPERBEE_LIMITER) ! SUPERBEE, Roe (1986) - !$OMP DO SCHEDULE(STATIC) PRIVATE(DU_LOC) - DO K=K1,K2 - DO J=J1,J2 - DO I=I1,I2 - DU_LOC = U(I+IP1,J+JP1,K+KP1)-U(I,J,K) - IF (A(I,J,K)>0._EB) THEN - F(I,J,K) = U(I,J,K) + 0.5_EB*BF(I,J,K)*DU_LOC - ELSE - F(I,J,K) = U(I+IP1,J+JP1,K+KP1) - 0.5_EB*BF(I,J,K)*DU_LOC - ENDIF - ENDDO - ENDDO - ENDDO - !$OMP END DO - CASE(CHARM_LIMITER) ! CHARM - !$OMP DO SCHEDULE(STATIC) PRIVATE(DU_UP) - DO K=K1,K2 - DO J=J1,J2 - DO I=I1,I2 - IF (A(I,J,K)>0._EB) THEN - DU_UP = U(I,J,K)-U(I+IM1,J+JM1,K+KM1) - F(I,J,K) = U(I,J,K) + 0.5_EB*BF(I,J,K)*DU_UP - ELSE - DU_UP = U(I+IP2,J+JP2,K+KP2)-U(I+IP1,J+JP1,K+KP1) - F(I,J,K) = U(I+IP1,J+JP1,K+KP1) - 0.5_EB*BF(I,J,K)*DU_UP - ENDIF - ENDDO - ENDDO - ENDDO - !$OMP END DO -END SELECT -!$OMP END PARALLEL - -END SUBROUTINE GET_SCALAR_FACE_VALUE_NEW - - !> \brief Random fluctuation, Theta'(t+dt) = R^2*Theta'(t) + Normal(0,sqrt(1-R^2)*SIGMA) ; R = exp(-dt/TAU) !> \param SIGMA Standard deviation of time series !> \param TAU Time scale or period of the time function (s) diff --git a/Source/wall.f90 b/Source/wall.f90 index b482fffa5b..5bba8080fd 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -31,7 +31,7 @@ SUBROUTINE WALL_BC(T,DT,NM) REAL(EB) :: TNOW REAL(EB), INTENT(IN) :: T,DT INTEGER, INTENT(IN) :: NM -LOGICAL :: CALL_HT_1D,SECOND_ORDER_INTERPOLATED_BOUNDARY +LOGICAL :: CALL_HT_1D REAL(EB) :: DT_BC,SLIP_COEF INTEGER :: IW,IP,ICF,ITW TYPE(WALL_TYPE), POINTER :: WC @@ -94,7 +94,7 @@ SUBROUTINE WALL_BC(T,DT,NM) WC=>WALL(IW) BC => BOUNDARY_COORD(WC%BC_INDEX) B1 => BOUNDARY_PROP1(WC%B1_INDEX) - IF (IW<=N_EXTERNAL_WALL_CELLS) CALL ASSIGN_GHOST_VALUE(IW,BC,B1,SECOND_ORDER_INTERPOLATED_BOUNDARY) + IF (IW<=N_EXTERNAL_WALL_CELLS) CALL ASSIGN_GHOST_VALUE(IW,BC,B1) IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_CELL_LOOP_0 SF => SURFACE(WC%SURF_INDEX) CALL NEAR_SURFACE_GAS_VARIABLES(T,SF,BC,B1,WALL_INDEX=IW) @@ -135,7 +135,7 @@ SUBROUTINE WALL_BC(T,DT,NM) SF => SURFACE(WC%SURF_INDEX) IF (.NOT.SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN - CALL SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX=IW,SOIB=SECOND_ORDER_INTERPOLATED_BOUNDARY) + CALL SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX=IW) ELSEIF (CALL_HT_1D) THEN CALL SOLID_HEAT_TRANSFER(NM,T,SF%HT_DIM*DT_BC,WALL_INDEX=IW) ENDIF @@ -262,16 +262,15 @@ SUBROUTINE WALL_BC(T,DT,NM) END SUBROUTINE WALL_BC -SUBROUTINE ASSIGN_GHOST_VALUE(WALL_INDEX,BC,B1,SECOND_ORDER_INTERPOLATED_BOUNDARY) +SUBROUTINE ASSIGN_GHOST_VALUE(WALL_INDEX,BC,B1) USE PHYSICAL_FUNCTIONS, ONLY : GET_SPECIFIC_GAS_CONSTANT USE COMPLEX_GEOMETRY, ONLY : CC_CGSC, CC_SOLID INTEGER, INTENT(IN) :: WALL_INDEX -LOGICAL, INTENT(OUT) :: SECOND_ORDER_INTERPOLATED_BOUNDARY REAL(EB) :: ARO,ZZ_GET(1:N_TRACKED_SPECIES),RHO_OTHER,RHO_ZZ_OTHER(1:N_TOTAL_SCALARS),& RSUM_TMP,RHO_OTHER_2,RHO_ZZ_OTHER_2(1:N_TOTAL_SCALARS),PBAR_P_2 INTEGER :: IIO,JJO,KKO,II2,JJ2,KK2,ICG,ICO -LOGICAL :: CC_SOLID_FLAG +LOGICAL :: CC_SOLID_FLAG,SECOND_ORDER_INTERPOLATED_BOUNDARY REAL(EB), POINTER, DIMENSION(:,:,:,:) :: OM_ZZP REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_RHOP TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC @@ -474,25 +473,24 @@ END SUBROUTINE NEAR_SURFACE_GAS_VARIABLES !> \param PARTICLE_INDEX Optional Lagrangian particle index !> \details Calculate the surface temperature TMP_F of either a WALL cell, immersed CFACE cell, or a Lagrangian particle. -SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX,SOIB) +SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX) -USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP,INTERPOLATE1D_UNIFORM,GET_SCALAR_FACE_VALUE +USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP,INTERPOLATE1D_UNIFORM, GET_SCALAR_FACE_VALUE +USE COMPLEX_GEOMETRY, ONLY : CC_CGSC, CC_SOLID USE PHYSICAL_FUNCTIONS, ONLY : GET_SPECIFIC_GAS_CONSTANT,GET_VISCOSITY,GET_MOLECULAR_WEIGHT USE DEVICE_VARIABLES, ONLY : PROPERTY,PROPERTY_TYPE REAL(EB), INTENT(IN) :: T INTEGER, INTENT(IN) :: NM INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX -LOGICAL, INTENT(IN), OPTIONAL :: SOIB !Second Order Interpolated Boundary -REAL(EB) :: ARO,QNET,RAMP_FACTOR,RSUM_F,PBAR_F,TSI,UN, & +REAL(EB) :: ARO,QNET,RAMP_FACTOR,RHO_G_2,RSUM_F,PBAR_F,TSI,UN, & RHO_ZZ_F(1:N_TOTAL_SCALARS),ZZ_GET(1:N_TRACKED_SPECIES),ZZ_GET_OTHER(1:N_TRACKED_SPECIES),& - RHO_OTHER, PBAR_OTHER,D_Z_N(0:I_MAX_TEMP),D_Z_G,D_Z_OTHER,TMP_OTHER,DDO,DENOM,PBAR_G, & - MU_DNS_G,MU_DNS_OTHER,MU_OTHER,RHO_D,RHO_D_TURB,RHO_D_DZDN,RHO_D_DZDN_OTHER,RSUM_OTHER, & - BBB,CCC,PPP,QQQ,RRR,UUU,YYY,WWW,HTC_OLD,RSC_LOC,DTMP,RSUM_G,MU_G,MW_G,MW_G_2,MW_OTHER,MW_OTHER_2, & - RHO_G_2,RHO_OTHER_2,RHO_ZZ_G,RHO_ZZ_G_2,RHO_ZZ_OTHER_2,RHO_ZZ_OTHER(1:N_TOTAL_SCALARS) -LOGICAL :: ATMOSPHERIC_INTERPOLATION -INTEGER :: IIO,JJO,KKO,N,ADCOUNT -REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_RHOP + RHO_OTHER,RHO_OTHER_2,RHO_ZZ_OTHER(1:N_TOTAL_SCALARS),RHO_ZZ_OTHER_2,RHO_ZZ_G,RHO_ZZ_G_2, & + DDO,PBAR_G,PBAR_OTHER,DENOM,D_Z_N(0:I_MAX_TEMP),D_Z_G,D_Z_OTHER,TMP_OTHER, & + MU_DNS_G,MU_DNS_OTHER,MU_OTHER,RHO_D,RHO_D_TURB,RHO_D_DZDN,RHO_D_DZDN_OTHER,RSUM_OTHER,& + BBB,CCC,PPP,QQQ,RRR,UUU,YYY,WWW,HTC_OLD,RSC_LOC,DTMP,RSUM_G,MU_G,MW_G,MW_G_2,MW_OTHER,MW_OTHER_2 +LOGICAL :: SECOND_ORDER_INTERPOLATED_BOUNDARY,ATMOSPHERIC_INTERPOLATION,CC_SOLID_FLAG +INTEGER :: IIO,JJO,KKO,N,ADCOUNT,ICG,ICO REAL(EB), POINTER, DIMENSION(:,:,:,:) :: OM_ZZP REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_MUP REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP @@ -508,6 +506,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC +REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_RHOP TYPE(PROPERTY_TYPE), POINTER :: PY IF (PRESENT(WALL_INDEX)) THEN @@ -582,9 +581,9 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I ! Ghost cell values - TMP(BC%II,BC%JJ,BC%KK) = B1%TMP_F + TMP(BC%II,BC%JJ,BC%KK) = B1%TMP_F TMP(BC%II2,BC%JJ2,BC%KK2) = B1%TMP_F - ZZP(BC%II,BC%JJ,BC%KK,1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES) + ZZP(BC%II,BC%JJ,BC%KK,1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES) ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES) ZZ_GET(1:N_TRACKED_SPECIES) = MAX(0._EB,ZZP(BC%II,BC%JJ,BC%KK,1:N_TRACKED_SPECIES)) CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM(BC%II,BC%JJ,BC%KK)) @@ -697,7 +696,6 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I EWC => EXTERNAL_WALL(WALL_INDEX) OM => OMESH(EWC%NOM) - MM => MESHES(EWC%NOM) IF (PREDICTOR) THEN OM_RHOP => OM%RHOS OM_ZZP => OM%ZZS @@ -706,9 +704,36 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I OM_ZZP => OM%ZZ ENDIF - ! interp or extrap RHO_OTHER for jump in vertical grid resolution, linear in temperature to match heat flux in divg + MM => MESHES(EWC%NOM) + RHO_G_2 = B1%RHO_G ! first-order + RHO_OTHER = RHOP(BC%II,BC%JJ,BC%KK) + RHO_OTHER_2 = RHOP(BC%II,BC%JJ,BC%KK) + RHO_ZZ_OTHER(1:N_TOTAL_SCALARS) = ZZP(BC%II,BC%JJ,BC%KK,1:N_TOTAL_SCALARS)*RHO_OTHER + + ! Determine if there are 4 equally sized cells spanning the interpolated boundary + + SECOND_ORDER_INTERPOLATED_BOUNDARY = .FALSE. + CC_SOLID_FLAG = .FALSE. + IF (ABS(EWC%AREA_RATIO-1._EB)<0.01_EB) THEN + IIO = EWC%IIO_MIN + JJO = EWC%JJO_MIN + KKO = EWC%KKO_MIN + SELECT CASE(BC%IOR) + CASE( 1) ; ICG = CELL_INDEX(BC%IIG+1,BC%JJG,BC%KKG) ; ICO = MM%CELL_INDEX(IIO-1,JJO,KKO) + CASE(-1) ; ICG = CELL_INDEX(BC%IIG-1,BC%JJG,BC%KKG) ; ICO = MM%CELL_INDEX(IIO+1,JJO,KKO) + CASE( 2) ; ICG = CELL_INDEX(BC%IIG,BC%JJG+1,BC%KKG) ; ICO = MM%CELL_INDEX(IIO,JJO-1,KKO) + CASE(-2) ; ICG = CELL_INDEX(BC%IIG,BC%JJG-1,BC%KKG) ; ICO = MM%CELL_INDEX(IIO,JJO+1,KKO) + CASE( 3) ; ICG = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG+1) ; ICO = MM%CELL_INDEX(IIO,JJO,KKO-1) + CASE(-3) ; ICG = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG-1) ; ICO = MM%CELL_INDEX(IIO,JJO,KKO+1) + END SELECT + IF (CC_IBM) THEN ! Test if one of surrounding cells is CC_SOLID. + IF(CCVAR(BC%IIG,BC%JJG,BC%KKG,CC_CGSC)==CC_SOLID .OR. CCVAR(BC%II,BC%JJ,BC%KK,CC_CGSC)==CC_SOLID) CC_SOLID_FLAG = .TRUE. + ENDIF + IF (.NOT.CELL(ICG)%SOLID.AND..NOT.MM%CELL(ICO)%SOLID.AND..NOT.CC_SOLID_FLAG) SECOND_ORDER_INTERPOLATED_BOUNDARY = .TRUE. + ENDIF + + ! Density - RHO_OTHER = RHOP(BC%II,BC%JJ,BC%KK) ATMOSPHERIC_INTERPOLATION = .FALSE. DDO = 1._EB KKO = EWC%KKO_MIN @@ -720,6 +745,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I ATMOSPHERIC_INTERPOLATION = .TRUE. IF (ATMOSPHERIC_INTERPOLATION) THEN + ! interp or extrap RHO_OTHER for jump in vertical grid resolution, linear in temperature to match heat flux in divg PBAR_G = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) PBAR_OTHER = EVALUATE_RAMP(MM%ZC(EWC%KKO_MIN),I_RAMP_P0_Z) ZZ_GET(1:N_TRACKED_SPECIES) = MAX(0._EB,MIN(1._EB,ZZP(BC%II,BC%JJ,BC%KK,1:N_TOTAL_SCALARS))) @@ -729,23 +755,12 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I TMP(BC%II,BC%JJ,BC%KK) = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/(RSUM(BC%II,BC%JJ,BC%KK)*RHOP(BC%II,BC%JJ,BC%KK)) ENDIF - ! Density, Species and temperature - ! Correction for variable molecular weights IF (FLUX_LIMITER_MW_CORRECTION) THEN - - RHO_G_2 = B1%RHO_G ! first-order - RHO_OTHER = RHOP(BC%II,BC%JJ,BC%KK) - RHO_OTHER_2 = RHOP(BC%II,BC%JJ,BC%KK) - RHO_ZZ_OTHER(1:N_TOTAL_SCALARS) = ZZP(BC%II,BC%JJ,BC%KK,1:N_TOTAL_SCALARS)*RHO_OTHER - - IIO = EWC%IIO_MIN - JJO = EWC%JJO_MIN - ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%IIG,BC%JJG,BC%KKG,1:N_TRACKED_SPECIES); CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_G) ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%II,BC%JJ,BC%KK,1:N_TRACKED_SPECIES) ; CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_OTHER) - IF (SOIB) THEN + IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN SELECT CASE(BC%IOR) CASE( 1) ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%IIG+1,BC%JJG,BC%KKG,1:N_TRACKED_SPECIES) @@ -776,7 +791,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I SELECT CASE(BC%IOR) CASE( 1) - IF (SOIB) THEN + IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN RHO_G_2 = RHOP(BC%IIG+1,BC%JJG,BC%KKG) RHO_OTHER_2 = OM_RHOP(IIO-1,JJO,KKO) ENDIF @@ -786,7 +801,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) B1%RHO_F = F_TEMP(1,1,1) CASE(-1) - IF (SOIB) THEN + IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN RHO_G_2 = RHOP(BC%IIG-1,BC%JJG,BC%KKG) RHO_OTHER_2 = OM_RHOP(IIO+1,JJO,KKO) ENDIF @@ -796,7 +811,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) B1%RHO_F = F_TEMP(1,1,1) CASE( 2) - IF (SOIB) THEN + IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN RHO_G_2 = RHOP(BC%IIG,BC%JJG+1,BC%KKG) RHO_OTHER_2 = OM_RHOP(IIO,JJO-1,KKO) ENDIF @@ -806,7 +821,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) B1%RHO_F = F_TEMP(1,1,1) CASE(-2) - IF (SOIB) THEN + IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN RHO_G_2 = RHOP(BC%IIG,BC%JJG-1,BC%KKG) RHO_OTHER_2 = OM_RHOP(IIO,JJO+1,KKO) ENDIF @@ -816,7 +831,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) B1%RHO_F = F_TEMP(1,1,1) CASE( 3) - IF (SOIB) THEN + IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN RHO_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG+1) RHO_OTHER_2 = OM_RHOP(IIO,JJO,KKO-1) ENDIF @@ -826,7 +841,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) B1%RHO_F = F_TEMP(1,1,1) CASE(-3) - IF (SOIB) THEN + IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN RHO_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG-1) RHO_OTHER_2 = OM_RHOP(IIO,JJO,KKO+1) ENDIF @@ -848,7 +863,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I SELECT CASE(BC%IOR) CASE( 1) - IF (SOIB) THEN + IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN RHO_ZZ_G_2 = RHOP(BC%IIG+1,BC%JJG,BC%KKG)*ZZP(BC%IIG+1,BC%JJG,BC%KKG,N) RHO_ZZ_OTHER_2 = OM_RHOP(IIO-1,JJO,KKO)*OM_ZZP(IIO-1,JJO,KKO,N) ENDIF @@ -858,7 +873,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I RHO_ZZ_F(N) = F_TEMP(1,1,1) PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) CASE(-1) - IF (SOIB) THEN + IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN RHO_ZZ_G_2 = RHOP(BC%IIG-1,BC%JJG,BC%KKG)*ZZP(BC%IIG-1,BC%JJG,BC%KKG,N) RHO_ZZ_OTHER_2 = OM_RHOP(IIO+1,JJO,KKO)*OM_ZZP(IIO+1,JJO,KKO,N) ENDIF @@ -868,7 +883,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I RHO_ZZ_F(N) = F_TEMP(1,1,1) PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) CASE( 2) - IF (SOIB) THEN + IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG+1,BC%KKG)*ZZP(BC%IIG,BC%JJG+1,BC%KKG,N) RHO_ZZ_OTHER_2 = OM_RHOP(IIO,JJO-1,KKO)*OM_ZZP(IIO,JJO-1,KKO,N) ENDIF @@ -878,7 +893,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I RHO_ZZ_F(N) = F_TEMP(1,1,1) PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) CASE(-2) - IF (SOIB) THEN + IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG-1,BC%KKG)*ZZP(BC%IIG,BC%JJG-1,BC%KKG,N) RHO_ZZ_OTHER_2 = OM_RHOP(IIO,JJO+1,KKO)*OM_ZZP(IIO,JJO+1,KKO,N) ENDIF @@ -888,7 +903,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I RHO_ZZ_F(N) = F_TEMP(1,1,1) PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) CASE( 3) - IF (SOIB) THEN + IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG+1)*ZZP(BC%IIG,BC%JJG,BC%KKG+1,N) RHO_ZZ_OTHER_2 = OM_RHOP(IIO,JJO,KKO-1)*OM_ZZP(IIO,JJO,KKO-1,N) ENDIF @@ -899,7 +914,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I PBAR_F = (PBAR_P(BC%KK,B1%PRESSURE_ZONE)*DZ(BC%KKG) + PBAR_P(BC%KKG,B1%PRESSURE_ZONE)*DZ(BC%KK)) / & (DZ(BC%KK)+DZ(BC%KKG)) CASE(-3) - IF (SOIB) THEN + IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG-1)*ZZP(BC%IIG,BC%JJG,BC%KKG-1,N) RHO_ZZ_OTHER_2 = OM_RHOP(IIO,JJO,KKO+1)*OM_ZZP(IIO,JJO,KKO+1,N) ENDIF @@ -926,16 +941,24 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I ! face value of temperature IF (ATMOSPHERIC_INTERPOLATION) THEN - B1%TMP_F = (TMP(BC%II,BC%JJ,BC%KK)*DZ(BC%KKG) + TMP(BC%IIG,BC%JJG,BC%KKG)*DZ(BC%KK)) / (DZ(BC%KK)+DZ(BC%KKG)) B1%ZZ_F(1:N_TOTAL_SCALARS) = (ZZP(BC%II,BC%JJ,BC%KK,1:N_TOTAL_SCALARS)*DZ(BC%KKG) + & ZZP(BC%IIG,BC%JJG,BC%KKG,1:N_TOTAL_SCALARS)*DZ(BC%KK)) / (DZ(BC%KK)+DZ(BC%KKG)) ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES) CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_F) - B1%RHO_F = PBAR_F/(RSUM_F*B1%TMP_F) + SELECT CASE(ABS(BC%IOR)) + CASE DEFAULT + PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) + CASE(3) + ! for very large cells with stratification the pressure variation in z can be significant + PBAR_F = (PBAR_P(BC%KK,B1%PRESSURE_ZONE)*DZ(BC%KKG) + PBAR_P(BC%KKG,B1%PRESSURE_ZONE)*DZ(BC%KK)) / & + (DZ(BC%KK)+DZ(BC%KKG)) + END SELECT + B1%TMP_F = PBAR_F/(RSUM_F*B1%RHO_F) ELSE B1%ZZ_F(1:N_TOTAL_SCALARS) = MAX(0._EB,MIN(1._EB,RHO_ZZ_F(1:N_TOTAL_SCALARS)/B1%RHO_F)) ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES) CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_F) + PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) B1%TMP_F = PBAR_F/(RSUM_F*B1%RHO_F) ENDIF