Radiation to solid sphere; energy conservation

[mcgratta]

I have a metal sphere weighing 0.005 kg in a 1 m cube box, sealed, adiabatic walls. The ambient temperature is 800 C, and the sphere is initially 20 C and has a specific heat of c=1 kJ/kg/K. The box is filled with N2 with a specific heat of 1 kJ/kg/K. If I turn off radiation and set the convective heat transfer coefficient to 100 W/m/K, the ball heats up to 783 C, gaining 3.81 kJ in internal energy: (T_f-T_i)mc. The N2 loses dE = dH - Vdp the same amount of internal energy.

Now, if I turn off convection by setting HEAT_TRANSFER_COEFFICIENT to 0 and turn on RADIATION, I do not get the same result. In both cases, the heat flux to the ball is about the same.

&HEAD CHID='particle_rad', TITLE='Test of energy conservation' /

&TIME T_END=900, DT=0.1, WALL_INCREMENT=1 /

&MISC TMPA=800. /

&MESH IJK=10,10,10 XB=0.0,1.0,0.0,1.0,0.0,1.0 /

&SPEC ID='NITROGEN', SPECIFIC_HEAT=1., BACKGROUND=T /

&SURF ID='WALL', ADIABATIC=T, DEFAULT=T /

&SURF ID                   = 'ball'
      MATL_ID              = 'metal'
      THICKNESS            = 0.006203505
      HEAT_TRANSFER_COEFFICIENT=100.
      TMP_INNER            = 20.
      GEOMETRY             = 'SPHERICAL' /

&RADI RADIATION=F /

&PART ID='ball'
      SURF_ID='ball'
      STATIC=.TRUE.  /

&INIT ID                   = 'ball',
      PART_ID              = 'ball'
      XYZ                  = 0.45,0.45,0.45
      N_PARTICLES          = 1 /

&MATL ID                 = 'metal'
      DENSITY            = 5000.
      CONDUCTIVITY       = 50.
      SPECIFIC_HEAT      = 1.0 /

&DUMP MASS_FILE=T, DT_DEVC=0.0000001 /

&DEVC XB=0,1,0,1,0,1, QUANTITY='MPUV', PART_ID='ball', ID='mass', STATISTICS='VOLUME INTEGRAL'  /
&DEVC ID='T_surf',INIT_ID='ball', QUANTITY='WALL TEMPERATURE' /
&DEVC XB=0,1,0,1,0,1, QUANTITY='ENTHALPY', SPATIAL_STATISTIC='VOLUME INTEGRAL', ID='Delta h_gas', RELATIVE=.TRUE. /
&DEVC XB=0,1,0,1,0,1, QUANTITY='PRESSURE', ID='Delta p', SPATIAL_STATISTIC='VOLUME MEAN', CONVERSION_FACTOR=0.001, UNITS='kPa' /
&DEVC XB=0,1,0,1,0,1, QUANTITY='TEMPERATURE', ID='mean temp', SPATIAL_STATISTIC='VOLUME MEAN' /

&SLCF PBY=0.45, QUANTITY='TEMPERATURE', CELL_CENTERED=T /
&SLCF PBY=0.45, QUANTITY='PRESSURE', CELL_CENTERED=T /
&SLCF PBY=0.45, QUANTITY='ENTHALPY', CELL_CENTERED=T /
&SLCF PBY=0.45, QUANTITY='PARTICLE RADIATION LOSS', CELL_CENTERED=T /
&SLCF PBY=0.45, QUANTITY='DIVERGENCE', CELL_CENTERED=T /

&TAIL /

[shostikk]

Which one goes wrong - gas enthalpy or sphere?

[mcgratta]

I am not sure what you mean. In the convection case, the gas temperature drops from 800 C to about 783 C. In the radiation case, the temperature only drops to about 797 C.

[mcgratta]

The source term looks correct; that is, the QR_W term is A*sigma*(Tw^4-Tb^4)/V_cell. Initially, Tw=1074 K, Tb=293 K, V_cell = 0.001 m3, A=4.83e-4 m2, and QR_W=36 kW/m3 in the cell with the particle, as seen in the slice file. So this might be a problem with the convection at the outer wall.

[rmcdermo]

