Modelling of two phase flow in packed bed using conservation equations

[Chestermiller]

Hi Chet, I had a separate report to finish today so I'm starting this now.

Yes to all this, except the solution vector is $h1,TS1,m1,etc$ where m1 is the flowrate. The way I've set up the mass holdup and mass flow does look suspect to me. I don't think this affects the solution but I could remove this from the solution vector in case I'm doing something wrong there in the setup

Yes to all of this, as far as I know I'm following this exact scheme unless the error is that I've coded this scheme incorrectly

Instead of this I'm just taking the solid energy balance, and generating the fluid temperature vector afterwards. The Ts array is known so I generate the dTs/dt array also and get the fluid temperature array from this. This image shows the flow for that:

View attachment 293952

I'll go for removing m from the solution array first

You are definitely not going to be integrating $dm_i/dt$ with respect to time. Once $dh_i/dt$ is known, $dm_i/dt$ can be evaluated and used to get $\dot{m}_i$ from $$\dot{m}_i=\dot{m}_{i-1}-\frac{dm_i}{dt}$$I guess that's what you do.

 

[Chestermiller]

Hi Chet, I had a separate report to finish today so I'm starting this now.

Yes to all this, except the solution vector is $h1,TS1,m1,etc$ where m1 is the flowrate. The way I've set up the mass holdup and mass flow does look suspect to me. I don't think this affects the solution but I could remove this from the solution vector in case I'm doing something wrong there in the setup

Yes to all of this, as far as I know I'm following this exact scheme unless the error is that I've coded this scheme incorrectly

Instead of this I'm just taking the solid energy balance, and generating the fluid temperature vector afterwards. The Ts array is known so I generate the dTs/dt array also and get the fluid temperature array from this. This image shows the flow for that:

View attachment 293952

I'll go for removing m from the solution array first

The fluid temperature and bed temperature must be solved simultaneously, not one and then the other.

 

[casualguitar]

I have a feeling that the thing that is causing your numerical problems is the overall mass balance for each tank.

Got it. Its the way I'm regenerating the derivatives. Solved now!:

So this leaves us at the point of having a working spatial variation model. We had discussed your $M_i$ model, however I haven't implemented this. Was this #M_i## model made purely to create a stepping stone between the single and multi tank model, or is there other reason to implement it?

If so I can implement this first thing tomorrow

 

[casualguitar]

The fluid temperature and bed temperature must be solved simultaneously, not one and then the other.

To answer this also, yes I am actually solving them simultaneously, but the way I've set up solve_ivp I don't have T_f in the solution array so I can't easily extract it. Thats why I regenerate it after

 

[casualguitar]

You are definitely not going to be integrating $dm_i/dt$ with respect to time. Once $dh_i/dt$ is known, $dm_i/dt$ can be evaluated and used to get $\dot{m}_i$ from $$\dot{m}_i=\dot{m}_{i-1}-\frac{dm_i}{dt}$$I guess that's what you do.

Missed this post. Yes correct I don't integrate it its just evaluated as part of the loop

 

[Chestermiller]

Got it. Its the way I'm regenerating the derivatives. Solved now!:
View attachment 293956

Excellent. Now how about plotting the temperature vs time for all three tanks, and also showing on the same plot the result from the single tank model for the overall bed.

So this leaves us at the point of having a working spatial variation model. We had discussed your $M_i$ model, however I haven't implemented this. Was this #M_i## model made purely to create a stepping stone between the single and multi tank model, or is there other reason to implement it?

No. It is mathematically equivalent to the multi-tank model, and I was offering it as a simpler alternative to program. Now that the existing multi-tank model is working, it would only be redundant and unnecessary.

How about next doing the liquid-only calculation with many more tanks and doing the 3 tank version with the full temperature range with inlet at 300 K?

 

[casualguitar]

Now how about plotting the temperature vs time for all three tanks, and also showing on the same plot the result from the single tank model for the overall bed.

On this now

How about next doing the liquid-only calculation with many more tanks and doing the 3 tank version with the full temperature range with inlet at 300 K?

Got it so the first model just involves increasing n (to 32 for example)?

The second model just involves changing the inlet temperature to 300. Actually I will make another small change - right now I'm inefficiently manually calculating the inlet enthalpy so I'll just create a function for that

I'll try find a good way to blend plots from separate files into the same plot, because this sounds like something I'll need again (rather than going the other route of copying/pasting into the same file and changing variable names)

 

[casualguitar]

How about next doing the liquid-only calculation with many more tanks

Got caught up with some report writing this evening. Back in action now

Heres the liquid-only plot with more tanks (n=32). Plotting the temperature profile for the first and last tank here:

I checked this against the single tank plot for a 32x flowrate and they do match

and doing the 3 tank version with the full temperature range with inlet at 300 K

I was about to post and noticed this plot was incorrect. I took out some gas functionality to simplify the code so I could spot the error. Its a quick job to put this back in so I'll do this first thing tomorrow

Now how about plotting the temperature vs time for all three tanks, and also showing on the same plot the result from the single tank model for the overall bed.

Looking into a good way to do this overlay plot still. For now here it is as two separate plots:

Temperature vs time for all three tanks (liquid phase only):

So now we seem to have a working model for a single and multi-tank setup. I can think roughly about what the future models would hopefully look like, however to get there, what do you think are the next best steps?

