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.