Merge pull request #12484 from drjfloyd/master

FDS Source: Issue #12481 Use Ferrari's method for the wall temperatu…
firemodels · Feb 15, 2024 · 81f1675 · 81f1675
2 parents 515d6e5 + 2e34eaf
commit 81f1675
Show file tree

Hide file tree

Showing 2 changed files with 21 additions and 18 deletions.
diff --git a/Source/part.f90 b/Source/part.f90
@@ -3927,7 +3927,7 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM)
          B2%LP_TEMP(LPC%ARRAY_INDEX) = B2%LP_TEMP(LPC%ARRAY_INDEX) + A_DROP*0.5_EB*(TMP_DROP+TMP_DROP_NEW)
          B2%LP_CPUA(LPC%ARRAY_INDEX) = B2%LP_CPUA(LPC%ARRAY_INDEX) + &
                                        (1._EB-LPC%RUNNING_AVERAGE_FACTOR_WALL)*WGT*Q_CON_WALL/(B1%AREA*DT)
-         B1%Q_RAD_IN = (B1%AREA*DT*B1%Q_RAD_IN - WGT*DT*Q_DOT_RAD) / (B1%AREA*DT)
+         B1%Q_RAD_IN = (B1%AREA*DT*B1%Q_RAD_IN - WGT*DT*Q_DOT_RAD*B1%Q_RAD_IN/(B1%Q_RAD_IN+B1%Q_RAD_OUT)) / (B1%AREA*DT)
       ENDIF
 
    ENDDO PARTICLE_LOOP

diff --git a/Source/wall.f90 b/Source/wall.f90
@@ -267,7 +267,8 @@ SUBROUTINE CALCULATE_TMP_F(WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX)
             RHO_ZZ_F(1:N_TOTAL_SCALARS),ZZ_GET(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
+            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
 
 LOGICAL :: SECOND_ORDER_INTERPOLATED_BOUNDARY,SOLID_OTHER,ATMOSPHERIC_INTERPOLATION,CC_SOLID_FLAG
 INTEGER :: II,JJ,KK,IIG,JJG,KKG,IOR,IIO,JJO,KKO,N,ADCOUNT,ICG,ICO
@@ -388,11 +389,11 @@ SUBROUTINE CALCULATE_TMP_F(WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX)
       ELSE
          TSI = T - B1%T_IGN
       ENDIF
-      TMP_OTHER = B1%TMP_F
+      HTC_OLD = B1%HEAT_TRANS_COEF
       RAMP_FACTOR = EVALUATE_RAMP(TSI,SF%RAMP(TIME_HEAT)%INDEX,TAU=SF%RAMP(TIME_HEAT)%TAU)
       IF (SF%SET_H) THEN
          B1%Q_CON_F = -RAMP_FACTOR*SF%CONVECTIVE_HEAT_FLUX*B1%AREA_ADJUST
-         B1%HEAT_TRANS_COEF = B1%Q_CON_F/(B1%TMP_G-TMP_OTHER+TWO_EPSILON_EB)
+         B1%HEAT_TRANS_COEF = B1%Q_CON_F/(B1%TMP_G-B1%TMP_F+TWO_EPSILON_EB)
       ELSE
          IF (SF%THERMAL_BC_INDEX==NET_FLUX_BC) THEN
             QNET = -RAMP_FACTOR*SF%NET_HEAT_FLUX*B1%AREA_ADJUST
@@ -410,7 +411,7 @@ SUBROUTINE CALCULATE_TMP_F(WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX)
          ADCOUNT = 0
          ADLOOP: DO
             ADCOUNT = ADCOUNT + 1
-            DTMP = B1%TMP_G - TMP_OTHER
+            DTMP = B1%TMP_G - B1%TMP_F
             IF (ABS(QNET) > 0._EB .AND. ABS(DTMP) <TWO_EPSILON_EB) DTMP=1._EB
             IF (PRESENT(WALL_INDEX)) THEN
                B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC,SURF_INDEX,WALL_INDEX_IN=WALL_INDEX)
@@ -419,22 +420,24 @@ SUBROUTINE CALCULATE_TMP_F(WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX)
             ELSEIF (PRESENT(PARTICLE_INDEX)) THEN
                B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC,SURF_INDEX,PARTICLE_INDEX_IN=PARTICLE_INDEX)
             ENDIF
-            IF (RADIATION .AND. SF%THERMAL_BC_INDEX/=CONVECTIVE_FLUX_BC) THEN
-               QEXTRA =  B1%HEAT_TRANS_COEF*DTMP + B1%Q_RAD_IN - B1%EMISSIVITY * SIGMA * TMP_OTHER ** 4 - QNET
-               FDERIV = -B1%HEAT_TRANS_COEF -  4._EB * B1%EMISSIVITY * SIGMA * TMP_OTHER ** 3
+            ! Use Ferrari's method
+            BBB = B1%HEAT_TRANS_COEF / (B1%EMISSIVITY * SIGMA)
+            CCC = -(B1%Q_RAD_IN+B1%HEAT_TRANS_COEF*B1%TMP_G-QNET)/(B1%EMISSIVITY * SIGMA)
+            PPP = -CCC
+            QQQ = -0.125_EB*BBB**2
+            RRR = -0.5_EB*QQQ+SQRT(0.25_EB*QQQ**2+PPP**3/27._EB)
+            UUU = RRR**ONTH
+            IF (UUU < TWO_EPSILON_EB) THEN
+               YYY = -QQQ**ONTH
             ELSE
-               QEXTRA = B1%HEAT_TRANS_COEF*DTMP - QNET
-               FDERIV = -B1%HEAT_TRANS_COEF
-            ENDIF
-            IF (ABS(FDERIV) > TWO_EPSILON_EB) TMP_OTHER = TMP_OTHER - QEXTRA / FDERIV
-            IF (ABS(QEXTRA) < 1.E-10_EB .OR. ADCOUNT > 20) THEN
-               B1%TMP_F = MAX(TMPMIN,MIN(TMPMAX,TMP_OTHER))
-               EXIT ADLOOP
-            ELSE
-               B1%TMP_F = MAX(TMPMIN,MIN(TMPMAX,TMP_OTHER))
-               CYCLE ADLOOP
+               YYY = UUU-ONTH*PPP/UUU
             ENDIF
+            WWW = SQRT(2._EB*YYY)
+            B1%TMP_F = 0.5_EB*(-WWW+SQRT(-2._EB*(YYY-BBB/WWW)))
+            HTC_OLD = 0.2_EB*HTC_OLD+0.8_EB*B1%HEAT_TRANS_COEF
+            IF (ABS(HTC_OLD-B1%HEAT_TRANS_COEF) < 1.E-6_EB .OR. ADCOUNT > 20) EXIT ADLOOP
          ENDDO ADLOOP
+         B1%HEAT_TRANS_COEF = HTC_OLD
          B1%Q_CON_F = B1%HEAT_TRANS_COEF*DTMP
       ENDIF
       IF (RADIATION) B1%Q_RAD_OUT = SIGMA*B1%EMISSIVITY*B1%TMP_F**4