Modelling of two phase flow in packed bed using conservation equations

[Chestermiller]

The idea is clear, and I can validation check by using the values of the other parameters (h,m,Ts).

One question - I'm not clear on the correct value of n. My thoughts are that temperature is evaluated at the 'centre' of a tank. n=0 is at the inlet and n=n is at the outlet. So to have 3 'tanks', we would need n=5?

Either way none of n=3,4,5 give exactly the same solution as the analytical one, so I need to do some testing regardless

The parameter n is the total number of tanks being considered, not a counting index. The counting index “i” is 1 for the first tank, 2 for the second tank, and 3 for the 3rd tank. The location of the center of the i’th tank is at $x=\frac{(2i-1)}{2}\Delta x$. The left side of the i'th tank is at $x=(i-1)\Delta x$. The right side of the i'th tank is at $x=i \Delta x$.

 

[casualguitar]

The parameter n is the total number of tanks being considered, not a counting index. The counting index “i” is 1 for the first tank, 2 for the second tank, and 3 for the 3rd tank. The location of the center of the i’th tank is at $x=\frac{(2i-1)}{2}\Delta x$. The left side of the i'th tank is at $x=(i-1)\Delta x$. The right side of the i'th tank is at $x=i \Delta x$.

Understood now. I read the original comments you had deriving the spatial variation model and I follow now.

Ok then sounds like n=3 and I should run some diagnostics

 

[casualguitar]

Understood now. I read the original comments you had deriving the spatial variation model and I follow now.

Ok then sounds like n=3 and I should run some diagnostics

So assuming the single tank model is correct (for now), given it matches the Mathematica model, I'm working with on getting the spatial variation model to match the single tank model. Specifically, getting tank index 1 with min = 0.05kg/s to match the single tank model, for min = 0.15kg/s
Some data for n=3, tank 1:

The mass holdup versus time:

This checks out because its just 8kg/3

dm/dt:

Zero as expected

dp/dh is also zero for all times

The differential equations all look ok to me, so I reckon the error is in how I have set up the initial and boundary conditions. Checking this now

 

[Chestermiller]

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

Are you still using an automatic integrator, and calculating the time derivative of the entire solution vector all at once for the integrator to use? If so, do you have the solution vector as the column vector ordered as $h_1,T_{S,1}, h_2, T_{S,2},h_3, T_{S,3},...$?

If you are using this scheme, you calculate the time derivative of $h_1$ for the first tank using $\dot{m}_0$ which does not change with time. Then you calculate $\dot{m}_1$ exiting tank 1 and entering tank 2 using the mass balance equation with the time derivative of $h_1$ that you just calculated for the first tank. Then you use this value of $\dot{m}_1$ entering tank 2 to calculate the time derivative of $h_2$ for the second tank. Then you calculate $\dot{m}_2$ exiting tank 2 and entering tank 3 using the mass balance equation with the time derivative of $h_2$ that you just calculated for the second tank. Then you use this value of $\dot{m}_2$ entering tank 3 to calculate the time derivative of $h_3$ for the third tank. etc. In the end, you will have all the time derivatives of the h's in all the tanks, which gets loaded into the time derivative vector. And, of course, in each of these steps you are also calculating the time derivative of the bed temperature for the specific tank, and also loading this into the derivative vector.

Is this what you are doing? If not, what is the scheme that you are applying?

 

[casualguitar]

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

Are you still using an automatic integrator, and calculating the time derivative of the entire solution vector all at once for the integrator to use? If so, do you have the solution vector as the column vector ordered as h1,TS,1,h2,TS,2,h3,TS,3,...?

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

If you are using this scheme, you calculate the time derivative of h1 for the first tank using m˙0 which does not change with time. Then you calculate m˙1 exiting tank 1 and entering tank 2 using the mass balance equation with the time derivative of h1 that you just calculated for the first tank. Then you use this value of m˙1 entering tank 2 to calculate the time derivative of h2 for the second tank. Then you calculate m˙2 exiting tank 2 and entering tank 3 using the mass balance equation with the time derivative of h2 that you just calculated for the second tank. Then you use this value of m˙2 entering tank 3 to calculate the time derivative of h3 for the third tank. etc. In the end, you will have all the time derivatives of the h's in all the tanks, which gets loaded into the time derivative vector.

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

And, of course, in each of these steps you are also calculating the time derivative of the bed temperature for the specific tank, and also loading this into the derivative vector.

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:

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

 

[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