Modelling of two phase flow in packed bed using conservation equations

[casualguitar]

I've spent the last few hours reading papers that include axial dispersion. I understand now why this isn't trivial. There are lots of correlations we could use to find a thermal axial dispersion coefficient (as a function of thermal conductivity, porosity etc).

This paper is quite similar to what we're doing, and it uses an axial dispersion coefficient correlation: https://www.sciencedirect.com/science/article/abs/pii/0038092X82902067

So essentially it seems to be 'lumped in' with the conduction term, by using a pseudoproperty $K_Z$ rather than the actual thermal conductivity.

Is this what you had in mind?

Also while reading I found this paper, which is sort of close to what we're doing: https://www.sciencedirect.com/science/article/pii/000925097480059X?via=ihub

I will take a close read of both of these tomorrow morning to prepare for this discussion

 

[Chestermiller]

I've spent the last few hours reading papers that include axial dispersion. I understand now why this isn't trivial. There are lots of correlations we could use to find a thermal axial dispersion coefficient (as a function of thermal conductivity, porosity etc).

This paper is quite similar to what we're doing, and it uses an axial dispersion coefficient correlation: https://www.sciencedirect.com/science/article/abs/pii/0038092X82902067

So essentially it seems to be 'lumped in' with the conduction term, by using a pseudoproperty $K_Z$ rather than the actual thermal conductivity.

Is this what you had in mind?

Yes. You may have to make this an adjustable parameter when you tune your model to the experimental data.

Also while reading I found this paper, which is sort of close to what we're doing: https://www.sciencedirect.com/science/article/pii/000925097480059X?via=ihub

I will take a close read of both of these tomorrow morning to prepare for this discussion

This paper doesn't seem to include a 2nd fluid phase.

 

[casualguitar]

Yes. You may have to make this an adjustable parameter when you tune your model to the experimental data.

Understood

This paper doesn't seem to include a 2nd fluid phase.

Yes you're absolutely right I linked the wrong one. This is the correct one: https://www.sciencedirect.com/science/article/pii/001793109090255S

This paper does assume axial dispersion is negligible so it won't be of use for that. It does model two phase flow though

 

[Chestermiller]

Understood

Yes you're absolutely right I linked the wrong one. This is the correct one: https://www.sciencedirect.com/science/article/pii/001793109090255S

This paper does assume axial dispersion is negligible so it won't be of use for that. It does model two phase flow though

I have access to only the abstract.

It looks like they are taking into account pressure variations, while we are assuming constant pressure.

In our analysis, a complexity in writing down the finite difference expressions for the spatial direction is that the density and mass rate of flow are both functions of time and spatial position; this needs to be accounted for in the mass balance equation and in the energy balance equation. This is what I've been working on.

 

[casualguitar]

I've set up a WeTransfer (allows sending of the PDF) link if you want to look at the paper: https://we.tl/t-9ebfBtGiuR

Yes correct they are. The library I am using can do temperature and pressure dependent parameters, this is an overview of the documentation if you want to take a look at the functionality: https://thermo.readthedocs.io

In our analysis, a complexity in writing down the finite difference expressions for the spatial direction is that the density and mass rate of flow are both functions of time and spatial position; this needs to be accounted for in the mass balance equation and in the energy balance equation. This is what I've been working on.

Ok so are you saying we will not just have two energy balance equations, but also a mass balance to the fluid equation, meaning our 'system' will involve solving three ODEs?

 

[Chestermiller]

I've set up a WeTransfer (allows sending of the PDF) link if you want to look at the paper: https://we.tl/t-9ebfBtGiuR

Yes correct they are. The library I am using can do temperature and pressure dependent parameters, this is an overview of the documentation if you want to take a look at the functionality: https://thermo.readthedocs.ioOk so are you saying we will not just have two energy balance equations, but also a mass balance to the fluid equation, meaning our 'system' will involve solving three ODEs?

Not exactly. I'll get to that soon.

But, before getting to the full problem, I'd like you to try something, please.

Try modifying the stirred tank version of the numerical model by solving the fluid heat balance equation $m\frac{dh}{dt}=\dot{m}_{in}(h_{in}-h)-UA(T(h)-T_S)$ (from post #12) in terms of the enthalpy h rather than the temperate T. (The heat balance for the solid can be left in terms of T and TS). The reason I'm asking this is that I think there would be an advantage to working in terms of h in the full model because the heat capacity is discontinuous, which would lead to logistical complexities in evaluating spatial derivatives numerically; such complexities could be avoided by working in terms of h (which is continuous). If this doesn't make a lot of sense, please take me at my word. Once the enthalpy is known, the temperature follows directly (in our scheme). Please make comparison plots so we can see how the results for the temperature dependence compare. Thanks.

 

[casualguitar]

Try modifying the stirred tank version of the numerical model by solving the fluid heat balance equation mdhdt=m˙in(hin−h)−UA(T(h)−TS) (from post #12) in terms of the enthalpy h rather than the temperate T. (The heat balance for the solid can be left in terms of T and TS).

Will do. Just some questions to clarify:
1) Am I to solve this for a single phase only, or for two phase flow?
2) The equation you supplied: $m\frac{dh}{dt}=\dot{m}_{in}(h_{in}-h)-UA(T(h)-T_S)$ seems to be in a useful form already. So if we are working on a single phase model then I actually only need an equation of state here I think, H(T) and T(H), and then we can solve using the ODE solver as usual