Kevin, These output lines make the devc file pretty easy to compare the energy diffs.

&DEVC XB=0,1,0,1,0,1, QUANTITY='INTERNAL ENERGY', SPATIAL_STATISTIC='VOLUME INTEGRAL', ID='Delta e_gas', RELATIVE=T /
&DEVC XB=0,1,0,1,0,1, QUANTITY='WALL ENTHALPY', SURF_ID='ball', INIT_ID='ball', ID='Delta e_part', RELATIVE=T /

[mcgratta]

Nice -- I wasn't even aware of INTERNAL ENERGY, but I see it's in the Guide. These lines

&DEVC XB=0,1,0,1,0,1, QUANTITY='PARTICLE RADIATION LOSS', ID='part loss',
SPATIAL_STATISTIC='VOLUME INTEGRAL' /
&DEVC XB=0,1,0,1,0,1, QUANTITY='NET HEAT FLUX', ID='wall loss',
SURF_ID='WALL', SPATIAL_STATISTIC='SURFACE INTEGRAL' /

seem to indicate that something is going wrong with the calculation of heat flux at the exterior wall. These two quantities are initially the same in absolute value and then are not for reasons that should become clear soon.

[rmcdermo]

I did not need the SURF_ID above. Should have been:

&DEVC XB=0,1,0,1,0,1, QUANTITY='WALL ENTHALPY', INIT_ID='ball', ID='Delta e_part', RELATIVE=T /

[mcgratta]

If I solve the radiation transport at every time step and all angles, the solution is better, but still not right. This is some odd behavior of the algorithm to find the wall temperature to attain no net heat flux. I think the radiation solver is working properly, as is the particle conduction, etc.

[drjfloyd]

If instead of ADIABATIC=T, you set HEAT_TRANSFER_COEFFICIENT to 0 and EMISSIVITY to 0 does that solve the problem?

[mcgratta]

I don't think so, because we need convective heat transfer between the gas and exterior wall. That is what actually cools off the cube of gas. There appears to be something funky going on with the simple Newton iteration to find the wall temperature such that the net heat flux is zero.

[rmcdermo]

Is the linearization getting us?

[drjfloyd]

ADIABATIC=T means that convection = - radiation at the wall and there is no net heat transfer at the wall to cool the cube of gas. This means we may redistribute energy since highly radiative areas can increase the wall temperature which is then convected to the gas, but total internal energy should remain the same. Setting HTC=0 and E=0 should just stop this resdistribution. If this noticeable improves the results then we need to consider if we need changes to how we set the wall temperature for adiabatic

[mcgratta]

Not sure. It's most likely the fact that we are making really small changes in wall temperature, and the iteration method might be too crude.

[marcosvanella]

Is the heat flux in the ball surface kept constant or (linearly varying) through the integration timestep? If not, and solid and gas sides use different integration schemes this might lead to the energy difference.

[MortezaHaghighi]

Can CHECK_HT=T help in any way?

[mcgratta]

CHECK_HT does lower the time step, but that doesn't seem to help.

[rmcdermo]

How close are you on this now? My working hypothesis is that we run into truncation error in the TMP_F calculation at the WALL. I cannot get a BNDF of 'NET HEAT FLUX' to be 0 on the WALL. The closest I get is something like 1E-4 (but there seems to be a problem with SMV on this so hard for me to tell for sure).

I modified the wall.f90 to do this (at about line 353):

         ! IF (ABS(TMP_OTHER - ONE_D%TMP_F) / ONE_D%TMP_F < 1.E-6_EB .OR. ADCOUNT > 20) THEN
         IF (ABS(QEXTRA) < 1.E-10_EB .OR. ADCOUNT > 20) THEN
            ONE_D%TMP_F = MIN(TMPMAX,TMP_OTHER)
            EXIT ADLOOP

The resulting energy balance gives:

s,kJ,kJ
Time,"Delta e_gas","Delta e_part"
...
9.0000000E+002,-3.7369159E+000, 3.8164733E+000

for this input file: particle_rad.fds.txt

Am I making progress or are you already this close?

Are we on the same page with not getting the BNDF of 'NET HEAT FLUX' to be zero for an ADIABATIC boundary?

[mcgratta]

