diff --git a/Source/mass.f90 b/Source/mass.f90 index 9c66336ef02..09747e2239f 100644 --- a/Source/mass.f90 +++ b/Source/mass.f90 @@ -178,7 +178,7 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM) ENDDO SPECIES_LOOP -FACE_CORRECTION_IF: IF (TEST_FLUX_LIMITER_FACE_CORRECTION .AND. N_TRACKED_SPECIES>2) THEN +FACE_CORRECTION_IF: IF (TEST_FLUX_LIMITER_FACE_CORRECTION) THEN ! Repeat the above for DENSITY diff --git a/Source/wall.f90 b/Source/wall.f90 index 7c4fd128d87..11177d048d9 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -358,18 +358,18 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I 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 +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 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), & + RHO_ZZ_F(1:N_TOTAL_SCALARS),ZZ_GET(1:N_TRACKED_SPECIES),ZZ_GET_OTHER(1:N_TRACKED_SPECIES),& 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,SF_HTC,& - BBB,CCC,PPP,QQQ,RRR,UUU,YYY,WWW,HTC_OLD,RSC_LOC,DTMP,RSUM_G,MU_G + 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 @@ -651,6 +651,40 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I RHOP(BC%II,BC%JJ,BC%KK) = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/RSUM(BC%II,BC%JJ,BC%KK)/DENOM ENDIF + ! Correction for variable molecular weights + + IF (TEST_FLUX_LIMITER_FACE_CORRECTION) THEN + 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 (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) + ZZ_GET_OTHER(1:N_TRACKED_SPECIES) = OM_ZZP(IIO-1,JJO,KKO,1:N_TRACKED_SPECIES) + CASE(-1) + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%IIG-1,BC%JJG,BC%KKG,1:N_TRACKED_SPECIES) + ZZ_GET_OTHER(1:N_TRACKED_SPECIES) = OM_ZZP(IIO+1,JJO,KKO,1:N_TRACKED_SPECIES) + CASE( 2) + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%IIG,BC%JJG+1,BC%KKG,1:N_TRACKED_SPECIES) + ZZ_GET_OTHER(1:N_TRACKED_SPECIES) = OM_ZZP(IIO,JJO-1,KKO,1:N_TRACKED_SPECIES) + CASE(-2) + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%IIG,BC%JJG-1,BC%KKG,1:N_TRACKED_SPECIES) + ZZ_GET_OTHER(1:N_TRACKED_SPECIES) = OM_ZZP(IIO,JJO+1,KKO,1:N_TRACKED_SPECIES) + CASE( 3) + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%IIG,BC%JJG,BC%KKG+1,1:N_TRACKED_SPECIES) + ZZ_GET_OTHER(1:N_TRACKED_SPECIES) = OM_ZZP(IIO,JJO,KKO-1,1:N_TRACKED_SPECIES) + CASE(-3) + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%IIG,BC%JJG,BC%KKG-1,1:N_TRACKED_SPECIES) + ZZ_GET_OTHER(1:N_TRACKED_SPECIES) = OM_ZZP(IIO,JJO,KKO+1,1:N_TRACKED_SPECIES) + END SELECT + CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_G_2) + CALL GET_MOLECULAR_WEIGHT(ZZ_GET_OTHER,MW_OTHER_2) + ELSE + MW_G_2 = MW_G + MW_OTHER_2 = MW_OTHER + ENDIF + ENDIF + SELECT CASE(BC%IOR) CASE( 1) IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN @@ -658,6 +692,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I RHO_OTHER_2 = OM_RHOP(IIO-1,JJO,KKO) ENDIF Z_TEMP(0:3,1,1) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/) + IF (TEST_FLUX_LIMITER_FACE_CORRECTION) Z_TEMP(0:3,1,1) = Z_TEMP(0:3,1,1)/(/MW_OTHER_2,MW_OTHER,MW_G,MW_G_2/) U_TEMP(1,1,1) = UU(BC%II,BC%JJ,BC%KK) 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) @@ -667,6 +702,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I RHO_OTHER_2 = OM_RHOP(IIO+1,JJO,KKO) ENDIF Z_TEMP(0:3,1,1) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/) + IF (TEST_FLUX_LIMITER_FACE_CORRECTION) Z_TEMP(0:3,1,1) = Z_TEMP(0:3,1,1)/(/MW_G_2,MW_G,MW_OTHER,MW_OTHER_2/) U_TEMP(1,1,1) = UU(BC%II-1,BC%JJ,BC%KK) 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) @@ -676,6 +712,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I RHO_OTHER_2 = OM_RHOP(IIO,JJO-1,KKO) ENDIF Z_TEMP(1,0:3,1) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/) + IF (TEST_FLUX_LIMITER_FACE_CORRECTION) Z_TEMP(1,0:3,1) = Z_TEMP(1,0:3,1)/(/MW_OTHER_2,MW_OTHER,MW_G,MW_G_2/) U_TEMP(1,1,1) = VV(BC%II,BC%JJ,BC%KK) 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) @@ -685,6 +722,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I RHO_OTHER_2 = OM_RHOP(IIO,JJO+1,KKO) ENDIF Z_TEMP(1,0:3,1) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/) + IF (TEST_FLUX_LIMITER_FACE_CORRECTION) Z_TEMP(1,0:3,1) = Z_TEMP(1,0:3,1)/(/MW_G_2,MW_G,MW_OTHER,MW_OTHER_2/) U_TEMP(1,1,1) = VV(BC%II,BC%JJ-1,BC%KK) 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) @@ -694,6 +732,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I RHO_OTHER_2 = OM_RHOP(IIO,JJO,KKO-1) ENDIF Z_TEMP(1,1,0:3) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/) + IF (TEST_FLUX_LIMITER_FACE_CORRECTION) Z_TEMP(1,1,0:3) = Z_TEMP(1,1,0:3)/(/MW_OTHER_2,MW_OTHER,MW_G,MW_G_2/) U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK) 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) @@ -703,6 +742,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I RHO_OTHER_2 = OM_RHOP(IIO,JJO,KKO+1) ENDIF Z_TEMP(1,1,0:3) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/) + IF (TEST_FLUX_LIMITER_FACE_CORRECTION) Z_TEMP(1,1,0:3) = Z_TEMP(1,1,0:3)/(/MW_G_2,MW_G,MW_OTHER,MW_OTHER_2/) U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK-1) 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) @@ -783,6 +823,17 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I END SELECT ENDDO SPECIES_LOOP + ! Correction for variable molecular weights + + IF (TEST_FLUX_LIMITER_FACE_CORRECTION) THEN + N=MAXLOC(RHO_ZZ_F(1:N_TRACKED_SPECIES),1) + MW_G = SPECIES_MIXTURE(N)%MW + RHO_ZZ_F(N) = MW_G*MAX( 0._EB, B1%RHO_F & + - SUM(RHO_ZZ_F(1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) & + - SUM(RHO_ZZ_F((N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) ) + B1%RHO_F = SUM(RHO_ZZ_F(1:N_TRACKED_SPECIES)) + ENDIF + ! face value of temperature IF (ATMOSPHERIC_INTERPOLATION) THEN