Modelling of two phase flow in packed bed using conservation equations

[casualguitar]

What does $\dot{m}$ vs tank number at various fixed times look like?

Here's the same for n=32. Theres a bit of chance involved in capturing good time increments given how short they are. n=32 does give a better sense of the 'wave' moving along the bed. Here it is:

The trend seems to be that the flow downstream of the saturation zone is much higher than the gas flow upstream of the saturation zone, and that this increased flow becomes less significant as we move downstream

 

[Chestermiller]

Given the spikes in mass flow are narrow, its hard to pick suitable time intervals that show this. However, this is a plot of the first 11 seconds in intervals of 1 second:
View attachment 294981

Its a poor quality plot (as I say difficult to plot when the saturation last <1 second at each position), so I'm working on presenting this clearly.

But yes for now this shows that after 1 second the mass flow in the tank is at a maximum of 2kg/s (much larger than the inlet flow), and this value matches the previous plot of time vs flowrate.

Also at approx 3 seconds we can see from the previous plot the flow is almost 0.25kg/s around the middle to the end of the tank (still much larger than the inlet flow), and this matches the position vs flow plot also

So yes while this plot isn't great it does line top with the last one

Your results do show the expected behavior that the mass flow rate increases monotonically with position through the bed. We ought to be able to get an idea of the amount of time it takes for the liquid phase to be fully cleared from the bed.

I'm still concerned about spikes in the results.

 

[Chestermiller]

Here's the same for n=32. Theres a bit of chance involved in capturing good time increments given how short they are. n=32 does give a better sense of the 'wave' moving along the bed. Here it is:
View attachment 294984
The trend seems to be that the flow downstream of the saturation zone is much higher than the gas flow upstream of the saturation zone, and that this increased flow becomes less significant as we move downstream

This is more like what we expected. It looks like it takes only about 11 seconds to fully clear the liquid from the bed. After that, all you have is single phase vapor.

 

[Chestermiller]

Question: In your numerical scheme, when do you update the $\dot{m}$ variation?
(a) for all tanks simultaneously, at the end of each time step (before the next time step is initiated)
(b) for each tank in sequence, after the time step for the tank under consideration is completed and the same time step for the subsequent tank is to be implemented.

If you do case (b), do you use the average $\dot{m}$ over the time step (exiting the previous tank and entering the subsequent tank) as the input for this subsequent tank.

 

[casualguitar]

Question: In your numerical scheme, when do you update the $\dot{m}$ variation?
(a) for all tanks simultaneously, at the end of each time step (before the next time step is initiated)
(b) for each tank in sequence, after the time step for the tank under consideration is completed and the same time step for the subsequent tank is to be implemented.

If you do case (b), do you use the average $\dot{m}$ over the time step (exiting the previous tank and entering the subsequent tank) as the input for this subsequent tank.

(b), if I understand the differences between the two correctly.

After a given time step, equations 4b and 4a are used to generate the new $\dot{m}$ values for each tank in sequence.

No I use the value entering the subsequent tank, and don't take the average.

To do this, I guess I could use effectively the same calculation flow as I'm currently doing i.e. calculating mdot for each tank in sequence. Except, instead of taking the flow entering the subsequent tank as input to tank n+1, I can take the numeric average (min-mout)/2 as input to tank n+1. Correct? If so this is an easy change

Is it correct to say the reason for taking the numeric average here is because I'm currently evaluating the flux at the boundary, rather than at the tank centre? Taking the average would in effect be evaluating the flux at the centre

 

[Chestermiller]

(b), if I understand the differences between the two correctly.

After a given time step, equations 4b and 4a are used to generate the new $\dot{m}$ values for each tank in sequence.

No I use the value entering the subsequent tank, and don't take the average.

To do this, I guess I could use effectively the same calculation flow as I'm currently doing i.e. calculating mdot for each tank in sequence. Except, instead of taking the flow entering the subsequent tank as input to tank n+1, I can take the numeric average (min-mout)/2 as input to tank n+1. Correct? If so this is an easy change

Is it correct to say the reason for taking the numeric average here is because I'm currently evaluating the flux at the boundary, rather than at the tank centre? Taking the average would in effect be evaluating the flux at the centre

No. For better accuracy, I'm recommending using the time average of the mass flow rate at the inlet to the tank (tank j) over the time interval, not the spatial average from inlet to outlet. $$\dot{m}_{j-1}(t+\Delta t /2)=\frac{(\dot{m}_{j-1}(t+\Delta t)+\dot{m}_{j-1}(t))}{2}$$

 

[casualguitar]

No. For better accuracy, I'm recommending using the time average of the mass flow rate at the inlet to the tank (tank j) over the time interval, not the spatial average from inlet to outlet. $$\dot{m}_{j-1}(t+\Delta t /2)=\frac{(\dot{m}_{j-1}(t+\Delta t)+\dot{m}_{j-1}(t))}{2}$$

Aha I see, time average not spatial. Got it

So to confirm - let's say the mass flow rate array is generated for t=1 (assuming t=0 is the I.C.), are you saying take the average the mass flow rate arrays for t=0 and t=1, and use this new array in the t=2 calculation, rather than the t=1 mass flow rate array. So essentially a 't=0.5' array would be used to calculate t=2. Then the 't=1.5' array (the average of the t=0.5 and t=2 arrays) would be used for t=3, etc.

Is this correct? This is a fast change if so

 

[Chestermiller]

Aha I see, time average not spatial. Got it