So some simplifications we've made:
- The enthalpy/mass holdup relations. In previous posts you mentioned we would stick to correlations like these. If its an option then we can also use H(T,P) methods already made in other Python libraries i.e. just use a property calculation library to calculate enthalpy
- The heat of vaporisation. Right now we're leaving this as a constant. I have methods again from thermo libraries that treat this as a function of temperatures so we have the option to use this
- The saturation conditions. Right now we've got Tsat defined at the top as a constant. I have Tsat(P) functionality so I could use this no bother to get the saturation temperature
- The quality. Right now we get around the quality with the saturation zone relations
- The fluid. So assuming a pure fluid let's us have one saturation temperature. Air is a mixture so we'll have a bubble/dew point. This definitely sounds difficult to implement. That said I have a lot of useful functionality from the thermo property libraries like dew/bubble point temperature calculations for mixtures, Pressure-Enthalpy flash functionality, really anything we would need. So the important bit here is the actual algorithm, the functionality is available
- Lastly obviously we have constant heat capacity, etc. This can easily be switched out for temperature/pressure dependent values later on but yes its probably best to leave them as constant for now to avoid cluttering etc
- The U value. I've put a value of 1 (out of thin air). We could improve on that

From the above it seems that it might be best to stick with a pure fluid until we get to a late stage model, and then switch to a mixture. Developing the pure fluid model to late-stage would also give us the ability to test pure fluids in the packed bed

So yes, in terms of the next step or two, what do you think is best?

 

[Chestermiller]

Got caught up with some report writing this evening. Back in action now

Heres the liquid-only plot with more tanks (n=32). Plotting the temperature profile for the first and last tank here:
View attachment 293994
I checked this against the single tank plot for a 32x flowrate and they do match

What happened to tank 32? This seems to be only tanks 1 and 31.

I was about to post and noticed this plot was incorrect. I took out some gas functionality to simplify the code so I could spot the error. Its a quick job to put this back in so I'll do this first thing tomorrow

Looking into a good way to do this overlay plot still. For now here it is as two separate plots:

Temperature vs time for all three tanks (liquid phase only):
View attachment 293997

What happened to the 3rd tank? This looks like only tanks 1 and 2. Also, I was hoping you would show on this same plot for comparison the results from the single tank model for 1X the throughput rate.

View attachment 293998

So now we seem to have a working model for a single and multi-tank setup. I can think roughly about what the future models would hopefully look like, however to get there, what do you think are the next best steps?

So some simplifications we've made:
- The enthalpy/mass holdup relations. In previous posts you mentioned we would stick to correlations like these. If its an option then we can also use H(T,P) methods already made in other Python libraries i.e. just use a property calculation library to calculate enthalpy
- The heat of vaporisation. Right now we're leaving this as a constant. I have methods again from thermo libraries that treat this as a function of temperatures so we have the option to use this

Are you saying that you want to be more precise on the representation of density and temperature as a function of enthalpy at a given pressure? Or, are you saying that you want to allow the pressure to vary?

- The saturation conditions. Right now we've got Tsat defined at the top as a constant. I have Tsat(P) functionality so I could use this no bother to get the saturation temperature
- The quality. Right now we get around the quality with the saturation zone relations
- The fluid. So assuming a pure fluid let's us have one saturation temperature. Air is a mixture so we'll have a bubble/dew point. This definitely sounds difficult to implement. That said I have a lot of useful functionality from the thermo property libraries like dew/bubble point temperature calculations for mixtures, Pressure-Enthalpy flash functionality, really anything we would need. So the important bit here is the actual algorithm, the functionality is available

This is very easy to change so that, given a specified pressure, we can express density and temperature of the mixture as a function of enthalpy.

- Lastly obviously we have constant heat capacity, etc. This can easily be switched out for temperature/pressure dependent values later on but yes its probably best to leave them as constant for now to avoid cluttering etc
- The U value. I've put a value of 1 (out of thin air). We could improve on that

From the results so far, my guess is that UA is much higher than we are using.

From the above it seems that it might be best to stick with a pure fluid until we get to a late stage model, and then switch to a mixture. Developing the pure fluid model to late-stage would also give us the ability to test pure fluids in the packed bed

So yes, in terms of the next step or two, what do you think is best?

I agree. I would like to see you put the pure fluid model through its paces. Things I suggest are,

1. in addition to temperature vs time at the various spatial positions, also prepare plots of temperature vs position at a selection of times. The latter will illustrate how the wave is traveling through the bed.

2. corresponding plots for the solid bed

3. Varying the parameters for the model like mass flow rate to the bed and UA (higher UA would be particularly interesting)

Here is a calculation that I would like to see you do to get your feet wet on the mixture model: I have 1 mole of air (80% nitrogen, 20% oxygen) at 2 bars. What is the dew point and what is the bubble point? At a temperature half-way between the dew point and bubble point, what is the molar split of liquid and vapor, the mole fractions in the liquid phase, the mole fractions in the vapor phase, the density of the liquid phase, and the density of the vapor phase?

 

[casualguitar]

What happened to tank 32? This seems to be only tanks 1 and 31.

What happened to the 3rd tank? This looks like only tanks 1 and 2.

Yes because its zero indexed its $tank number = position number +1$ so that is tank 32 and tank 3 respectively

Also, I was hoping you would show on this same plot for comparison the results from the single tank model for 1X the throughput rate.

Working on a function to merge two plots now. Ok so we're saying the 3 tank plot for 0.05kg/s and the single tank plot for 0.05kg/s on the same plot? So we're expecting the last tank of the three tank plot should match the single tank plot

Are you saying that you want to be more precise on the representation of density and temperature as a function of enthalpy at a given pressure? Or, are you saying that you want to allow the pressure to vary?