I'm coming to the conclusion that the problem is ill-posed. Yes, I did run into issues of accuracy, but no matter how small the time step, and other code hacks, I could not get it.

I'm trying the modified case where we have a thin layer of insulated metal around the box. At least this problem is well-posed.

[rmcdermo]

So you are getting the same numbers as me for the energy balance? My numbers are off by 2% relative error. Are you in this ball park?

[mcgratta]

Bravo! Yes, I'm getting those numbers. Weird because I tried several variants of that tolerance statement. Clearly not the right one, and then went down many rabbit holes. I'll play around with this a bit more. BTW, the thin metal with insulated backing case is essentially the same thing as this one. That's the theory of the plate thermometer.

[rmcdermo]

Right now I’m trying a variation with a small OBST instead of a particle.

On Mon, Jun 22, 2020 at 7:30 PM Kevin McGrattan notifications@github.com wrote:

Bravo! Yes, I'm getting those numbers. Weird because I tried several variants of that tolerance statement. Clearly not the right one, and then went down many rabbit holes. I'll play around with this a bit more. BTW, the thin metal with insulated backing case is essentially the same thing as this one. That's the theory of the plate thermometer.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/firemodels/fds/issues/8474#issuecomment-647819826, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABBWXQNPGXEWZBPOKXZKO73RX7SRDANCNFSM4OESBINA .

-- Sent from my iPhone

[MortezaHaghighi]

@rmcdermo if you still have access to the verification cases for HT3D&PYRO3D I've sent to you a couple of months ago, you may see similar deviations there as well. For example, this one was a gas-filled domain (without any OBST) heated up from the exterior:

next, just inserting an adiabatic OBST in the middle of the domain:

and lastly, only activating the conduction (eps=0) for the OBST:

weird that the red and black lines match in the third one, but not in the other ones.

[rmcdermo]

Using OBST did not solve the problem.

I'm becoming more and more convinced that we are being bitten by the time accuracy of the pressure update. CHECK_HT did not affect the time step much because (a) it does not yet account for radiative flux and (b) it has not yet been applied to thermally thick particles. If you do this, you see much smaller time steps. The energy absorbed by the particle does not change, while the gas phase energy changes drastically.

In this sense, the sealed box case my be sending us in circles. On the other hand, coming to grips with this pressure change may have many other benefits (e.g., tunnels).

[mcgratta]

Before this issue gets more complicated, I'm going to create verification cases for my little metal sphere, one with convection only and one with radiation only. We can then discuss how the radiation only case might be tightened up a bit.

[mcgratta]

My initial problem is solved. We can create a new issue to deal with similar issues. I'll transfer this issue to Randy, and Simo you're off the hook ;). The radiation solver is still rock-solid.

Jason, you can weigh in if you're interested, but I think things are OK for now. Oddly, I cannot unassign you from this issue.

[drjfloyd]

I just took the original case from the issue opening and ran it once with SPEC=NITROGEN (not setting CP) and once with SPEC=SOOT. This means one case where all the radiation from the particle is essentially dealt with in that immediate grid cell since pure soot is extremely optically thick and one case where all the radiation interacts with adiabatic surfaces since N2 doesn't abosrb anything. I added outputs for INTERNAL ENERGY and WALL ENTHALPY. For both cases there was a loss of ~2 W of total energy or a relative error of ~1E-5. I agree that it looks like the issue is solved.

[mcgratta]

What quantities did you compare and what was the difference?

[drjfloyd]

Added these device:

&DEVC XB=0,1,0,1,0,1, QUANTITY='INTERNAL ENERGY', ID='E', SPATIAL_STATISTIC='VOLUME INTEGRAL'/ &DEVC ID='E_surf',INIT_ID='ball', QUANTITY='WALL ENTHALPY' /

The sum is the total energy in the box which was 38.2 kJ for SOOT as the background and 172 kJ for NITROGEN as the background (I didn't set SPECIFIC_HEAT for either). At the end in both cases there was a drop of 0.002 kJ.

[mcgratta]

Pages from FDS_Verification_Guide.pdf

[rmcdermo]

Where do we stand on this issue? Firebot is passing. If there is a verification case that demonstrates a problem, please put it in the verification suite. If not, let's close the issue.