Is the above correct? i.e. your equation is in a useful form and I just need an EOS (thermodynamic property library) to solve?

 

[Chestermiller]

Will do. Just some questions to clarify:
1) Am I to solve this for a single phase only, or for two phase flow?

Two phase.

2) The equation you supplied: $m\frac{dh}{dt}=\dot{m}_{in}(h_{in}-h)-UA(T(h)-T_S)$ seems to be in a useful form already. So if we are working on a single phase model then I actually only need an equation of state here I think, H(T) and T(H), and then we can solve using the ODE solver as usual

Is the above correct? i.e. your equation is in a useful form and I just need an EOS (thermodynamic property library) to solve?

You use the property equations for M and h that we have written for the 2 phase model. It's easiest to express T and M as functions of h over the three intervals.

 

[casualguitar]

You use the property equations for M and H that we have written for the 2 phase model.

Are we assuming $C_{PL}$ and $C_{PG}$ are constant?

So the differences in terms of calculations here would be:
1) We are first calculating the inlet enthalpy using the previously developed property equations
2) The mass holdup will now be a function of enthalpy instead of temperature, meaning that we should solve for $T_{sat}$ in each of the three enthalpy property equations
3) At each time step we solve for the fluid $H$, and then convert to $H$ to $T$ for input into the solid equation
4) Also we would need the enthalpy at $T1$ and $T2$ rather than the temperature bounds

I will start on this solution this evening however I'm not sure why this enthalpy model would output different results to the temperature model. Why would it?

I'll do the necessary equation adjusting and comment with a clear algorithm before I start coding

 

[Chestermiller]

Are we assuming $C_{PL}$ and $C_{PG}$ are constant?

So the differences in terms of calculations here would be:
1) We are first calculating the inlet enthalpy using the previously developed property equations
2) The mass holdup will now be a function of enthalpy instead of temperature, meaning that we should solve for $T_{sat}$ in each of the three enthalpy property equations
3) At each time step we solve for the fluid $H$, and then convert to $H$ to $T$ for input into the solid equation
4) Also we would need the enthalpy at $T1$ and $T2$ rather than the temperature bounds

I will start on this solution this evening however I'm not sure why this enthalpy model would output different results to the temperature model. Why would it?

I'll do the necessary equation adjusting and comment with a clear algorithm before I start coding

OK. Let's do the interval thing with this approach.

1. If $h\leq h_{sat,L}$ , then $$T=T_{sat}+(h-h_{sat,L})/C_{PL}$$ and $$m=\rho_L V$$

2. If $h_{sat,L}\leq h \leq h_{sat,V}$, then $$T=T_{sat}$$and $$m=\frac{V}{\left[\frac{1}{\rho_L}+\left(\frac{RT_{sat}}{PM}-\frac{1}{\rho_L}\right)\left(\frac{(h-h_{sat,L})}{(h_{sat,V}-h_{sat,L})}\right)\right]}$$

3. If $h_{sat,V}\leq h$, then $$T=T_{sat}+(h-h_{sat,V})/C_{PV}$$ and $$m=\frac{PM}{RT} V$$ with T given by the previous equation.

Hmmm...This method gets us out of the artificial approach of using T1 and T2, and is "exact" for the relations between enthalpy, temperature. and mass.

We expect that the results of this calculation will match the results from the temperature-version, so it should provide a good validation of the model. Plus, it will show the effect of the temperature remaining exactly constant with time at the saturation temperature while the fluid is converting from vapor to liquid.

Note finally, that, with no loss of generality, we can simplify further by taking the reference temperature for zero enthalpy of the liquid as $T_{sat}$, such that $h_{sat,L}=0$ and $h_{sat,V}=\Delta h_{vaporization}$.

 

[casualguitar]

Note finally, that, with no loss of generality, we can simplify further by taking the reference temperature for zero enthalpy of the liquid as Tsat, such that hsat,L=0 and hsat,V=Δhvaporization.

Got it, so just adding that in:

1. If $h\leq h_{sat,L}$, then
$$T=T_{sat}+h/C_{PL}$$
$$m=\rho_L V$$

2. If $h_{sat,L}\leq h \leq h_{sat,V}$ then
$$T=T_{sat}$$
$$m=\frac{V}{\left[\frac{1}{\rho_L}+\left(\frac{RT_{sat}}{PM}-\frac{1}{\rho_L}\right)\left(\frac{h}{\Delta h_{vap}}\right)\right]}$$

