diff --git a/Source/mass.f90 b/Source/mass.f90 index b2d4ce0aab3..9c66336ef02 100644 --- a/Source/mass.f90 +++ b/Source/mass.f90 @@ -20,15 +20,16 @@ MODULE MASS SUBROUTINE MASS_FINITE_DIFFERENCES(NM) USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE +USE PHYSICAL_FUNCTIONS, ONLY: GET_MOLECULAR_WEIGHT INTEGER, INTENT(IN) :: NM -REAL(EB) :: TNOW +REAL(EB) :: TNOW,MW_F,MW_G,ZZ_GET(1:N_TRACKED_SPECIES) REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,Z_TEMP,F_TEMP REAL(EB), PARAMETER :: DUMMY=0._EB INTEGER :: I,J,K,N,IOR,IW,IIG,JJG,KKG,II,JJ,KK,IC REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,RHOP REAL(EB), POINTER, DIMENSION(:,:,:,:) :: ZZP -REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P +REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P,RHO_RMW TYPE(WALL_TYPE), POINTER :: WC TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 @@ -181,11 +182,25 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM) ! Repeat the above for DENSITY - CALL GET_SCALAR_FACE_VALUE(UU,RHOP,FX(:,:,:,0),1,IBM1,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_VALUE(VV,RHOP,FY(:,:,:,0),1,IBAR,1,JBM1,1,KBAR,2,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_VALUE(WW,RHOP,FZ(:,:,:,0),1,IBAR,1,JBAR,1,KBM1,3,I_FLUX_LIMITER) + RHO_RMW=>WORK1 - !$OMP PARALLEL DO PRIVATE(IW,WC,BC,B1,II,JJ,KK,IIG,JJG,KKG,IOR,IC,U_TEMP,Z_TEMP,F_TEMP) + !$OMP PARALLEL DO PRIVATE(ZZ_GET,MW_G) SCHEDULE(STATIC) + DO K=0,KBP1 + DO J=0,JBP1 + DO I=0,IBP1 + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(I,J,K,1:N_TRACKED_SPECIES) + CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_G) + RHO_RMW(I,J,K) = RHOP(I,J,K)/MW_G + ENDDO + ENDDO + ENDDO + !$OMP END PARALLEL DO + + CALL GET_SCALAR_FACE_VALUE(UU,RHO_RMW,FX(:,:,:,0),1,IBM1,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) + CALL GET_SCALAR_FACE_VALUE(VV,RHO_RMW,FY(:,:,:,0),1,IBAR,1,JBM1,1,KBAR,2,I_FLUX_LIMITER) + CALL GET_SCALAR_FACE_VALUE(WW,RHO_RMW,FZ(:,:,:,0),1,IBAR,1,JBAR,1,KBM1,3,I_FLUX_LIMITER) + + !$OMP PARALLEL DO PRIVATE(IW,WC,BC,B1,II,JJ,KK,IIG,JJG,KKG,IOR,IC,U_TEMP,Z_TEMP,F_TEMP,ZZ_GET,MW_F) WALL_LOOP_3: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS WC=>WALL(IW) IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_LOOP_3 @@ -211,13 +226,15 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM) CASE(-3); FZ(IIG,JJG,KKG,0) = 0._EB END SELECT ELSE + ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES) + CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_F) SELECT CASE(IOR) - CASE( 1); FX(IIG-1,JJG,KKG,0) = B1%RHO_F - CASE(-1); FX(IIG,JJG,KKG,0) = B1%RHO_F - CASE( 2); FY(IIG,JJG-1,KKG,0) = B1%RHO_F - CASE(-2); FY(IIG,JJG,KKG,0) = B1%RHO_F - CASE( 3); FZ(IIG,JJG,KKG-1,0) = B1%RHO_F - CASE(-3); FZ(IIG,JJG,KKG,0) = B1%RHO_F + CASE( 1); FX(IIG-1,JJG,KKG,0) = B1%RHO_F/MW_F + CASE(-1); FX(IIG,JJG,KKG,0) = B1%RHO_F/MW_F + CASE( 2); FY(IIG,JJG-1,KKG,0) = B1%RHO_F/MW_F + CASE(-2); FY(IIG,JJG,KKG,0) = B1%RHO_F/MW_F + CASE( 3); FZ(IIG,JJG,KKG-1,0) = B1%RHO_F/MW_F + CASE(-3); FZ(IIG,JJG,KKG,0) = B1%RHO_F/MW_F END SELECT ENDIF @@ -231,7 +248,7 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM) ! /// II /// II+1 | II+2 | ... ! ^ WALL_INDEX(II+1,+1) IF ((UU(II+1,JJ,KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(II+1,JJ,KK))%WALL_INDEX(+1)>0)) THEN - Z_TEMP(0:3,1,1) = (/RHOP(II+1,JJ,KK),RHOP(II+1:II+2,JJ,KK),DUMMY/) + Z_TEMP(0:3,1,1) = (/RHO_RMW(II+1,JJ,KK),RHO_RMW(II+1:II+2,JJ,KK),DUMMY/) U_TEMP(1,1,1) = UU(II+1,JJ,KK) CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) FX(II+1,JJ,KK,0) = F_TEMP(1,1,1) @@ -241,35 +258,35 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM) ! ... | II-2 | II-1 /// II /// ! ^ WALL_INDEX(II-1,-1) IF ((UU(II-2,JJ,KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(II-1,JJ,KK))%WALL_INDEX(-1)>0)) THEN - Z_TEMP(0:3,1,1) = (/DUMMY,RHOP(II-2:II-1,JJ,KK),RHOP(II-1,JJ,KK)/) + Z_TEMP(0:3,1,1) = (/DUMMY,RHO_RMW(II-2:II-1,JJ,KK),RHO_RMW(II-1,JJ,KK)/) U_TEMP(1,1,1) = UU(II-2,JJ,KK) CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) FX(II-2,JJ,KK,0) = F_TEMP(1,1,1) ENDIF CASE( 2) OFF_WALL_SELECT_3 IF ((VV(II,JJ+1,KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ+1,KK))%WALL_INDEX(+2)>0)) THEN - Z_TEMP(1,0:3,1) = (/RHOP(II,JJ+1,KK),RHOP(II,JJ+1:JJ+2,KK),DUMMY/) + Z_TEMP(1,0:3,1) = (/RHO_RMW(II,JJ+1,KK),RHO_RMW(II,JJ+1:JJ+2,KK),DUMMY/) U_TEMP(1,1,1) = VV(II,JJ+1,KK) CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) FY(II,JJ+1,KK,0) = F_TEMP(1,1,1) ENDIF CASE(-2) OFF_WALL_SELECT_3 IF ((VV(II,JJ-2,KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ-1,KK))%WALL_INDEX(-2)>0)) THEN - Z_TEMP(1,0:3,1) = (/DUMMY,RHOP(II,JJ-2:JJ-1,KK),RHOP(II,JJ-1,KK)/) + Z_TEMP(1,0:3,1) = (/DUMMY,RHO_RMW(II,JJ-2:JJ-1,KK),RHO_RMW(II,JJ-1,KK)/) U_TEMP(1,1,1) = VV(II,JJ-2,KK) CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) FY(II,JJ-2,KK,0) = F_TEMP(1,1,1) ENDIF CASE( 3) OFF_WALL_SELECT_3 IF ((WW(II,JJ,KK+1)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ,KK+1))%WALL_INDEX(+3)>0)) THEN - Z_TEMP(1,1,0:3) = (/RHOP(II,JJ,KK+1),RHOP(II,JJ,KK+1:KK+2),DUMMY/) + Z_TEMP(1,1,0:3) = (/RHO_RMW(II,JJ,KK+1),RHO_RMW(II,JJ,KK+1:KK+2),DUMMY/) U_TEMP(1,1,1) = WW(II,JJ,KK+1) CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) FZ(II,JJ,KK+1,0) = F_TEMP(1,1,1) ENDIF CASE(-3) OFF_WALL_SELECT_3 IF ((WW(II,JJ,KK-2)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ,KK-1))%WALL_INDEX(-3)>0)) THEN - Z_TEMP(1,1,0:3) = (/DUMMY,RHOP(II,JJ,KK-2:KK-1),RHOP(II,JJ,KK-1)/) + Z_TEMP(1,1,0:3) = (/DUMMY,RHO_RMW(II,JJ,KK-2:KK-1),RHO_RMW(II,JJ,KK-1)/) U_TEMP(1,1,1) = WW(II,JJ,KK-2) CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) FZ(II,JJ,KK-2,0) = F_TEMP(1,1,1) @@ -281,20 +298,30 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM) ENDDO WALL_LOOP_3 !$OMP END PARALLEL DO - ! Now correct the face value of (RHO*ZZ) such that SUM(RHO*ZZ)_FACE = RHO_FACE + ! Now correct the max face value of (RHO*ZZ) such that SUM(RHO*ZZ/MW)_FACE = RHO_FACE/MW_FACE + ! (necessary condition to preserve isothermal flow) - !$OMP PARALLEL DO PRIVATE(N) SCHEDULE(STATIC) + !$OMP PARALLEL DO PRIVATE(N,MW_G) SCHEDULE(STATIC) DO K=0,KBAR DO J=0,JBAR DO I=0,IBAR N=MAXLOC(FX(I,J,K,1:N_TRACKED_SPECIES),1) - FX(I,J,K,N) = MAX( 0._EB, FX(I,J,K,0) - SUM(FX(I,J,K,1:(N-1))) - SUM(FX(I,J,K,(N+1):N_TRACKED_SPECIES)) ) + MW_G = SPECIES_MIXTURE(N)%MW + FX(I,J,K,N) = MW_G*MAX( 0._EB, FX(I,J,K,0) & + - SUM(FX(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) & + - SUM(FX(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) ) N=MAXLOC(FY(I,J,K,1:N_TRACKED_SPECIES),1) - FY(I,J,K,N) = MAX( 0._EB, FY(I,J,K,0) - SUM(FY(I,J,K,1:(N-1))) - SUM(FY(I,J,K,(N+1):N_TRACKED_SPECIES)) ) + MW_G = SPECIES_MIXTURE(N)%MW + FY(I,J,K,N) = MW_G*MAX( 0._EB, FY(I,J,K,0) & + - SUM(FY(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) & + - SUM(FY(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) ) N=MAXLOC(FZ(I,J,K,1:N_TRACKED_SPECIES),1) - FZ(I,J,K,N) = MAX( 0._EB, FZ(I,J,K,0) - SUM(FZ(I,J,K,1:(N-1))) - SUM(FZ(I,J,K,(N+1):N_TRACKED_SPECIES)) ) + MW_G = SPECIES_MIXTURE(N)%MW + FZ(I,J,K,N) = MW_G*MAX( 0._EB, FZ(I,J,K,0) & + - SUM(FZ(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) & + - SUM(FZ(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) ) ENDDO ENDDO ENDDO