Modelling of two phase flow in packed bed using conservation equations

[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

 

[Chestermiller]

You seem intent on doing this. For an ideal liquid mixture and vapor mixture, the VLE equations are going to be $$\frac{xP_{N2}^*(T)}{P}+\frac{(1-x)P_{O2}^*(T)}{P}=1$$$$y=\frac{xP_{N2}^*(T)}{P}$$$$L+V=1$$and $$Lx+Vy=z$$or$$L(z-x)=V(y-z)$$where P is the total pressure, x is the mole fraction N2 in the liquid, y is the mole fraction N2 in the vapor, z is the mole fraction N2 in the combination of vapor and liquid (i.e., 0.8), L is the molar fraction liquid, V is the molar fraction vapor, and starred quantities are the equilibrium vapor pressures of the pure components. OK so far?

 

[casualguitar]

You seem intent on doing this. For an ideal liquid mixture and vapor mixture, the VLE equations are going to be $$\frac{xP_{N2}^*(T)}{P}+\frac{(1-x)P_{O2}^*(T)}{P}=1$$$$y=\frac{xP_{N2}^*(T)}{P}$$$$L+V=1$$and $$Lx+Vy=z$$or$$L(z-x)=V(y-z)$$where P is the total pressure, x is the mole fraction N2 in the liquid, y is the mole fraction N2 in the vapor, z is the mole fraction N2 in the combination of vapor and liquid (i.e., 0.8), L is the molar fraction liquid, V is the molar fraction vapor, and starred quantities are the equilibrium vapor pressures of the pure components. OK so far?

Good so far!

 

[Chestermiller]

Good so far!

Once the overall mole fraction of N2 is specified (0.8), and the temperature and pressure are known, the mole fractions of N2 and O2 in the liquid and vapor are known, and so is the liquid vapor molar split L and V. This would then allow us to determine the molar enthalpy and density of the liquid vapor mixture. So for our mixture of z = 0.8, we can determine molar H(T,P) and molar $\rho(T,P)$. This would also enable us to determine T(H,P) and $\rho(H,P)$ in the 2 phase region. So, for the 2 phase region, we could parameterize temperature and density as functions of molar enthalpy and pressure, say, using a quadratic fit. This would give us all we need for the 2 phase region.

 

[casualguitar]

Once the overall mole fraction of N2 is specified (0.8), and the temperature and pressure are known, the mole fractions of N2 and O2 in the liquid and vapor are known, and so is the liquid vapor molar split L and V.

I fully follow up to this point

This would then allow us to determine the molar enthalpy and density of the liquid vapor mixture. So for our mixture of z = 0.8, we can determine molar H(T,P) and molar ρ(T,P).

I agree that knowing the temperature, pressure, the composition of the liquid/vapour phases, and the liquid vapour molar split is enough information to find the molar enthalpy/density of the mixture. So yes I also agree that we can get $H(T,P)$ and $\rho(T,P)$.

I have a point of confusion in that its not clear to me what use $H(T,P)$ is, if we already have the fluid energy/mass balances though. I think if we continue on with the algorithm discussion this will become clear though

So, for the 2 phase region, we could parameterize temperature and density as functions of molar enthalpy and pressure, say, using a quadratic fit. This would give us all we need for the 2 phase region.

I follow this if you are saying we can get existing $T(H,P)$ and $\rho(T,P)$ fits from literature data?

Also a note: The above is for an ideal liquid/vapour mixture. I think it will be useful to code this manually once I understand it, and then later use existing functionality to implement this for real mixtures. What do you think?

 

[Chestermiller]

I fully follow up to this point

I agree that knowing the temperature, pressure, the composition of the liquid/vapour phases, and the liquid vapour molar split is enough information to find the molar enthalpy/density of the mixture. So yes I also agree that we can get $H(T,P)$ and $\rho(T,P)$.

I have a point of confusion in that its not clear to me what use $H(T,P)$ is, if we already have the fluid energy/mass balances though. I think if we continue on with the algorithm discussion this will become clear though

These functionalities are independent of the energy and mass balances. They are the physical properties of N2-O2 mixtures from the thermodynamics of ideal solutions.

I follow this if you are saying we can get existing $T(H,P)$ and $\rho(T,P)$ fits from literature data?

No. We are assuming we are dealing with ideal solutions. Otherwise, we could not apply Raoult's law. The only literature data we need to do this are the physical properties of pure N2 and O2 as functions of T and P. Are familiar with the thermodynamics of mixtures/solutions, the concept of partial molar properties, and the behavior of these properties for ideal solutions? If not, see Chapter 10 of Smith and Van Ness, Introduction to Chemical Engineering Thermodynamics.

Also a note: The above is for an ideal liquid/vapour mixture. I think it will be useful to code this manually once I understand it, and then later use existing functionality to implement this for real mixtures. What do you think?

To do this for real mixtures, you need to know the activity coefficients (and their variation with concentration) of N2 and O2 in the liquid. Do you have data on this?

I still think you are focusing too much on a portion of the model with much less uncertainty (inaccuracy) than other parts of the model. I don't think it's work your time.

 

[casualguitar]

The only literature data we need to do this are the physical properties of pure N2 and O2 as functions of T and P.

This answers the confusion. Its clear to me now. I would like to avoid coding this data myself if possible, as its already been coded up by others

Are familiar with the thermodynamics of mixtures/solutions, the concept of partial molar properties, and the behavior of these properties for ideal solutions? If not, see Chapter 10 of Smith and Van Ness, Introduction to Chemical Engineering Thermodynamics.

I am familiar with these concepts yes. I took at look at chapter 10 to confirm this

To do this for real mixtures, you need to know the activity coefficients (and their variation with concentration) of N2 and O2 in the liquid. Do you have data on this?

I do have data on this. This link is to the table of contents of the 'Thermo' library I am using. The table of contents gives a good idea of the functionality available: https://thermo.readthedocs.io

It has basically anything we would need, even including things like heat of fusion correlations for CO2 etc

I still think you are focusing too much on a portion of the model with much less uncertainty (inaccuracy) than other parts of the model. I don't think it's work your time.

Which part(s) of the model have more variation than this one? By this one, I mean using constant heat capacity/hVap, and using out temperature/density correlations, rather than using real data

I am happy to work on any part of the model. I don't want to direct the modelling that much at all for now because I don't want to take away from your intuition. I would like to take any guidance you can give (I'll challenge it where I can but really your intuition for model progression is what should drive the changes I think).

But yes I am 100% open to suggestions on how best to approach improving the model from here, so if leaving these correlations/constants as they are is fine then I'm happy to do that. At some point I'll make those changes so I can sleep easy!

 

[Chestermiller]

This answers the confusion. Its clear to me now. I would like to avoid coding this data myself if possible, as its already been coded up by others

I am familiar with these concepts yes. I took at look at chapter 10 to confirm this

I do have data on this. This link is to the table of contents of the 'Thermo' library I am using. The table of contents gives a good idea of the functionality available: https://thermo.readthedocs.io

It has basically anything we would need, even including things like heat of fusion correlations for CO2 etc

Which part(s) of the model have more variation than this one? By this one, I mean using constant heat capacity/hVap, and using out temperature/density correlations, rather than using real data

I am happy to work on any part of the model. I don't want to direct the modelling that much at all for now because I don't want to take away from your intuition. I would like to take any guidance you can give (I'll challenge it where I can but really your intuition for model progression is what should drive the changes I think).

But yes I am 100% open to suggestions on how best to approach improving the model from here, so if leaving these correlations/constants as they are is fine then I'm happy to do that. At some point I'll make those changes so I can sleep easy!

If it were me, I would be looking at how to quantify the heat transfer coefficient for the 3 regions and I would be looking at the pressure variation in the bed. For the latter, the first thing I would consider is how the mass flow rate of liquid would look if the bed were full of liquid at the initial condition and there was atmospheric pressure (or other constant pressure) at the inlet and outlet of the bed. In other words, free flow of liquid through the bed under gravity. How would this compare with the imposed mass flow rate (0.05 kg/sec). This would provide an idea, at least initially, on whether the pressure varies significantly in the bed. I would also look at the number of tanks in the bed, according to present calculations, in which there are two phases present in the tanks (to determine whether this middle region needs to be seriously addressed in terms of fluid flow).

 

[casualguitar]

Ok I'm going to head out for a walk and think about this for a bit

If it were me, I would be looking at how to quantify the heat transfer coefficient for the 3 regions and I would be looking at the pressure variation in the bed.

First two things that come to mind here are the Gnielinski and Ergun correlations

For the latter, the first thing I would consider is how the mass flow rate of liquid would look if the bed were full of liquid at the initial condition and there was atmospheric pressure (or other constant pressure) at the inlet and outlet of the bed. In other words, free flow of liquid through the bed under gravity. How would this compare with the imposed mass flow rate (0.05 kg/sec). This would provide an idea, at least initially, on whether the pressure varies significantly in the bed.

I can do this. Although, would it not be as easy to implement the Ergun equation directly? Actually yes I'll do both. Similar to previous models we can use one to validate the other

I would also look at the number of tanks in the bed, according to present calculations, in which there are two phases present in the tanks (to determine whether this middle region needs to be seriously addressed in terms of fluid flow).

It isn't clear to me if our 'tank' length equates to a real distance in metres. I guess the answer to this will be in your original spatial variation model derivation, so I'll take a look

I would also look at the number of tanks in the bed, according to present calculations, in which there are two phases present in the tanks (to determine whether this middle region needs to be seriously addressed in terms of fluid flow).

Got it

And yes since the inaccuracies caused by using constant Cp/Hvap etc are small, I won't put time into modelling this manually. I'll just make the change using the thermo library. The change is as tiny as changing the code from cp = 1, to cp = air.cp(T), for example. Similar for hVap. I still maintain though just for my own model clarity i.e. following the model in my own head, I'd like to replace the temperature(H) and mass(H) functions we have currently with the library correlations at some point, to get it out of my head and into the model, even if it isn't a source of large inaccuracy.

So yes I'll think about this and propose a flow to the above tasks this evening, if this is suitable for you

 

[casualguitar]

How does this flow sound to you?

Pressure Gradient:
1) Independently implement the Ergun equation to get a sense of the expected pressure drop for liquid only and gas only flow across the packed bed i.e. don't include it in the spatial variation model initially.

2) Make a function that allows calculation of the pressure at any point along the bed. I guess this will just be the Ergun equation again except instead of taking the total bed length as input, this function would take a distance from the inlet as input. This would give the $\Delta P$ between the inlet and that point. Then we just take that $\Delta P$ from the inlet pressure to get the pressure at that point.

3) This function could then go into the spatial variation model. It just affects the density and $\frac{d\rho}{dh}$ functions currently so its easy to slot in