3. If $h_{sat,V}\leq h$ then
$$T=T_{sat}+(h-\Delta h_{vap})/C_{PV}$$
$$m=\frac{PM}{RT} V$$

We expect that the results of this calculation will match the results from the temperature-version, so it should provide a good validation of the model.

The post #38 model yes?

This method gets us out of the artificial approach of using T1 and T2, and is "exact" for the relations between enthalpy, temperature. and mass.

It does, I guess dealing with a mixture means that we should get an average value for $T_{sat}$ here (N2 boils at 77K, O2 boils at 90K) so maybe a Tsat of 73.5K is suitable here? The same for any other constant I could take the average of the N2 and O2 values.

I think I follow this model fully so, I will aim to implement this tomorrow morning!

 

[Chestermiller]

Got it, so just adding that in:

More like

1. If $h\leq 0$, then
$$T=T_{sat}+h/C_{PL}$$
$$m=\rho_L V$$

2. If $0\leq h \leq \Delta h_{vap}$ then
$$T=T_{sat}$$
$$m=\frac{V}{\left[\frac{1}{\rho_L}+\left(\frac{RT_{sat}}{PM}-\frac{1}{\rho_L}\right)\left(\frac{h}{\Delta h_{vap}}\right)\right]}$$

3. If $\Delta h_{vap}\leq h$ then
$$T=T_{sat}+(h-\Delta h_{vap})/C_{PV}$$
$$m=\frac{PM}{RT} V$$

The post #38 model yes?

Yes.

It does, I guess dealing with a mixture means that we should get an average value for $T_{sat}$ here (N2 boils at 77K, O2 boils at 90K) so maybe a Tsat of 73.5K is suitable here? The same for any other constant I could take the average of the N2 and O2 values.

Please use your best judgment for this. I guess at some point the model can be modified to more the precise VLE behavior. I think the 73.5 is a typo, and you meant 83.5, right?

 

[Chestermiller]

I think I can solve for the case with U = 0 analytically.

 

[casualguitar]

Working on the above model now I should have results for this in an hour or so. I'll change the inputs of the model we're comparing to so that they're close to equivalent. What would a U = 0 model show us? The effect of having zero heat transfer between fluid and solid?

 

[Chestermiller]

Working on the above model now I should have results for this in an hour or so. I'll change the inputs of the model we're comparing to so that they're close to equivalent. What would a U = 0 model show us? The effect of having zero heat transfer between fluid and solid?

It would tell us the effect of purging liquid initially at a low temperature from the lumped mixed tank with a vapor at a high temperature.

 

[casualguitar]

It would tell us the effect of purging liquid initially at a low temperature from the lumped mixed tank with a vapor at a high temperature.

So I made the above model, and ran it for ranges of input parameters etc. I've printed the output and I can see the expected increase of temperature to the saturation temperature, a (short) constant temperature zone, and then an increase up to the fluid inlet temperature. I've checked that the correct enthalpy/mass equations are being used for the right enthalpy ranges. I've also made this models inputs as close as possible to the old model inputs. This one produces results that are quite a bit different, in that the temperature of the fluid and gas increase at a much slower rate. I'm attaching the output to model 1 and 2, and also the code to the new model below. Note I plotted over 10000s as this was the time required to see a levelling out of the temperature curve in the new model. I'm assuming for now there is a bug in my code somewhere so I'll hopefully find it this evening!

New model:

Old model:

Enthalpy model code:

 

[Chestermiller]

I don’t see anything like a constant temperature zone in either plot. What were the parameter values you used, including the mass flow rate in, and the solid mass and heat capacity?

 

[casualguitar]

I don’t see anything like a constant temperature zone in either plot. What were the parameter values you used, including the mass flow rate in, and the solid mass and heat capacity?

There is a constant temperature zone in the new plot, however it’s very small (I think it spans only a few seconds). I’m away from my laptop currently so I will zoom in on this in the plot later this evening. In regards to mass flow in, solid mass and heat capacity, these are minL, mS and cpS in the first code image. I’ll also change these parameters this evening to get a longer constant temperature zone

 

[Chestermiller]

There is a constant temperature zone in the new plot, however it’s very small (I think it spans only a few seconds). I’m away from my laptop currently so I will zoom in on this in the plot later this evening. In regards to mass flow in, solid mass and heat capacity, these are minL, mS and cpS in the first code image. I’ll also change these parameters this evening to get a longer constant temperature zone

You omitted multiplying by the bed volume when you calculated the mass of fluid in the bed for pure vapor. This is an indication that you really need to spend some time testing and checking your code.

What are the units of the physical parameters in your code?

Heat transfer coefficeint
Heat transfer area
Heat capacity of bed
mass of bed

 

