Modelling of two phase flow in packed bed using conservation equations

[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

 

[casualguitar]

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

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

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

 

[Chestermiller]

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

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

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

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

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

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

I'm still concerned about spikes in the results.

 

[Chestermiller]

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

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

 

[Chestermiller]

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

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

 

[casualguitar]

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

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

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

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

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

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

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

 

[Chestermiller]

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

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

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

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

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

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

 

[casualguitar]

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

Aha I see, time average not spatial. Got it

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

Is this correct? This is a fast change if so

 

[Chestermiller]

Aha I see, time average not spatial. Got it

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

Is this correct? This is a fast change if so

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

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

 

[Chestermiller]

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

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

 

[casualguitar]

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

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

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

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

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

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

 

[Chestermiller]

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

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

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

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

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

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

 

[casualguitar]

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

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

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

Time to return to the list of remaining tasks?

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

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

 

[Chestermiller]

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

Time to return to the list of remaining tasks?

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

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

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

 

[casualguitar]

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

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

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

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

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

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

 

[Chestermiller]

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

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

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

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

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

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

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