4) Later, this function can be used to get the input pressure for real mixture heat capacity, density, etc. Not priority as you say because it won't affect the accuracy that much, but just to clear my heat and for model completion, I'll add it in

The author of the thermo library I mentioned also has some other nice libraries. The 'fluids' library has packed bed pressure drop functionality, here: https://fluids.readthedocs.io/fluids.packed_bed.html?highlight=ergun
This will give an immediate sense of the pressure drop

The real system:
I've chosen non-exact volume, mass, area, solid property values, etc. I know the dimensions of the actual tank being used and any information about the material (diameter, $C_p$, etc). I'll add in these so that model output resembles the real system:
- tank volume
- solid area
- inlet fluid mass flow
- solid particle diameter
- solid density
- solid $C_p$
- voidage
- tank length
- tank diameter

Heat transfer coefficient:
I don't know much about this. I know there are a few correlations used to get the Nusselt number, like Gnielinski and Achenbach, and we can use this with the Reynolds number I think to get the heat transfer coefficient.

This library (another from the same author) called 'heat transfer' calculates the Nusselt number for a packed bed given the particle diameter, voidage, superficial velocity, fluid density, fluid viscosity and the Prandtl number: https://ht.readthedocs.io/en/release/ht.conv_packed_bed.html?highlight=packed bed
So again we can get an immediate sense of the heat transfer coefficient using this (I am keen to use these libraries where possible, mostly out of interest in the functionality available)

