firemodels / fds

Fire Dynamics Simulator
https://pages.nist.gov/fds-smv/
Other
664 stars 622 forks source link

Curious ENTHALPY output #5248

Closed mcgratta closed 7 years ago

mcgratta commented 7 years ago

In the following case:

&HEAD CHID = 'divergence_demo' TITLE = 'Demonstration of adding gas to a sealed compartment' /

&MESH IJK=9,1,9, XB=0.0,1.0,-0.05,0.05,0.0,1.0 /

&TIME T_END=1.0, DT=0.01 /

&MISC STRATIFICATION=.FALSE., NOISE=.FALSE. /

&PART ID='ball', SURF_ID='blower', STATIC=.TRUE., PROP_ID='ball object' /
&SURF ID='blower', MASS_FLUX(1)=10.0, SPEC_ID(1)='AIR', TAU_MF(1)=0.0, RADIUS=0.005, GEOMETRY='SPHERICAL', HEAT_TRANSFER_COEFFICIENT=0. /
&INIT PART_ID='ball', XYZ=0.5,0.0,0.5, N_PARTICLES=1 /
&PROP ID='ball object', SMOKEVIEW_ID='SPHERE', SMOKEVIEW_PARAMETERS(1)='D=0.01' /

&SURF ID='WALL', ADIABATIC=.TRUE., DEFAULT=.TRUE. /

&SLCF PBY=0., QUANTITY = 'TEMPERATURE', CELL_CENTERED = .TRUE., VECTOR=.TRUE. /
&SLCF PBY=0., QUANTITY = 'DIVERGENCE', CELL_CENTERED = .TRUE. /
&SLCF PBY=0., QUANTITY = 'DENSITY', CELL_CENTERED = .TRUE. /
&SLCF PBY=0., QUANTITY = 'PRESSURE', CELL_CENTERED = .TRUE. /
&SLCF PBY=0., QUANTITY = 'ENTHALPY', CELL_CENTERED = .TRUE. /

&DEVC ID='Delta m', XB=0,1,-.05,0.05,0,1, QUANTITY='DENSITY', STATISTICS='VOLUME INTEGRAL', RELATIVE=.TRUE. /
&DEVC ID='Delta T', XB=0,1,-.05,0.05,0,1, QUANTITY='TEMPERATURE', STATISTICS='VOLUME MEAN', RELATIVE=.TRUE. /
&DEVC ID='Delta h', XB=0,1,-.05,0.05,0,1, QUANTITY='ENTHALPY', STATISTICS='VOLUME INTEGRAL', RELATIVE=.TRUE. /
&DEVC ID='Delta p', XB=0,1,-.05,0.05,0,1, QUANTITY='PRESSURE', STATISTICS='VOLUME MEAN', RELATIVE=.TRUE. /

&TAIL /

I see an unusual pattern in the ENTHALPY device output. I do not understand why it is not smooth like the other quantities. The final result also does not satisfy delta e = delta h - V delta p

mcgratta commented 7 years ago

Glenn -- there is a minor Smokeview issue here. Run the case I pasted above. Bring up SMV and look at a vector slice of DIVERGENCE. You should see red lines emanating from the center cell, and blue elsewhere, because the center cell is red and all others are blue. But in this mode, SMV is coloring the lines with the cell value. It should color with the average of the two cells on each side of the vector line. In other words, velocity components "live" at cell faces; scalars at cell centers. So the color for ui should be 0.5*(s(i+1)+s_i)

gforney commented 7 years ago

I'll take a look at it.

gforney commented 7 years ago

think I've figured out a fix. There a couple of tricky issues.

  1. my fix will involve additional subscript computations. The tricky part will be to make sure I don't have subscript out of bound errors when a slice (or parts of a slice) is on the edge of a domain. I may need another test case to make sure I'm doing this correctly. we can talk about this when I'm ready.

  2. the 2nd issue may be an opportunity to learn a little git! I'm in the middle of some big changes to improve boundary and slice menus for geometry files (want files that are "logically" slice files to appear in the slice menu - now all geometry data files appear in the boundary file menu. Anyway I've got a lot of changes not committed.

So, the question is can I create a topic branch for this issue, do the fix then merge my fix into the branch that contains my changes for geometry file menus. (I had tried this earlier and it told me I needed to stash my changes )

gforney commented 7 years ago

got 2 in previous comment working. was able to do a git stash to save my uncommitted files then merged changes from a topic branch I created. When I did a git stash apply (to get my uncommited files back) my working files were merged not copied back i.e. the change I made in my topic branch was applied to my working/uncommitted files. pretty cool.

drjfloyd commented 7 years ago

This may be because do not interpolate for temperature but do NINT(T) when doing property lookups. (interpolation was done originally but was removed to save a little cost).

CALL GET_AVERAGE_SPECIFIC_HEAT(ZZ_GET,CPBAR,TMP(II,JJ,KK)) ... IF (IND==47) GAS_PHASE_OUTPUT_RES = RHO(II,JJ,KK)CPBARTMP(II,JJ,KK)*0.001_EB

rmcdermo commented 7 years ago

On the SMV "issue", CELL_CENTERED=T was intended to as a way to see raw scalar values. The application of CC to VECTOR=T was a hack, intended for velocity only.

