Modelling of two phase flow in packed bed using conservation equations

[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?

 

[Chestermiller]

Got it

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$

I would call it the net rate at which mass leaves

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 would call it "rate of increase of enthalpy per unit volume" and "the net rate at which enthalpy per unit volume leaves the control volume in the x direction"

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

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?

Yes.

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

As we said before, it represents the net advection of enthalpy per unit volume out of the control volume.

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?

Like you said, it is a catch all for the net of all these effects. It's probably not going to be worthwhile trying to quantify it to a more detailed extent, and, in my judgment, it will be preferable to treat the dispersion length l as an adjustable parameter that the experimental results will be used to calibrate.

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?

No. The heat transfer coefficient will depend on the fluid parameters, velocities, and bed geometry, so, in reality, it can be a function of location and time. I too will be difficult to quantify, and I recommend, at least initially, treating it as a constant and, like the dispersibity, calibrating it to the experimental results. (This is how we treated it in the lumped model) Later, if necessary, we can "sharpen our pencils."

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?

No. The solid energy balance is still used. We have just had the add the mass balance because the density is a function of x and t (which needs to be accounted for, since the density varies by a factor of 1000 from liquid to vapor).

 

[casualguitar]

I would call it the net rate at which mass leaves

I would call it "rate of increase of enthalpy per unit volume" and "the net rate at which enthalpy per unit volume leaves the control volume in the x direction"

It is really the microscopic energy balance expressed in terms of enthalpy, assuming negligible viscous heating. See Eqn. S in Table 11.4-1 of Transport Phenomena by Bird, et al.

As we said before, it represents the net advection of enthalpy per unit volume out of the control volume.

Like you said, it is a catch all for the net of all these effects. It's probably not going to be worthwhile trying to quantify it to a more detailed extent, and, in my judgment, it will be preferable to treat the dispersion length l as an adjustable parameter that the experimental results will be used to calibrate.

No. The heat transfer coefficient will depend on the fluid parameters, velocities, and bed geometry, so, in reality, it can be a function of location and time. I too will be difficult to quantify, and I recommend, at least initially, treating it as a constant and, like the dispersibity, calibrating it to the experimental results. (This is how we treated it in the lumped model) Later, if necessary, we can "sharpen our pencils."

All understood. Ok so we're leaving the dispersion length and the heat transfer coefficient as calibration parameters

No. The solid energy balance is still used. We have just had the add the mass balance because the density is a function of x and t (which needs to be accounted for, since the density varies by a factor of 1000 from liquid to vapor).

Understood yes this is what I was trying to communicate also

I follow all of that. No further questions on any of the above

I do have two short questions on the assumptions we will make here:
1) Are we assuming here that saturation occurs at a specific temperature (rather than find the temperature with bubble/dew point correlations)
2) Should I prepare to start using correlations for any parameters previously kept constant (heat capacity, enthalpy, etc), or will we be using the equations used previously

 

[Chestermiller]

All understood. Ok so we're leaving the dispersion length and the heat transfer coefficient as calibration parameters

Understood yes this is what I was trying to communicate also

I follow all of that. No further questions on any of the above

I do have two short questions on the assumptions we will make here:
1) Are we assuming here that saturation occurs at a specific temperature (rather than find the temperature with bubble/dew point correlations)
2) Should I prepare to start using correlations for any parameters previously kept constant (heat capacity, enthalpy, etc), or will we be using the equations used previously

A foundation of success in doing modeling is to start with simple models and gradually build to more advanced versions incorporating additional features of scomplexity until one arrives at a final version that satisfies the modeler (and matches whatever experimental data has been developed to characterize the system and test the model). So, typically, the data development and the model development are done in tandem. So you should have in mind going to more advanced version of the model as needed.

Depending on your judgment, it may be necessary to employ more accurate correlations of physical property parameters as the model evolves to future versions.

Understand that it is a very bad idea to try to include everything in the model all at once. The concept is that, if you can't solve the more simple versions, you'll never be able to solve the most complex version. Moreover, for each version along the way, you will already have model results to show for your effort (which can be compared with the previous results to quantify the effect of the added complexity and to compare with experimental data).

Right now, as initially specified by you, the model assumes a pure component phase change. But there is no reason you can't start doing separate calculations by hand or by computer to evaluate the VLE behavior at a specified pressure of the effect of temperature in the 2-phase region. Start with one mole of an air mixture, and, using Raout's law, as a function of temperature, calculate the molar split of vapor and liquid, the specific volume of the vapor, the specific volume of the liquid, the average specific volume of the combination of vapor and liquid, the specific enthalpy of the liquid, the specific enthalpy of the vapor, the average specific enthalpy of the combination of vapor and liquid. Get an idea of what this all looks like. Maybe it can be incorporated pretty easily into the model.

