From 549044f4082f59416d69d7f62ff8a43b58c928c9 Mon Sep 17 00:00:00 2001 From: mcgratta Date: Fri, 25 Aug 2023 17:04:12 -0400 Subject: [PATCH] FDS Source: Introduce the concept of a LINING for HT3D --- .../FDS_Verification_Guide.tex | 7 +- Source/init.f90 | 102 +++++++++++++++--- Source/read.f90 | 5 +- Source/type.f90 | 1 + .../FDS_verification_dataplot_inputs.csv | 1 + Verification/FDS_Cases.sh | 1 + .../ht3d_energy_conservation_6.fds | 52 +++++++++ 7 files changed, 149 insertions(+), 20 deletions(-) create mode 100644 Verification/Heat_Transfer/ht3d_energy_conservation_6.fds diff --git a/Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex b/Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex index 67509b5e248..3f77296d6da 100644 --- a/Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex +++ b/Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex @@ -4424,14 +4424,15 @@ \subsection{Energy Conservation in a 3-D Solid (\texorpdfstring{\textct{ht3d\_en \label{fig:ht3d_energy_conservation_2} \end{figure} -Figure~\ref{fig:ht3d_energy_conservation_4} compares the integrated net heat flux versus the internal enthalpy in a case where a small hot block is placed above a larger block for which 3-D heat conduction is invoked. In one case (left plot), the block is made of one material with a thermal conductivity of 0.1~W/m/K, typical of insulation. The internal cells cluster near all six faces of the block as would be done in a 1-D case. The calculation of the internal enthalpy can be performed in any coordinate direction, both positive and negative. In this case, the positive $x$ direction is chosen, whereas the hot block is placed above the larger block. This ensures that the energy is properly transferred in the lateral directions. In a second case (right plot), the solid block consists of a thin, 1~cm thick layer of insulation over a solid block of steel. +Figure~\ref{fig:ht3d_energy_conservation_4} compares the integrated net heat flux versus the internal enthalpy in a case where a small hot block is placed above a larger block for which 3-D heat conduction is invoked. In one case (upper left plot), the block is made of one material with a thermal conductivity of 0.1~W/m/K, typical of insulation. The internal cells cluster near all six faces of the block as would be done in a 1-D case. The calculation of the internal enthalpy can be performed in any coordinate direction, both positive and negative. In this case, the positive $x$ direction is chosen, whereas the hot block is placed above the larger block. This ensures that the energy is properly transferred in the lateral directions. In a second case (upper right plot), the solid block consists of a thin, 1~cm thick layer of insulation over a solid block of steel. The insulation layer is specified using a one-cell thick obstruction on top of the obstruction representing the steel. In a third case (lower plot), the solid block of steel is covered by two layers of insulation of thickness 1.5~cm and 1.0~cm that are specified using the conventional {\ct SURF} line with arrays of {\ct MATL\_ID} and {\ct THICKNESS}. Each layer of insulation is composed of two different material components. \begin{figure}[ht] \begin{tabular*}{\textwidth}{l@{\extracolsep{\fill}}r} \includegraphics[height=2.2in]{SCRIPT_FIGURES/ht3d_energy_conservation_4} & -\includegraphics[height=2.2in]{SCRIPT_FIGURES/ht3d_energy_conservation_5} +\includegraphics[height=2.2in]{SCRIPT_FIGURES/ht3d_energy_conservation_5} \\ +\multicolumn{2}{c}{\includegraphics[height=2.2in]{SCRIPT_FIGURES/ht3d_energy_conservation_6}} \end{tabular*} -\caption[Additional \textct{ht3d\_energy\_conservation} test cases, 4 and 5]{Comparison of the integrated net heat flux versus the internal enthalpy of a homogenous solid block with relativley low thermal conductivity (left), and one with multiple layers (right).} +\caption[Additional \textct{ht3d\_energy\_conservation} test cases, 4, 5, and 6]{Comparison of the integrated net heat flux versus the internal enthalpy for a homogenous solid block of insulation material (upper left), a block of steel with a single layer of insulation (upper right), and a block with multiple layers of multi-component insulation (bottom).} \label{fig:ht3d_energy_conservation_4} \end{figure} diff --git a/Source/init.f90 b/Source/init.f90 index e6ac13fbf7e..35200acc222 100644 --- a/Source/init.f90 +++ b/Source/init.f90 @@ -998,10 +998,11 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_2(NM) INTEGER :: N,I,J,K,II,JJ,KK,IPTS,JPTS,KPTS,N_EDGES_DIM,IW,IC,ICG,IOR,IERR,IPZ,NOM,ITER,IZERO,IIG,JJG,KKG,ICF,NSLICE,ITW,IEC,ITW2 INTEGER, INTENT(IN) :: NM REAL(EB) :: ZZ_GET(1:N_TRACKED_SPECIES),VC,RTRM,CP,XXC,YYC,ZZC -INTEGER :: IBP1,JBP1,KBP1,IBAR,JBAR,KBAR,OBST_INDEX,OBST_INDEX_PREVIOUS,NN,NNN,NL,MATL_INDEX(MAX_MATERIALS),MI -REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: MATL_MASS_FRACTION -REAL(EB), DIMENSION(MAX_LAYERS) :: LAYER_THICKNESS -REAL(EB) :: XS,XF,YS,YF,ZS,ZF,THICKNESS,OLD_THICKNESS +INTEGER :: IBP1,JBP1,KBP1,IBAR,JBAR,KBAR,OBST_INDEX,OBST_INDEX_PREVIOUS,NN,NNN,NL +INTEGER, DIMENSION(MAX_MATERIALS) :: MATL_INDEX,MATL_INDEX_NEW +REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: MATL_MASS_FRACTION,MATL_MASS_FRACTION_NEW +REAL(EB), DIMENSION(MAX_LAYERS) :: LAYER_THICKNESS,LAYER_THICKNESS_NEW +REAL(EB) :: XS,XF,YS,YF,ZS,ZF,THICKNESS,OLD_THICKNESS,TOTAL_THICKNESS,UPPER_BOUND,LOWER_BOUND,LINING_THICKNESS,BACK_LINING_THICKNESS TYPE (MESH_TYPE), POINTER :: M,M4 TYPE (OMESH_TYPE), POINTER :: M3 TYPE (WALL_TYPE), POINTER :: WC @@ -1013,10 +1014,10 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_2(NM) TYPE (MESH_TYPE), POINTER :: OM,OM_PREV TYPE (VENTS_TYPE), POINTER :: VT TYPE (OBSTRUCTION_TYPE), POINTER :: OB,OB_PREV -TYPE (SURFACE_TYPE), POINTER :: SF +TYPE (SURFACE_TYPE), POINTER :: SF,SF_BACK TYPE (STORAGE_TYPE), POINTER :: OS LOGICAL :: SOLID_CELL -INTEGER :: WSC,WSCS,N_LAYERS,N_MATLS +INTEGER :: WSC,WSCS,N_LAYERS,N_MATLS,N_LAYERS_NEW,N_MATLS_NEW TYPE WALL_SAVE_TYPE INTEGER :: COUNTER=0 INTEGER :: START=1 @@ -1252,15 +1253,7 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_2(NM) LAYER_THICKNESS(N_LAYERS) = LAYER_THICKNESS(N_LAYERS) + THICKNESS - OLD_THICKNESS - MATL_LOOP: DO NN=1,MAX_MATERIALS - MI = OB%MATL_INDEX(NN) - IF (MI<1) EXIT - DO NNN=1,N_MATLS - IF (MI==MATL_INDEX(NNN)) EXIT MATL_LOOP - ENDDO - N_MATLS = N_MATLS + 1 - MATL_INDEX(N_MATLS) = MI - ENDDO MATL_LOOP + CALL ADD_MATERIAL(N_MATLS,OB%MATL_INDEX,MATL_INDEX) IF (OBST_INDEX/=OBST_INDEX_PREVIOUS .AND. OBST_INDEX_PREVIOUS>0 .AND. OBST_INDEX>0) THEN OB_PREV => OM_PREV%OBSTRUCTION(OBST_INDEX_PREVIOUS) @@ -1305,6 +1298,66 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_2(NM) ENDDO FIND_BACK_WALL_CELL + ! If the user has specified LINING materials (HT3D or HT1D with SURF MATLs and THICKNESS), add this information to + ! existing lists of layers and materials. + ! The new arrays of thickness and material information are temporarily stored in arrays with _NEW suffix. + + IF (SF%LINING) THEN + SF_BACK => SURFACE(WC%BACK_SURF) + LINING_THICKNESS = SUM(SF%LAYER_THICKNESS(1:SF%N_LAYERS)) + BACK_LINING_THICKNESS = SUM(SF_BACK%LAYER_THICKNESS(1:SF_BACK%N_LAYERS)) + ! Copy the front face SURF layers into the NEW holding arrays + N_LAYERS_NEW = SF%N_LAYERS + LAYER_THICKNESS_NEW(1:N_LAYERS_NEW) = SF%LAYER_THICKNESS(1:SF%N_LAYERS) + N_MATLS_NEW = N_MATLS + MATL_MASS_FRACTION_NEW = 0._EB + MATL_INDEX_NEW = -1 + MATL_INDEX_NEW(1:N_MATLS_NEW) = MATL_INDEX(1:N_MATLS) ! MATL_INDEX is taken from the OBSTs that make up the solid + CALL ADD_MATERIAL(N_MATLS_NEW,SF%MATL_INDEX,MATL_INDEX_NEW) ! Add new materials from the front surface lining + CALL ADD_MATERIAL(N_MATLS_NEW,SF_BACK%MATL_INDEX,MATL_INDEX_NEW) ! Add new materials from the back surface lining + TOTAL_THICKNESS = SUM(LAYER_THICKNESS(1:N_LAYERS)) ! Thickness of the solid made up of OBSTs + DO NL=1,SF%N_LAYERS + DO NN=1,SF%N_LAYER_MATL(NL) + DO NNN=1,N_MATLS_NEW + IF (SF%LAYER_MATL_INDEX(NL,NN)==MATL_INDEX_NEW(NNN)) MATL_MASS_FRACTION_NEW(NL,NNN) = SF%MATL_MASS_FRACTION(NL,NN) + ENDDO + ENDDO + ENDDO + ! Add layers made made up of the OBSTs + DO NL=1,N_LAYERS + IF (NL>1) THEN ; LOWER_BOUND=SUM(LAYER_THICKNESS(1:NL-1)) ; ELSE ; LOWER_BOUND=0._EB ; ENDIF + UPPER_BOUND = SUM(LAYER_THICKNESS(1:NL)) + IF (UPPER_BOUND=TOTAL_THICKNESS-BACK_LINING_THICKNESS) CYCLE ! Lining covers layer + N_LAYERS_NEW = N_LAYERS_NEW + 1 + LAYER_THICKNESS_NEW(N_LAYERS_NEW) = MIN(TOTAL_THICKNESS-BACK_LINING_THICKNESS,UPPER_BOUND) - & + MAX(SUM(LAYER_THICKNESS_NEW(1:N_LAYERS_NEW-1)),LOWER_BOUND) + DO NN=1,N_MATLS + DO NNN=1,N_MATLS_NEW + IF (MATL_INDEX(NN)==MATL_INDEX_NEW(NNN)) MATL_MASS_FRACTION_NEW(N_LAYERS_NEW,NNN) = MATL_MASS_FRACTION(NL,NN) + ENDDO + ENDDO + ENDDO + ! Add layers from the back surface lining + DO NL=1,SF_BACK%N_LAYERS + N_LAYERS_NEW = N_LAYERS_NEW + 1 + LAYER_THICKNESS_NEW(N_LAYERS_NEW) = SF_BACK%LAYER_THICKNESS(SF_BACK%N_LAYERS-NL+1) + DO NN=1,SF_BACK%N_LAYER_MATL(NL) + DO NNN=1,N_MATLS_NEW + IF (SF_BACK%LAYER_MATL_INDEX(SF_BACK%N_LAYERS-NL+1,NN)==MATL_INDEX_NEW(NNN)) & + MATL_MASS_FRACTION_NEW(N_LAYERS_NEW,NNN) = SF_BACK%MATL_MASS_FRACTION(SF_BACK%N_LAYERS-NL+1,NN) + ENDDO + ENDDO + ENDDO + ! Copy the temporary _NEW arrays back into the original holding arrays + N_MATLS = N_MATLS_NEW + N_LAYERS = N_LAYERS_NEW + LAYER_THICKNESS(1:N_LAYERS) = LAYER_THICKNESS_NEW(1:N_LAYERS_NEW) + MATL_MASS_FRACTION(1:N_LAYERS,1:N_MATLS) = MATL_MASS_FRACTION_NEW(1:N_LAYERS_NEW,1:N_MATLS_NEW) + MATL_INDEX(1:N_MATLS) = MATL_INDEX_NEW(1:N_MATLS_NEW) + ENDIF + + ! If HT1D or HT3D, reallocate ONE_D arrays holding layer and material info + IF (SF%HT1D .OR. SF%HT_DIM>1) THEN ONE_D%N_LAYERS = N_LAYERS ONE_D%N_MATL = N_MATLS @@ -1564,6 +1617,25 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_2(NM) CONTAINS +SUBROUTINE ADD_MATERIAL(N_MATLS_X,MATL_INDEX_SEARCH,MATL_INDEX_X) + +INTEGER, INTENT(INOUT) :: N_MATLS_X +INTEGER :: JJ,MI,JJJ +INTEGER, INTENT(INOUT), DIMENSION(MAX_MATERIALS) :: MATL_INDEX_SEARCH,MATL_INDEX_X + +MATL_LOOP: DO JJ=1,MAX_MATERIALS + MI = MATL_INDEX_SEARCH(JJ) + IF (MI<1) EXIT MATL_LOOP + DO JJJ=1,N_MATLS_X + IF (MI==MATL_INDEX_X(JJJ)) CYCLE MATL_LOOP + ENDDO + N_MATLS_X = N_MATLS_X + 1 + MATL_INDEX_X(N_MATLS_X) = MI +ENDDO MATL_LOOP + +END SUBROUTINE ADD_MATERIAL + + SUBROUTINE REALLOCATE_WALL_SAVE INTEGER, ALLOCATABLE, DIMENSION(:) :: DUMMY diff --git a/Source/read.f90 b/Source/read.f90 index bb78a8f97da..78a4790ede6 100644 --- a/Source/read.f90 +++ b/Source/read.f90 @@ -7234,8 +7234,9 @@ SUBROUTINE READ_SURF ! Set up a dummy surface for HT1D and HT3D. The properties will be changed later. - If (HT1D .OR. HT3D .AND. THICKNESS(1)TWO_EPSILON_EB .AND. MATL_ID(1,1)/='null') SF%LINING = .TRUE. + If ((HT1D .OR. HT3D) .AND. THICKNESS(1),0,100000,,0,100000,-1.00E+09,1.00E+09,20,Heat_Transfer/ht3d_ibeam_devc.csv,2,3,Time,TS_x195-40|TS_x145-30|TS_x095-20|TS_x025-40|TS_x195-01|TS_x025-01,FDS 1|FDS 2|FDS 3|FDS 4|FDS 5|FDS 6,r-|k-|b-|g-|m-|c-,0,100000,,0,100000,-1.00E+09,1.00E+09,20,HT3D I-beam Surface Temperature (ht3d\_ibeam),Time (s),Temperature (°C),0,3600,1,0,1000,1,no,0.05 0.90,EastOutside,,1.2,linear,FDS_Verification_Guide/SCRIPT_FIGURES/ht3d_ibeam_TS,Relative Error,end,8.00E-02,Heat Transfer,r^,r,TeX d,ht3d_mass_conservation,Heat_Transfer/ht3d_mass_conservation_git.txt,Heat_Transfer/ht3d_mass_conservation.csv,1,2,Time,Mass,Exact (Mass),ko,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Heat_Transfer/ht3d_mass_conservation_mass.csv,2,3,Time,WOOD MOISTURE,FDS (WOOD MOISTURE),k-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Mass Balance (ht3d\_mass\_conservation),Time (s),Mass (kg),0,180,1,0,0.3,1,no,0.05 0.90,SouthEast,,1,linear,FDS_Verification_Guide/SCRIPT_FIGURES/ht3d_mass_conservation,Relative Error,end,1.00E-02,Heat Transfer,r^,r,TeX d,ht3d_mass_conservation,Heat_Transfer/ht3d_mass_conservation_2_git.txt,Heat_Transfer/ht3d_mass_conservation.csv,1,2,Time,Mass,Exact (Mass),ko,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Heat_Transfer/ht3d_mass_conservation_2_mass.csv,2,3,Time,WOOD MOISTURE,FDS (WOOD MOISTURE),k-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Mass Balance (ht3d\_mass\_conservation\_2),Time (s),Mass (kg),0,180,1,0,0.3,1,no,0.05 0.90,SouthEast,,1,linear,FDS_Verification_Guide/SCRIPT_FIGURES/ht3d_mass_conservation_2,Relative Error,end,1.00E-02,Heat Transfer,r^,r,TeX diff --git a/Verification/FDS_Cases.sh b/Verification/FDS_Cases.sh index 4f6acc5291c..3bf2e6862cf 100755 --- a/Verification/FDS_Cases.sh +++ b/Verification/FDS_Cases.sh @@ -252,6 +252,7 @@ $QFDS -d Heat_Transfer ht3d_energy_conservation_2.fds $QFDS -d Heat_Transfer ht3d_energy_conservation_3.fds $QFDS -p 8 -d Heat_Transfer ht3d_energy_conservation_4.fds $QFDS -p 8 -d Heat_Transfer ht3d_energy_conservation_5.fds +$QFDS -p 8 -d Heat_Transfer ht3d_energy_conservation_6.fds $QFDS -d Heat_Transfer ht3d_ibeam.fds $QFDS -d Heat_Transfer ht3d_mass_conservation.fds $QFDS -d Heat_Transfer ht3d_mass_conservation_2.fds diff --git a/Verification/Heat_Transfer/ht3d_energy_conservation_6.fds b/Verification/Heat_Transfer/ht3d_energy_conservation_6.fds new file mode 100644 index 00000000000..af29a615f41 --- /dev/null +++ b/Verification/Heat_Transfer/ht3d_energy_conservation_6.fds @@ -0,0 +1,52 @@ +&HEAD CHID='ht3d_energy_conservation_6' / + +&MESH IJK=25,25,25, XB=-0.25,0.00,-0.25,0.00,-0.25,0.00, MULT_ID='mesh' / +&MULT ID='mesh', DX=0.25, DY=0.25, DZ=0.25, I_UPPER=1, J_UPPER=1, K_UPPER=1 / + +&TIME T_END=100 / + +&VENT DB='XMIN', SURF_ID='OPEN' / +&VENT DB='XMAX', SURF_ID='OPEN' / +&VENT DB='YMIN', SURF_ID='OPEN' / +&VENT DB='YMAX', SURF_ID='OPEN' / +&VENT DB='ZMIN', SURF_ID='OPEN' / +&VENT DB='ZMAX', SURF_ID='OPEN' / + +&OBST XB=-0.20, 0.20,-0.20, 0.20,-0.05, 0.05, SURF_ID='SLAB', MATL_ID='STEEL' / + +&OBST XB=-0.01, 0.01,-0.01, 0.01, 0.06, 0.07, SURF_ID='HOT' / + +&SURF ID='SLAB', HT3D=T, COLOR='BEIGE', THICKNESS=0.015,0.01, + MATL_ID(1,1:2)='STUFF A','STUFF B', MATL_MASS_FRACTION(1,1:2)=0.3,0.7, + MATL_ID(2,1:2)='STUFF B','STUFF C', MATL_MASS_FRACTION(2,1:2)=0.2,0.8 / + +&MATL ID='STUFF A', DENSITY= 50, SPECIFIC_HEAT=1.0, CONDUCTIVITY=0.1 / +&MATL ID='STUFF B', DENSITY= 80, SPECIFIC_HEAT=1.5, CONDUCTIVITY=0.2 / +&MATL ID='STUFF C', DENSITY= 30, SPECIFIC_HEAT=2.5, CONDUCTIVITY=0.3 / +&MATL ID='STEEL', DENSITY=7500, SPECIFIC_HEAT=0.5, CONDUCTIVITY=50. / + +&SURF ID='HOT', TMP_FRONT=1000, COLOR='RED' / + +&BNDF QUANTITY='WALL TEMPERATURE', CELL_CENTERED=T / + +&SLCF PBY=0.001, QUANTITY='TEMPERATURE', CELL_CENTERED=T / + +&DUMP DT_PROF=5., DT_DEVC=4. / + +&PROF ID='top', QUANTITY='TEMPERATURE', XYZ=0.005,0.005,0.050, IOR= 3 / +&PROF ID='side', QUANTITY='TEMPERATURE', XYZ=0.200,0.005,0.040, IOR= 1 / +&PROF ID='bottom', QUANTITY='TEMPERATURE', XYZ=0.005,0.005,-.050, IOR=-3 / + +'WALL ENTHALPY' is the energy (kJ) of the volume of solid bounded by the surface cell. The CONVERSION_FACTOR is intended to +cancel out the cell area 0.01 m x 0.01 m + +&DEVC XB=-0.24,0.24,-0.24,0.24,-0.24,0.24, QUANTITY='WALL ENTHALPY', SPATIAL_STATISTIC='SURFACE INTEGRAL', ID='H1', IOR=-1, TIME_AVERAGED=F, RELATIVE=T, CONVERSION_FACTOR=10000, SURF_ID='SLAB' / +&DEVC XB=-0.24,0.24,-0.24,0.24,-0.24,0.24, QUANTITY='WALL ENTHALPY', SPATIAL_STATISTIC='SURFACE INTEGRAL', ID='H2', IOR= 1, TIME_AVERAGED=F, RELATIVE=T, CONVERSION_FACTOR=10000, SURF_ID='SLAB' / +&DEVC XB=-0.24,0.24,-0.24,0.24,-0.24,0.24, QUANTITY='WALL ENTHALPY', SPATIAL_STATISTIC='SURFACE INTEGRAL', ID='H3', IOR=-2, TIME_AVERAGED=F, RELATIVE=T, CONVERSION_FACTOR=10000, SURF_ID='SLAB' / +&DEVC XB=-0.24,0.24,-0.24,0.24,-0.24,0.24, QUANTITY='WALL ENTHALPY', SPATIAL_STATISTIC='SURFACE INTEGRAL', ID='H4', IOR= 2, TIME_AVERAGED=F, RELATIVE=T, CONVERSION_FACTOR=10000, SURF_ID='SLAB' / +&DEVC XB=-0.24,0.24,-0.24,0.24,-0.24,0.24, QUANTITY='WALL ENTHALPY', SPATIAL_STATISTIC='SURFACE INTEGRAL', ID='H5', IOR=-3, TIME_AVERAGED=F, RELATIVE=T, CONVERSION_FACTOR=10000, SURF_ID='SLAB' / +&DEVC XB=-0.24,0.24,-0.24,0.24,-0.24,0.24, QUANTITY='WALL ENTHALPY', SPATIAL_STATISTIC='SURFACE INTEGRAL', ID='H6', IOR= 3, TIME_AVERAGED=F, RELATIVE=T, CONVERSION_FACTOR=10000, SURF_ID='SLAB' / + +&DEVC XB=-0.24,0.24,-0.24,0.24,-0.24,0.24, QUANTITY='NET HEAT FLUX', SPATIAL_STATISTIC='SURFACE INTEGRAL', TEMPORAL_STATISTIC='TIME INTEGRAL', ID='Q_net', SURF_ID='SLAB' / + +&TAIL /