Yes to the pressure variation. Maybe not the most important next step but it would be nice to implement. Yes the functionality is there to have T(H,P) or something similar, and for representing density, there is also rho(T,P) or something like rho(H(T,P)). I'm not sure on which to use but we can use whichever is most convenient because I have property calculators for this

This is very easy to change so that, given a specified pressure, we can express density and temperature of the mixture as a function of enthalpy.

Interesting so I can see how specifying pressure will get us Tsat and quality. I guess you're saying this for a pure fluid though? Actually though that's right I think knowing P and H at a point (and pure fluid properties) is enough to use a PH Flash calculation so we could get the composition of the liquid/gas, assuming P and H are known. This PH flash calculation functionality is also available and I have actually used it. However seeing your suggestion below to implement a simpler version of this I think is a good idea because I only have a vague idea of how that works

1. in addition to temperature vs time at the various spatial positions, also prepare plots of temperature vs position at a selection of times. The latter will illustrate how the wave is traveling through the bed.

Can do

2. corresponding plots for the solid bed

Can do

3. Varying the parameters for the model like mass flow rate to the bed and UA (higher UA would be particularly interesting)

I'll plot a range of UA values on the same plot

Here is a calculation that I would like to see you do to get your feet wet on the mixture model: I have 1 mole of air (80% nitrogen, 20% oxygen) at 2 bars. What is the dew point and what is the bubble point? At a temperature half-way between the dew point and bubble point, what is the molar split of liquid and vapor, the mole fractions in the liquid phase, the mole fractions in the vapor phase, the density of the liquid phase, and the density of the vapor phase?

I can do this very quickly in code as I have the functionality to set up an 80/20 N2/O2 mixture. I also have dew and bubble point methods set up. I can try to do this manually though because I don't know how it works. Is there any resource in particular that would give pointers for calculating dew/bubble point temperatures for mixtures? For the second bit I can use Raoults law

 

[Chestermiller]

Yes because its zero indexed its $tank number = position number +1$ so that is tank 32 and tank 3 respectively

On your plot, position 0 looks like nothing happened.

Working on a function to merge two plots now. Ok so we're saying the 3 tank plot for 0.05kg/s and the single tank plot for 0.05kg/s on the same plot?

Yes.

So we're expecting the last tank of the three tank plot should match the single tank plot

Not exactly. There should be some similarities, although the 1 tank is inherently more disperse than the 3rd tank of the three tank.

Yes to the pressure variation. Maybe not the most important next step but it would be nice to implement. Yes the functionality is there to have T(H,P) or something similar, and for representing density, there is also rho(T,P) or something like rho(H(T,P)). I'm not sure on which to use but we can use whichever is most convenient because I have property calculators for this

Interesting so I can see how specifying pressure will get us Tsat and quality. I guess you're saying this for a pure fluid though? Actually though that's right I think knowing P and H at a point (and pure fluid properties) is enough to use a PH Flash calculation so we could get the composition of the liquid/gas, assuming P and H are known. This PH flash calculation functionality is also available and I have actually used it. However seeing your suggestion below to implement a simpler version of this I think is a good idea because I only have a vague idea of how that works

Can do

Can do

I'll plot a range of UA values on the same plot

I can do this very quickly in code as I have the functionality to set up an 80/20 N2/O2 mixture. I also have dew and bubble point methods set up. I can try to do this manually though because I don't know how it works. Is there any resource in particular that would give pointers for calculating dew/bubble point temperatures for mixtures? For the second bit I can use Raoults law

I would start out by using Raoult's Law on all cases, unless you know more details about the non-iddalalities of the liquid or if you have actual VLE data on air.

 

[casualguitar]

On your plot, position 0 looks like nothing happened.

Yes.

Not exactly. There should be some similarities, although the 1 tank is inherently more disperse than the 3rd tank of the three tank.

I would start out by using Raoult's Law on all cases, unless you know more details about the non-iddalalities of the liquid or if you have actual VLE data on air.

Plot merging function is finished now (was a bit involved but it will be very useful), so I'll put together the plots mentioned in the posts above first thing tomorrow. All plots for both liquid/gas phase. These are the plots, unless I've forgotten any:

I'll actually do that last calculation in code rather than as a hand calculation first

 

[casualguitar]

Three tank plot for 0.05kg/s and single tank plot for 0.05kg/s:

Fluid temperature versus position for a range of times (0,200,400,600,800,1000):

Solid temperature versus position for a range of times (0,200,400,600,800,1000):

The above are for n=3. Might be more visually interesting to see more spatial plots, so the same plots for n=32:
Fluid:

Solid:

Note: The time for the n=32 simulations to run to completion is about 2 minutes

 

[casualguitar]

Hmm there is some code editing required to run a range of UA values. What I'll do for now is plot U = 1, 1.5 and 2 for the single tank:
View attachment 294079

Setting up the code for the bubble/dew point/mole fraction calculation now
Edit: These UA plots don't make much sense to me yet. I thought increasing U would decrease the time required to reach 300K. So I need to think about this a bit

 

[Chestermiller]

Three tank plot for 0.05kg/s and single tank plot for 0.05kg/s:
View attachment 294069
Fluid temperature versus position for a range of times (0,200,400,600,800,1000):
View attachment 294073
Solid temperature versus position for a range of times (0,200,400,600,800,1000):
View attachment 294074
The above are for n=3. Might be more visually interesting to see more spatial plots, so the same plots for n=32:
Fluid:
View attachment 294075
Solid:
View attachment 294076
Note: The time for the n=32 simulations to run to completion is about 2 minutes

