Skip to content

Commit 66c41e1

Browse files
drjfloyddrjfloyd
committed
FDS Source: Fix Kn number calculation and associated verification Issue # 13926
1 parent 287cc5e commit 66c41e1

13 files changed

+3913
-54
lines changed

‎Manuals/FDS_Technical_Reference_Guide/Aerosol_Chapter.tex‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -111,11 +111,11 @@ \section{Aerosol Agglomeration}
111111
\be
112112
\frac{\d N_i}{\d t} = \sum_{k=1}^{M} \sum_{j=k}^{M} \left( 1- \frac{1}{2} \delta_{j,k} \right) \eta \Phi(j,k) N_j N_k - N_i \sum_{k=1}^M \Phi(j,k) N_k - R_i + S_i
113113
\ee
114-
where $\eta$, shown below, is a function for apportioning mass between adjacent particle bins.
114+
where $\eta$, shown below, is a function for apportioning the new particle mass, $m_{j+k}=x_j+x_k$, between adjacent particle bins.
115115
\be
116116
\eta = \begin{cases}
117-
\frac{x_{i+1}-m_i}{x_{i+1}-x_i} & x_i< m_i \leq x_{i+1} \\
118-
\frac{m_i-x_{i-1}}{x_i-x_{i-1}} & x_{i-1}< m_i \leq x_i
117+
\frac{x_{i+1}-m_{j+k}}{x_{i+1}-x_i} & x_i< m_{j+k} \leq x_{i+1} \\
118+
\frac{m_{j+k}-x_{i-1}}{x_i-x_{i-1}} & x_{i-1}< m_{j+k} \leq x_i
119119
\end{cases}
120120
\ee
121121
The agglomeration kernel is given by the sum of the Brownian (Br), gravitational (Gr), shear (Sh), and inertial (In) agglomeration terms.

‎Source/func.f90‎

Lines changed: 3850 additions & 1 deletion
Large diffs are not rendered by default.

‎Source/read.f90‎

Lines changed: 15 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1732,8 +1732,7 @@ SUBROUTINE READ_MISC
17321732
CFL_MAX,CFL_MIN,CFL_VELOCITY_NORM,CHECK_FO,CHECK_HT,CHECK_VN, &
17331733
CNF_CUTOFF,CONSTANT_SPECIFIC_HEAT_RATIO,&
17341734
C_SMAGORINSKY,C_VREMAN,C_WALE,DEPOSITION,EXTERNAL_FILENAME,&
1735-
FIXED_LES_FILTER_WIDTH,FLUX_LIMITER,FLUX_LIMITER_MW_CORRECTION,&
1736-
FREEZE_VELOCITY,FYI,GAMMA,GRAVITATIONAL_DEPOSITION,&
1735+
FIXED_LES_FILTER_WIDTH,FLUX_LIMITER,FREEZE_VELOCITY,FYI,GAMMA,GRAVITATIONAL_DEPOSITION,&
17371736
GRAVITATIONAL_SETTLING,GVEC,H_F_REFERENCE_TEMPERATURE,&
17381737
HUMIDITY,HVAC_LOCAL_PRESSURE,HVAC_MASS_TRANSPORT_CELL_L,HVAC_PRES_RELAX,HVAC_QFAN,IBLANK_SMV,I_MAX_TEMP,&
17391738
LES_FILTER_TYPE,LEVEL_SET_ELLIPSE,LEVEL_SET_MODE,&
@@ -1746,7 +1745,7 @@ SUBROUTINE READ_MISC
17461745
RAMP_UX,RAMP_UY,RAMP_UZ,RAMP_VX,RAMP_VY,RAMP_VZ,RAMP_WX,RAMP_WY,RAMP_WZ,&
17471746
RESTART,RESTART_CHID,SC,&
17481747
RND_SEED,SIMULATION_MODE,SMOKE3D_16,SMOKE_ALBEDO,SOLID_PHASE_ONLY,SOOT_DENSITY,SOOT_OXIDATION,&
1749-
TAU_DEFAULT,TENSOR_DIFFUSIVITY,TERRAIN_IMAGE,TEST_NEW_CHAR_MODEL,TEXTURE_ORIGIN,&
1748+
TAU_DEFAULT,TENSOR_DIFFUSIVITY,TERRAIN_IMAGE,FLUX_LIMITER_MW_CORRECTION,TEXTURE_ORIGIN,&
17501749
THERMOPHORETIC_DEPOSITION,THERMOPHORETIC_SETTLING,THICKEN_OBSTRUCTIONS,&
17511750
TMPA,TURBULENCE_MODEL,TURBULENT_DEPOSITION,UVW_FILE,&
17521751
VERBOSE,VISIBILITY_FACTOR,VN_MAX,VN_MIN,Y_CO2_INFTY,Y_O2_INFTY,&
@@ -7432,9 +7431,6 @@ SUBROUTINE PROC_MATL
74327431
H_ADJUST = ML%REFERENCE_ENTHALPY - ANS
74337432
ML%H = ML%H + H_ADJUST
74347433