So to confirm - let's say the mass flow rate array is generated for t=1 (assuming t=0 is the I.C.), are you saying take the average the mass flow rate arrays for t=0 and t=1, and use this new array in the t=2 calculation, rather than the t=1 mass flow rate array. So essentially a 't=0.5' array would be used to calculate t=2. Then the 't=1.5' array (the average of the t=0.5 and t=2 arrays) would be used for t=3, etc.

Is this correct? This is a fast change if so

No. Suppose you are looking at the time interval between t = 1 and t = 2, and suppose you have just solved tank j-1 for h at time t = 2; you then use the mass balance relationship to calculate m-dot exiting tank j-1 at time t = 2. This is also the mass flow into tank j at time t = 2. For the heat balance in tank j over the interval t = 1 to t = 2, you then use as the inflow to tank j over the time interval t = 1 to t = 2 the m-dot at t = 1 and j-1 averaged with the m-dot at t = 2 and j-1 to solve for h in tank j at time t = 2.

I hope that this makes sense. I feel that this will reduce the spikes problem and give you more accurate results.

 

[Chestermiller]

I'm having second thoughts about switching to this ad hoc integration approach in an effort to remove the unphysical spikes from the numerical solution. I think we should table this ad hod approach for now and stick with the automatic integrator. There are things we can test within the framework of the automatic integrator to make the solution more physically realistic and accurate. I recommend testing the following:
1. Limit the size of the time step that the integrator is allowed to build
2. Limit the order of the integration that the integrator is allowed to build (to nothing higher than 1st order)
3. Specify tighter error control on the integrator.

Please consider trying a few tests with these to see if any or all of them can improved the spike issue.

 

[casualguitar]

I'm having second thoughts about switching to this ad hoc integration approach in an effort to remove the unphysical spikes from the numerical solution. I think we should table this ad hod approach for now and stick with the automatic integrator. There are things we can test within the framework of the automatic integrator to make the solution more physically realistic and accurate. I recommend testing the following:
1. Limit the size of the time step that the integrator is allowed to build
2. Limit the order of the integration that the integrator is allowed to build (to nothing higher than 1st order)
3. Specify tighter error control on the integrator.

Please consider trying a few tests with these to see if any or all of them can improved the spike issue.

So the time step and error limit adjustments did not change the plot/spike issue. I tested some very low time steps and low tolerances but no it did not visibly change the plot. So they don't seem to be the issue if there is one

Limiting the order of the integrator to 1st order doesn't seem to be readily possible with solve_ivp. Looking at the documentation, there are a number of possible integrators (RK45, RK23, LSODA, etc) and the order is specified for all of them. One of these integrators called 'BDF' is of variable order, in that the integrator chooses the most suitable order. It is not possible to manually set it though. That said, I tried all of the available integrators (explicit and implicit) available in solve_ivp and none changed the output, so it seems the integrator is not the issue either.

Is there definitely an issue with spiking? Or maybe this is correct?

Also, what would be the reasoning for specifying first order? Surely we would want higher order, or is this incorrect?

 

[Chestermiller]

So the time step and error limit adjustments did not change the plot/spike issue. I tested some very low time steps and low tolerances but no it did not visibly change the plot. So they don't seem to be the issue if there is one

Limiting the order of the integrator to 1st order doesn't seem to be readily possible with solve_ivp. Looking at the documentation, there are a number of possible integrators (RK45, RK23, LSODA, etc) and the order is specified for all of them. One of these integrators called 'BDF' is of variable order, in that the integrator chooses the most suitable order. It is not possible to manually set it though. That said, I tried all of the available integrators (explicit and implicit) available in solve_ivp and none changed the output, so it seems the integrator is not the issue either.

Is there definitely an issue with spiking? Or maybe this is correct?

Also, what would be the reasoning for specifying first order? Surely we would want higher order, or is this incorrect?

I guess that this means that the spikes are a real feature of the solution?

As far as the order is concerned, specifying lower order would make a solution like this one with discontinuities in the derivatives of enthalpy vs temperature and enthalpy vs density behave better.

 

[casualguitar]

I guess that this means that the spikes are a real feature of the solution?

As far as the order is concerned, specifying lower order would make a solution like this one with discontinuities in the derivatives of enthalpy vs temperature and enthalpy vs density behave better.

We could go the route of manually creating a 1st order integrator as a check, but I think the odds are low that it would help the behaviour, given the number of other integrators I tried

Time to return to the list of remaining tasks?

Theres the addition of temperature dependent parameters, and the U value work, as far as I can see

What do you think? We can ignore both of these for now if there is further core model work to do

 

[Chestermiller]

We could go the route of manually creating a 1st order integrator as a check, but I think the odds are low that it would help the behaviour, given the number of other integrators I tried

Time to return to the list of remaining tasks?

Theres the addition of temperature dependent parameters, and the U value work, as far as I can see

What do you think? We can ignore both of these for now if there is further core model work to do

Now that you have shown that the mass flow rate of liquid downstream of the 2 phase zone is much higher than 0.05 kg/sec., I think it would be worthwhile estimating the pressure variation spatially during the time that this effect is present (i.e., liquid is still present in the bed). Or, do you think that this is not too worthwhile, since it is only at times less than about 10 seconds that this occurs (according to your results)? Thoughts?

 

[casualguitar]

Now that you have shown that the mass flow rate of liquid downstream of the 2 phase zone is much higher than 0.05 kg/sec., I think it would be worthwhile estimating the pressure variation spatially during the time that this effect is present (i.e., liquid is still present in the bed). Or, do you think that this is not too worthwhile, since it is only at times less than about 10 seconds that this occurs (according to your results)? Thoughts?