I still have a problem with tank 1 in these which seems to show an instantaneous response of the temperature to the inlet temperature. It seems to me that this is just the inlet temperature to the bed. This part of the results just doesn't seem right to me. Maybe you can show the comparison for the liquid-only model for these same cases?

 

[casualguitar]

I still have a problem with tank 1 in these which seems to show an instantaneous response of the temperature to the inlet temperature

I agree here maybe by inlet conditions for the gas setup (300K) are incorrect

Maybe you can show the comparison for the liquid-only model for these same cases?

Yes sounds good I'll compare to the liquid plots

 

[casualguitar]

I still have a problem with tank 1 in these which seems to show an instantaneous response of the temperature to the inlet temperature. It seems to me that this is just the inlet temperature to the bed. This part of the results just doesn't seem right to me. Maybe you can show the comparison for the liquid-only model for these same cases?

Liquid plot for n=3. I thought the initial conditions might have been set up correctly but no, at position 1 (index 0) the initial temperature is 80 and it rapidly rises to 90K (similar behaviour in the gas plot)

The notation on the graph above is poor. 'Position 1' is at index 0. Position 1 is where the boundary conditions are applied. The way the code is set up, the bed/fluid temperatures for the entire bed are both set to 80K, and then the first position only is overwritten, and this is the boundary condition:

So using that approach it is expected that the index 0 position would go straight to the inlet fluid temperature. It does this in the model that matches the analytical solution also

Because it matches the analytical solution for the single tank, I think this seems ok

 

[casualguitar]

I'll actually do that last calculation in code rather than as a hand calculation first

Also this dew point/bubble point/composition calculation I will do this evening (have some report writing to finish out today)

 

[casualguitar]

Here is a calculation that I would like to see you do to get your feet wet on the mixture model: I have 1 mole of air (80% nitrogen, 20% oxygen) at 2 bars. What is the dew point and what is the bubble point? At a temperature half-way between the dew point and bubble point, what is the molar split of liquid and vapor, the mole fractions in the liquid phase, the mole fractions in the vapor phase, the density of the liquid phase, and the density of the vapor phase?

Here is a rough solution to the calculation you mentioned above. The code should be fairly readable. You can see the Mixture setup about half way down, where the components and the mole fractions are defined. From there once we have the mixture defined there are bubble and dew point methods, and a flash method to get the liquid/vapour fractions and compositions.

One note here is that this is an ideal flash calculation. I do have the functionality for non-ideal flash. Its in another Github repository so I need to do some file moving to get it. This also has density methods for the liquid and vapour phases (or anything else we would need). For now here are some preliminary results from he above simulation:

 

[Chestermiller]

Liquid plot for n=3. I thought the initial conditions might have been set up correctly but no, at position 1 (index 0) the initial temperature is 80 and it rapidly rises to 90K (similar behaviour in the gas plot)
View attachment 294099
The notation on the graph above is poor. 'Position 1' is at index 0. Position 1 is where the boundary conditions are applied. The way the code is set up, the bed/fluid temperatures for the entire bed are both set to 80K, and then the first position only is overwritten, and this is the boundary condition:
View attachment 294101
So using that approach it is expected that the index 0 position would go straight to the inlet fluid temperature. It does this in the model that matches the analytical solution also

Because it matches the analytical solution for the single tank, I think this seems ok

I'm totally confused about this. Are you saying that the overwrite the results for the first tank with the inlet condition, so it is not really the temperature vs time of the material in the first tank? If so, we never get to see the temperature vs time for the first tank? What about the subsequent tanks, are they somehow overwritten too? Why do you overwrite the first tank? That is of much more interest than the inlet condition. If this is what you do, please undo it.

 

[Chestermiller]

View attachment 294124
Here is a rough solution to the calculation you mentioned above. The code should be fairly readable. You can see the Mixture setup about half way down, where the components and the mole fractions are defined. From there once we have the mixture defined there are bubble and dew point methods, and a flash method to get the liquid/vapour fractions and compositions.

One note here is that this is an ideal flash calculation. I do have the functionality for non-ideal flash. Its in another Github repository so I need to do some file moving to get it. This also has density methods for the liquid and vapour phases (or anything else we would need). For now here are some preliminary results from he above simulation:
View attachment 294125
View attachment 294126

Does the flash calculation assume constant enthalpy? if so, the temperature changes, right? What is the final temperature. I was hoping for a hand calculation on your part. Do you know how to set up the equations for such a hand calculation, assuming ideal mixing and ideal gas behavior?

 

[casualguitar]

I'm totally confused about this. Are you saying that the overwrite the results for the first tank with the inlet condition, so it is not really the temperature vs time of the material in the first tank? If so, we never get to see the temperature vs time for the first tank? What about the subsequent tanks, are they somehow overwritten too? Why do you overwrite the first tank? That is of much more interest than the inlet condition. If this is what you do, please undo it.

There we go. You were right the results for tank 1 were being overwritten by the initial conditions. Obvious mistake looking back. The highest graph now shows the actual behaviour in tank 1, rather than the initial conditions.

Does the flash calculation assume constant enthalpy? if so, the temperature changes, right? What is the final temperature. I was hoping for a hand calculation on your part. Do you know how to set up the equations for such a hand calculation, assuming ideal mixing and ideal gas behavior?

I don't know actually. I took a look at the calculation (in code form) and its not clear to me. Yes I agree a hand calculation for the ideal flash would be useful to get a grip on how it works. I don't know how to set up the equations for this, I can figure this out though

