From 6589a426be75607d69f1e546977d0f03ef77d76f Mon Sep 17 00:00:00 2001 From: rmcdermo Date: Tue, 3 Feb 2026 17:17:14 -0500 Subject: [PATCH] FDS Source: fix a few things in TENSOR_DIFFUSIVITY_MODEL --- Source/divg.f90 | 2 +- Source/turb.f90 | 80 +++++++++++++++++++++---------------------------- Source/wall.f90 | 2 ++ 3 files changed, 37 insertions(+), 47 deletions(-) diff --git a/Source/divg.f90 b/Source/divg.f90 index f315634fbc..3809e1205e 100644 --- a/Source/divg.f90 +++ b/Source/divg.f90 @@ -236,7 +236,7 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) ! Ensure RHO_D terms sum to zero over all species. Gather error into largest mass fraction present. - IF (SIM_MODE==DNS_MODE .OR. SIM_MODE==LES_MODE) THEN + IF (SIM_MODE==DNS_MODE .OR. SIM_MODE==LES_MODE .OR. TENSOR_DIFFUSIVITY) THEN ! for VLES and SVLES modes, the diffusivity is the same for all species ! so, as long as ZZP is realizable, the sum of diffusive fluxes will be zero ! and the flux corrections below are not required diff --git a/Source/turb.f90 b/Source/turb.f90 index 93c09b7fb3..96d70b3c22 100644 --- a/Source/turb.f90 +++ b/Source/turb.f90 @@ -1611,7 +1611,7 @@ SUBROUTINE TENSOR_DIFFUSIVITY_MODEL(NM,OPT_N) INTEGER, INTENT(IN) :: NM INTEGER, INTENT(IN), OPTIONAL :: OPT_N INTEGER :: I,J,K,N -REAL(EB) :: DRHOZDX,DRHOZDY,DRHOZDZ,DUDX,DUDY,DUDZ,DVDX,DVDY,DVDZ,DWDX,DWDY,DWDZ,DTDX,DTDY,DTDZ +REAL(EB) :: DZDX,DZDY,DZDZ,DUDX,DUDY,DUDZ,DVDX,DVDY,DVDZ,DWDX,DWDY,DWDZ,DTDX,DTDY,DTDZ,RHOBAR REAL(EB), POINTER, DIMENSION(:,:,:,:) :: ZZP,RHO_D_DZDX,RHO_D_DZDY,RHO_D_DZDZ REAL(EB), POINTER, DIMENSION(:,:,:) :: RHOP,UU,VV,WW,KDTDX,KDTDY,KDTDZ REAL(EB), PARAMETER :: C_NL=0.083_EB ! C_NL=1/12, See Pope Exercise 13.28 @@ -1646,22 +1646,18 @@ SUBROUTINE TENSOR_DIFFUSIVITY_MODEL(NM,OPT_N) DO K=1,KBAR DO J=1,JBAR DO I=0,IBAR + RHOBAR = 0.5_EB*(RHOP(I,J,K)+RHOP(I+1,J,K)) - DUDX = (UU(I+1,J,K)-UU(I-1,J,K))/(DX(I)+DX(I+1)) - DUDY = (UU(I,J+1,K)-UU(I,J-1,K))/(DYN(J)+DYN(J+1)) - DUDZ = (UU(I,J,K+1)-UU(I,J,K-1))/(DZN(K)+DZN(K+1)) + DUDX = (UU(I+1,J,K)-UU(MAX(0,I-1),J,K))/(DX(I)+DX(I+1)) + DUDY = (UU(I,J+1,K)-UU(I,J-1,K))/(DYN(J-1)+DYN(J)) + DUDZ = (UU(I,J,K+1)-UU(I,J,K-1))/(DZN(K-1)+DZN(K)) - DRHOZDX = RDXN(I)*(RHOP(I+1,J,K)*ZZP(I+1,J,K,N)-RHOP(I,J,K)*ZZP(I,J,K,N)) - - DRHOZDY = 0.25_EB*RDY(J)*( RHOP(I,J+1,K)*ZZP(I,J+1,K,N) + RHOP(I+1,J+1,K)*ZZP(I+1,J+1,K,N) & - - RHOP(I,J-1,K)*ZZP(I,J-1,K,N) - RHOP(I+1,J-1,K)*ZZP(I+1,J-1,K,N) ) - - DRHOZDZ = 0.25_EB*RDZ(K)*( RHOP(I,J,K+1)*ZZP(I,J,K+1,N) + RHOP(I+1,J,K+1)*ZZP(I+1,J,K+1,N) & - - RHOP(I,J,K-1)*ZZP(I,J,K-1,N) - RHOP(I+1,J,K-1)*ZZP(I+1,J,K-1,N) ) + DZDX = RDXN(I)*(ZZP(I+1,J,K,N)-ZZP(I,J,K,N)) + DZDY = 0.25_EB*RDY(J)*( ZZP(I,J+1,K,N) + ZZP(I+1,J+1,K,N) - ZZP(I,J-1,K,N) - ZZP(I+1,J-1,K,N) ) + DZDZ = 0.25_EB*RDZ(K)*( ZZP(I,J,K+1,N) + ZZP(I+1,J,K+1,N) - ZZP(I,J,K-1,N) - ZZP(I+1,J,K-1,N) ) RHO_D_DZDX(I,J,K,N) = RHO_D_DZDX(I,J,K,N) & - + C_NL*(DX(I)**2*DUDX*DRHOZDX+DY(J)**2*DUDY*DRHOZDY+DZ(K)**2*DUDZ*DRHOZDZ) - + + RHOBAR*C_NL*( DXN(I)**2*DZDX*DUDX + DY(J)**2*DZDY*DUDY + DZ(K)**2*DZDZ*DUDZ ) ENDDO ENDDO ENDDO @@ -1669,22 +1665,18 @@ SUBROUTINE TENSOR_DIFFUSIVITY_MODEL(NM,OPT_N) DO K=1,KBAR DO J=0,JBAR DO I=1,IBAR + RHOBAR = 0.5_EB*(RHOP(I,J,K)+RHOP(I,J+1,K)) - DVDX = (VV(I+1,J,K)-VV(I-1,J,K))/(DXN(I)+DXN(I+1)) - DVDY = (VV(I,J+1,K)-VV(I,J-1,K))/(DY(J)+DY(J+1)) - DVDZ = (VV(I,J,K+1)-VV(I,J,K-1))/(DZN(K)+DZN(K+1)) - - DRHOZDX = 0.25_EB*RDX(I)*( RHOP(I+1,J,K)*ZZP(I+1,J,K,N) + RHOP(I+1,J+1,K)*ZZP(I+1,J+1,K,N) & - - RHOP(I-1,J,K)*ZZP(I-1,J,K,N) - RHOP(I-1,J+1,K)*ZZP(I-1,J+1,K,N) ) - - DRHOZDY = RDYN(J)*(RHOP(I,J+1,K)*ZZP(I,J+1,K,N)-RHOP(I,J,K)*ZZP(I,J,K,N)) + DVDX = (VV(I+1,J,K)-VV(I-1,J,K))/(DXN(I-1)+DXN(I)) + DVDY = (VV(I,J+1,K)-VV(I,MAX(0,J-1),K))/(DY(J)+DY(J+1)) + DVDZ = (VV(I,J,K+1)-VV(I,J,K-1))/(DZN(K-1)+DZN(K)) - DRHOZDZ = 0.25_EB*RDZ(K)*( RHOP(I,J,K+1)*ZZP(I,J,K+1,N) + RHOP(I,J+1,K+1)*ZZP(I,J+1,K+1,N) & - - RHOP(I,J,K-1)*ZZP(I,J,K-1,N) - RHOP(I,J+1,K-1)*ZZP(I,J+1,K-1,N) ) + DZDX = 0.25_EB*RDX(I)*( ZZP(I+1,J,K,N) + ZZP(I+1,J+1,K,N) - ZZP(I-1,J,K,N) - ZZP(I-1,J+1,K,N) ) + DZDY = RDYN(J)*(ZZP(I,J+1,K,N)-ZZP(I,J,K,N)) + DZDZ = 0.25_EB*RDZ(K)*( ZZP(I,J,K+1,N) + ZZP(I,J+1,K+1,N) - ZZP(I,J,K-1,N) - ZZP(I,J+1,K-1,N) ) RHO_D_DZDY(I,J,K,N) = RHO_D_DZDY(I,J,K,N) & - + C_NL*(DX(I)**2*DVDX*DRHOZDX+DY(J)**2*DVDY*DRHOZDY+DZ(K)**2*DVDZ*DRHOZDZ) - + + RHOBAR*C_NL*( DX(I)**2*DZDX*DVDX + DY(J)**2*DZDY*DVDY + DZ(K)**2*DZDZ*DVDZ ) ENDDO ENDDO ENDDO @@ -1692,22 +1684,18 @@ SUBROUTINE TENSOR_DIFFUSIVITY_MODEL(NM,OPT_N) DO K=0,KBAR DO J=1,JBAR DO I=1,IBAR + RHOBAR = 0.5_EB*(RHOP(I,J,K)+RHOP(I,J,K+1)) - DWDX = (WW(I+1,J,K)-WW(I-1,J,K))/(DXN(I+1)+DXN(I) ) - DWDY = (WW(I,J+1,K)-WW(I,J-1,K))/(DYN(J+1)+DYN(J) ) - DWDZ = (WW(I,J,K+1)-WW(I,J,K-1))/(DZ(K) +DZ(K+1)) - - DRHOZDX = 0.25_EB*RDX(I)*( RHOP(I+1,J,K)*ZZP(I+1,J,K,N) + RHOP(I+1,J,K+1)*ZZP(I+1,J,K+1,N) & - - RHOP(I-1,J,K)*ZZP(I-1,J,K,N) - RHOP(I-1,J,K+1)*ZZP(I-1,J,K+1,N) ) + DWDX = (WW(I+1,J,K)-WW(I-1,J,K))/(DXN(I-1)+DXN(I)) + DWDY = (WW(I,J+1,K)-WW(I,J-1,K))/(DYN(J-1)+DYN(J)) + DWDZ = (WW(I,J,K+1)-WW(I,J,MAX(0,K-1)))/(DZ(K)+DZ(K+1)) - DRHOZDY = 0.25_EB*RDY(J)*( RHOP(I,J+1,K)*ZZP(I,J+1,K,N) + RHOP(I,J+1,K+1)*ZZP(I,J+1,K+1,N) & - - RHOP(I,J-1,K)*ZZP(I,J-1,K,N) - RHOP(I,J-1,K+1)*ZZP(I,J-1,K+1,N) ) - - DRHOZDZ = RDZN(K)*(RHOP(I,J,K+1)*ZZP(I,J,K+1,N)-RHOP(I,J,K)*ZZP(I,J,K,N)) + DZDX = 0.25_EB*RDX(I)*( ZZP(I+1,J,K,N) + ZZP(I+1,J,K+1,N) - ZZP(I-1,J,K,N) - ZZP(I-1,J,K+1,N) ) + DZDY = 0.25_EB*RDY(J)*( ZZP(I,J+1,K,N) + ZZP(I,J+1,K+1,N) - ZZP(I,J-1,K,N) - ZZP(I,J-1,K+1,N) ) + DZDZ = RDZN(K)*(ZZP(I,J,K+1,N)-ZZP(I,J,K,N)) RHO_D_DZDZ(I,J,K,N) = RHO_D_DZDZ(I,J,K,N) & - + C_NL*(DX(I)**2*DWDX*DRHOZDX+DY(J)**2*DWDY*DRHOZDY+DZ(K)**2*DWDZ*DRHOZDZ) - + + RHOBAR*C_NL*( DX(I)**2*DZDX*DWDX + DY(J)**2*DZDY*DWDY + DZ(K)**2*DZDZ*DWDZ ) ENDDO ENDDO ENDDO @@ -1735,9 +1723,9 @@ SUBROUTINE TENSOR_DIFFUSIVITY_MODEL(NM,OPT_N) DO J=1,JBAR DO I=0,IBAR - DUDX = (UU(I+1,J,K)-UU(I-1,J,K))/(DX(I)+DX(I+1)) - DUDY = (UU(I,J+1,K)-UU(I,J-1,K))/(DYN(J)+DYN(J+1)) - DUDZ = (UU(I,J,K+1)-UU(I,J,K-1))/(DZN(K)+DZN(K+1)) + DUDX = (UU(I+1,J,K)-UU(MAX(0,I-1),J,K))/(DX(I)+DX(I+1)) + DUDY = (UU(I,J+1,K)-UU(I,J-1,K))/(DYN(J-1)+DYN(J)) + DUDZ = (UU(I,J,K+1)-UU(I,J,K-1))/(DZN(K-1)+DZN(K)) DTDX = RDXN(I)*(TMP(I+1,J,K)-TMP(I,J,K)) DTDY = 0.25_EB*RDY(J)*( TMP(I,J+1,K) + TMP(I+1,J+1,K) - TMP(I,J-1,K) - TMP(I+1,J-1,K) ) @@ -1753,9 +1741,9 @@ SUBROUTINE TENSOR_DIFFUSIVITY_MODEL(NM,OPT_N) DO J=0,JBAR DO I=1,IBAR - DVDX = (VV(I+1,J,K)-VV(I-1,J,K))/(DXN(I)+DXN(I+1)) - DVDY = (VV(I,J+1,K)-VV(I,J-1,K))/(DY(J)+DY(J+1)) - DVDZ = (VV(I,J,K+1)-VV(I,J,K-1))/(DZN(K)+DZN(K+1)) + DVDX = (VV(I+1,J,K)-VV(I-1,J,K))/(DXN(I-1)+DXN(I)) + DVDY = (VV(I,J+1,K)-VV(I,MAX(0,J-1),K))/(DY(J)+DY(J+1)) + DVDZ = (VV(I,J,K+1)-VV(I,J,K-1))/(DZN(K-1)+DZN(K)) DTDX = 0.25_EB*RDX(I)*( TMP(I+1,J,K) + TMP(I+1,J+1,K) - TMP(I-1,J,K) - TMP(I-1,J+1,K) ) DTDY = RDYN(J)*(TMP(I,J+1,K)-TMP(I,J,K)) @@ -1771,9 +1759,9 @@ SUBROUTINE TENSOR_DIFFUSIVITY_MODEL(NM,OPT_N) DO J=1,JBAR DO I=1,IBAR - DWDX = (WW(I+1,J,K)-WW(I-1,J,K))/(DXN(I+1)+DXN(I) ) - DWDY = (WW(I,J+1,K)-WW(I,J-1,K))/(DYN(J+1)+DYN(J) ) - DWDZ = (WW(I,J,K+1)-WW(I,J,K-1))/(DZ(K) +DZ(K+1)) + DWDX = (WW(I+1,J,K)-WW(I-1,J,K))/(DXN(I-1)+DXN(I)) + DWDY = (WW(I,J+1,K)-WW(I,J-1,K))/(DYN(J-1)+DYN(J)) + DWDZ = (WW(I,J,K+1)-WW(I,J,MAX(0,K-1)))/(DZ(K)+DZ(K+1)) DTDX = 0.25_EB*RDX(I)*( TMP(I+1,J,K) + TMP(I+1,J,K+1) - TMP(I-1,J,K) - TMP(I-1,J,K+1) ) DTDY = 0.25_EB*RDY(J)*( TMP(I,J+1,K) + TMP(I,J+1,K+1) - TMP(I,J-1,K) - TMP(I,J-1,K+1) ) diff --git a/Source/wall.f90 b/Source/wall.f90 index bddea4ec28..ee0378c8a1 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -944,6 +944,8 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I ! store for use in divg B1%RHO_D_DZDN_F(N) = RHO_D_DZDN_OTHER ENDDO SPECIES_LOOP_2 + ! In DNS and LES, the diffusivity is not a constant due to the addition of the molecular value. + ! The following correction assures the sum of the diffusive fluxes is zero. IF (SIM_MODE==DNS_MODE .OR. SIM_MODE==LES_MODE) THEN N=MAXLOC(B1%ZZ_F(1:N_TRACKED_SPECIES),1) B1%RHO_D_DZDN_F(N) = -(SUM(B1%RHO_D_DZDN_F(1:N_TRACKED_SPECIES))-B1%RHO_D_DZDN_F(N))