water_evaporation_3

[drjfloyd]

water_evaporation_3 is the only case in that series not meeting its tolerance. A little investigation shows that as we go from complete evaporation to remaining liquid that error increases. I modified the case to have just a little bit of water that evaporates in one timestep and the new gas enthalpy is being correctly computed as is the new gas temperature. I then went and made a very simple test case:

&HEAD CHID='we', TITLE='Water Evaporation Test, SVN $Revision$' / &MESH IJK=4,1,4, XB=0.0,1.301073075,0.0,1.301073075,0.0,1.301073075 / &TIME T_END=.01 DT=0.001/ &MISC TMPA=499.85,GVEC=0,0,0,STRATIFICATION=.FALSE.,NOISE=.FALSE., FREEZE_VELOCITY=.TRUE./ &RADI RADIATION=.FALSE./ &SPEC ID='BACK',MW=28.8,SPECIFIC_HEAT=1.,BACKGROUND=.TRUE./ &SPEC ID='LIQUID', MW=28.8, SPECIFIC_HEAT=1., SPECIFIC_HEAT_LIQUID=4., HEAT_OF_VAPORIZATION=2000., H_V_REFERENCE_TEMPERATURE=-0.15, MELTING_TEMPERATURE=-0.15 VAPORIZATION_TEMPERATURE=99.85, DENSITY_LIQUID=1000./ &INIT PART_ID='water drops', MASS_PER_VOLUME=0.000454041, N_PARTICLES_PER_CELL=1 / &PART ID='water drops', SPEC_ID='LIQUID', STATIC=.TRUE., QUANTITIES(1:2)='PARTICLE DIAMETER','PARTICLE TEMPERATURE', DIAMETER=2., INITIAL_TEMPERATURE=19.85, SAMPLING_FACTOR=5, MONODISPERSE=.TRUE./

At t= 0 • 1 kg of gas and 0.001 kg of liquid • 1 kg x 773 K x 1000 kJ/kg/K = 773 kJ of enthalpy in the gas. • 0.001 kg x (293 K x 4000 kJ/kg/K -2819 kJ/kg [offset to give the correct H_V at the reference temperature]) = -1.647 kJ of enthalpy in the liquid

After one time step the liquid evaporates completely. Therefore, the total energy of the system (gas plus pressure work) must be 773 kJ – 1.647 kJ = 771.353 kJ. Molecular weights are the same and gas enthalpies are the same; therefore, we can easily solve for P and T to give the correct total energy: (P/1000 – 101.325) V + m c T = E. Which gives P = 101157.5 Pa and T = 770.95101 K.

We have a sealed box with no movement, homogenous conditions, and adiabatic surfaces; therefore, there is no velocity at the wall, and the only divergence component is D_LAGRANGIAN.

• dP/dt = -167.525113/0.001 = -167525.113 = D_LAGRANGIAN / PSUM. • PSUM = 1/PBAR - RTRM = 1/101325 + V/771353 = 7.01393E-6. • D_LAGRANGIAN must be -1.17500969

In our case in each grid cell we have 0.001/16 kg of liquid and 1/16 kg of gas. D_LAGRANGIAN(II,JJ,KK) = (MW_RATIO_M_VAP/M_GAS + (M_VAP_DELTA_H_G - Q_CON_GAS)/H_G_OLD)* WGT / DT

MW_RATIO_M_VAP_WGT /M_GAS /DT = 28.8/28.8*0.001/16 / (1/16)/0.001 = 1

M_VAP_WGT_DELTA_H_G/H_G_OLD/DT = 0.001/16 * (1000(773-293))/(1/16 \ 773) / 0.001 = -0.620957

Since the drop completely evaporates this means that the convective loss is equal to energy of vaporization.

-Q_CON_GAS WGT / H_G_OLD / DT = -(0.001/16 * 1940000 [H_V])/(1/16 * 773) / 0.001 = -2.50970246

D_LAGRANGIAN = 1 -0.620957 -2.50970246 = -2.1306605 /= -1.17500969

Something seems to be wrong with our D_LAGRANGIAN term, but I am not sure what it is.

I feel like I must be missing something but not sure what it is.

[rmcdermo]

Your analysis so far looks like a forward Euler step. FDS will use RK2 and the PBAR will be different each substep. Not sure it is significant off hand.

[rmcdermo]

Jason, Quick question: Haver you ever noticed the kink in the plot of relative humidity? It may not be related to the tolerance problem but I cannot think of any reason for it to be there. I tried locking the time step, but that did not fix it. water_evaporation_3_humidity.pdf

[rmcdermo]

Interesting, the image attached. I wonder whether a plain text file will attach... nope. What if I change the extension...

[rmcdermo]

Yes, this works. So, we can upload input files by just adding .png to the filename and retrieve with right-click Save Link As... and renaming the hash that gets put in your Downloads folder.