After that hand calculation I'll redo the code with the real flash (and real mixture) methods. I guess this would be too involved to do by hand, however knowing how ideal flash works will give some insight into real flash calculations

I have a report due that I will finish this evening so I'll be back on the ideal flash hand calculation first thing tomorrow morning

 

[casualguitar]

Here is a calculation that I would like to see you do to get your feet wet on the mixture model: I have 1 mole of air (80% nitrogen, 20% oxygen) at 2 bars. What is the dew point and what is the bubble point? At a temperature half-way between the dew point and bubble point, what is the molar split of liquid and vapor, the mole fractions in the liquid phase, the mole fractions in the vapor phase, the density of the liquid phase, and the density of the vapor phase?

The bubble/dew point:

So we can use the Antoine equation and Raoult's law iteratively (bisection method) here to solve for both the de and bubble point temperatures

Algorithm for this could be as follows:
1) Set a temperature range that will capture the dew/bubble points (the boiling points of N2 and O2 at 2 bar for example)
2) Calculate the vapour pressure of N2 and O2 with the Antoine equation (using temperature)
3) Calculate the vapour composition of both components using Raoult's law
4) Sum both vapour compositions (O2+N2)

For the bubble point: Depending on if the sum of the vapour compositions is greater than or less than 1, update the temperature bounds accordingly and repeat

For the dew point: Depending on if the sum of the vapour compositions is greater than or less than 1, update the temperature bounds accordingly and repeat

Now that I follow the algorithm, this is the algorithms in code (I don't think there's value in writing this out):

Air dew and bubble point temperatures at 2 bar:

The density of the liquid and vapour phases should be fine. I could use two correlations (one for vapour and and one for liquid density) and get a weighted average based on the composition of each phase

Getting the phase molar split and phase compositions seems more involved. Reading up on this now

 

[Chestermiller]

The bubble/dew point:

So we can use the Antoine equation and Raoult's law iteratively (bisection method) here to solve for both the de and bubble point temperatures

Algorithm for this could be as follows:
1) Set a temperature range that will capture the dew/bubble points (the boiling points of N2 and O2 at 2 bar for example)
2) Calculate the vapour pressure of N2 and O2 with the Antoine equation (using temperature)
3) Calculate the vapour composition of both components using Raoult's law
4) Sum both vapour compositions (O2+N2)

For the bubble point: Depending on if the sum of the vapour compositions is greater than or less than 1, update the temperature bounds accordingly and repeat

For the dew point: Depending on if the sum of the vapour compositions is greater than or less than 1, update the temperature bounds accordingly and repeat

Now that I follow the algorithm, this is the algorithms in code (I don't think there's value in writing this out):
View attachment 294188
Air dew and bubble point temperatures at 2 bar:
View attachment 294189

The density of the liquid and vapour phases should be fine. I could use two correlations (one for vapour and and one for liquid density) and get a weighted average based on the composition of each phase

Getting the phase molar split and phase compositions seems more involved. Reading up on this now

Given that there is only a 3 degree range between the bubble point and dew point, it won't make a significant difference (considering all the other uncertainties in the model), if you just use the existing model with Tsat equal to 86.7 K, and the change in enthalpy of vaporization equal to the molar weighted average of the molar heats of vaporization. What do you think?

 

[casualguitar]

Given that there is only a 3 degree range between the bubble point and dew point, it won't make a significant difference (considering all the other uncertainties in the model), if you just use the existing model with Tsat equal to 86.7 K, and the change in enthalpy of vaporization equal to the molar weighted average of the molar heats of vaporization. What do you think?

This could be very useful yes I like this idea also. It would also simplify the next stage of modelling fairly significantly so let's do it. I will change my Tsat constant (in the existing model) to a dew/bubble temperature point calculation, as the pressure may vary from 2 bar. This fits in well. For the enthalpy of vaporisation I have temperature dependent correlations so these could fit well here instead of a weighted average also. This would mean it would no longer be constant but would update based on fluid temperature.

Because of the above change I want to mention another model goal that I had not mentioned previously now. My initial thought was that this model would come after the air liquefaction model, however given your message above to remove the dew/bubble point and use a single Tsat, I want to mention it now instead, in case it does have a short term affect on the direction this model takes.

So this cryogenic packed bed currently takes liquid air and vaporises it, leaving a cold packed bed. Or vice versa. The bed will at some point be required to take a CO2 enriched air stream (CO2/N2/O2), and separate the CO2 via CO2 solidification (or liquefaction). If we are taking the air dew/bubble points as a single temperature (say the average of the two), then this might be a good time to add CO2 to the model, because we could also assume this has its own liquefaction temperature i.e. we would have two condensation temperatures, one for CO2 and one for 'Air'. I wanted to mention this now in case it has any effect on the direction of our discussion

So the idea here is to take a CO2/N2/O2 stream as input, and to freeze (or potentially just liquefy) out the CO2 in the stream, to separate the CO2 from the air. It is essentially giving the packed bed both the capability to store energy and to capture CO2.

As for modelling capability, I have any functionality we would need to set up flue gas mixtures of varying compositions, methods to calculate densities of each phase, heat capacity correlations etc

So looking at the above, maybe first incorporating the enthalpy of vaporisation to the air only model might be a good next step. I can do either the weighted average approach like you mentioned, or I can use the temperature dependent correlations for latent heat.