I think the heat transfer coefficient modelling sounds like something that would slot into the model best after the above two sections are done (they shouldn't take long given the available functionality).

Question: Does '$n$', the number of tanks, correspond to a real length?

So yes in condensed form:
- setup of Ergun
- implementation of Ergun
- addition of $T$ and $P$ dependent parameters where possible
- editing tank/solid values to match exact real system values
- heat transfer coefficient modelling (I guess this will not be short so I won't aim to plan this yet)

Does this above flow look reasonable to you?

 

[Chestermiller]

How does this flow sound to you?

Pressure Gradient:
1) Independently implement the Ergun equation to get a sense of the expected pressure drop for liquid only and gas only flow across the packed bed i.e. don't include it in the spatial variation model initially.

2) Make a function that allows calculation of the pressure at any point along the bed. I guess this will just be the Ergun equation again except instead of taking the total bed length as input, this function would take a distance from the inlet as input. This would give the $\Delta P$ between the inlet and that point. Then we just take that $\Delta P$ from the inlet pressure to get the pressure at that point.

3) This function could then go into the spatial variation model. It just affects the density and $\frac{d\rho}{dh}$ functions currently so its easy to slot in

4) Later, this function can be used to get the input pressure for real mixture heat capacity, density, etc. Not priority as you say because it won't affect the accuracy that much, but just to clear my heat and for model completion, I'll add it in