Ah yes, I think this is worth doing. I can probably get an immediate sense of the max/min pressure variation by using the max/min flows in the tank during that time?

I still think its worth doing even given the short timespan, because some other parameters (the real U value for example) might increase this time.

Also, this is the first configuration, in which we've got liquid in the tank with a higher temperature gas entering. There is also the reverse case of a gas in the tank and a liquid inlet

Anyway yes I think we could get an initial sense of the range of possible pressure drops with the Ergun equation pretty quickly. I could for example use the max flow (about 2kg/s) and assume gas flow (this wouldn't happen I know because of the gas density but it would give an idea of an overly pessimistic pressure drop). If the pressure drop for this case is still insignificant then we can possibly say the pressure drop for all cases is insignificant?

If the above is reasonable, I'll try this Ergun case

 

[Chestermiller]

Ah yes, I think this is worth doing. I can probably get an immediate sense of the max/min pressure variation by using the max/min flows in the tank during that time?

I still think its worth doing even given the short timespan, because some other parameters (the real U value for example) might increase this time.

Also, this is the first configuration, in which we've got liquid in the tank with a higher temperature gas entering. There is also the reverse case of a gas in the tank and a liquid inlet

Anyway yes I think we could get an initial sense of the range of possible pressure drops with the Ergun equation pretty quickly. I could for example use the max flow (about 2kg/s) and assume gas flow (this wouldn't happen I know because of the gas density but it would give an idea of an overly pessimistic pressure drop). If the pressure drop for this case is still insignificant then we can possibly say the pressure drop for all cases is insignificant?

If the above is reasonable, I'll try this Ergun case

I don't exactly follow what you are proposing. You know where you have vapor and where you have liquid, and you can get the pressure drop for each of these regions using the Ergun equation. The only big uncertainty is the region where vapor and liquid are both present, but this is only 1 tank (according to what you have said previously). So tentatively neglect it.

Of course, when liquid is coming in replacing vapor, that needs to still be addressed.

 

[casualguitar]

You know where you have vapor and where you have liquid, and you can get the pressure drop for each of these regions using the Ergun equation. The only big uncertainty is the region where vapor and liquid are both present, but this is only 1 tank (according to what you have said previously). So tentatively neglect it.

I follow, so yes I can get the pressure drop in both regions. Regarding the big uncertainty about where both vapour and liquid are present - will it not be the case though that this region will be present at every position in the bed i.e. moving from inlet to outlet. So, getting the pressure drop for the all vapour case, and the all liquid case should give us the max and min pressure drop that can occur? The pressure would move from the initial pressure drop when the tank is mostly liquid (the lowest pressure drop), up to the max pressure drop when it is all gas. So is it correct to say the saturated region can be ignored, and we can just consider all gas and all liquid?

Also - in the previous Ergun calculation I was able to use the actual length of the bed (about 1.4m). However I am not sure how 'long' each tank is. We know that $L=n\Delta x$, so knowing the total length, say 1.4m, we could decide a value of delta x and this would give us the value of n. Is it correct to say that the length of a tank 'n' in this model is not constant?

 

[Chestermiller]

I follow, so yes I can get the pressure drop in both regions. Regarding the big uncertainty about where both vapour and liquid are present - will it not be the case though that this region will be present at every position in the bed i.e. moving from inlet to outlet. So, getting the pressure drop for the all vapour case, and the all liquid case should give us the max and min pressure drop that can occur? The pressure would move from the initial pressure drop when the tank is mostly liquid (the lowest pressure drop), up to the max pressure drop when it is all gas. So is it correct to say the saturated region can be ignored, and we can just consider all gas and all liquid?

I would say that you are correct, and that these would bound the answer.

Also - in the previous Ergun calculation I was able to use the actual length of the bed (about 1.4m). However I am not sure how 'long' each tank is. We know that $L=n\Delta x$, so knowing the total length, say 1.4m, we could decide a value of delta x and this would give us the value of n. Is it correct to say that the length of a tank 'n' in this model is not constant?

The parameter n is the total number of tanks, so the length of each tank is the same: $$\Delta x=\frac{L}{n}$$Once you specify n, the length of each tank is determined.

 

[casualguitar]

I would say that you are correct, and that these would bound the answer.

Great this seems straightforward then. So, checking these two cases:
1) Where liquid flows through the tank at 2kg/s, for the entire length of the tank
2) Where gas flows through the tank at 0.05kg/s, for the entire length of the tank

The max pressure drop occurs in the liquid case, and returns a value of about 4 bar. In the context of a 30 bar operating pressure, this still seems minimal. Also given that this situation doesn't occur in the model, and cannot occur in reality (the model doesn't show the flow increase to 2 kg/s for long, just a quick spike), the real pressure drop is likely still much smaller.

Safe to say pressure drop is negligible for now?

Temperature dependent parameters and the U value are two options for continuation
Also replacing the temperature(H) and mass(H) functions with real EoSs is an option. Open to ideas

 

[Chestermiller]

Great this seems straightforward then. So, checking these two cases:
1) Where liquid flows through the tank at 2kg/s, for the entire length of the tank
2) Where gas flows through the tank at 0.05kg/s, for the entire length of the tank

The max pressure drop occurs in the liquid case, and returns a value of about 4 bar. In the context of a 30 bar operating pressure, this still seems minimal. Also given that this situation doesn't occur in the model, and cannot occur in reality (the model doesn't show the flow increase to 2 kg/s for long, just a quick spike), the real pressure drop is likely still much smaller.