One thing we haven't discussed yet is the possibility that the liquid will be traveling through the bed at a different velocity than the vapor (because of viscous drag on the liquid). You should also be thinking about the possibility of pressure variations due to drag, and how that might come into play for a two phase mixture in a packed bed (as well as for the single phase regions). (You have chosen a very complicated problem to work on). All these things should be in the back of your mind, and you should be starting to think about whether they need to be included in later versions of your model (or, at least, how to estimate whether to include them).

 

[Chestermiller]

RECOMMENDED FINITE DIFFERENCE SCHEME

The finite difference scheme I am recommending employs a finite difference grid with spacing $\Delta x$, grid cell boundaries at $(x=0,\ \Delta x,\ 2\Delta x...)$, and grid cell centers at $(x=\Delta x/2,\ 3\Delta x/2,\ 5\Delta x/2,\ ...$). The mass flux into the bed $\phi_0$ is known at the cell boundary x = 0, and we calculate the mass flux at all the other cell boundaries. On the other hand, we calculate the fluid temperatures, enthalpies, and densities at the cell centers.

Based on this description, and employing the method of lines, the 2nd order accurate finite difference representations of the mass balance and heat balance equations for the fluid are expressed as:
$$\frac{d\rho_x}{dt}=\frac{(\phi_{x-\Delta x/2}-\phi_{x+\Delta x/2})}{\Delta x}$$and$$\frac{d(\rho_xh_x)}{dt}=\rho_x\frac{dh_x}{dt}+h_x\frac{d\rho_x}{dt}$$$$=\frac{\phi_{x-\Delta x/2}\left(\frac{h_{x-\Delta x}+h_x}{2}\right)-\phi_{x+\Delta x/2}\left(\frac{h_{x+\Delta x}+h_x}{2}\right)}{\Delta x}$$$$+l\frac{\phi_{x+\Delta x/2}(h_{x+\Delta x}-h_x)-\phi_{x-\Delta x/2}(h_x-h_{x-\Delta x})}{(\Delta x)^2}$$$$+\frac{Ua(T_{S,x}-T_x)}{\epsilon}$$
OK so far?

Are you familiar with the finite differencing concepts of upwind differencing and numerical dispersion associated with modeling advection in transport models?

 

[casualguitar]

A foundation of success in doing modeling is to start with simple models and gradually build to more advanced versions incorporating additional features of complexity until one arrives at a final version that satisfies the modeler

Might have to add this as my laptop wallpaper

Understand that it is a very bad idea to try to include everything in the model all at once.

Got it. From my own limited experience I don't think I've never actually been able to implement the full version of any model first try. I always end up fighting the urge start simply though, as it feels unintuitive and as though it would be slower overall than tackling the entire thing first go I suppose. I'll remember to avoid this approach from now on, and start simply

But there is no reason you can't start doing separate calculations by hand or by computer to evaluate the VLE behavior at a specified pressure of the effect of temperature in the 2-phase region. Start with one mole of an air mixture, and, using Raout's law, as a function of temperature, calculate the molar split of vapor and liquid, the specific volume of the vapor, the specific volume of the liquid, the average specific volume of the combination of vapor and liquid, the specific enthalpy of the liquid, the specific enthalpy of the vapor, the average specific enthalpy of the combination of vapor and liquid. Get an idea of what this all looks like. Maybe it can be incorporated pretty easily into the model.

Understood. I've actually got some dew/bubble point functionality coded up, and also flash functionality (TP,PH, etc) for pure fluids and mixtures, so really in including two components instead of 1, a boiling range instead of a single temperature etc, the difficult bit for me is conceptualising the model. The functionality in code is ready to be used

One thing we haven't discussed yet is the possibility that the liquid will be traveling through the bed at a different velocity than the vapor (because of viscous drag on the liquid). You should also be thinking about the possibility of pressure variations due to drag, and how that might come into play for a two phase mixture in a packed bed (as well as for the single phase regions). (You have chosen a very complicated problem to work on)

Got it. I suppose with the approach of gradually building in complexity, this is something that doesn't have to be included now but can be looked at later (once other things, like the boiling range, have been added)

Actually one note I want to add is that I haven't yet described the existing functionality I've got available in code. Effectively any property/basic functionality like property calculations for pure fluids/mixtures, temperature and/or pressure dependent properties, flash for mixtures to calculate liquid fractions, or the composition of the gas/liquid at given conditions, etc is all possible. So as I mentioned above the most difficult part is likely the conceptualisation of the model. I mention this mainly because I don't want to limit the model capability (where possible) because of lack of access to property/process functionality like the above