What are your thoughts on that as a next step? Apologies I know dropping CO2 into the discussion is potentially a big deviation. If we can continue on developing the pure air model, forget about CO2 for now, and add CO2 once we've put the pure model through its paces that suits perfectly also. As I say I just wanted to mention it to put it in the background of the discussion

 

[casualguitar]

So some simplifications we've made:
1) The enthalpy/mass holdup relations. In previous posts you mentioned we would stick to correlations like these. If its an option then we can also use H(T,P) methods already made in other Python libraries i.e. just use a property calculation library to calculate enthalpy
2) The heat of vaporisation. Right now we're leaving this as a constant. I have methods again from thermo libraries that treat this as a function of temperatures so we have the option to use this
3) The saturation conditions. Right now we've got Tsat defined at the top as a constant. I have Tsat(P) functionality so I could use this no bother to get the saturation temperature
4) The quality. Right now we get around the quality with the saturation zone relations
5) The fluid. So assuming a pure fluid let's us have one saturation temperature. Air is a mixture so we'll have a bubble/dew point. This definitely sounds difficult to implement. That said I have a lot of useful functionality from the thermo property libraries like dew/bubble point temperature calculations for mixtures, Pressure-Enthalpy flash functionality, really anything we would need. So the important bit here is the actual algorithm, the functionality is available
6) We have constant heat capacity, etc. This can easily be switched out for temperature/pressure dependent values later on but yes its probably best to leave them as constant for now to avoid cluttering etc
7) The U value. I've put a value of 1 (out of thin air). We could improve on that

Going back to this list above, I wanted to comment on each item to get a sense for what the next best steps might be (also I'm mostly forgetting about the above CO2 comment for now unless you think we should do otherwise):

Commenting on each point in the list above:
1) Forget about this for now

2) Use a single value for both the heat of vaporisation of O2 and N2 i.e. do not have a temperature dependent heat of vaporisation, and use the weighted value as you mentioned i.e. $\Delta H_{AIR} = x\Delta H_{N2} + (1-x)\Delta H_{O2}$ where $x$ is the mole fraction of $N2$ in the feed

3) As mentioned by you we can simplify the dew/bubble point model. I will add a dew/bubble point calculation to the model dewT(P), bubbleT(P), then take the average of the two values and use this as Tsat. This would be nice as we would have a constant temperature zone that would allow us to see the change in vapour/liquid fraction across this zone.

4) Similar to above, given the new assumption of a single saturation temperature we should be able to model the quality in the constant temperature zone

5) As mentioned above, simplifying to one Tsat rather than dew/bubble points. Those flash libraries will still be useful here (possibly for the quality calculations)

6) I can sub in temperature/pressure dependent heat capacity, heat of vaporisation, density at any point with the thermo property libraries. Maybe this is best left until nearer the end of this model to avoid cluttering and to simplify any debugging

7) I could use some correlations here to get a better estimate however it seems like this value will be a tuning parameter anyway

So right now I could:
1) Add the weighted heat of vaporisation
2) Add the Tsat calculation (average of dew and bubble T at a given pressure)

This is straightforward so I'll do this right now

 

[Chestermiller]

Going back to this list above, I wanted to comment on each item to get a sense for what the next best steps might be (also I'm mostly forgetting about the above CO2 comment for now unless you think we should do otherwise):

Commenting on each point in the list above:
1) Forget about this for now

2) Use a single value for both the heat of vaporisation of O2 and N2 i.e. do not have a temperature dependent heat of vaporisation, and use the weighted value as you mentioned i.e. $\Delta H_{AIR} = x\Delta H_{N2} + (1-x)\Delta H_{O2}$ where $x$ is the mole fraction of $N2$ in the feed

For the model, you would, of course, dividing by the MW to get the heat of vaporization per unit mass.

3) As mentioned by you we can simplify the dew/bubble point model. I will add a dew/bubble point calculation to the model dewT(P), bubbleT(P), then take the average of the two values and use this as Tsat. This would be nice as we would have a constant temperature zone that would allow us to see the change in vapour/liquid fraction across this zone.

4) Similar to above, given the new assumption of a single saturation temperature we should be able to model the quality in the constant temperature zone

5) As mentioned above, simplifying to one Tsat rather than dew/bubble points. Those flash libraries will still be useful here (possibly for the quality calculations)

6) I can sub in temperature/pressure dependent heat capacity, heat of vaporisation, density at any point with the thermo property libraries. Maybe this is best left until nearer the end of this model to avoid cluttering and to simplify any debugging

7) I could use some correlations here to get a better estimate however it seems like this value will be a tuning parameter anyway

So right now I could:
1) Add the weighted heat of vaporisation
2) Add the Tsat calculation (average of dew and bubble T at a given pressure)

This is straightforward so I'll do this right now

If you wanted, you could also replace constant Tsat zone with a range between the dew point and bubble point. In this range, you would parameterize T and $\rho$ as functions of h. So the model would still have 3 zones, but T would not be constant in the middle zone. I still contend that this minor change would not be worth the effort.

I would still like to see a comparison between the position model and the single tank model using the full temperature range, and comparing the 32-tank results from the position model with for the first tank with those using the single tank model with 32x the mass flow rate.

Regarding systems involving air and CO2, you can start figuring out how to do stand-alone phase change calculations involving lowering the temperature at constant pressure until a liquid or solid starts forming. Check out the thermodynamic diagrams for CO2 to get an idea what happens first for anticipated the concentrations in your mixtures.

 

[casualguitar]