[casualguitar]

You omitted multiplying by the bed volume when you calculated the mass of fluid in the bed for pure vapor.

Thanks!

This is an indication that you really need to spend some time testing and checking your code.

Yes I agree completely. I did spend over an hour looking, however this is really a bug I should have been able to find if I'd spent some more time. I'll take more care to find these myself from now on

What are the units of the physical parameters in your code?

Units included in the comments here:

Updated graph below. Shows the constant temperature zone and levels out at the same rate:

I'm going to make a graph of this one over the previous model graph to see how well they actually overlay. However from visual inspection of the individual graphs they look ok

 

[Chestermiller]

It seems to me what is happening in this calculation is that the heat transfer resistance between the gas an the bed is very low (equivalent to very high heat transfer coefficient) so that the gas always comes very close to the bed temperature. Moreover, once the temperature rises above $T_{sat}$, there is all gas in the tank, and the thermal inertia of this gas is very low compared to the thermal inertia of the bed. So it can be neglected. So the rate at which the incoming gas gives up heat to the bed is $\dot{m}C_{PV}(T_{in}-T_S)$ and the heat balance on the bed in this range of temperatures simplifies to $$m_SC_{PS}\frac{dT_S}{dt}=\dot{m}C_{PV}(T_{in}-T_S)$$
This same result can be obtained by (a) eliminating the heat transfer term between the gas heat balance and solid heat balance, then (b) setting the instantaneous gas temperature in the bed equal to the bed temperature and then (c) neglecting the thermal inertia of the gas.

The characteristic time constant for the system in this range of temperatures is then $$\tau=\frac{m_SC_{PS}}{\dot{m}C_{PV}}$$What is the value that you calculate for this time constant for your model inputs?

 

[casualguitar]

What is the value that you calculate for this time constant for your model inputs?

Just adding the comparison graph here as mentioned before

Moreover, once the temperature rises above Tsat, there is all gas in the tank, and the thermal inertia of this gas is very low compared to the thermal inertia of the bed. So it can be neglected.

Understood

This same result can be obtained by (a) eliminating the heat transfer term between the gas heat balance and solid heat balance, then (b) setting the instantaneous gas temperature in the bed equal to the bed temperature and then (c) neglecting the thermal inertia of the gas.

Understood also.

The characteristic time constant for the system in this range of temperatures is then τ=mSCPSm˙CPVWhat is the value that you calculate for this time constant for your model inputs?

$M_S$ = 1, $C_{PS}$ = 10, $mass flow (kg/s) = 0.05$, $C_{PV}$ = 1, which returns the time constant:
$$\tau=200$$

I haven't looked into what a time constant of 200s means yet (quite late here). My initial sense is that it is very large. I will do this first thing tomorrow

I've also spotted an error in my previous model. The heat capacity of the solid is higher than any reasonable value for materials of its type (typical values are between 0.8 and 2). I'll rerun with a reasonable value tomorrow morning.

Looking at the time constant equation, a solid heat capacity of 1, for example, reduces the time constant to 20s which seems a lot more reasonable.

As I say I'll make these changes in the morning

 

[Chestermiller]

View attachment 292463
Just adding the comparison graph here as mentioned before

UnderstoodUnderstood also.
$M_S$ = 1, $C_{PS}$ = 10, $mass flow (kg/s) = 0.05$, $C_{PV}$ = 1, which returns the time constant:
$$\tau=200$$

I haven't looked into what a time constant of 200s means yet (quite late here). My initial sense is that it is very large. I will do this first thing tomorrow

I've also spotted an error in my previous model. The heat capacity of the solid is higher than any reasonable value for materials of its type (typical values are between 0.8 and 2). I'll rerun with a reasonable value tomorrow morning.

Looking at the time constant equation, a solid heat capacity of 1, for example, reduces the time constant to 20s which seems a lot more reasonable.

As I say I'll make these changes in the morning

Nice. Please make a semi-log plot of $T_{in}-T$ as a function of t, and let's see what you get.

I should also mention that the R in the equations I wrote is the universal gas constant, not the customized value for air. This has caused the masses of vapor you have used in the calculations to be a factor of 28 too high.

 

[casualguitar]

Nice. Please make a semi-log plot of $T_{in}-T$ as a function of t, and let's see what you get.

I should also mention that the R in the equations I wrote is the universal gas constant, not the customized value for air. This has caused the masses of vapor you have used in the calculations to be a factor of 28 too high.

Updating the overlay plot here with a solid heat capacity of $1.1 kJ/kg.K$. As expected the temperature levels off at a faster rate

I should also mention that the R in the equations I wrote is the universal gas constant, not the customized value for air.

Got it. So I should divide my value by the molecular weight, meaning that the new units for $R$ are $J/mol.K$?

Making the semi-log plot now

 

[casualguitar]

Nice. Please make a semi-log plot of $T_{in}-T$ as a function of t, and let's see what you get.