Safe to say pressure drop is negligible for now?

All the liquid downstream is at 2 kg/sec, not just the spike. And all the vapor upstream is at the higher pressure. How does this effect the VLE? How are the flow and pressure boundary pressures imposed on this bed?

Temperature dependent parameters and the U value are two options for continuation
Also replacing the temperature(H) and mass(H) functions with real EoSs is an option. Open to ideas

If I were you, my next focus would be on U. Added accuracy in the property relationships is secondary (if needed at all) in my judgment.

 

[casualguitar]

All the liquid downstream is at 2 kg/sec, not just the spike. And all the vapor upstream is at the higher pressure. How does this effect the VLE? How are the flow and pressure boundary pressures imposed on this bed?

Whoops, yes I read the plot incorrectly. So yes the max pressure drop across the system could be close to 4 bar.

The VLE: How does the higher pressure vapour affect the VLE, or how does the pressure drop (about 4 bar) affect the VLE in the lower pressure regions? When you say how are the flow and pressure boundary pressures imposed on the bed, what do you mean by this?

If I were you, my next focus would be on U. Added accuracy in the property relationships is secondary (if needed at all) in my judgment.

Sounds good. U can be the next point of focus. I'll do some reading today and get a general plan for this

 

[Chestermiller]

Whoops, yes I read the plot incorrectly. So yes the max pressure drop across the system could be close to 4 bar.

The VLE: How does the higher pressure vapour affect the VLE, or how does the pressure drop (about 4 bar) affect the VLE in the lower pressure regions? When you say how are the flow and pressure boundary pressures imposed on the bed, what do you mean by this?

When the pressure changes in the 2 phase region, this affects the VLE.

I looked over the pressure-enthalpy diagrams of N2 and O2 available online, and can now see what you are driving at with respect to the complexity of the thermodynamic behavior. At a nominal 30 bars, the vapor behavior will definitely be non-ideal. I assume that the data parameterizations you have available apply to the vapor region only? Do you have any information on activity coefficients (or excess free energy) of liquid mixtures of N2 and O2? What is understanding of how close the liquid behavior is to an ideal solution?

With regard to boundary conditions, I assume you are controlling the exit pressure from the bed by some sort of valve arrangement and controlling the vapor flow rate to the bed by continuously adjusting the inlet pressure (since the vapor is compressible).

 

[casualguitar]

I looked over the pressure-enthalpy diagrams of N2 and O2 available online, and can now see what you are driving at with respect to the complexity of the thermodynamic behavior. At a nominal 30 bars, the vapor behavior will definitely be non-ideal.

Yes I was actually just looking at the heat capacity/thermal conductivity graphs of N2/O2 vs temperature and they do vary quite a bit also across our temperature range.

I assume that the data parameterizations you have available apply to the vapor region only?

The data parameterisations apply to both vapour and liquid. There is also flash functionality that can calculate the heat capacity etc in the saturated zone as an average

Do you have any information on activity coefficients (or excess free energy) of liquid mixtures of N2 and O2?

I think so actually. There is functionality in the thermo library to return the liquid activity coefficients:

What is understanding of how close the liquid behavior is to an ideal solution?

Effectively none. I'm not sure how to quantify this. What I could do, is take some things like heat capacity etc and plot the percentage variation between the ideal liquid heat capacity and the real liquid heat capacity versus temperature. Something like this?

 

[Chestermiller]

Yes I was actually just looking at the heat capacity/thermal conductivity graphs of N2/O2 vs temperature and they do vary quite a bit also across our temperature range.

Does this include the effect of pressure on the heat capacity, or is it just the ideal gas heat capacity (which is just the molar average at a given temperature).

The data parameterisations apply to both vapour and liquid. There is also flash functionality that can calculate the heat capacity etc in the saturated zone as an average

I think so actually. There is functionality in the thermo library to return the liquid activity coefficients:

View attachment 295181

Effectively none. I'm not sure how to quantify this. What I could do, is take some things like heat capacity etc and plot the percentage variation between the ideal liquid heat capacity and the real liquid heat capacity versus temperature. Something like this?

I was wondering what the effect of composition on the liquid heat capacity was. There is also the question of the VLE deviation from Raoult's Law (as a result of liquid phase non-ideality).

Check out the chapters in Introduction to Chemical Engineering Thermodynamics by Smith and Van Ness (Chapters 10 and beyond).

 

[casualguitar]

I was wondering what the effect of composition on the liquid heat capacity was. There is also the question of the VLE deviation from Raoult's Law (as a result of liquid phase non-ideality).

Yes I actually started Denbigh yesterday and it seems to have similar stuff, so I might continue with this (the effect of composition on liquid heat capacity, deviation from Raoults law). Is this reasonable?

Does this include the effect of pressure on the heat capacity, or is it just the ideal gas heat capacity (which is just the molar average at a given temperature).

The latter, just the ideal gas heat capacity (function of temperature only). I think I can check the variation of liquid heat capacity with temperature for a given pressure, though.

So the plan is loosely to tackle the non ideality of the vapour and liquid phases? And yes I agree that this will require textbook reading on my part, for me to be able to discuss this (two days or so of reading)

 

[casualguitar]

Check out the chapters in Introduction to Chemical Engineering Thermodynamics by Smith and Van Ness (Chapters 10 and beyond).

Two other thoughts for model development:
- Effect of adding radial heat transfer to the system
- Inclusion of a conduction term (seems to be common in packed bed models), rather than having this encompassed in U

But yes for now I'm just reading

 

[Chestermiller]