If you wanted, you could also replace constant Tsat zone with a range between the dew point and bubble point. In this range, you would parameterize T and ρ as functions of h. So the model would still have 3 zones, but T would not be constant in the middle zone. I still contend that this minor change would not be worth the effort.

Ah I follow, like in one of our previous models with T1>T<T2. Yes I will note this as a minor change that I can make at a later point but yes I agree the constant Tsat does sound like a nice simplification to work with for now

I would still like to see a comparison between the position model and the single tank model using the full temperature range, and comparing the 32-tank results from the position model with for the first tank with those using the single tank model with 32x the mass flow rate.

Got it, yes I haven't remade these after the initial condition fix. I'll regenerate these using the correct initial conditions now

Regarding systems involving air and CO2, you can start figuring out how to do stand-alone phase change calculations involving lowering the temperature at constant pressure until a liquid or solid starts forming. Check out the thermodynamic diagrams for CO2 to get an idea what happens first for anticipated the concentrations in your mixtures.

Will do. It sounds like you're saying leave this until the air model is more developed, so yes I'll just take a look at the thermodynamic diagrams for CO2 and a few papers that model this kind of thing for now

Anyway that Tsat calculation is finished (yet to do the Hvap one but this is even simpler). I'll get the single tank/position model plots done now

 

[casualguitar]

I would still like to see a comparison between the position model and the single tank model using the full temperature range, and comparing the 32-tank results from the position model with for the first tank with those using the single tank model with 32x the mass flow rate.

Position versus temperature for a range of times.

Time versus temperature for the n=32 model, for a number of tanks in the spatial variation model, and the single tank model. The first tank matches up with the single tank model almost exactly:

The single tank is hard to see but its there behind position 1. One thing that stands out to me here is the shape of the temperature distribution for each point. Its close to right angled. So the temperature at each point stays at about 80K, until the 'wave' reaches the tank, then the temperature jumps quickly, and after the wave passes the temperature just gradually increases again. This would suggest very low axial dispersion I think? How does the above look to you?

Will add those small Tsat/Hvap changes this evening

 

[Chestermiller]

Position versus temperature for a range of times.
View attachment 294269
Time versus temperature for the n=32 model, for a number of tanks in the spatial variation model, and the single tank model. The first tank matches up with the single tank model almost exactly:
View attachment 294276
The single tank is hard to see but its there behind position 1. One thing that stands out to me here is the shape of the temperature distribution for each point. Its close to right angled. So the temperature at each point stays at about 80K, until the 'wave' reaches the tank, then the temperature jumps quickly, and after the wave passes the temperature just gradually increases again. This would suggest very low axial dispersion I think? How does the above look to you?

I'm wondering why there is no temporary leveling off at Tsat, and why there are sharp changes in slope with time at various positions. Also, what are the solid bed temperature variations with time at the various positions? (I did expect the wave to travel through the bed gradually at the lowest temperature as the graph shows)

 

[casualguitar]

I'm wondering why there is no temporary leveling off at Tsat, and why there are sharp changes in slope with time at various positions.

The sharp change I think is because that's the 3 tank model. This is the same plot for n=32, so these changes are smoothed out:

The model takes 2-3 minutes to run now each time so I'll look into speeding this up

The temporary levelling off at Tsat does happen, however since the system is increasing in temperature quite quickly this zone is short:

Zooming in on it (a lot) you can see at 105K (Tsat) there is a levelling off of the curve to a constant saturation temperature of 105K. Its just a few seconds of saturation and then single phase again. Note this graph is n=32 (incorrect title however I left it as it takes a while to regenerate)

Also, what are the solid bed temperature variations with time at the various positions?

Solid bed temperature variations at various positions here:

If the above looks ok to you I'll add the new Tsat (based on pressure) and the new heat of vaporisation weighted average (and check that the output remains about the same)

 

[casualguitar]

If the above looks ok to you I'll add the new Tsat (based on pressure) and the new heat of vaporisation weighted average (and check that the output remains about the same)

Added both the Tsat and Hvap changes. So now Tsat is the average of the dew and bubble points at a given pressure, and the heat of vaporisation is the weighted average of the O2 and N2 heats of vaporisation, and I've taken them both as constants for now. I can add temperature dependent Hvap later easily:

Position versus temperature for a range of times:

Time versus temperature for a range of positions, for n=32:

Zooming into show the levelling off at Tsat, 87K or so as expected

Going back to these points:

Commenting on each point in the list above:
1) Forget about this for now

2) Use a single value for both the heat of vaporisation of O2 and N2 i.e. do not have a temperature dependent heat of vaporisation, and use the weighted value as you mentioned i.e. ΔHAIR=xΔHN2+(1−x)ΔHO2 where x is the mole fraction of N2 in the feed

3) As mentioned by you we can simplify the dew/bubble point model. I will add a dew/bubble point calculation to the model dewT(P), bubbleT(P), then take the average of the two values and use this as Tsat. This would be nice as we would have a constant temperature zone that would allow us to see the change in vapour/liquid fraction across this zone.

4) Similar to above, given the new assumption of a single saturation temperature we should be able to model the quality in the constant temperature zone

5) As mentioned above, simplifying to one Tsat rather than dew/bubble points. Those flash libraries will still be useful here (possibly for the quality calculations)

6) I can sub in temperature/pressure dependent heat capacity, heat of vaporisation, density at any point with the thermo property libraries. Maybe this is best left until nearer the end of this model to avoid cluttering and to simplify any debugging

7) I could use some correlations here to get a better estimate however it seems like this value will be a tuning parameter anyway