7435-
! Store index of char material
7436-
IF (TRIM(ML%ID)=='CHAR') CHAR_INDEX=N
7437-
74387434
ENDDO PROC_MATL_LOOP
74397435

74407436
! Construct and solve linear system to adjust MATL enthalpies to conserve energy. For reaction i for MATL m the equation is
@@ -7678,6 +7674,7 @@ SUBROUTINE PROC_MATL
76787674
ENDIF
76797675
ENDDO
76807676

7677+
76817678
END SUBROUTINE PROC_MATL
76827679

76837680

@@ -9101,10 +9098,10 @@ END SUBROUTINE PROC_SURF_1
91019098

91029099
SUBROUTINE PROC_SURF_2
91039100
USE MATH_FUNCTIONS, ONLY: INTERPOLATE1D_UNIFORM
9104-
USE PHYSICAL_FUNCTIONS, ONLY: Q_REF_FIT
9101+
USE PHYSICAL_FUNCTIONS, ONLY: Q_REF_FIT,Q_REF_FIT_2
91059102
INTEGER :: ILPC,N,NN,NNN,NL,N_LIST,NLPC,NR,NS
91069103
LOGICAL :: INDEX_LIST(MAX_LPC)
9107-
REAL(EB) :: R_L(0:MAX_LAYERS), FUEL_MF,HRRPUA,E
9104+
REAL(EB) :: R_L(0:MAX_LAYERS), FUEL_MF,HRRPUA,E,MF_SUM,HR_SUM
91089105
REAL(EB), ALLOCATABLE, DIMENSION(:) :: Q_REF_INT
91099106
INTEGER :: I_GRAD,I_CONE_RAMP,E_PTR
91109107
INTEGER :: PROCESSED(N_RAMP)
@@ -9324,13 +9321,20 @@ SUBROUTINE PROC_SURF_2
93249321
ALLOCATE(SF%E2T_INDEX(SF%N_QDOTPP_REF))
93259322
SF%HOC_EFF = SF%HRRPUA / SUM(SF%MASS_FLUX)
93269323
SF%Y_S_EFF = 0._EB
9324+
SF%CHI_R_EFF = 0._EB
9325+
MF_SUM = TWO_EPSILON_EB
9326+
HR_SUM = TWO_EPSILON_EB
93279327
IF (SOOT_INDEX > 0) THEN
93289328
DO NR=1,N_REACTIONS
93299329
RN => REACTION(NR)
93309330
IF (SF%MASS_FLUX(RN%FUEL_SMIX_INDEX) <= 0._EB) CYCLE
93319331
SF%Y_S_EFF = SF%Y_S_EFF + RN%NU_SPECIES(SOOT_INDEX) * MW_SOOT / RN%MW_FUEL * SF%MASS_FLUX(RN%FUEL_SMIX_INDEX)
9332+
SF%CHI_R_EFF = SF%CHI_R_EFF + RN%CHI_R * SF%MASS_FLUX(RN%FUEL_SMIX_INDEX) * RN%HOC_COMPLETE
9333+
MF_SUM = MF_SUM + SF%MASS_FLUX(RN%FUEL_SMIX_INDEX)
9334+
HR_SUM = HR_SUM + SF%MASS_FLUX(RN%FUEL_SMIX_INDEX) * RN%HOC_COMPLETE
93329335
ENDDO
9333-
SF%Y_S_EFF = SF%Y_S_EFF/SUM(SF%MASS_FLUX)
9336+
SF%Y_S_EFF = SF%Y_S_EFF/MF_SUM
9337+
SF%CHI_R_EFF = SF%CHI_R_EFF/HR_SUM
93349338
ENDIF
93359339
N_QDOTPP_DO: DO NN=1,SF%N_QDOTPP_REF
93369340
IF (PROCESSED(SF%RAMP(TIME_HEAT-NN+1)%INDEX)>0) THEN
@@ -9361,7 +9365,8 @@ SUBROUTINE PROC_SURF_2
93619365
DO NNN=0,RP%NUMBER_INTERPOLATION_POINTS
93629366
IF (.NOT. SF%INERT_Q_REF) THEN
93639367
HRRPUA = SF%HRRPUA*RP%INTERPOLATED_DATA(NNN)
9364-
RP_QREF%INTERPOLATED_DATA(NNN) = Q_REF_FIT(HRRPUA,SF%HOC_EFF,SF%Y_S_EFF,SF%REFERENCE_HEAT_FLUX(NN))
9368+
!RP_QREF%INTERPOLATED_DATA(NNN) = Q_REF_FIT(HRRPUA,SF%HOC_EFF,SF%Y_S_EFF,SF%REFERENCE_HEAT_FLUX(NN))
9369+
RP_QREF%INTERPOLATED_DATA(NNN) = Q_REF_FIT_2(HRRPUA,SF%HOC_EFF,SF%Y_S_EFF,SF%CHI_R_EFF,SF%REFERENCE_HEAT_FLUX(NN))
93659370
ELSE
93669371
RP_QREF%INTERPOLATED_DATA(NNN) = SF%REFERENCE_HEAT_FLUX(NN)
93679372
ENDIF