The author of the thermo library I mentioned also has some other nice libraries. The 'fluids' library has packed bed pressure drop functionality, here: https://fluids.readthedocs.io/fluids.packed_bed.html?highlight=ergun
This will give an immediate sense of the pressure drop

The real system:
I've chosen non-exact volume, mass, area, solid property values, etc. I know the dimensions of the actual tank being used and any information about the material (diameter, $C_p$, etc). I'll add in these so that model output resembles the real system:
- tank volume
- solid area
- inlet fluid mass flow
- solid particle diameter
- solid density
- solid $C_p$
- voidage
- tank length
- tank diameter

Heat transfer coefficient:
I don't know much about this. I know there are a few correlations used to get the Nusselt number, like Gnielinski and Achenbach, and we can use this with the Reynolds number I think to get the heat transfer coefficient.

This library (another from the same author) called 'heat transfer' calculates the Nusselt number for a packed bed given the particle diameter, voidage, superficial velocity, fluid density, fluid viscosity and the Prandtl number: https://ht.readthedocs.io/en/release/ht.conv_packed_bed.html?highlight=packed bed
So again we can get an immediate sense of the heat transfer coefficient using this (I am keen to use these libraries where possible, mostly out of interest in the functionality available)

I think the heat transfer coefficient modelling sounds like something that would slot into the model best after the above two sections are done (they shouldn't take long given the available functionality).

Question: Does '$n$', the number of tanks, correspond to a real length?

The number of tanks is related to the total bed length and the grid spacing by $$L=n\Delta x$$

So yes in condensed form:
- setup of Ergun
- implementation of Ergun
- addition of $T$ and $P$ dependent parameters where possible
- editing tank/solid values to match exact real system values
- heat transfer coefficient modelling (I guess this will not be short so I won't aim to plan this yet)

Does this above flow look reasonable to you?

Yes

 

[casualguitar]

The number of tanks is related to the total bed length and the grid spacing by $$L=n\Delta x$$

Yes

Hi Chet, just letting you know I have other work to do today. I'll aim to start the above changes later this evening and update tomorrow. Thanks

 

[casualguitar]

The number of tanks is related to the total bed length and the grid spacing by $$L=n\Delta x$$

Yes

Hi Chet, back on this (took some time away for Christmas). Hope your enjoying the holidays

I'll comment the above stages I mentioned as I do them.

For the expected pressure drop across the bed, I first collected the actual values of flow/pressure/bed dimensions that will be used in the lab. I've modeled the two extreme cases, all vapour and all liquid. The actual pressure drop will be somewhere in the middle.

For reference:
Inlet pressure = 30bar
Mass flow = 0.05 kg/s
Bed length = 1.365 m
Bed diameter = 0.08 m
Voidage = 0.28
Particle diameter = 5mm

Heres the code for the expected pressure drop across the bed:

Output:

So for an inlet pressure of 30 bar, the expected pressure drop is at most 0.5bar (or in this ballpark anyway). This small drop indicates that the pressure drop won't have a large effect on the system.

Going off the list above:

So yes in condensed form:
- setup of Ergun
- implementation of Ergun
- addition of T and P dependent parameters where possible
- editing tank/solid values to match exact real system values
- heat transfer coefficient modelling (I guess this will not be short so I won't aim to plan this yet)

Ergun has been set up so that's task 1 done. I haven't implemented it in the spatial variation model yet so I'll do this now, even though it should have minimal impact. As part of this I'll need to switch the thermo parameters to T/P dependent. All of this will have a small effect on the model output.

This will just leave the heat transfer coefficient. Will aim to be ready to talk about the heat transfer coefficient model in a day or two

 

[Chestermiller]

Hi Chet, back on this (took some time away for Christmas). Hope your enjoying the holidays

I'll comment the above stages I mentioned as I do them.

For the expected pressure drop across the bed, I first collected the actual values of flow/pressure/bed dimensions that will be used in the lab. I've modeled the two extreme cases, all vapour and all liquid. The actual pressure drop will be somewhere in the middle.

For reference:
Inlet pressure = 30bar
Mass flow = 0.05 kg/s
Bed length = 1.365 m
Bed diameter = 0.08 m
Voidage = 0.28
Particle diameter = 5mm

Heres the code for the expected pressure drop across the bed:
View attachment 294774

Output:
View attachment 294775

So for an inlet pressure of 30 bar, the expected pressure drop is at most 0.5bar (or in this ballpark anyway). This small drop indicates that the pressure drop won't have a large effect on the system.

Going off the list above:

Ergun has been set up so that's task 1 done. I haven't implemented it in the spatial variation model yet so I'll do this now, even though it should have minimal impact. As part of this I'll need to switch the thermo parameters to T/P dependent. All of this will have a small effect on the model output.

This will just leave the heat transfer coefficient. Will aim to be ready to talk about the heat transfer coefficient model in a day or two

In previous calculations, you showed that, in tanks where the phase change is taking place, $\dot{m}$ exiting the tank is much larger than $\dot{m}_{in}=0.05\ kg/sec$. How do you think that affects the mass flow rates in and out of the tanks downstream of this? Doesn't it also cause them to have that same high value?

At 30 bars, what is the dew point and bubble point of the air?

 

[casualguitar]

At 30 bars, what is the dew point and bubble point of the air?

Just adding the bubble and dew points at 30bar:

In previous calculations, you showed that, in tanks where the phase change is taking place, m˙ exiting the tank is much larger than m˙in=0.05 kg/sec.

Yes I did show this I have the plots of the behaviour. I don't understand yet why this happens (or have the answer to the subsequent questions). Looking into this now

 

[Chestermiller]

Just adding the bubble and dew points at 30bar:
View attachment 294875Yes I did show this I have the plots of the behaviour. I don't understand yet why this happens (or have the answer to the subsequent questions). Looking into this now

It happens because the fluid is expanding from a liquid to a vapor. This causes all the liquid downstream to move with the same increased vapor velocity and mass flow rate. So the downstream liquid has a much greater mass flow rate than 0.05 kg/sec.

 

[casualguitar]

It happens because the fluid is expanding from a liquid to a vapor. This causes all the liquid downstream to move with the same increased vapor velocity and mass flow rate. So the downstream liquid has a much greater mass flow rate than 0.05 kg/sec.

Ok so in tanks where there is evaporation taking place, we will have a mixture of liquid and vapour. The liquid will therefore travel at the faster velocity (the vapour velocity). As the liquid density is greater than the vapour density, this will result in an increase in the mass flow rate.

Question on that - why would the liquid take the same value as the velocity of the gas in the saturated zone, and not at an 'averaged' speed i.e. $v_{mix} = (1-X)v_l + Xv_g$

How do you think that affects the mass flow rates in and out of the tanks downstream of this? Doesn't it also cause them to have that same high value?

Hmm well yes if the mass flow out of the saturated zone is increased then you would expect this value to propagate downstream.

However I don't see why it would be increased in our case. I would have thought the mass flow would take on this higher value initially, and then actually decrease back down to a lower mass flow rate, because we've got a low density gas downstream of the saturated zone, and a high density liquid upstream. If the volumetric flow is constant across the bed, then surely the mass flow downstream of the saturated zone will be lower than the upstream flow rates, due to the lower density?

So what I thought was - in the saturated zone we would see a mass flow rate higher than the inlet liquid mass flow, and as evaporation takes place this mass flow rate would gradually decrease back down to the inlet flow, and possibly even lower due to the lower vapour density. Incorrect?

 

[Chestermiller]

Ok so in tanks where there is evaporation taking place, we will have a mixture of liquid and vapour. The liquid will therefore travel at the faster velocity (the vapour velocity). As the liquid density is greater than the vapour density, this will result in an increase in the mass flow rate.

Question on that - why would the liquid take the same value as the velocity of the gas in the saturated zone, and not at an 'averaged' speed i.e. $v_{mix} = (1-X)v_l + Xv_g$

Hmm well yes if the mass flow out of the saturated zone is increased then you would expect this value to propagate downstream.

However I don't see why it would be increased in our case. I would have thought the mass flow would take on this higher value initially, and then actually decrease back down to a lower mass flow rate, because we've got a low density gas downstream of the saturated zone, and a high density liquid upstream. If the volumetric flow is constant across the bed, then surely the mass flow downstream of the saturated zone will be lower than the upstream flow rates, due to the lower density?

So what I thought was - in the saturated zone we would see a mass flow rate higher than the inlet liquid mass flow, and as evaporation takes place this mass flow rate would gradually decrease back down to the inlet flow, and possibly even lower due to the lower vapour density. Incorrect?

I think you have downstream and upstream switched. What does your model predict about $\dot{m}$ leaving the various tanks as a function of the tank number (i.e., spatially) at times when there is still pure liquid remaining in some of the (downstream) tanks?

 

[casualguitar]

I think you have downstream and upstream switched. What does your model predict about $\dot{m}$ leaving the various tanks as a function of the tank number (i.e., spatially) at times when there is still pure liquid remaining in some of the (downstream) tanks?

Apologies for the delay I found some script errors that I needed to fix.

And yes you were right I did have them switched. Ok yes now get the idea I think. The high speed gas upstream should 'push' the liquid downstream of it. The gas and liquid will have to travel at the same velocity, and therefore increase the mass flow rate of liquid to a value above 0.05kg/s, due to its higher density.

Here is what the spatial model predicts about $\dot{m}$.

Flow exiting each tank as a function of time. I've used 5 tanks:

and dm/dt:

So I understand in theory how we would see an increased flow downstream of the saturation zone, but I don't understand how (or if) we have incorporated that properly in our model. I'll think about this.

I can clearly see though there is a spike in mass flow in each tank, and that spike gets progressively smaller as we move through the tanks. Hmm I guess this has to do with there being less liquid downstream in the later tanks, but again how we have incorporated that in the model I'm not sure. I will think about this before the morning

 

[Chestermiller]

Apologies for the delay I found some script errors that I needed to fix.

And yes you were right I did have them switched. Ok yes now get the idea I think. The high speed gas upstream should 'push' the liquid downstream of it. The gas and liquid will have to travel at the same velocity, and therefore increase the mass flow rate of liquid to a value above 0.05kg/s, due to its higher density.

Here is what the spatial model predicts about $\dot{m}$.

Flow exiting each tank as a function of time. I've used 5 tanks:
View attachment 294928and dm/dt:
View attachment 294927

So I understand in theory how we would see an increased flow downstream of the saturation zone, but I don't understand how (or if) we have incorporated that properly in our model. I'll think about this.

I can clearly see though there is a spike in mass flow in each tank, and that spike gets progressively smaller as we move through the tanks. Hmm I guess this has to do with there being less liquid downstream in the later tanks, but again how we have incorporated that in the model I'm not sure. I will think about this before the morning

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

 

[casualguitar]

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

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

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

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

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

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