So right now I could:
1) Add the weighted heat of vaporisation
2) Add the Tsat calculation (average of dew and bubble T at a given pressure)

Points 2,3,5 are now implemented. Forgetting point 1 for now, we've got 4,6,7 left

4) Modelling the quality.
6) Using temperature dependent heat capacity, heat of vaporisation, density
7) A correlation for U

Leaving 6 until the end, and leaving 7 until the end also (or leaving it out completely) seems like a good approach. So this leaves us with 4) modelling the quality

My thoughts here are that we've got Tsat, and we know fluid enthalpy across the whole phase change zone. So if we got a correlation for saturated liquid enthalpy, and saturated vapour enthalpy, we can get X (liquid quality) at any time with: $$h = xh_{SAT,L} + (1-x)h_{SAT,V}$$ So we really don't need any core model adjustments to model quality, we just need those two saturated liquid/vapour enthalpy correlations. We can get a correlation or I can use the thermo property library to get these. Does this sound correct, and does this sound like a reasonable next step to you?

 

[Chestermiller]

Added both the Tsat and Hvap changes. So now Tsat is the average of the dew and bubble points at a given pressure, and the heat of vaporisation is the weighted average of the O2 and N2 heats of vaporisation, and I've taken them both as constants for now. I can add temperature dependent Hvap later easily:

Position versus temperature for a range of times:
View attachment 294403

Time versus temperature for a range of positions, for n=32:
View attachment 294405Zooming into show the levelling off at Tsat, 87K or so as expected
View attachment 294404

Going back to these points:

Points 2,3,5 are now implemented. Forgetting point 1 for now, we've got 4,6,7 left

4) Modelling the quality.
6) Using temperature dependent heat capacity, heat of vaporisation, density
7) A correlation for U

Leaving 6 until the end, and leaving 7 until the end also (or leaving it out completely) seems like a good approach. So this leaves us with 4) modelling the quality

My thoughts here are that we've got Tsat, and we know fluid enthalpy across the whole phase change zone. So if we got a correlation for saturated liquid enthalpy, and saturated vapour enthalpy, we can get X (liquid quality) at any time with: $$h = xh_{SAT,L} + (1-x)h_{SAT,V}$$ So we really don't need any core model adjustments to model quality, we just need those two saturated liquid/vapour enthalpy correlations. We can get a correlation or I can use the thermo property library to get these. Does this sound correct, and does this sound like a reasonable next step to you?

This can be evaluated directly from the model, since we are taking the enthalpy of the saturated liquid as zero in the model and the enthalpy of the saturated vapor as the heat of vaporization (our region 2 equations). So the vapor quality from the model is just $h/\Delta h_{vap}$.

 

[casualguitar]

since we are taking the enthalpy of the saturated liquid as zero in the model and the enthalpy of the saturated vapor as the heat of vaporization (our region 2 equations)

Ah understood. I guess we'll also need a conditional statement for the quality array saying 'if H/Hvap is < 0 then vapour quality is 0, or H>hVap then quality is 1'.

Here is time vs H and T for tank 1 (3 tank system). The liquid and saturated zones do seem to allow for rapid increase of temperature (the graph is almost vertical in these zones). This slows down a lot for the vapour phase, and it looks a lot more natural.

Zooming in on the saturation zone:

So as enthalpy reaches 200, the temperature also does the large jump from approx 87K to 210K, meaning that our saturation zone correlations predict very fast increases in temperature across this zone

Anyway here's the plot of vapour quality for each of the 3 tanks:

So the saturation zone lasts 1 to 2 seconds approx in each tank

And here's the code I used to create the vapour quality array:

Hmm, looking at the remaining model additions, maybe the temperature/pressure dependent properties (heat capacity, density, heat of vaporisation, dp/dh) might be a reasonable next addition? This would mean replacing the temperature/mass functions also, with correlations available in the thermo library

Do you think there are other changes we should make before this?

Also just as a note - here's the bubble/dew point code just to get an idea of how it looks to set up:

 

[casualguitar]

Hmm, looking at the remaining model additions, maybe the temperature/pressure dependent properties (heat capacity, density, heat of vaporisation, dp/dh) might be a reasonable next addition? This would mean replacing the temperature/mass functions also, with correlations available in the thermo library

Just as a side note I've imported the required functionality for temperature dependent heat capacity and heat of vaporisation (not implemented, just ready to implement).

The heat capacity functionality calculates $Cp$ at a given temperature. The heat of vaporisation function returns the heat of vaporisation for 'all elements in the mixture at the given temperature' so in our case it will return $hVap$ for $N2$ and $O2$. I can use the weighted average function I made so we can use one value for $hVap$ of air

For density, I've also got this functionality imported and ready to implement if you think this is a good next step. It calculates the weighted average of the O2 and N2 densities at a given temperature. If the phases of $O2$ and $N2$ are different (this occurs for a very small range of temperatures) then it will still calculate a weighted average.

However, actually looking at the code for the model, $Cp$ and $hVap$ only occur in the temperature(H), mass(H) and dpdh(H) functions, so it may be better so sub out these entire functions rather than just $Cp$ and $hVap$

A possible substitution for the first two functions might be:
temperature(H):
- We know $H$ and $P$, and also the composition of the mixture, so I think we have enough to do a flash calculation here to calculate the temperature at these conditions of $H$ and $P$

mass(H):
- if vessel volume is known (it is), then we could use the density functionality from the thermo library to get density, which allows us to calculate mass holdup

dpdh(H):
- haven't found a solution to this one yet. I guess this functionality exists also. but I haven't found it yet