[drjfloyd]

I haven't noticed the kink before. Not to say it wasn't there before. Is there a kink in the 6.2 release?

Did the very simple case give you correct results? For me I got the wrong pressure and temperature.

[rmcdermo]

"Very simple" for you means is complex for me. I have an old script that was set up to look at evap 3 in Matlab. That is what I am working with at the moment to try to understand what the correct solution should be. I got pretty confused going back in time when the enthalpy outputs jumped from 10 kJ to 1000 kJ, so it is hard to even pinpoint which commit broke the case.

On Thu, Aug 13, 2015 at 10:44 AM, Jason Floyd notifications@github.com wrote:

I haven't noticed the kink before. Not to say it wasn't there before. Is there a kink in the 6.2 release?

Did the very simple case give you correct results? For me I got the wrong pressure and temperature.

— Reply to this email directly or view it on GitHub https://github.com/firemodels/fds-smv/issues/2610#issuecomment-130710917 .

[drjfloyd]

The jump was from renorming the enthalpy so H(298.15 K) = H_F to make us consistent with the combustion codes.

[rmcdermo]

Looks like the kink was there in the 6.2.0 documentation. But h_gas was not very tight. Looks like we never had a great handle on this one.

[rmcdermo]

I'm no closer to understanding this. Every time I look at evap stuff I get dizzy. But I have an observation and a question. The observation is that if I set SPECIFIC_HEAT=1 for both WATER VAPOR and AIR then I get saturations above 100% relative humidity for water_evaporation_3, something like 104%. Either I've done something wrong in my setup or we still are not out of the woods with over saturation.

The question has to do with the functions WATER_VAPOR_MASS_FRACTION and RELATIVE_HUMIDITY, which seem to be inverses of each other. Yet, they use different definitions for H_V. And this definition seems to be different than the one in the liquid JANAF table, which is just a constant value of H_V. Should these H_V definitions be made the same?

[drjfloyd]