‎Source/soot.f90‎

Lines changed: 25 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -67,9 +67,10 @@ SUBROUTINE SETTLING_VELOCITY(NM)
6767
PRES_G=0.5_EB*(PBAR(K,PRESSURE_ZONE(I,J,K))+PBAR(K,PRESSURE_ZONE(I+1,J,K)))
6868
DTDN = -(TMP(I+1,J,K)-TMP(I,J,K))*RDXN(I)
6969
CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP_G)
70-
KN_FAC = SQRT(0.5_EB*PI/(PRES_G*RHO_G))*MU_G
70+
! Kn=2 lambda/d, lambda has sqrt(1/2). kn_fac has 2*sqrt(1/2)=sqrt(4/2)=sqrt(2)
71+
KN_FAC = SQRT(2._EB*PI/(PRES_G*RHO_G))*MU_G
7172
IF (THERMOPHORETIC_SETTLING) THEN
72-
KN = KN_FAC/SPECIES_MIXTURE(N)%THERMOPHORETIC_DIAMETER
73+
KN = KN_FAC/SPECIES_MIXTURE(N)%THERMOPHORETIC_DIAMETER
7374
CN = CUNNINGHAM(KN)
7475
CALL GET_CONDUCTIVITY(ZZ_GET,K_G,TMP_G)
7576
ALPHA = K_G/SPECIES_MIXTURE(N)%CONDUCTIVITY_SOLID
@@ -89,7 +90,7 @@ SUBROUTINE SETTLING_VELOCITY(NM)
8990
PRES_G=0.5_EB*(PBAR(K,PRESSURE_ZONE(I,J,K))+PBAR(K,PRESSURE_ZONE(I,J+1,K)))
9091
DTDN = -(TMP(I,J+1,K)-TMP(I,J,K))*RDYN(J)
9192
CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP_G)
92-
KN_FAC = SQRT(0.5_EB*PI/(PRES_G*RHO_G))*MU_G
93+
KN_FAC = SQRT(2._EB*PI/(PRES_G*RHO_G))*MU_G
9394
IF (THERMOPHORETIC_SETTLING) THEN
9495
KN = KN_FAC/SPECIES_MIXTURE(N)%THERMOPHORETIC_DIAMETER
9596
CN = CUNNINGHAM(KN)
@@ -111,7 +112,7 @@ SUBROUTINE SETTLING_VELOCITY(NM)
111112
PRES_G=0.5_EB*(PBAR(K,PRESSURE_ZONE(I,J,K))+PBAR(K,PRESSURE_ZONE(I,J,K+1)))
112113
DTDN = -(TMP(I,J,K+1)-TMP(I,J,K))*RDZN(K)
113114
CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP_G)
114-
KN_FAC = SQRT(0.5_EB*PI/(PRES_G*RHO_G))*MU_G
115+
KN_FAC = SQRT(2._EB*PI/(PRES_G*RHO_G))*MU_G
115116
IF (THERMOPHORETIC_SETTLING) THEN
116117
KN = KN_FAC/SPECIES_MIXTURE(N)%THERMOPHORETIC_DIAMETER
117118
CN = CUNNINGHAM(KN)
@@ -177,9 +178,9 @@ END SUBROUTINE SETTLING_VELOCITY
177178