mcgratta commented 7 years ago

Understood, but there's no reason the CC and VECTOR cannot be combined to indicate the value of the velocity component times the scalar at the cell face. It's just u_i0.5(si+s(i+1))

rmcdermo commented 7 years ago

Also understood... but this is getting pretty confusing. My suggestion is that we create a namelist parameter called FACE_CENTERED=T for this purpose.

mcgratta commented 7 years ago

I don't think we need that. If SMV draws the velocity component at the cell face (V=T and CC=T), it should just take the average of the scalars in the adjacent cells. It has this info already.

gforney commented 7 years ago

so when drawing a cell centered vector slice file, I'll

  1. draw values in each grid cell colored using actual values
  2. overlay these colored grid cells with vectors colored black (vectors centered on cell face - where they are now)
mcgratta commented 7 years ago

Yes

Jason -- as for the ENTHALPY issue. Here's where I'm stuck:

I start with 0.12 kg of air, and add about 0.003 kg. I assume the specific heat of air is about 1 kJ/kg/K. The temperature increases about 2.5 K. I would think that the enthalpy rise would be (m+delta m)cp(T+delta T) - mcpT = 1.2 kJ, but I am getting about 0.12 kJ in my simulation, even when I account for the rounding issue. I'm not sure where I'm off.

mcgratta commented 7 years ago

If I set the specific heat to 1, I get results consistent with my simple estimates. If I do not, I get a decrease in the internal energy when I add air to the box.

drjfloyd commented 7 years ago

Enthalpy of added mass = -79.5 kJ/kg 0.003 kg = -0.25 kJ Enthalpy due to temperature rise = 3 K 1 kJ/kg/K * 0.12 K = 0.37 kJ Total enthalpy rise = -0.25 kJ + 0.37 kJ = 0.12 kJ

mcgratta commented 7 years ago

So why do the values change so much if I set the specific heat to 1?

drjfloyd commented 7 years ago

As is the case is DT=2.98 K and DP=3720 Pa If I set CONSTANT_SPECIFIC_HEAT_RATIO I get DT=3.00 K and DP=3720 Pa

So the end conditions are the same (the CP is not quite the same so a little difference is expected).

With constant specific heat Enthalpy of added mass = 1.01 kJ/kg/K 293.15 K 0.003 kg = 0.93 kJ Enthalpy due to temperature rise = 3 K 1 kJ/kg/K 0.12 kg = 0.37 kJ Total enthalpy rise = 0.93 kJ + 0.37 kJ = 1.3 kJ

mcgratta commented 7 years ago

So if I do not set cp=1, what are the values of de = dh - V dp? If dh=0.12 kJ and V dp=0.37 kJ, then my change in internal energy, de, is negative. How can that be, if I'm adding air to the sealed box?

gforney commented 7 years ago

I changed smokeview so that 1. smokeview colors vectors in cell centered vector files black and 2. by default draws these vectors overlayed on the corresponding slice file. You can uncheck an option in the bounds slice dialog box to just draw the black vectors. This checkbox option was there previously - my recent update is to let you select this option for node and cell centered vector slice file separately.

try smokeview here and let me know if you see any problems. https://drive.google.com/drive/folders/0B_wB1pJL2bFQc1F4cjJWY2duWTA?usp=sharing

drjfloyd commented 7 years ago

U = H - PV H = h * m

With temperature dependent specific heat U0 = h0 m0 + P0 V = -79.5 kJ/kg 0.1196 kg - 101325 Pa 0.1 m^3 = -10142 kJ U1 = h1 m1 + P1 V h1 = h0 + cp dT = -79.6 kJ/kg +1.015 kJ/kg/K 2.98 K = -76.5 kJ/kg P1 = 105042.3 Pa U1 = -76.5 kJ/kg 0.1227 kg - 105042.3 Pa * 0.1 m^3 = -10512 kJ dU = U1 - U0 = -371.604

With cp = 1 U0 = h0 m0 + P0 V = cp T0 m0 + P0 V = -10097.4 P1 = 105065.3 Pa dT = 3.05 K U1 = h1 m1 +P1 V = cp T1 m1 + P1 V = -10470.2 dU = -372.739

mcgratta commented 7 years ago

Pa * m^3 is not a kJ

Pa = kg/(m*s^2)

Pa * m^3 = kg m^2/s^2 = J

Also, do you mean U0 = h0 m0 - P0 V?

drjfloyd commented 7 years ago

Thanks, forgot the 1000 factor, and yes I typed plus and not minus in the formula (the math is correct way)

Temperature dependent: U0=-19.64 kJ U1 = -19.89 kJ dU = -0.25 kJ h dm = -79.5 * 0.00314 = -0.25 kJ = dU

cp=1 U0=24.92 kJ U1=25.84 kJ dU= 0.92 kJ h dm = 293.15 * 0.00314 = 0.92 kJ = dU

In both cases the internal energy change computed using the beginning and end states is equal to the enthalpy addition of the gas.

mcgratta commented 7 years ago

Here's a way to work it

&SPEC ID='AIR', BACKGROUND=.TRUE., REFERENCE_TEMPERATURE=-273.15 /

Now I get results that are consistent with simplified calcs.

Glenn -- Smokeview works too.

Thanks all.