The value in data.f90 is not a constant. It is the value at the H_V reference temperature (see 14.3.1 in the User's Guide). This fixed point is used to adjust the data processing of the thermophysical properties to try and ensure consistency between the gas and liquid phases.

R_H and W_V_M_F are different from data.f90 due to how we process inputs. We allow the user to define HUMIDITY which is used to define the water vapor mass fraction in the AIR lumped species. At the point in time where we are constructing the species arrays we haven't processed the thermophysical property tables to have H_V as a function of temperature but we need to be able deal with the HUMIDITY input in order to know what mass fraction of water is needed. Hence a fit of the H_V that results from data.f90. I suppose ideally we should before reading misc, create arrays outside of SPEC to hold water vapor data so we can use the same fit (i.e. a duplication of the PROC_SPEC and PROC_PART routines somewhere before calling READ_MISC).

WATER_VAPOR_MASS_FRACTION should be the same as RELATIVE_HUMIDTY, that was an error on my part when I updated RELATIVE_HUMIDTY but not WATER_VAPOR_MASS_FRACTION.

If you set only the SPECIFIC_HEAT for water vapor but no other properties, then you will not have consistency between the liquid and gas phases. By definition H_V = H_G - H_L. We peg H_L(0) to get the H_V at the reference temperature. If you have changed H_G(T) from the NASA fit to H_G(T)=T but done nothing to adjust the default values H_V or H_L, then H_V(T) may have odd behavior. This could be causing garbage in-garbage out issues in the evaporation routine. This is why I created the simple input at the start of the issue where I set both liquid and gas properties for the liquid.

[rmcdermo]

OK, thanks, that answers another question I had. I get drastically different enthalpies if I set CPs on the SPEC line. I think we need to automatically make the liquid values consistent if SPECIFIC_HEAT is set. I would not expect to have to do this as a user.

[drjfloyd]

I don't think we should be doing that. How is FDS to know what assumptions to make to force inconsistent properties to be consistent?

I am planning to expand the writeup in the user's guide to provide more detail on specifying properties and how those properties are used as well as try and capture input errors better (water_evaporation_3 sidetracked me from that).

[rmcdermo]

Took a look at the equation of state near line 2387 in part. The EOS is not exactly satisfied. But it's interesting, if we reset TMP (and TMP_G_NEW) using the current PBAR then the solution does not change much. But we we leave TMP and instead set PBAR then the solution changes quite a bit.

[drjfloyd]

Does this suggest that a solution is we track the updated temperature in a WORK array in case there are multiple drops in a cell but otherwise pass a Q term like we do for combustion?

[rmcdermo]

I think the first step is for us to carefully write down what qdot_b''' and mdot_b''' are supposed to be in the FDS Tech Guide. I think this my fault from when I reworked the divergence section. We say that these terms are spelled out in Chapter 8, but they really are not.

I have a suspicion that something to do with h_ref is going to show up and that we are not really accounting for this correctly.

If you look at the current form of D_LAGRANGIAN and compare this with Eq. (3.22) in the tech guide, my interpretation is that DELTA_H_G accounts for a term in q_b''' and also the -h_s,alpha term that multiplies mdot_b'''. But we are gathering the enthalpies with "sensible enthalpy". I worry that h_ref is not consistent somehow.

[drjfloyd]

Thinking in parallel to combustion - we have a term for the molecular weight change (which I think we have correct) and a term for the change in enthalpy. For water evaporation the net change in enthalpy is the vapor enthalpy at the initial drop temperature (as this is the H_V we are using) minus any convection to the drop. I'll give this a try and see.

[rmcdermo]

Jason, I have another clue. Hopefully it points to something and is not a diversion. I am bothered by the DELTA_H_G term in D_LAGRANGIAN because this is not how it is written in the formulation of the divergence expression in the tech guide. See Eq. (3.22), for example. DELTA_HG should be (-h{s,alpha}) in (3.22). To me, that means DELTA_H_G should be -H_S_B. (There should be no H_S there, plus the sign should be reversed). If I make this change, then evap_3 is fixed. See attached. Unfortunately, now a few other evap plots are broken, for example, press and temp (but not enthalpy or density) in evap_1. The worst is evap_6. Well, maybe these other solutions need to be revisited(?). I have not made changes to the code because of these other cases. But I think this is the term, DELTA_H_G, that needs to be addressed. Gotta go, now I'll be on vacation for a week. Good luck. R

On Fri, Aug 14, 2015 at 11:42 AM, Jason Floyd notifications@github.com wrote:

Thinking in parallel to combustion - we have a term for the molecular weight change (which I think we have correct) and a term for the change in enthalpy. For water evaporation the net change in enthalpy is the vapor enthalpy at the initial drop temperature (as this is the H_V we are using) minus any convection to the drop. I'll give this a try and see.

— Reply to this email directly or view it on GitHub https://github.com/firemodels/fds-smv/issues/2610#issuecomment-131154818 .

[drjfloyd]

haha, I was just about to post that I changed DELTA_H_G to +H_S_B (yes + not -) in the simple case above and it worked for the simple case but did not for the water evap cases. I agree that this indicates the problem is in the DELTA_H_G term.

[drjfloyd]

I looked at two extremes of the simple case above. One where I made the liquid enthalpy and H_V very low (delta H term in D_LAGRANGIAN goes to 0) and one where I made the liquid MW very high (molar expansion term in D_LAGRANGIAN goes to 0). Neither case obtains the correct end state. In the first case, reduce c and H_V has no effect on the current calculation of the molar expansion term suggesting we are missing something in our concept of what the delta H term is.

[rmcdermo]

I think the question is whether or not h_v should be included in h_s,w. In my mind it should be, because h_s,w is the enthalpy of the gas phase vapor crossing the boundary into the gas phase domain. q_b should account for the heat transfer that vaporizes the water.

On Tuesday, August 18, 2015, Jason Floyd notifications@github.com wrote:

I looked at two extremes of the simple case above. One where I made the liquid enthalpy and H_V very low (delta H term in D_LAGRANGIAN goes to 0) and one where I made the liquid MW very high (molar expansion term in D_LAGRANGIAN goes to 0). Neither case obtains the correct end state. In the first case, reduce c and H_V has no effect on the current calculation of the molar expansion term suggesting we are missing something in our concept of what the delta H term is.

— Reply to this email directly or view it on GitHub https://github.com/firemodels/fds-smv/issues/2610#issuecomment-132321354 .

Sent from my iPhone

[drjfloyd]

Yes, I agree. I am just having a difficult time getting some formulation where I get a D_LAGRANGIAN that is what it needs to be to get the right pressure. I note that the insentropic cases show similar issues. isentropic gets the correct answer, but isentropic2 which has a delta H term has a much larger error (within tolerance but the tolerance is somewhat large).

In the simple example above if I force D_LAGRANGIAN to be what I compute by hand, all is good. If I assume the molar expansion term is correct, I have yet to find the combination of enthalpy variables that gives the remaining value needed.

[rmcdermo]

Here is an idea. I can't get to this right now, but it might be interesting to derive the divergence constraint starting from the internal energy transport equation instead of enthalpy. If you can find a copy of Poinsot and Veynante, Theoretical and Numerical Combustion, I think this equation is spelled out.

Kevin, this book is on my shelf. For the simple case, if we are missing a term it should show up there.

On Wednesday, August 19, 2015, Jason Floyd notifications@github.com wrote:

Yes, I agree. I am just having a difficult time getting some formulation where I get a D_LAGRANGIAN that is what it needs to be to get the right pressure. I note that the insentropic cases show similar issues. isentropic gets the correct answer, but isentropic2 which has a delta H term has a much larger error (within tolerance but the tolerance is somewhat large).

In the simple example above if I force D_LAGRANGIAN to be what I compute by hand, all is good. If I assume the molar expansion term is correct, I have yet to find the combination of enthalpy variables that gives the remaining value needed.

— Reply to this email directly or view it on GitHub https://github.com/firemodels/fds-smv/issues/2610#issuecomment-132594135 .

Sent from my iPhone

[mcgratta]

I don't see a divergence expression spelled out. I see various forms of the energy conservation equations. Do you know where in the book they discuss this?

[rmcdermo]

They give the internal energy transport equation. You have to factor out the divergence. I'll have a look Monday.

On Friday, August 21, 2015, Kevin McGrattan notifications@github.com wrote:

I don't see a divergence expression spelled out. I see various forms of the energy conservation equations. Do you know where in the book they discuss this?

— Reply to this email directly or view it on GitHub https://github.com/firemodels/fds-smv/issues/2610#issuecomment-133504881 .

Sent from my iPhone

[rmcdermo]

Testing markdown support...

code snippet

[rmcdermo]

Trying a big chunk of code...

            ! Limit gas temperature change

            IF (ABS(TMP_G_NEW/TMP_G - 1._EB) > 0.05_EB) THEN
               DT_SUBSTEP = DT_SUBSTEP * 0.5_EB
               N_SUBSTEPS = NINT(DT/DT_SUBSTEP)
               IF (DT_SUBSTEP <= 0.00001_EB*DT) THEN
                  CALL SHUTDOWN('Numerical instability in particle energy transport, TMP_G')
               ENDIF
               CYCLE TIME_ITERATION_LOOP
            ENDIF

            ! Update gas cell density, temperature, and mass fractions

            RHO(II,JJ,KK) = M_GAS_NEW*RVC
            ZZ(II,JJ,KK,1:N_TRACKED_SPECIES) = ZZ_GET2(1:N_TRACKED_SPECIES)
            CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET2,RSUM(II,JJ,KK))
            TMP(II,JJ,KK) = TMP_G_NEW

            ! Compute contribution to the divergence

            CALL GET_SPECIFIC_HEAT(ZZ_GET,CP,TMP_G)
            H_G_OLD = CP * TMP_G * M_GAS
            ZZ_GET = 0._EB
            ZZ_GET(Z_INDEX) = 1._EB
            CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,H_S_B,TMP_DROP)
            ! CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,H_S,TMP_G)
            ! DELTA_H_G = H_S_B - H_S
            DELTA_H_G = -H_S_B
            D_LAGRANGIAN(II,JJ,KK) = D_LAGRANGIAN(II,JJ,KK) &
                                   + (MW_RATIO*M_VAP/M_GAS + (M_VAP*DELTA_H_G - Q_CON_GAS)/H_G_OLD) * WGT / DT
            CALC_D_LAGRANGIAN = .TRUE.

            ! Add energy losses and gains to overall energy budget array

            Q_DOT(7,NM) = Q_DOT(7,NM) - (Q_CON_GAS + Q_CON_WALL + Q_RAD)*WGT/DT  ! Q_PART
            Q_DOT(3,NM) = Q_DOT(3,NM) + M_VAP*H_S_B*WGT/DT                       ! Q_CONV
            Q_DOT(2,NM) = Q_DOT(2,NM) + Q_RAD*WGT/DT                             ! Q_RADI
            Q_DOT(4,NM) = Q_DOT(4,NM) + Q_CON_WALL*WGT/DT                        ! Q_COND

            IF (LPC%Z_INDEX==I_FUEL .AND. I_FUEL>0) M_DOT(2,NM) = M_DOT(2,NM) + WGT*M_VAP/DT/LPC%ADJUST_EVAPORATION
            M_DOT(4,NM) = M_DOT(4,NM) + WGT*M_VAP/DT  ! Total mass loss rate

            ! Keep track of total mass evaporated in cell

            MVAP_TOT(II,JJ,KK) = MVAP_TOT(II,JJ,KK) + WGT*M_VAP

            ! Update PARTICLE quantities

            LP%ONE_D%X(1)   = (M_DROP/FTPR)**ONTH
            LP%ONE_D%LAYER_THICKNESS(1) = LP%ONE_D%X(1)
            LP%ONE_D%TMP(1) = TMP_DROP_NEW
            LP%ONE_D%TMP_F  = TMP_DROP_NEW
            LP%MASS = M_DROP

[rmcdermo]

We are trying to figure out how to force a vertical scroll bar to that it may be easier to add FDS input files to Issues. So far, no luck.

[drjfloyd]

now working again (my error in doing the analytic solution)