178179

179180
SUBROUTINE INITIALIZE_AGGLOMERATION
181+
USE GLOBAL_CONSTANTS, ONLY: GRAVITATIONAL_SETTLING
180182
INTEGER :: N,I,II,III
181183
REAL(EB) :: E_PK,MIN_PARTICLE_MASS,MAX_PARTICLE_MASS
182-
183184
ALLOCATE(BIN_S(N_AGGLOMERATION_SPECIES))
184185
ALLOCATE(BIN_M(N_AGGLOMERATION_SPECIES,0:MAXVAL(N_PARTICLE_BINS)))
185186
ALLOCATE(BIN_X(N_AGGLOMERATION_SPECIES,1:MAXVAL(N_PARTICLE_BINS)))
@@ -211,8 +212,6 @@ SUBROUTINE INITIALIZE_AGGLOMERATION
211212
BIN_X(N,I) = 2._EB*BIN_M(N,I)/(1._EB+BIN_S(N))
212213
ENDDO
213214

214-
215-
216215
DO I=1,N_PARTICLE_BINS(N)
217216
PARTICLE_RADIUS(N,I) = (BIN_X(N,I) / FOTHPI / SPECIES(AGGLOMERATION_SPEC_INDEX(N))%DENSITY_SOLID)**ONTH
218217
MOBILITY_FAC(N,I) = 1._EB/(6._EB*PI*PARTICLE_RADIUS(N,I))
@@ -222,7 +221,11 @@ SUBROUTINE INITIALIZE_AGGLOMERATION
222221
DO II=1,N_PARTICLE_BINS(N)
223222
PARTICLE_MASS(N,I,II) = BIN_X(N,I) + BIN_X(N,II)
224223
E_PK = MIN(PARTICLE_RADIUS(N,I),PARTICLE_RADIUS(N,II))**2/(2._EB*(PARTICLE_RADIUS(N,I)+PARTICLE_RADIUS(N,II))**2)
225-
PHI_G_FAC(N,I,II) = E_PK*(PARTICLE_RADIUS(N,I)+PARTICLE_RADIUS(N,II))**2
224+
IF (GRAVITATIONAL_SETTLING) THEN
225+
PHI_G_FAC(N,I,II) = E_PK*(PARTICLE_RADIUS(N,I)+PARTICLE_RADIUS(N,II))**2
226+
ELSE
227+
PHI_G_FAC(N,I,II) = 0._EB
228+
ENDIF
226229
PHI_B_FAC(N,I,II)= 4._EB*PI*K_BOLTZMANN*(PARTICLE_RADIUS(N,I)+PARTICLE_RADIUS(N,II))
227230
!Check what Re and r are for PHI_I and _S
228231
PHI_S_FAC(N,I,II) = E_PK*(PARTICLE_RADIUS(N,I)+PARTICLE_RADIUS(N,II))**3*SQRT(8._EB*PI/15._EB)
@@ -234,21 +237,21 @@ SUBROUTINE INITIALIZE_AGGLOMERATION
234237
BINDO:DO III=2,N_PARTICLE_BINS(N)
235238
IF (PARTICLE_MASS(N,I,II) > BIN_X(N,N_PARTICLE_BINS(N))) THEN
236239
BIN_ETA_INDEX(N,I,II,:) = N_PARTICLE_BINS(N)
237-
BIN_ETA(N,I,II,:) = 0.5_EB*BIN_X(N,N_PARTICLE_BINS(N))/PARTICLE_MASS(N,I,II)
240+
BIN_ETA(N,I,II,:) = 0.5_EB*PARTICLE_MASS(N,I,II)/BIN_X(N,N_PARTICLE_BINS(N))
238241
EXIT BINDO
239242
ELSE
240243
IF (PARTICLE_MASS(N,I,II) > BIN_X(N,III-1) .AND. PARTICLE_MASS(N,I,II) < BIN_X(N,III)) THEN
241244
BIN_ETA_INDEX(N,I,II,1) = III-1
242245
BIN_ETA(N,I,II,1) = (BIN_X(N,III)-PARTICLE_MASS(N,I,II))/(BIN_X(N,III)-BIN_X(N,III-1))
243246
BIN_ETA_INDEX(N,I,II,2) = III
244247
BIN_ETA(N,I,II,2) = (PARTICLE_MASS(N,I,II)-BIN_X(N,III-1))/(BIN_X(N,III)-BIN_X(N,III-1))
245-
IF (I==II) BIN_ETA(N,I,II,:) = BIN_ETA(N,I,II,:) *0.5_EB
246248
EXIT BINDO
247249
ENDIF
248250
ENDIF
249251
ENDDO BINDO
250252
ENDDO
251253
ENDDO
254+
252255
BIN_ETA(N,N_PARTICLE_BINS(N),N_PARTICLE_BINS(N),1) = 1._EB
253256
BIN_ETA_INDEX(N,N_PARTICLE_BINS(N),N_PARTICLE_BINS(N),1) = N_PARTICLE_BINS(N)
254257
BIN_ETA(N,N_PARTICLE_BINS(N),N_PARTICLE_BINS(N),2) = 0._EB
@@ -265,8 +268,8 @@ SUBROUTINE CALC_AGGLOMERATION(DT,NM)
265268
INTEGER, INTENT(IN) :: NM
266269
REAL(EB), INTENT(IN) :: DT
267270
REAL(EB) :: DUDX,DVDY,DWDZ,ONTHDIV,S11,S22,S33,DUDY,DUDZ,DVDX,DVDZ,DWDX,DWDY,S12,S23,S13,STRAIN_RATE,DISSIPATION_RATE
268-
REAL(EB) :: KN,MFP,N_I(MAXVAL(N_PARTICLE_BINS)),N0(MAXVAL(N_PARTICLE_BINS)),N1(MAXVAL(N_PARTICLE_BINS)),&
269-
N2(MAXVAL(N_PARTICLE_BINS)),N3(MAXVAL(N_PARTICLE_BINS)),&
271+
REAL(EB) :: KN,KN_FAC,N_I(MAXVAL(N_PARTICLE_BINS)),N0(MAXVAL(N_PARTICLE_BINS)),N1(MAXVAL(N_PARTICLE_BINS)),&
272+
N2(MAXVAL(N_PARTICLE_BINS)),N3(MAXVAL(N_PARTICLE_BINS)),DN,&
270273
RHO_G,TMP_G,MUG,TERMINAL(MAXVAL(N_PARTICLE_BINS)),&
271274
FU,MOBILITY(MAXVAL(N_PARTICLE_BINS)),ZZ_GET(1:N_TRACKED_SPECIES),AM,AMT(MAXVAL(N_PARTICLE_BINS)),&
272275
PHI_B,PHI_S,PHI_G,PHI_I,PHI(MAXVAL(N_PARTICLE_BINS),MAXVAL(N_PARTICLE_BINS)),VREL,FU1,FU2,DT_SUBSTEP,DT_SUM,TOL
@@ -287,7 +290,7 @@ SUBROUTINE CALC_AGGLOMERATION(DT,NM)
287290
ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(I,J,K,1:N_TRACKED_SPECIES)
288291
TMP_G = TMP(I,J,K)
289292
CALL GET_VISCOSITY(ZZ_GET,MUG,TMP_G)
290-
MFP = MUG*SQRT(0.5_EB*PI/(PBAR(K,PRESSURE_ZONE(I,J,K))*RHO_G))
293+
KN_FAC = MUG*SQRT(PI/(2._EB*PBAR(K,PRESSURE_ZONE(I,J,K))*RHO_G))
291294
IM1 = MAX(0,I-1)
292295
JM1 = MAX(0,J-1)
293296
KM1 = MAX(0,K-1)
@@ -316,8 +319,9 @@ SUBROUTINE CALC_AGGLOMERATION(DT,NM)
316319
STRAIN_RATE = 2._EB*(S11**2 + S22**2 + S33**2 + 2._EB*(S12**2 + S13**2 + S23**2))
317320
DISSIPATION_RATE = MU(I,J,K)/RHO_G*STRAIN_RATE
318321
N_I = N0
322+
319323
DO N=1,N_PARTICLE_BINS(NS)
320-
KN=0.5_EB*MFP/PARTICLE_RADIUS(NS,N)
324+
KN=KN_FAC/PARTICLE_RADIUS(NS,N)
321325
!Verify CN
322326
MOBILITY(N) = CUNNINGHAM(KN)*MOBILITY_FAC(NS,N)/MUG
323327
TERMINAL(N) = MOBILITY(N)*GRAV*BIN_X(NS,N)
@@ -326,6 +330,7 @@ SUBROUTINE CALC_AGGLOMERATION(DT,NM)
326330
(3._EB*PARTICLE_RADIUS(NS,N)*AM)-PARTICLE_RADIUS(NS,N)
327331
ENDDO
328332
PHI = 0._EB
333+
329334
DO N=1,N_PARTICLE_BINS(NS)
330335
DO NN=1,N_PARTICLE_BINS(NS)
331336
IF (NN<N) CYCLE
@@ -356,15 +361,15 @@ SUBROUTINE CALC_AGGLOMERATION(DT,NM)
356361
DO NN=N,N_PARTICLE_BINS(NS)
357362
IF (N0(N)<MIN_AGGLOMERATION .OR. N0(NN)<MIN_AGGLOMERATION) CYCLE
358363
!Remove particles that agglomerate
359-
N1(N)=N1(N)-PHI(NN,N)*N0(N)*N0(NN)*DT_SUBSTEP
360-
IF (NN/=N) N1(NN)=N1(NN)-PHI(NN,N)*N0(N)*N0(NN)*DT_SUBSTEP
364+
DN = PHI(NN,N)*N0(N)*N0(NN)*DT_SUBSTEP
365+
N1(N) = N1(N) - DN
366+
N1(NN) = N1(NN) - DN
361367
! Create new particles from agglomeration
362-
N2(BIN_ETA_INDEX(NS,N,NN,1)) = N2(BIN_ETA_INDEX(NS,N,NN,1)) + &
363-
BIN_ETA(NS,N,NN,1)*PHI(N,NN)*N0(N)*N0(NN)*DT_SUBSTEP
364-
N2(BIN_ETA_INDEX(NS,N,NN,2)) = N2(BIN_ETA_INDEX(NS,N,NN,2)) + &
365-
BIN_ETA(NS,N,NN,2)*PHI(N,NN)*N0(N)*N0(NN)*DT_SUBSTEP
368+
N2(BIN_ETA_INDEX(NS,N,NN,1)) = N2(BIN_ETA_INDEX(NS,N,NN,1)) + BIN_ETA(NS,N,NN,1) * DN
369+
N2(BIN_ETA_INDEX(NS,N,NN,2)) = N2(BIN_ETA_INDEX(NS,N,NN,2)) + BIN_ETA(NS,N,NN,2) * DN
366370
ENDDO
367371
ENDDO AGGLOMERATE_LOOP
372+
368373
N3 = N1 + N2
369374
N3 = N3 + N0
370375
TOL = MAXVAL((N0-N3)/(N0+TINY_EB))