$T_{in} - T$ plot for the enthalpy model:

To be clear, what I did here was took the final temperature arrays for the fluid and solid, and then created two new arrays $(T_{in}-T_{fluid})$ and $(T_{in}-T_{solid})$ and plotted these versus $t$

 

[Chestermiller]

Updating the overlay plot here with a solid heat capacity of $1.1 kJ/kg.K$. As expected the temperature levels off at a faster rate
View attachment 292479Got it. So I should divide my value by the molecular weight, meaning that the new units for $R$ are $J/mol.K$?

Making the semi-log plot now

No. You multiply your value by 28 kg/(k-mole) to get $R=8314\ J/(k-mole\ K)$

 

[Chestermiller]

$T_{in} - T$ plot for the enthalpy model:
View attachment 292480
To be clear, what I did here was took the final temperature arrays for the fluid and solid, and then created two new arrays $(T_{in}-T_{fluid})$ and $(T_{in}-T_{solid})$ and plotted these versus $t$

Good. Now have your plotting software change the ordinate scale to logarithmic.

 

[casualguitar]

No. You multiply your value by 28 kg/(k-mole) to get $R=8314\ J/(k-mole\ K)$

Ok, new plot below with the updated R value. The effect of the updated R value is that the time taken to reach steady state is approximately cut in half:

Good. Now have your plotting software change the ordinate scale to logarithmic.

Both time and temperature axes to log scale?

 

[casualguitar]

Good. Now have your plotting software change the ordinate scale to logarithmic.

Here is the plot when both x and y are plotted in log scale:

What is the value of this plot? It shows us that at the later times (in the gas phase) the temperatures rises more rapidly towards the inlet gas temperature, meaning that (as you said) the thermal inertia of the gas is very low compared to the thermal inertia of the bed?

 

[Chestermiller]

Ok, new plot below with the updated R value. The effect of the updated R value is that the time taken to reach steady state is approximately cut in half:
View attachment 292481

Both time and temperature axes to log scale?

no. Just temperature.

 

[Chestermiller]

Here is the plot when both x and y are plotted in log scale:
View attachment 292482
What is the value of this plot? It shows us that at the later times (in the gas phase) the temperatures rises more rapidly towards the inlet gas temperature, meaning that (as you said) the thermal inertia of the gas is very low compared to the thermal inertia of the bed?

just log scale on temperature, not time.

 

[casualguitar]

no. Just temperature.

Here is the plot when only temperature is plotted in log scale:

Actually, I extended the range out to 10k seconds and this happens:

I'm going to see if any values have unexpected behaviour at about 5k seconds.

Actually, has this got something to do with stability, given that the derivatives will approach zero around this range?

 

[Chestermiller]

Here is the plot when only temperature is plotted in log scale:
View attachment 292485

Actually, I extended the range out to 10k seconds and this happens:
View attachment 292486
I'm going to see if any values have unexpected behaviour at about 5k seconds.

Actually, has this got something to do with stability, given that the derivatives will approach zero around this range?

Don't you think that the issue here is one of precision (round-off error). Do you really feel that it is necessary to take the temperature to within 0.001 degrees of the inlet temperature? In my judgement, it is adequate to show the results out to about 1500 sec. If you want better precision than that, go to double precision.

 

[casualguitar]

Don't you think that the issue here is one of precision (round-off error). Do you really feel that it is necessary to take the temperature to within 0.001 degrees of the inlet temperature? In my judgement, it is adequate to show the results out to about 1500 sec. If you want better precision than that, go to double precision.

Yes exactly that's what I mean also. No its not necessary to go to within 0.001K of the inlet temperature I just thought it was interesting to see.

Plot to 1500s:

I'm going to collect each of these models into separate scripts now (to make a 'trail' of the model progression). Is there anything else important that I should do before we talk about introducing an axial temperature gradient?

 

[Chestermiller]

Yes exactly that's what I mean also. No its not necessary to go to within 0.001K of the inlet temperature I just thought it was interesting to see.

Plot to 1500s:
View attachment 292488

I'm going to collect each of these models into separate scripts now (to make a 'trail' of the model progression). Is there anything else important that I should do before we talk about introducing an axial temperature gradient?

Is there a reason that the graph doesn't show a constant gas temperature during the phase change? Have you incorporated the change to the gas constant noted in post #76 yet?

Hold your horses. What is the slope of the straight line on the graph at times > 200 seconds, in terms of $\Delta {lnT}/\Delta t$? The reciprocal of this is the time constant for the response in the vapor region. What do you get for the time constant? How does this compare with our expectations?

 

[casualguitar]

Is there a reason that the graph doesn't show a constant gas temperature during the phase change? Have you incorporated the change to the gas constant noted in post #76 yet?