Yes I actually started Denbigh yesterday and it seems to have similar stuff, so I might continue with this (the effect of composition on liquid heat capacity, deviation from Raoults law). Is this reasonable?The latter, just the ideal gas heat capacity (function of temperature only). I think I can check the variation of liquid heat capacity with temperature for a given pressure, though.

So the plan is loosely to tackle the non ideality of the vapour and liquid phases? And yes I agree that this will require textbook reading on my part, for me to be able to discuss this (two days or so of reading)

The heat capacity is a function of pressure also (for non-ideal gas behavior). You really need to include the effect of pressure on gas enthalpy.

 

[casualguitar]

The heat capacity is a function of pressure also (for non-ideal gas behavior). You really need to include the effect of pressure on gas enthalpy.

HI Chet, yes there is h(T,P) functionality in the thermo library so I can use that

In regards to the non-ideality of air, yes as you said the vapour is definitely non-ideal, and actually it seems at lower temperatures air is increasingly non-ideal i.e. liquid phase is even more non-ideal. Increasing the pressure results in an increase in non-ideality also. So yes its safe to say vapour and liquid phase air are both non-ideal

The functionality is absolutely there to model this. There is a bit of setup work required in python to get the actual activity coefficients. I'll do that to get an exact idea of how non-ideal air is at various temperatures.

So in terms of implementing non-ideal air, how can we do this in stages? Are these reasonable as a first two?

1) Implement the Ergun equation in the ideal model, to show the pressure decreasing across the system (wont affect any parameter except mass as this is the only place we use P). I have this implemented in a separate script, just need to edit it to find the pressure drop at any point in the system
2) add non ideal correlations for density and heat capacity (as a function of T and P)

 

[Chestermiller]

HI Chet, yes there is h(T,P) functionality in the thermo library so I can use that

In regards to the non-ideality of air, yes as you said the vapour is definitely non-ideal, and actually it seems at lower temperatures air is increasingly non-ideal i.e. liquid phase is even more non-ideal. Increasing the pressure results in an increase in non-ideality also. So yes its safe to say vapour and liquid phase air are both non-ideal

I'm not so sure it is unreasonable to treat the liquid phase as ideal. This would certainly simplify things greatly. I suggest checking the literature relating to the VLE behavior of O2-N2 liquid mixtures and the activity coefficients for such mixtures.

The functionality is absolutely there to model this. There is a bit of setup work required in python to get the actual activity coefficients. I'll do that to get an exact idea of how non-ideal air is at various temperatures.

So in terms of implementing non-ideal air, how can we do this in stages? Are these reasonable as a first two?

There is certainly going to be data on the activity coefficients for O2 and N2 in the vapor phase. Also, maybe the Lewis-Randle rule applies to the vapor. Spend some time digging in the literature.