‎Source/type.f90‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -907,6 +907,7 @@ MODULE TYPES
907907
REAL(EB) :: M_DOT_G_PP_ADJUST_FAC=1._EB !< For MLRPUA, scales the gas production if gas H_o_C /= solid H_o_C
908908
REAL(EB) :: HOC_EFF !< Effective heat of combustion for S_pyro
909909
REAL(EB) :: Y_S_EFF !< Effective soot yield for S_pyro
910+
REAL(EB) :: CHI_R_EFF !< Effective radiative fraction for S_pyro
910911
REAL(EB) :: TIME_STEP_FACTOR=10._EB !< Maximum amount to reduce solid phase conduction time step
911912
REAL(EB) :: REMESH_RATIO=0.05 !< Fraction change in wall node DX to trigger a remesh
912913

‎Source/wall.f90‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1160,8 +1160,8 @@ SUBROUTINE CALCULATE_ZZ_F(T,DT,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX)
11601160
DT_SPYRO(1:SF%N_THICK_REF) = DT * SF%SPYRO_TH_FACTOR(1:SF%N_THICK_REF)
11611161
TSI = MIN(T-B1%T_IGN+DT, SF%REFERENCE_HEAT_FLUX_TIME_INTERVAL+DT)
11621162
IF (SOLID_PHASE_ONLY .AND. .NOT. SF%INERT_Q_REF) THEN
1163-
B1%Q_IN_SMOOTH = (B1%Q_IN_SMOOTH *(TSI-DT) + &
1164-
DT*Q_REF_FIT(SUM(B1%M_DOT_G_PP_ACTUAL)*SF%HOC_EFF,SF%HOC_EFF,SF%Y_S_EFF,B1%Q_RAD_IN/B1%EMISSIVITY))/TSI
1163+
B1%Q_IN_SMOOTH = (B1%Q_IN_SMOOTH *(TSI-DT) + DT * &
1164+
Q_REF_FIT(SUM(B1%M_DOT_G_PP_ACTUAL)*SF%HOC_EFF,SF%HOC_EFF,SF%Y_S_EFF,B1%Q_RAD_IN/B1%EMISSIVITY))/TSI
11651165
ELSE
11661166
IF (B1%EMISSIVITY > 0._EB) THEN
11671167
B1%Q_IN_SMOOTH = (B1%Q_IN_SMOOTH *(TSI-DT) + DT*(B1%Q_CON_F+B1%Q_RAD_IN/B1%EMISSIVITY))/TSI
@@ -1623,7 +1623,8 @@ SUBROUTINE CALC_DEPOSITION(DT,BC,B1,B2,WALL_INDEX,CFACE_INDEX)
16231623
TMP_FILM = 0.5_EB*(B1%TMP_G+B1%TMP_F)
16241624
CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP_FILM)
16251625
CALL GET_CONDUCTIVITY(ZZ_GET,K_G,TMP_FILM)
1626-
KN_FAC = MU_G*SQRT(0.5_EB*PI/(PBAR(BC%KKG,PRESSURE_ZONE(BC%IIG,BC%JJG,BC%KKG))*B1%RHO_G))
1626+
! Kn=2 lambda/d, lambda has sqrt(1/2). kn_fac has 2*sqrt(1/2)=sqrt(4/2)=sqrt(2)
1627+
KN_FAC = MU_G*SQRT(2._EB*PI/(PBAR(BC%KKG,PRESSURE_ZONE(BC%IIG,BC%JJG,BC%KKG))*B1%RHO_G))
16271628
ALPHA = K_G/SM%CONDUCTIVITY_SOLID
16281629
DTMPDX = B1%HEAT_TRANS_COEF*(B1%TMP_G-B1%TMP_F)/K_G
16291630
U_THERM = 0._EB
@@ -1634,7 +1635,6 @@ SUBROUTINE CALC_DEPOSITION(DT,BC,B1,B2,WALL_INDEX,CFACE_INDEX)
16341635
KN = KN_FAC/SM%THERMOPHORETIC_DIAMETER
16351636
U_THERM = CS2*(ALPHA+CT*KN)*CUNNINGHAM(KN)/((1._EB+CM3*KN)*(1+2*ALPHA+CT2*KN)) * MU_G/(B1%TMP_G*B1%RHO_G)*DTMPDX
16361637
ENDIF
1637-
16381638
IF (GRAVITATIONAL_DEPOSITION) THEN
16391639
KN = KN_FAC/SM%MEAN_DIAMETER
16401640
U_GRAV = - DOT_PRODUCT(GVEC,BC%NVEC)*CUNNINGHAM(KN)*SM%MEAN_DIAMETER**2*SM%DENSITY_SOLID/(18._EB*MU_G)
@@ -3894,4 +3894,3 @@ SUBROUTINE TGA_ANALYSIS(NM)
38943894
END SUBROUTINE TGA_ANALYSIS
38953895