Yes it seems to be because I was plotting 100 points only, and the phase change zone is quite short so it mostly bypassed it in the context of a 1500s plot. A higher resolution plot below shows the constant temperature phase change zone. Its quite small:

And zooming in on the phase change zone:

Hold your horses. What is the slope of the straight line on the graph at times > 200 seconds, in terms of ΔlnT/Δt? The reciprocal of this is the time constant for the response in the vapor region. What do you get for the time constant? How does this compare with our expectations?

The slope for t>200:
$$\Delta {lnT}/\Delta t = -0.16$$

The time constant (taking the absolute value of the slope):
$$\tau = 6.09s$$

Hmm I'm not sure how this compares. So, it is smaller than the previously calculated value of 15s using this equation:
$$\tau=\frac{m_SC_{PS}}{\dot{m}C_{PV}}$$

Is this value in the realm of 'close enough for the purpose of this model'?

One thing that does concern me is that there is a sharp change in slope at around 26s in the vapour region. Why would this be? We would probably get a higher slope value in the vapour region if this was an error in the model, meaning that the estimated tau and actual tau would also be closer

 

[Chestermiller]

Yes it seems to be because I was plotting 100 points only, and the phase change zone is quite short so it mostly bypassed it in the context of a 1500s plot. A higher resolution plot below shows the constant temperature phase change zone. Its quite small:
View attachment 292493
And zooming in on the phase change zone:
View attachment 292494

Thanks. Much better. How come, not that you've lowered the heat capacity of the bed, the bed temperature does not more closely approach the gas saturation temperature in this zone. It did in earlier plots before the bed heat capacity was lowered.

The slope for t>200:
$$\Delta {lnT}/\Delta t = -0.16$$

The time constant (taking the absolute value of the slope):
$$\tau = 6.09s$$

I get a temperature difference change by a factor of 10 every 700 seconds. I guess I made an error in specifying the slope I wanted. What I really wanted was $\Delta \ln{(T_{in}-T)}/\Delta t$. This gives a time constant of about 300 seconds.

Hmm I'm not sure how this compares. So, it is smaller than the previously calculated value of 15s using this equation:
$$\tau=\frac{m_SC_{PS}}{\dot{m}C_{PV}}$$

Is this value in the realm of 'close enough for the purpose of this model'?

Now we are comparing 300 seconds with 15 seconds. This seems too large. Please provide the latest version of your enthalpy code.

One thing that does concern me is that there is a sharp change in slope at around 26s in the vapour region. Why would this be? We would probably get a higher slope value in the vapour region if this was an error in the model, meaning that the estimated tau and actual tau would also be closer

I am concerned about the sharp change too. Needs more study.

 

[casualguitar]

How come, not that you've lowered the heat capacity of the bed, the bed temperature does not more closely approach the gas saturation temperature in this zone.

It seems that decreasing the value of R (back to my incorrect value which was 28 times smaller) has the effect of making the approach temperature value smaller. This is the plot with my incorrect R value:

I am concerned about the sharp change too. Needs more study.

I'll look into this now

Please provide the latest version of your enthalpy code.

I am looking for the reason for the sharp change now. I guess this will be different to the reason the time constant value is different

 

[Chestermiller]

It seems that decreasing the value of R (back to my incorrect value which was 28 times smaller) has the effect of making the approach temperature value smaller.

Im not sure I know what you mean by the "approach temperature." The lower value of R makes the vapor density higher so that it has more thermal inertia. This would be roughly equivalent to the effect of increasing the heat capacity of the vapor.

This is the plot with my incorrect R value:
View attachment 292530

I'll look into this now

Run some diagnostic calculations where you print out things like h and m vs t to see if they are behaving the way you expect them to behave.

View attachment 292531
View attachment 292532

I am looking for the reason for the sharp change now. I guess this will be different to the reason the time constant value is different

The diagnostic calculations should provide some insight into this as well.

It looks like your latest code is using the high value for CS again. Is that what you intended?

 

[casualguitar]

Im not sure I know what you mean by the "approach temperature." The lower value of R makes the vapor density higher so that it has more thermal inertia. This would be roughly equivalent to the effect of increasing the heat capacity of the vapor.

You asked previously why "the bed temperature does not more closely approach the gas saturation temperature in this zone", that is what I mean by the approach temperature i.e. the difference between the bed temperature and the gas saturation temperature

Run some diagnostic calculations where you print out things like h and m vs t to see if they are behaving the way you expect them to behave.

Got it, working on this now

 

[casualguitar]

Run some diagnostic calculations where you print out things like h and m vs t to see if they are behaving the way you expect them to behave.

I ran some tests to see if this is expected behaviour, and it seems to be. These plots are of $T_{fluid}$, $T_{solid}$, $H_{fluid}$ and mass holdup versus time. I plotted it like this to see what other variables changed significantly around the time temperature does the weird jump:

Zoomed in version around the saturation zone:

The 'resolution' of my plot was too low. I was calculating one data point for each time value, nothing in between. So the temperature seemed to jump up randomly in this case. However when I increased the number of points (10 points between every second), the jump is really a curve and it happens exactly when the enthalpy goes above the heat of vaporisation (h>200):

This seems to answer the question about the sudden temperature jump. Looking into the time constant difference now

 

[Chestermiller]

I ran some tests to see if this is expected behaviour, and it seems to be. These plots are of $T_{fluid}$, $T_{solid}$, $H_{fluid}$ and mass holdup versus time. I plotted it like this to see what other variables changed significantly around the time temperature does the weird jump:
View attachment 292572
Zoomed in version around the saturation zone:
View attachment 292573
The 'resolution' of my plot was too low. I was calculating one data point for each time value, nothing in between. So the temperature seemed to jump up randomly in this case. However when I increased the number of points (10 points between every second), the jump is really a curve and it happens exactly when the enthalpy goes above the heat of vaporisation (h>200):

View attachment 292576

This seems to answer the question about the sudden temperature jump. Looking into the time constant difference now

Very nice. Excellent

For the mass holdup, it would more instructive to examine it on a logarithmic scale.

 

[casualguitar]

For the mass holdup, it would more instructive to examine it on a logarithmic scale.

Understood yes the values are tiny in comparison to Temperature/Enthalpy. I'll change the scale for the mass holdup. Although it seems to do what is expected i.e. start at 8kg (V*rhoL = 0.01*800), gradually reduce across the phase change zone and then level out in the vapour zone

For the time constant, are we assuming this equation produces an accurate time constant?:
$$\tau=\frac{m_SC_{PS}}{\dot{m}C_{PV}}$$

And that the model should produce a value close to this if its correct? i.e. close to 200

I'll check the slope now with the higher resolution model

 

[casualguitar]

I get a temperature difference change by a factor of 10 every 700 seconds. I guess I made an error in specifying the slope I wanted. What I really wanted was Δln⁡(Tin−T)/Δt. This gives a time constant of about 300 seconds.

So using this equation:
$$\tau=\frac{m_SC_{PS}}{\dot{m}C_{PV}}$$

And the updated $C_s$ value of 1.1 (I was actually still using the old value up until recently as you spotted):
$$\tau=22s$$
And with these equations:
$$ slope = \Delta \ln{(T_{in}-T)}/\Delta t$$
$$\tau = 1/slope$$
I get:
$$\tau=18s$$
So there is relatively good similarity now it seems

Heres the code snippet:

The questions I would have here are:
1) Where do the two equations for tau come from?
2) Is it ok that I've only calculated the slope in the vapour region? i.e. I haven't found an 'averaged' slope across all regions?

 

[Chestermiller]

So using this equation:
$$\tau=\frac{m_SC_{PS}}{\dot{m}C_{PV}}$$

And the updated $C_s$ value of 1.1 (I was actually still using the old value up until recently as you spotted):
$$\tau=22s$$
And with these equations:
$$ slope = \Delta \ln{(T_{in}-T)}/\Delta t$$
$$\tau = 1/slope$$
I get:
$$\tau=18s$$
So there is relatively good similarity now it seems

Heres the code snippet:
View attachment 292577

The questions I would have here are:
1) Where do the two equations for tau come from?

See post #71. The following is an approximation for the vapor region only, based on the approximation in post #71:
$$m_SC_{PS}\frac{dT_S}{dt}=\dot{m}C_{PV}(T_{in}-T_S)$$The solution to this equation is $$\frac{T_{in}-T_S}{T_{in}-T_{sat}}=e^{-t/\tau}$$where t is measured from the time that we have all vapor, with $$\tau=\frac{m_SC_{PS}}{\dot{m}C_{PV}}$$

In this approximation, the vapor temperature T is approximately equal to TS.

2) Is it ok that I've only calculated the slope in the vapour region? i.e. I haven't found an 'averaged' slope across all regions?

Yes. The approximation applies only to the pure vapor.

 

[casualguitar]

See post #71. The following is an approximation for the vapor region only, based on the approximation in post #71:
$$m_SC_{PS}\frac{dT_S}{dt}=\dot{m}C_{PV}(T_{in}-T_S)$$The solution to this equation is $$\frac{T_{in}-T_S}{T_{in}-T_{sat}}=e^{-t/\tau}$$where t is measured from the time that we have all vapor, with $$\tau=\frac{m_SC_{PS}}{\dot{m}C_{PV}}$$

In this approximation, the vapor temperature T is approximately equal to TS.Yes. The approximation applies only to the pure vapor.

Ah I understand, and slope here is $\frac{-1}{\tau}$?

I'll clean up the code so far and split it into the various models done so far. Are there other things I should add before we spatial variation?

Actually I'm also going to map this timeline of model development into a powerpoint

 

[Chestermiller]