1) Implement the Ergun equation in the ideal model, to show the pressure decreasing across the system (wont affect any parameter except mass as this is the only place we use P). I have this implemented in a separate script, just need to edit it to find the pressure drop at any point in the system
2) add non ideal correlations for density and heat capacity (as a function of T and P

Once you have the non-ideal PVT behavior of the vapor mixtures, you can, knowing the ideal gas heat capacities of N2 and O2 vapor as a function of temperature, get the non-ideal variations in enthalpy with temperature and pressure.

 

[casualguitar]

I'm not so sure it is unreasonable to treat the liquid phase as ideal. This would certainly simplify things greatly. I suggest checking the literature relating to the VLE behavior of O2-N2 liquid mixtures and the activity coefficients for such mixtures.

Checking the literature for the O2-N2 liquid mixture VLE behaviour doesn't seem to return much regarding the activity coefficients.

To possibly simplify further I have one question - why not use an existing EOS model rather than the activity coefficients to model the vapour and liquid phases? If we used an EOS (one already written in code) it seems we could calculate any T or P dependent property easily, and we could assume non-ideality for the liquid phase because the work it would add would not cause a significant increase in modelling difficulty

That said, although I very much prefer the idea of using existing functionality rather than hard coding my own, I don't want to move too fast on that and lose understanding myself

What are your thoughts on using a prewritten EOS instead of an activity coefficient model?

Also, maybe the Lewis-Randle rule applies to the vapor.

I did read about this rule actually in Denbigh. I'll read over this

Once you have the non-ideal PVT behavior of the vapor mixtures, you can, knowing the ideal gas heat capacities of N2 and O2 vapor as a function of temperature, get the non-ideal variations in enthalpy with temperature and pressure.

So yes back the EOS point, the thermo library in Python has existing functionality to calculate both ideal and non-ideal heat capacities of pure O2/N2 and mixtures of these. It can also get the non-ideal enthalpy(T,P)

https://thermo.readthedocs.io/thermo.phases.html?highlight=mixture enthalpy#module-thermo.phases

I think we could save lots of time, and add more complex functionality if we use this. Do you think it is reasonable to use this functionality rather than hard code our own?

EDIT: Also, I'm working on getting the activity coefficients for air plotted as a function of temperature using this thermo library. I will post this plot once its working

 

[Chestermiller]

Checking the literature for the O2-N2 liquid mixture VLE behaviour doesn't seem to return much regarding the activity coefficients.

To possibly simplify further I have one question - why not use an existing EOS model rather than the activity coefficients to model the vapour and liquid phases? If we used an EOS (one already written in code) it seems we could calculate any T or P dependent property easily, and we could assume non-ideality for the liquid phase because the work it would add would not cause a significant increase in modelling difficulty

That said, although I very much prefer the idea of using existing functionality rather than hard coding my own, I don't want to move too fast on that and lose understanding myself

What are your thoughts on using a prewritten EOS instead of an activity coefficient model?I did read about this rule actually in Denbigh. I'll read over this

So yes back the EOS point, the thermo library in Python has existing functionality to calculate both ideal and non-ideal heat capacities of pure O2/N2 and mixtures of these. It can also get the non-ideal enthalpy(T,P)

https://thermo.readthedocs.io/thermo.phases.html?highlight=mixture enthalpy#module-thermo.phases

I think we could save lots of time, and add more complex functionality if we use this. Do you think it is reasonable to use this functionality rather than hard code our own?

EDIT: Also, I'm working on getting the activity coefficients for air plotted as a function of temperature using this thermo library. I will post this plot once its working

I looked at the reference, but I don't exactly see how they do this. I assume they have an appropriate data base for pure o2 and n2, and somehow mixture parameters. I guess I'll leave it up to you to figure how to work with this software to get the needed thermodynamic functionalities.

 

[casualguitar]

I looked at the reference, but I don't exactly see how they do this. I assume they have an appropriate data base for pure o2 and n2, and somehow mixture parameters. I guess I'll leave it up to you to figure how to work with this software to get the needed thermodynamic functionalities.

Exactly yes there is appropriate data taken from a number of sources (DIPPR, Perrys, etc) to allow for pure and mixture parameter calculations.

I can work with the software absolutely (I have used it for the dew/bubble and Ergun calculations).

I could start with say some plots of air heat capacity and density for our range of temperatures, for a few pressures? To get an idea of how we would expect the model output to change by adding these two in. I'll add them in separately then

 

[Chestermiller]

Exactly yes there is appropriate data taken from a number of sources (DIPPR, Perrys, etc) to allow for pure and mixture parameter calculations.

I can work with the software absolutely (I have used it for the dew/bubble and Ergun calculations).

I could start with say some plots of air heat capacity and density for our range of temperatures, for a few pressures? To get an idea of how we would expect the model output to change by adding these two in. I'll add them in separately then

What we really need is temperature and density (or specific volume) as functions of enthalpy and pressure for an overall mixture of 79% n2 and 21% o2 (including the two phase region). Starting from enthalpy and density as functions of temperature and pressure for such an overall mixture.

 

[casualguitar]

What we really need is temperature and density (or specific volume) as functions of enthalpy and pressure for an overall mixture of 79% n2 and 21% o2 (including the two phase region). Starting from enthalpy and density as functions of temperature and pressure for such an overall mixture.

To confirm, we are looking to get to this?
$$T(H,P)$$ $$\rho(T,P)$$

And we're assuming we will have this available to start:
$$H(T,P)$$ $$\rho(T,P)$$

So we would need to 'invert' these functions to get to the objective functions?

I can check to see if the goal functions already exist just in case

 

[casualguitar]

What we really need is temperature and density (or specific volume) as functions of enthalpy and pressure for an overall mixture of 79% n2 and 21% o2 (including the two phase region). Starting from enthalpy and density as functions of temperature and pressure for such an overall mixture.

Just as a side note - here are the plots mentioned earlier, 79%/21% air for a pressure of 30 bar:

Density:

Heat capacity:

Thermal conductivity (included this in case we eventually split U into components):

 

[casualguitar]

What we really need is temperature and density (or specific volume) as functions of enthalpy and pressure for an overall mixture of 79% n2 and 21% o2 (including the two phase region). Starting from enthalpy and density as functions of temperature and pressure for such an overall mixture.

Ahh ok it just clicked with me now why we need those specifically. To replace the temperature and mass functions we've got currently implemented. Got it

 

[casualguitar]

Ahh ok it just clicked with me now why we need those specifically. To replace the temperature and mass functions we've got currently implemented. Got it

Hi Chet, away for 1 week on annual leave. Will update here once I'm back

 

[casualguitar]

What we really need is temperature and density (or specific volume) as functions of enthalpy and pressure for an overall mixture of 79% n2 and 21% o2 (including the two phase region). Starting from enthalpy and density as functions of temperature and pressure for such an overall mixture.

Hi Chet, returned to work this morning! Looking at the last sentence, you seem to suggest that we can get to H,P formulation from T,P formulation? Is this correct?

I have yet to confirm if H,P formulation is available in the thermo library. Looking into this today

 

[casualguitar]

What we really need is temperature and density (or specific volume) as functions of enthalpy and pressure for an overall mixture of 79% n2 and 21% o2 (including the two phase region). Starting from enthalpy and density as functions of temperature and pressure for such an overall mixture.

Hi Chet, I don't think there are direct PH dependent properties available in thermo

However, if we know the pressure and enthalpy of the system, we could do a PH flash and that will solve the equations to figure out the temperature. We can then query for all the required properties at those specific temperature and pressure values.

PH flash functionality is definitely available in thermo

Is this suitable for us?

 

[Chestermiller]

I guess you can do that. Try it can see if it works.

 

[casualguitar]

I guess you can do that. Try it can see if it works.

Hi Chet, so the temp(H) and mass(H) functions have now been replaced by a flash function that returns the temperature and density at a given pressure and enthalpy (density is then multiplied by volume to get mass).

One point of confusion I have on this is in relation to the heat of vaporisation. Previously you showed me that we could use different correlations for the temperature/mass based on the enthalpy value i.e. if the enthalpy was lower than the heat of vaporisation then we would use the liquid correlation etc. This naturally allowed us to include the heat of vaporisation for the mixed phase and gas phase.

Now that I'm using a PH flash to calculate the temperature and density, I'm wondering if we have taken the heat of vaporisation into account, and if so, how have we done this?

I guess we have sort of 'rounded' this problem by using the flash calculation however I'm not sure. Does PH flash 'account' for the heat of vaporisation?

Some notes: pressure is assumed constant, liquid and gas heat capacity are assumed constant, U correlation has not been added. These are all minor changes that I will add this evening and tomorrow

 

[Chestermiller]

Hi Chet, so the temp(H) and mass(H) functions have now been replaced by a flash function that returns the temperature and density at a given pressure and enthalpy (density is then multiplied by volume to get mass).

One point of confusion I have on this is in relation to the heat of vaporisation. Previously you showed me that we could use different correlations for the temperature/mass based on the enthalpy value i.e. if the enthalpy was lower than the heat of vaporisation then we would use the liquid correlation etc. This naturally allowed us to include the heat of vaporisation for the mixed phase and gas phase.

Now that I'm using a PH flash to calculate the temperature and density, I'm wondering if we have taken the heat of vaporisation into account, and if so, how have we done this?

I guess we have sort of 'rounded' this problem by using the flash calculation however I'm not sure. Does PH flash 'account' for the heat of vaporisation?

Some notes: pressure is assumed constant, liquid and gas heat capacity are assumed constant, U correlation has not been added. These are all minor changes that I will add this evening and tomorrow

Pressure is not constant in a flash calculation, right? And, in a flash calculation, heat of vaporization is implicitly included.

Now that you are adopting this approach, you should do the calculations outside the bed model, and parameterize the relationships between temperature and density vs enthalpy and pressure using analytic fits to the results of the flash calculations. That seems like the easiest thing to do. Then you won't have to be doing flash calculations or VLE calculations within your main model.

 

[casualguitar]

Pressure is not constant in a flash calculation, right? And, in a flash calculation, heat of vaporization is implicitly included.

Yes you're right I think my terminology is poor. I assume when you say flash you mean going from P1 to P2 where P1>P2. What I meant was using the known P and H values to calculate T, density, etc. In the thermo library this is called 'flash', however yes I agree its not really 'flashing' in the real sense. Are we on the same page there? And what is the name for the calculation I'm describing, if not flash?

Now that you are adopting this approach, you should do the calculations outside the bed model, and parameterize the relationships between temperature and density vs enthalpy and pressure using analytic fits to the results of the flash calculations. That seems like the easiest thing to do. Then you won't have to be doing flash calculations or VLE calculations within your main model.

Ah I see. For reference, here is the line of code that I use to calculate the temperature and mass holdup of the air mixture given the enthalpy and pressure:

The setup of the air mixture is a bit more complex. However as you can see its very simple to use the flasher itself to calculate H and P dependent properties. We mentioned earlier getting analytic fits for T(H,P) and density(H,P). Unless I'm misunderstanding, the flash calculation gives us this fit indirectly (and easily) i.e. we in effect have T(H,P) and density(H,P) through the flash calculation. Do you agree with that?

If the above is true, then it seems like its slightly easier to use the above code, rather than extract the analytic relationships T(H,P) and density(H,P)?

I also think it would be useful to post the 'flow' of the code with the added flash calculation now. I'll do this today

 

[casualguitar]

If the above is ok, then I could also remove the $\frac{d\rho}{dH}$ function with an analytic alternative?

So the thermo library does provide property derivatives for mixtures. However, it does not provide $\frac{d\rho}{dH}$ or $\frac{dH}{d\rho}$. One idea I had was to use the chain rule i.e. to multiply two other property derivatives together to get to the one we want, $\frac{d\rho}{dH}$.

For reference, these are the available density derivatives:

And these are the available enthalpy derivatives:

So maybe we could choose two that multiply together to give the objective derivative. To use one of these derivatives, we would need to know what is 'constant' while density and enthalpy change. Do we have any constant properties in our system that are mentioned above? Can we assume that pressure is effectively constant, given the change across the system is small?

For later reference, the list of available derivatives is here: https://thermo.readthedocs.io/thermo.stream.html?highlight=stream#thermo.stream.Stream

 

[Chestermiller]

Yes you're right I think my terminology is poor. I assume when you say flash you mean going from P1 to P2 where P1>P2. What I meant was using the known P and H values to calculate T, density, etc. In the thermo library this is called 'flash', however yes I agree its not really 'flashing' in the real sense. Are we on the same page there? And what is the name for the calculation I'm describing, if not flash?

I think you need to ascertain exactly what this calculation does. Is it a VLE calculation or an actual flash calculation.

Ah I see. For reference, here is the line of code that I use to calculate the temperature and mass holdup of the air mixture given the enthalpy and pressure:
View attachment 296088
View attachment 296089

The setup of the air mixture is a bit more complex. However as you can see its very simple to use the flasher itself to calculate H and P dependent properties. We mentioned earlier getting analytic fits for T(H,P) and density(H,P). Unless I'm misunderstanding, the flash calculation gives us this fit indirectly (and easily) i.e. we in effect have T(H,P) and density(H,P) through the flash calculation. Do you agree with that?

OK, provided we are convinced that it is just a VLE calculation, where you specify the T and P and it calculates everything else, including split. Still, there is an iteration that is going to be involved in specifying H and P and extracting T and density. After all, in the model, you dependent variable is H, not T.

 

[casualguitar]

I think you need to ascertain exactly what this calculation does. Is it a VLE calculation or an actual flash calculation.

A VLE calculation, where P and H are known, and are used to find T and density yes

OK, provided we are convinced that it is just a VLE calculation, where you specify the T and P and it calculates everything else, including split. Still, there is an iteration that is going to be involved in specifying H and P and extracting T and density. After all, in the model, you dependent variable is H, not T.

Exactly, H and P are specified, and T and density are extracted (after some iteration etc). However what I'm saying is that the thermo library takes care of this iteration for us. So all we have to do is supply the library with P and H, and it will do the iteration required to solve for T and density. This is useful because it means there are no large blocks of iteration code to be seen in the main script

So what I can do is replace our previous temperature function here:

With this new function:

The advantage of the second one is that it is analytic, so that single line allows us to take our known H and P values and convert them to T values

The same thing can be done for the mass holdup function, as the VLE calculation done in the thermo library also allows us to get density. I can just multiply by volume to give mass holdup

The final function to be replaced by an analytic one is the $\frac{d\rho}{dH}$ one. As shown in the earlier message above, are there any properties we can assume as constant (while density and enthalpy change)? We can multiply these derivatives together (chain rule) to get to the objective derivative $\frac{d\rho}{dH}$ I think

EDIT: The T(H,P) and density(H,P) functionality is now in place so we have analytic values for T and density at each H and P value. I think replacing the $\frac{d\rho}{dH}$ derivative is a natural next step as it is the only remaining non analytic function

 

[Chestermiller]

A VLE calculation, where P and H are known, and are used to find T and density yes

Exactly, H and P are specified, and T and density are extracted (after some iteration etc). However what I'm saying is that the thermo library takes care of this iteration for us. So all we have to do is supply the library with P and H, and it will do the iteration required to solve for T and density. This is useful because it means there are no large blocks of iteration code to be seen in the main script

So what I can do is replace our previous temperature function here:
View attachment 296093

With this new function:
View attachment 296094

The advantage of the second one is that it is analytic, so that single line allows us to take our known H and P values and convert them to T values

The same thing can be done for the mass holdup function, as the VLE calculation done in the thermo library also allows us to get density. I can just multiply by volume to give mass holdup

The final function to be replaced by an analytic one is the $\frac{d\rho}{dH}$ one. As shown in the earlier message above, are there any properties we can assume as constant (while density and enthalpy change)? We can multiply these derivatives together (chain rule) to get to the objective derivative $\frac{d\rho}{dH}$ I think

EDIT: The T(H,P) and density(H,P) functionality is now in place so we have analytic values for T and density at each H and P value. I think replacing the $\frac{d\rho}{dH}$ derivative is a natural next step as it is the only remaining non analytic function

Maybe you can supply some graphs of T(H,P) and rho(H,P) to examine?

 

[casualguitar]

Maybe you can supply some graphs of T(H,P) and rho(H,P) to examine?

Yep can do

 

[casualguitar]

Maybe you can supply some graphs of T(H,P) and rho(H,P) to examine?

Hi Chet, T(H,P) and rho(H,P) plots as mentioned (at 30 bar for a range of enthalpies):

We can clearly see the phase change zone in the temperature plot (constant temperature).

These plots were produced assuming constant pressure of 30bar. They were also produced with the d(rho)/dH function you derived. As mentioned, we can replace this derivative with an analytic one, as they are available in the thermo library, however its not clear to me which derivatives to chain rule to get d(rho)/dH

What info were you hoping to examine in these plots?

Thanks

 

[Chestermiller]

Hi Chet, T(H,P) and rho(H,P) plots as mentioned (at 30 bar for a range of enthalpies):
View attachment 296198
View attachment 296199

We can clearly see the phase change zone in the temperature plot (constant temperature).

These plots were produced assuming constant pressure of 30bar. They were also produced with the d(rho)/dH function you derived. As mentioned, we can replace this derivative with an analytic one, as they are available in the thermo library, however its not clear to me which derivatives to chain rule to get d(rho)/dH

What info were you hoping to examine in these plots?

Thanks

The temperature should not be constant in the 2 phase region. It should be varying a little with enthalpy. Just for clarification, you are assuming in these calculations an overall mix of 80 mole % N2 and 20 mole % oxygen, right? Also, is it possible to provide the vapor-liquid split (molar or mass) vs enthalpy? Other pressures? The enthalpies presented in the plots are per mole of mixture or per kg of mixture? The densities are molar densities or mass densities?

 

[casualguitar]

The temperature should not be constant in the 2 phase region. It should be varying a little with enthalpy.

Apologies yes it does vary. Zoomed in plot of enthalpy versus temperature for the two phase region, showing the temperature increase of approx 2-3C :

Just for clarification, you are assuming in these calculations an overall mix of 80 mole % N2 and 20 mole % oxygen, right?

Effectively. The exact mole% breakdown is: N2 = 78.08, O2 = 20.95, Ar = 0.97%
The Argon boiling point falls between that of N2 and O2 so our assumption of a single boiling point rather than an envelope is valid still

Also, is it possible to provide the vapor-liquid split (molar or mass) vs enthalpy?

Vapour-liquid split (molar) vs enthalpy (J/mol)

The enthalpies presented in the plots are per mole of mixture or per kg of mixture? The densities are molar densities or mass densities?

Enthalpy units: J/mol
Density units: kg/m3

Yes the units are not of the same basis, however I just picked kg/m3 units for the density plot as its easiest to visualise. Are there preferred units? I can provide either easily

Note also these plots are not packed bed related model output plots just plots of various properties of air as a function of enthalpy and pressure

If there is any other enthalpy/pressure dependent info required at this point I can provide it