38963896
END MODULE WALL_ROUTINES
3897-

‎Verification/Aerosols/aerosol_agglomeration.csv‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,3 @@
11
Time,Exact Total,Exact Bin 1,Exact Bin 2
22
0,0.00001,0.00001,0
3-
60,0.00001,9.9989E-06,1.06E-09
3+
60,0.00001,1.00E-05,2.26E-09

‎Verification/Aerosols/aerosol_agglomeration.fds‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22

33
&MESH IJK=10,10,10, XB=0.0,0.01,0.0,0.01,0.00,0.01 /
44

5-
&TIME T_END=60., DT=0.01/
5+
&TIME T_END=60.,DT=0.004/
66

77
&DUMP FLUSH_FILE_BUFFERS=T, NFRAMES=20 /
88

@@ -13,7 +13,7 @@
1313
TURBULENT_DEPOSITION =.FALSE.
1414
STRATIFICATION =.FALSE.
1515
NOISE =.FALSE./
16-
16+
1717
&RADI RADIATION=.FALSE./
1818

1919
&SPEC ID='AIR',MW=28.8,CONDUCTIVITY=0.025,VISCOSITY=2.E-5,SPECIFIC_HEAT=1.,BACKGROUND=.TRUE./

‎Verification/Aerosols/aerosol_agglomeration_2.csv‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,3 @@
11
Time,Exact Total,Exact Bin 1,Exact Bin 2
22
0,1.00E-05,1.00E-05,0.00E+00
3-
60,1.00E-05,9.9989E-06,1.08E-09
3+
60,1.00E-05,1.00E-05,2.30E-09

‎Verification/Aerosols/aerosol_agglomeration_2.fds‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22

33
&MESH IJK=10,10,10, XB=0.0,0.01,0.0,0.01,0.00,0.01 /
44

5-
&TIME T_END=60./
5+
&TIME T_END=60.,DT=0.004/
66

77
&DUMP FLUSH_FILE_BUFFERS=T, NFRAMES=20 /
88

@@ -14,7 +14,7 @@
1414
STRATIFICATION =.FALSE.
1515
NOISE =.FALSE.
1616
SIMULATION_MODE='DNS'/
17-
17+
1818
&RADI RADIATION=.FALSE./
1919

2020
&SPEC ID='AIR',MW=28.8,CONDUCTIVITY=0.025,VISCOSITY=2.E-5,SPECIFIC_HEAT=1.,BACKGROUND=.TRUE./

0 commit comments

Comments
 (0)