Ah I understand, and slope here is $\frac{-1}{\tau}$?

I'll clean up the code so far and split it into the various models done so far. Are there other things I should add before we spatial variation?

Actually I'm also going to map this timeline of model development into a powerpoint

I think the next step is to move on to spatial variation. I'll write something up tomorrow.

 

[casualguitar]

I think the next step is to move on to spatial variation. I'll write something up tomorrow.

Looking forward to it! Again thank you Chet for walking through the modelling process with me, this is hugely appreciated.

 

[Chestermiller]

FORMULATION

Considering the fact that the fluid density is varying substantially with time and spatial position within this system, in my judgment, there is asignificant advantage to formulating and solving the heat balance equation for the fluid using the so-called divergence representation of the equation (see below); this also provides the added advantage of automatically converting mass not only in the PDE form of the equation, but also when it is expressed in finite difference form. The divergence form of the mass and energy balance for the fluid are thus written as: $$\frac{\partial \rho}{\partial t}+\frac{\partial (\rho u)}{\partial x}=0$$and$$\frac{\partial (\rho h)}{\partial t}+\frac{\partial (\rho uh)}{\partial x}=\frac{\partial}{\partial x}\left(\rho D\frac{\partial h}{\partial x}\right)+\frac{Ua}{\epsilon}(T_S-T)$$where D is the axial thermal dispersion coefficient (assumed to be much larger than the thermal diffusivity of the fluid at all typical operating conditions of the bed).

In these equations, if we substitute $\phi=\rho u$ for the axial mass flux, and express the thermal dispersion coefficient in the typical approximate form D = ul (where l is the characteristic thermal dispersion length parameter), these equations become:
$$\frac{\partial \rho}{\partial t}+\frac{\partial \phi}{\partial x}=0$$and$$\frac{\partial (\rho h)}{\partial t}+\frac{\partial (\phi h)}{\partial x}=\frac{\partial}{\partial x}\left(\phi l\frac{\partial h}{\partial x}\right)+\frac{Ua}{\epsilon}(T_S-T)$$
In lieu of any definite information on the effects of fluid density on the dispersivity parameter l, it will probably be necessary to assume that l is a constant depending only on packing geometry and to treat it as an adjustable 'tuning" parameter.

Thoughts so far?

 

[casualguitar]

in my judgment, there is asignificant advantage to formulating and solving the heat balance equation for the fluid using the so-called divergence representation of the equation (see below); this also provides the added advantage of automatically converting mass not only in the PDE form of the equation, but also when it is expressed in finite difference form.

Got it

∂ρ∂t+∂ϕ∂x=0and∂(ρh)∂t+∂(ϕh)∂x=∂∂x(ϕl∂h∂x)+Uaϵ(TS−T)

To describe the above equations I would say the following, where every term is per unit volume:
1) The rate of increase of mass plus the rate at which mass leaves in the $x$ direction $= 0$
2) The rate of increase of energy plus the rate at which energy leaves in the $x$ direction = the rate of axial dispersion in the x direction plus the heat transfer between the fluid and solid

(I notice we have arrived at the same equation I started the post with (!), just using enthalpy instead of temperature. Nice)

In these equations, if we substitute ϕ=ρu for the axial mass flux, and express the thermal dispersion coefficient in the typical approximate form D = ul (where l is the characteristic thermal dispersion length parameter), these equations become

Understood

In lieu of any definite information on the effects of fluid density on the dispersivity parameter l, it will probably be necessary to assume that l is a constant depending only on packing geometry and to treat it as an adjustable 'tuning" parameter.

Understood. I think I have seen some correlations for $D$ the axial dispersion coefficient but I guess this is of secondary importance for now.

The mass balance is straightforward, however I do have questions on the energy balance to clear up.

Points of confusion:
1) Are you using the term 'heat balance' and 'energy balance' interchangeably?
2) Is this an advection term: $\frac{\partial (\phi h)}{\partial x}$? I'm not familiar with advection, other than its definition as the transfer of heat through fluid flow, so I guess this is it
3) So the axial dispersion coefficient is a 'catch all' term for a number of energy transport processes as far as I know (diffusion and conduction within the fluid, convection by the fluid in the axial direction, axial and transverse mixing of the fluid, etc, as far as I know). What terms are included in this term in our case?
4) The final term seems to be the convection between the fluid and the solid term. Why does this term not have a time or space derivative, because it will surely depend on space and time? Is it because 'time/space dependence' is built in, in the sense that $T$ and $T_S$ will change with time and space, meaning that the convection term will also change with time and space indirectly?

In lieu of any definite information on the effects of fluid density on the dispersivity parameter l, it will probably be necessary to assume that l is a constant depending only on packing geometry and to treat it as an adjustable 'tuning" parameter.

Got it

So now instead of solving a fluid and solid energy balance, we will be solving the mass balance for the fluid, and the energy balance for the fluid/solid instead?