Merge pull request firemodels#12096 from mcgratta/natural_convection

Natural convection for downward facing hot plate

Source/turb.f90

@@ -1290,10 +1290,10 @@ END SUBROUTINE WALL_MODEL
 !> \param SURF_GEOMETRY_INDEX Indicator of the surface geometry
 !> \param IOR Index of the surface orientation
 
-SUBROUTINE NATURAL_CONVECTION_MODEL(NUSSELT,RA,SURF_INDEX,SURF_GEOMETRY_INDEX,IOR)
+SUBROUTINE NATURAL_CONVECTION_MODEL(NUSSELT,RA,SURF_INDEX,SURF_GEOMETRY_INDEX,IOR,DELTA_TMP)
 
 REAL(EB), INTENT(OUT) :: NUSSELT
-REAL(EB), INTENT(IN) :: RA
+REAL(EB), INTENT(IN) :: RA,DELTA_TMP
 INTEGER, INTENT(IN) :: SURF_INDEX,SURF_GEOMETRY_INDEX,IOR
 
 SELECT CASE(SURF_GEOMETRY_INDEX)
@@ -1302,10 +1302,14 @@ SUBROUTINE NATURAL_CONVECTION_MODEL(NUSSELT,RA,SURF_INDEX,SURF_GEOMETRY_INDEX,IO
          CASE(0:2)
             NUSSELT = (0.825_EB + 0.324_EB*RA**ONSI)**2  ! Incropera and DeWitt, 7th edition, Eq. 9.26
          CASE(3)
-            IF (RA<1.E7_EB) THEN
-               NUSSELT = 0.54_EB*RA**0.25_EB  ! Incropera and DeWitt, 7th edition, Eq. 9.30
+            IF ((IOR==3.AND.DELTA_TMP<0._EB) .OR. (IOR==-3.AND.DELTA_TMP>=0._EB)) THEN
+               IF (RA<1.E7_EB) THEN
+                  NUSSELT = 0.54_EB*RA**0.25_EB  ! Incropera and DeWitt, 7th edition, Eq. 9.30
+               ELSE
+                  NUSSELT = 0.15_EB*RA**ONTH     ! Incropera and DeWitt, 7th edition, Eq. 9.31
+               ENDIF
             ELSE
-               NUSSELT = 0.15_EB*RA**ONTH     ! Incropera and DeWitt, 7th edition, Eq. 9.31
+               NUSSELT = 0.52_EB*RA**0.2     ! Incropera and DeWitt, 7th edition, Eq. 9.32
             ENDIF
       END SELECT
    CASE (SURF_CYLINDRICAL)  ! Simplification of Eq. 9.34, Incropera and DeWitt, 7th edition

Source/wall.f90

@@ -3211,7 +3211,7 @@ REAL(EB) FUNCTION HEAT_TRANSFER_COEFFICIENT(DELTA_TMP,H_FIXED,SURF_INDEX_IN,WALL
                CALL FORCED_CONVECTION_MODEL(NUSSELT_FORCED,RE,PR_ONTH,SURF_GEOMETRY)
             ENDIF
             RA = GR*PR_AIR
-            CALL NATURAL_CONVECTION_MODEL(NUSSELT_FREE,RA,SURF_INDEX_IN,SFX%GEOMETRY,BCX%IOR)
+            CALL NATURAL_CONVECTION_MODEL(NUSSELT_FREE,RA,SURF_INDEX_IN,SFX%GEOMETRY,BCX%IOR,DELTA_TMP)
             IF (PRESENT(PARTICLE_INDEX_IN)) THEN
                HEAT_TRANSFER_COEFFICIENT = MAX(NUSSELT_FORCED,NUSSELT_FREE)*K_G/CONV_LENGTH
             ELSE

Utilities/Matlab/FDS_verification_script.m

@@ -104,6 +104,7 @@
 disp('natconh...');                       natconh
 disp('natconv...');                       natconv
 disp('freecon_sphere...');                freecon_sphere
+disp('nat_conv_hot_plate...');            nat_conv_hot_plate
 disp('tree_shapes...');                   tree_shapes
 
 display('verification scripts completed successfully!')

Utilities/Matlab/scripts/freecon_sphere.m

@@ -120,7 +120,7 @@
     end
 end
 
-lh=legend(marker_handle,'Yuge (1960), Nu = 2 + 0.43 Ra^{1/4}','Amato & Tien (1972), Nu = 2 + 0.5 Ra^{1/4}','FDS {\itD/\deltax}=8','FDS {\itD/\deltax}=16','Location','Northwest');
+lh=legend(marker_handle,'Yuge (1960)','Amato & Tien (1972)','FDS {\itD/\deltax}=8','FDS {\itD/\deltax}=16','Location','Northwest');
 set(lh,'FontName',Font_Name,'FontSize',Key_Font_Size)
 
 % title('Free Convection from a Sphere','FontSize',Label_Font_Size)

Utilities/Matlab/scripts/nat_conv_hot_plate.m

@@ -0,0 +1,104 @@
+% McGrattan
+% September 1, 2023
+% nat_conv_hot_plate.m
+
+close all
+clear all
+
+plot_style
+
+results_dir = ['../../../out/Convection/'];
+
+% calculations below were used for input file setup
+
+g  = 9.80665;
+S  = [1.0  0.1  5.0  10 .05];
+T1 = [503  503  503 323 323];
+T2 = 293;
+MW = 28.85476; % FDS 'LJ AIR'
+P0 = 101325;
+mu = 1.8216e-5;
+cp = 1000;
+k  = 0.018216; % for Pr=1 fluid
+
+RAYLEIGH_DOWN = logspace(4, 9,100);
+RAYLEIGH_UP   = logspace(4,11,100);
+for i=1:length(RAYLEIGH_DOWN)
+   NUSSELT_DOWN(i)=0.52*RAYLEIGH_DOWN(i)^0.2;  % Incropera and DeWitt, 7th, Eq. 9.32
+end
+for i=1:length(RAYLEIGH_UP)
+   if RAYLEIGH_UP(i)<1e7; 
+      NUSSELT_UP(i)=0.54*RAYLEIGH_UP(i)^0.25;  % Incropera and DeWitt, 7th, Eq. 9.30
+   else
+      NUSSELT_UP(i)=0.15*RAYLEIGH_UP(i)^0.333; % Incropera and DeWitt, 7th, Eq. 9.31
+   end
+end
+
+figure
+set(gcf,'Visible',Figure_Visibility);
+set(gca,'Units',Plot_Units)
+set(gca,'Position',[Plot_X Plot_Y Plot_Width Plot_Height])
+
+marker_handle(1)=loglog(RAYLEIGH_DOWN,NUSSELT_DOWN,'k-'); hold on
+marker_handle(2)=loglog(RAYLEIGH_UP  ,NUSSELT_UP  ,'k--'); hold on
+set(gca,'FontName',Font_Name)
+set(gca,'FontSize',Label_Font_Size)
+axis([1e0 1e15 1 1000])
+ytick = [1 2 5 10 20 50 100 200 500 1000];
+yticks(ytick);
+xlabel('Rayleigh Number','FontSize',Label_Font_Size)
+ylabel('Nusselt Number','FontSize',Label_Font_Size)
+
+% FDS results
+
+casename={...
+'nat_conv_hot_plate_1',...
+'nat_conv_hot_plate_2',...
+'nat_conv_hot_plate_3',...
+'nat_conv_hot_plate_4',...
+'nat_conv_hot_plate_5',...
+};
+
+marker_style = {'gsq','rsq','ksq','go','ro','ko'};
+res = {'8','16','32'};
+
+for j=1:length(res)
+    for i=1:length(S)
+
+        M = importdata([results_dir,casename{i},'_',res{j},'_devc.csv']);
+
+        % check for steady state
+        Qdot_down = 1000*M.data(end,2);
+        Qdot_up = 1000*M.data(end,3);
+        A = S(i)^2;
+        Tm = 0.5*(T1(i)+T2);
+        beta = 1/Tm;
+        rho = P0*MW/(8341.5*Tm);
+        alpha = k/(rho*cp);
+        nu = mu/rho;
+        L = S(i)/4;
+        Ra_FDS = (g*beta*(T1(i)-T2)*L^3)/(alpha*nu);
+        Nu_FDS_down = (Qdot_down/A)*(L/k)/(T1(i)-T2);
+        Nu_FDS_up   = (Qdot_up  /A)*(L/k)/(T1(i)-T2);
+
+        marker_handle(j+2)=plot(Ra_FDS,Nu_FDS_down,marker_style{j},'MarkerSize',8);
+        marker_handle(j+5)=plot(Ra_FDS,Nu_FDS_up,marker_style{j+3},'MarkerSize',8);
+    end
+end
+
+lh=legend(marker_handle,'Correlation; Down','Correlation; Up','Down ($S/\delta x=8$)','Down ($S/\delta x=16$)','Down ($S/\delta x=32$)',...
+  'Up ($S/\delta x=8$)','Up ($S/\delta x=16$)','Up ($S/\delta x=32$)','Location','Northwest','Interpreter','LaTeX');
+set(lh,'FontName',Font_Name,'FontSize',Key_Font_Size)
+
+% add Git if file is available
+
+Git_Filename = [results_dir,'nat_conv_hot_plate_1_8_git.txt'];
+addverstr(gca,Git_Filename,'loglog')
+
+set(gcf,'Visible',Figure_Visibility);
+set(gcf,'Units',Paper_Units);
+set(gcf,'PaperSize',[Paper_Width Paper_Height]);
+set(gcf,'Position',[0 0 Paper_Width Paper_Height]);
+print(gcf,'-dpdf','../../Manuals/FDS_Verification_Guide/SCRIPT_FIGURES/nat_conv_hot_plate');
+
+

Validation/Convection/FDS_Input_Files/nat_conv_hot_plate_1_16.fds

@@ -0,0 +1,39 @@
+&HEAD CHID='nat_conv_hot_plate_1_16', TITLE='Natural convection from downward facing hot plate' /
+
+&MULT DX=1, DY=1, I_UPPER=3, J_UPPER=1 ID='m1' /
+&MESH IJK=16,16,16, XB=-1,0.0,-1,0.0,0.0,1, MULT_ID='m1' /
+
+&TIME T_END=100 /
+
+&DUMP DT_DEVC=1. /
+
+&MISC STRATIFICATION=F /
+
+&RADI RADIATION=F /
+
+&SPEC ID='LJ AIR', SPECIFIC_HEAT=1., CONDUCTIVITY=0.018216, VISCOSITY=1.8216E-5, BACKGROUND=.TRUE./ ! Pr=1.00
+
+&OBST XB=-0.5,0.5,-0.5,0.5,0.45,0.55, SURF_IDS='WALL','WALL','HOT' /
+&OBST XB=1.5,2.5,-0.5,0.5,0.45,0.55, SURF_IDS='HOT','WALL','WALL' /
+
+&SURF ID='WALL', ADIABATIC=T, DEFAULT=T /
+&SURF ID='HOT', COLOR='RED', TMP_FRONT=230, CONVECTION_LENGTH_SCALE=0.25 /
+
+&VENT MB='XMIN', SURF_ID='OPEN' /
+&VENT MB='XMAX', SURF_ID='OPEN' /
+&VENT MB='YMIN', SURF_ID='OPEN' /
+&VENT MB='YMAX', SURF_ID='OPEN' /
+&VENT MB='ZMIN', SURF_ID='OPEN' /
+&VENT MB='ZMAX', SURF_ID='OPEN' /
+
+&BNDF QUANTITY='HEAT TRANSFER COEFFICIENT', CELL_CENTERED=T/
+&BNDF QUANTITY='CONVECTIVE HEAT FLUX', CELL_CENTERED=T/
+&BNDF QUANTITY='THERMAL WALL UNITS', CELL_CENTERED=T /
+&BNDF QUANTITY='VISCOUS WALL UNITS', CELL_CENTERED=T /
+
+&SLCF PBY=0.001, QUANTITY='TEMPERATURE', CELL_CENTERED=.TRUE. /
+
+&DEVC ID='qdot1', QUANTITY='CONVECTIVE HEAT FLUX', XB=-0.5,0.5,-0.5,0.5,0.0,1, SURF_ID='HOT', SPATIAL_STATISTIC='SURFACE INTEGRAL', CONVERSION_FACTOR=-1. /
+&DEVC ID='qdot2', QUANTITY='CONVECTIVE HEAT FLUX', XB=1.5,2.5,-0.5,0.5,0.0,1, SURF_ID='HOT', SPATIAL_STATISTIC='SURFACE INTEGRAL', CONVERSION_FACTOR=-1. /
+
+&TAIL /

Validation/Convection/FDS_Input_Files/nat_conv_hot_plate_1_32.fds

@@ -0,0 +1,39 @@
+&HEAD CHID='nat_conv_hot_plate_1_32', TITLE='Natural convection from downward facing hot plate' /
+
+&MULT DX=1, DY=1, I_UPPER=3, J_UPPER=1 ID='m1' /
+&MESH IJK=32,32,32, XB=-1,0.0,-1,0.0,0.0,1, MULT_ID='m1' /
+
+&TIME T_END=100 /
+
+&DUMP DT_DEVC=1. /
+
+&MISC STRATIFICATION=F /
+
+&RADI RADIATION=F /
+
+&SPEC ID='LJ AIR', SPECIFIC_HEAT=1., CONDUCTIVITY=0.018216, VISCOSITY=1.8216E-5, BACKGROUND=.TRUE./ ! Pr=1.00
+
+&OBST XB=-0.5,0.5,-0.5,0.5,0.45,0.55, SURF_IDS='WALL','WALL','HOT' /
+&OBST XB=1.5,2.5,-0.5,0.5,0.45,0.55, SURF_IDS='HOT','WALL','WALL' /
+
+&SURF ID='WALL', ADIABATIC=T, DEFAULT=T /
+&SURF ID='HOT', COLOR='RED', TMP_FRONT=230, CONVECTION_LENGTH_SCALE=0.25 /
+
+&VENT MB='XMIN', SURF_ID='OPEN' /
+&VENT MB='XMAX', SURF_ID='OPEN' /
+&VENT MB='YMIN', SURF_ID='OPEN' /
+&VENT MB='YMAX', SURF_ID='OPEN' /
+&VENT MB='ZMIN', SURF_ID='OPEN' /
+&VENT MB='ZMAX', SURF_ID='OPEN' /
+
+&BNDF QUANTITY='HEAT TRANSFER COEFFICIENT', CELL_CENTERED=T/
+&BNDF QUANTITY='CONVECTIVE HEAT FLUX', CELL_CENTERED=T/
+&BNDF QUANTITY='THERMAL WALL UNITS', CELL_CENTERED=T /
+&BNDF QUANTITY='VISCOUS WALL UNITS', CELL_CENTERED=T /
+
+&SLCF PBY=0.001, QUANTITY='TEMPERATURE', CELL_CENTERED=.TRUE. /
+
+&DEVC ID='qdot1', QUANTITY='CONVECTIVE HEAT FLUX', XB=-0.5,0.5,-0.5,0.5,0.0,1, SURF_ID='HOT', SPATIAL_STATISTIC='SURFACE INTEGRAL', CONVERSION_FACTOR=-1. /
+&DEVC ID='qdot2', QUANTITY='CONVECTIVE HEAT FLUX', XB=1.5,2.5,-0.5,0.5,0.0,1, SURF_ID='HOT', SPATIAL_STATISTIC='SURFACE INTEGRAL', CONVERSION_FACTOR=-1. /
+
+&TAIL /

Validation/Convection/FDS_Input_Files/nat_conv_hot_plate_1_8.fds

@@ -0,0 +1,39 @@
+&HEAD CHID='nat_conv_hot_plate_1_8', TITLE='Natural convection from downward facing hot plate' /
+
+&MULT DX=1, DY=1, I_UPPER=3, J_UPPER=1 ID='m1' /
+&MESH IJK=8,8,8, XB=-1,0.0,-1,0.0,0.0,1, MULT_ID='m1' /
+
+&TIME T_END=100 /
+
+&DUMP DT_DEVC=1. /
+
+&MISC STRATIFICATION=F /
+
+&RADI RADIATION=F /
+
+&SPEC ID='LJ AIR', SPECIFIC_HEAT=1., CONDUCTIVITY=0.018216, VISCOSITY=1.8216E-5, BACKGROUND=.TRUE./ ! Pr=1.00
+
+&OBST XB=-0.5,0.5,-0.5,0.5,0.45,0.55, SURF_IDS='WALL','WALL','HOT' /
+&OBST XB=1.5,2.5,-0.5,0.5,0.45,0.55, SURF_IDS='HOT','WALL','WALL' /
+
+&SURF ID='WALL', ADIABATIC=T, DEFAULT=T /
+&SURF ID='HOT', COLOR='RED', TMP_FRONT=230, CONVECTION_LENGTH_SCALE=0.25 /
+
+&VENT MB='XMIN', SURF_ID='OPEN' /
+&VENT MB='XMAX', SURF_ID='OPEN' /
+&VENT MB='YMIN', SURF_ID='OPEN' /
+&VENT MB='YMAX', SURF_ID='OPEN' /
+&VENT MB='ZMIN', SURF_ID='OPEN' /
+&VENT MB='ZMAX', SURF_ID='OPEN' /
+
+&BNDF QUANTITY='HEAT TRANSFER COEFFICIENT', CELL_CENTERED=T/
+&BNDF QUANTITY='CONVECTIVE HEAT FLUX', CELL_CENTERED=T/
+&BNDF QUANTITY='THERMAL WALL UNITS', CELL_CENTERED=T /
+&BNDF QUANTITY='VISCOUS WALL UNITS', CELL_CENTERED=T /
+
+&SLCF PBY=0.001, QUANTITY='TEMPERATURE', CELL_CENTERED=.TRUE. /
+
+&DEVC ID='qdot1', QUANTITY='CONVECTIVE HEAT FLUX', XB=-0.5,0.5,-0.5,0.5,0.0,1, SURF_ID='HOT', SPATIAL_STATISTIC='SURFACE INTEGRAL', CONVERSION_FACTOR=-1. /
+&DEVC ID='qdot2', QUANTITY='CONVECTIVE HEAT FLUX', XB=1.5,2.5,-0.5,0.5,0.0,1, SURF_ID='HOT', SPATIAL_STATISTIC='SURFACE INTEGRAL', CONVERSION_FACTOR=-1. /
+
+&TAIL /

Validation/Convection/FDS_Input_Files/nat_conv_hot_plate_2_16.fds

@@ -0,0 +1,39 @@
+&HEAD CHID='nat_conv_hot_plate_2_16', TITLE='Natural convection from downward facing hot plate' /
+
+&MULT DX=0.1, DY=0.1, I_UPPER=3, J_UPPER=1 ID='m1' /
+&MESH IJK=16,16,16, XB=-0.1,0.0,-0.1,0.0,0.0,0.1, MULT_ID='m1' /
+
+&TIME T_END=10 /
+
+&DUMP DT_DEVC=1. /
+
+&MISC STRATIFICATION=F /
+
+&RADI RADIATION=F /
+
+&SPEC ID='LJ AIR', SPECIFIC_HEAT=1., CONDUCTIVITY=0.018216, VISCOSITY=1.8216E-5, BACKGROUND=.TRUE./ ! Pr=1.00
+
+&OBST XB=-0.05,0.05,-0.05,0.05,0.045,0.055, SURF_IDS='WALL','WALL','HOT' /
+&OBST XB=0.15,0.25,-0.05,0.05,0.045,0.055, SURF_IDS='HOT','WALL','WALL' /
+
+&SURF ID='WALL', ADIABATIC=T, DEFAULT=T /
+&SURF ID='HOT', COLOR='RED', TMP_FRONT=230, CONVECTION_LENGTH_SCALE=0.025 /
+
+&VENT MB='XMIN', SURF_ID='OPEN' /
+&VENT MB='XMAX', SURF_ID='OPEN' /
+&VENT MB='YMIN', SURF_ID='OPEN' /
+&VENT MB='YMAX', SURF_ID='OPEN' /
+&VENT MB='ZMIN', SURF_ID='OPEN' /
+&VENT MB='ZMAX', SURF_ID='OPEN' /
+
+&BNDF QUANTITY='HEAT TRANSFER COEFFICIENT', CELL_CENTERED=T/
+&BNDF QUANTITY='CONVECTIVE HEAT FLUX', CELL_CENTERED=T/
+&BNDF QUANTITY='THERMAL WALL UNITS', CELL_CENTERED=T /
+&BNDF QUANTITY='VISCOUS WALL UNITS', CELL_CENTERED=T /
+
+&SLCF PBY=0.001, QUANTITY='TEMPERATURE', CELL_CENTERED=.TRUE. /
+
+&DEVC ID='qdot1', QUANTITY='CONVECTIVE HEAT FLUX', XB=-0.05,0.05,-0.05,0.05,0.0,0.1, SURF_ID='HOT', SPATIAL_STATISTIC='SURFACE INTEGRAL', CONVERSION_FACTOR=-1. /
+&DEVC ID='qdot2', QUANTITY='CONVECTIVE HEAT FLUX', XB=0.15,0.25,-0.05,0.05,0.0,0.1, SURF_ID='HOT', SPATIAL_STATISTIC='SURFACE INTEGRAL', CONVERSION_FACTOR=-1. /
+
+&TAIL /

Validation/Convection/FDS_Input_Files/nat_conv_hot_plate_2_32.fds

@@ -0,0 +1,39 @@
+&HEAD CHID='nat_conv_hot_plate_2_32', TITLE='Natural convection from downward facing hot plate' /
+
+&MULT DX=0.1, DY=0.1, I_UPPER=3, J_UPPER=1 ID='m1' /
+&MESH IJK=32,32,32, XB=-0.1,0.0,-0.1,0.0,0.0,0.1, MULT_ID='m1' /
+
+&TIME T_END=10 /
+
+&DUMP DT_DEVC=1. /
+
+&MISC STRATIFICATION=F /
+
+&RADI RADIATION=F /
+
+&SPEC ID='LJ AIR', SPECIFIC_HEAT=1., CONDUCTIVITY=0.018216, VISCOSITY=1.8216E-5, BACKGROUND=.TRUE./ ! Pr=1.00
+
+&OBST XB=-0.05,0.05,-0.05,0.05,0.045,0.055, SURF_IDS='WALL','WALL','HOT' /
+&OBST XB=0.15,0.25,-0.05,0.05,0.045,0.055, SURF_IDS='HOT','WALL','WALL' /
+
+&SURF ID='WALL', ADIABATIC=T, DEFAULT=T /
+&SURF ID='HOT', COLOR='RED', TMP_FRONT=230, CONVECTION_LENGTH_SCALE=0.025 /
+
+&VENT MB='XMIN', SURF_ID='OPEN' /
+&VENT MB='XMAX', SURF_ID='OPEN' /
+&VENT MB='YMIN', SURF_ID='OPEN' /
+&VENT MB='YMAX', SURF_ID='OPEN' /
+&VENT MB='ZMIN', SURF_ID='OPEN' /
+&VENT MB='ZMAX', SURF_ID='OPEN' /
+
+&BNDF QUANTITY='HEAT TRANSFER COEFFICIENT', CELL_CENTERED=T/
+&BNDF QUANTITY='CONVECTIVE HEAT FLUX', CELL_CENTERED=T/
+&BNDF QUANTITY='THERMAL WALL UNITS', CELL_CENTERED=T /
+&BNDF QUANTITY='VISCOUS WALL UNITS', CELL_CENTERED=T /
+
+&SLCF PBY=0.001, QUANTITY='TEMPERATURE', CELL_CENTERED=.TRUE. /
+
+&DEVC ID='qdot1', QUANTITY='CONVECTIVE HEAT FLUX', XB=-0.05,0.05,-0.05,0.05,0.0,0.1, SURF_ID='HOT', SPATIAL_STATISTIC='SURFACE INTEGRAL', CONVERSION_FACTOR=-1. /
+&DEVC ID='qdot2', QUANTITY='CONVECTIVE HEAT FLUX', XB=0.15,0.25,-0.05,0.05,0.0,0.1, SURF_ID='HOT', SPATIAL_STATISTIC='SURFACE INTEGRAL', CONVERSION_FACTOR=-1. /
+
+&TAIL /