Modelling of two phase flow in packed bed using conservation equations

[casualguitar]

The objective of my plot was to compare the two correlations on a common basis. I don't care about anything else.

Eqn. 29 is for use inside the bed/tank model. There is a different U value for each tank in the model. The U value for each tank should be calculated within the bed/tank model. Given P, T, and $G=\dot{m}/A$ for a given tank, you calculate Pr and Re for the tank, then calculate U for that tank. This value of U also varies with time. Use the value of U at the beginning of the time step (end of the previous time step.)

I fully follow this now and can model it. I can get a sense of how U will vary with temperature (at constant pressure) by doing something similar to what I did above except with the new correlations. U doesn't seem to vary all that much with temperature but I will redo these to be sure.

What is the purpose of the Re vs Nu plots? What information can you get from these? I suppose it tells you how convection dominated the heat transfer is, for a given flowrate? Meaning that for our model we would want Nu to be as high as possible, to avoid axial conduction along the bed?

I'll get on scripting U (outside the main model) to see what the approximate range of U values is

 

[Chestermiller]

I fully follow this now and can model it. I can get a sense of how U will vary with temperature (at constant pressure) by doing something similar to what I did above except with the new correlations. U doesn't seem to vary all that much with temperature but I will redo these to be sure.

It will vary very much with $\dot{m}_{j-1}$

What is the purpose of the Re vs Nu plots? What information can you get from these? I suppose it tells you how convection dominated the heat transfer is, for a given flowrate? Meaning that for our model we would want Nu to be as high as possible, to avoid axial conduction along the bed?

My only purpose was to see the comparison between the two correlations.

I'll get on scripting U (outside the main model) to see what the approximate range of U values is

 

[casualguitar]

It will vary very much with mj-1

Ok yes I was thinking about that. Given the range of mdot and T values it seems it might be useful to calculate U at some specific points like (max mdot, max T), or (min mdot, min T) etc to see how much it varies

There was a storm here today so we had no power. Back now so I'll get to implementing this

 

[casualguitar]

Ok yes I was thinking about that. Given the range of mdot and T values it seems it might be useful to calculate U at some specific points like (max mdot, max T), or (min mdot, min T) etc to see how much it varies

There was a storm here today so we had no power. Back now so I'll get to implementing this

Hi Chet, have I got both the intended computational flow, and correlations to be used, correct here? The only change from he previous version was to include the new Nu correlation from BSL rather than from Dixon et al.

You mentioned your version gives Nu directly, and the previous version I did does not do that. I don't see where that comes into play here. Is this correct?

 

[casualguitar]

If the above is correct, then for the approximate range of flow rates the heat transfer coefficient varies from approximately $50$ to $200\frac{W}{m2.K}$. I'm just providing that range to give an idea of what range of values U will take.

The U calculation is straightforward to implement in the main script. I'll make a function that calculates U(T, P, m_dot) tomorrow morning, and will add this to the main model.

After this, there are remaining items such as using a better approximation of pressure drop (using relative permeability). Do you see any critical remaining items here, or is this model mostly ready to be fitted to experimental data?

If it is largely ready, I would very much like to add some other points of discussion to this discussion regarding adding CO2-air mixtures to this model (the idea being to freeze out the CO2 during air liquefaction). I have already checked a number of mixture parameters (enthalpy, density, etc) against literature data for CO2-air mixtures.

 

[Chestermiller]

Hi Chet, have I got both the intended computational flow, and correlations to be used, correct here? The only change from he previous version was to include the new Nu correlation from BSL rather than from Dixon et al.

View attachment 297395

You mentioned your version gives Nu directly, and the previous version I did does not do that. I don't see where that comes into play here. Is this correct?

The BSL definitions of Nu and Re are different from these. So, if you use these equations, you will get the wrong answer for the convective heat transfer coefficient.

 

[casualguitar]

This is the approach using the BSL Re correlation (eq.2, or 14.5-4 in BSL). The Nu definition I used is the same though. This approach does allow you to get. Nu 'directly' as it cancels out as you mentioned. Have I used the correct Nu correlation here?

 

[Chestermiller]

This is the approach using the BSL Re correlation (eq.2, or 14.5-4 in BSL). The Nu definition I used is the same though. This approach does allow you to get. Nu 'directly' as it cancels out as you mentioned. Have I used the correct Nu correlation here?
View attachment 297457

This is the approach using the BSL Re correlation (eq.2, or 14.5-4 in BSL). The Nu definition I used is the same though. This approach does allow you to get. Nu 'directly' as it cancels out as you mentioned. Have I used the correct Nu correlation here?
View attachment 297457

This version is correct, but different from the incorrect version you provided in post #389.

 

[casualguitar]

This version is correct, but different from the incorrect version you provided in post #389.

Got it I see the difference. Great. I've implemented this in a separate script so I will add this U(T,P,m_dot) function into the main script now

You mentioned previously we should include the internal resistance of the solid particle in the lumped parameter model. Is this the bottom right term in eq.8 above?

 

[Chestermiller]

Got it I see the difference. Great. I've implemented this in a separate script so I will add this U(T,P,m_dot) function into the main script now

You mentioned previously we should include the internal resistance of the solid particle in the lumped parameter model. Is this the bottom right term in eq.8 above?

Sure

 

[casualguitar]

Sure

Ok this function is now integrated into the main model. I am getting some unusual output (not because of the new U function, I was getting it before that). My guess its a units thing so I'm going to do a full debug to see what is going on

 

[casualguitar]

Sure

Hi Chet, I've done a 1st debug of the script and these are the output plots, which are definitely incorrect, however I'm sharing them anyway to show where the model is at. The S shaped curves are not as 'nice' as I expected. I think the source of the error is that I am mixing up the ICs/BCs for m and m_dot somewhere, however that's just a guess for now

 

[casualguitar]

Hi Chet, apologies I'm actually still debugging. I think the bug is in the mass holdup and mass flow rate usage. I'm almost sure I'm not setting up the boundary and initial conditions correctly for these. The reason I think this is because 1) after thinking about it I'm actually not sure if I should use m or m_dot an an initial condition and 2) The plots above show n=1 as being shaped differently to the other positions

Its late here so I'll reframe this as both a question and a clear summary of how I am using m and mdot in the morning

 

[Chestermiller]

Hi Chet, apologies I'm actually still debugging. I think the bug is in the mass holdup and mass flow rate usage. I'm almost sure I'm not setting up the boundary and initial conditions correctly for these. The reason I think this is because 1) after thinking about it I'm actually not sure if I should use m or m_dot an an initial condition and 2) The plots above show n=1 as being shaped differently to the other positions

Its late here so I'll reframe this as both a question and a clear summary of how I am using m and mdot in the morning

By the initial condition, you mean time zero, right? Why do you need this mdot?

 

[casualguitar]

By the initial condition, you mean time zero, right? Why do you need this mdot?

The initial condition array is y0 = [h0,Ts0,m0, h1,Ts1,m1,...hn,Tsn,mn]
The solution array however is y = [h0,Ts0,mdot0, h1,Ts1,mdot1,...hn,Tsn,mdotn]

So I am supplying mass holdup as one of the initial conditions (mdot is only supplied as a boundary condition i.e. the inlet flow), and the solution array has mdot, not mass holdup

This is the code I am currently using:

Its a lengthy script, however it is just the my_system function where the algebraic solutions for each derivative at each position are set up that are important, and the few lines after this where the initial conditions are set up

I could put mass holdup in the solution array also, and then solve for mdot in the after processing stage. I don't know if this would affect the output though. I'll try this

 

[Chestermiller]

The initial condition array is y0 = [h0,Ts0,m0, h1,Ts1,m1,...hn,Tsn,mn]
The solution array however is y = [h0,Ts0,mdot0, h1,Ts1,mdot1,...hn,Tsn,mdotn]

So I am supplying mass holdup as one of the initial conditions (mdot is only supplied as a boundary condition i.e. the inlet flow), and the solution array has mdot, not mass holdup

This is the code I am currently using:

Its a lengthy script, however it is just the my_system function where the algebraic solutions for each derivative at each position are set up that are important, and the few lines after this where the initial conditions are set up

I could put mass holdup in the solution array also, and then solve for mdot in the after processing stage. I don't know if this would affect the output though. I'll try this

What parameters are being integrated with respect to time (I.e.,your dependent variables)?

 

[casualguitar]

What parameters are being integrated with respect to time (I.e.,your dependent variables)?

The parameters in the dydt list (in the script) are being integrated w.r.t. time and these are fluid enthalpy, solid temperature, and m_dot (mass flow rate out), which are equations 2, 3, and 4a in the image below (your original equations).

Now that you mention it, eq.2. and 3 have time derivatives, but 4a does not. We can solve directly for dm/dt and m_dot, so would I just need to integrate eq.2. and 3 here? I'm actually not sure why I assumed I would need to integrate 3 equations

 

[Chestermiller]

The parameters in the dydt list (in the script) are being integrated w.r.t. time and these are fluid enthalpy, solid temperature, and m_dot (mass flow rate out), which are equations 2, 3, and 4a in the image below (your original equations).

Now that you mention it, eq.2. and 3 have time derivatives, but 4a does not. We can solve directly for dm/dt and m_dot, so would I just need to integrate eq.2. and 3 here? I'm actually not sure why I assumed I would need to integrate 3 equations
View attachment 297588

Yes, just two. The other equation is spatial only.

 

[casualguitar]

Yes, just two. The other equation is spatial only.

Hi Chet, I have removed dm/dt from the dependent variables (required some small code adjustments to the after processing etc also). The simulation takes close to 10 minutes to run 1000s of simulation time, so this evening I am looking into optimising this code (replacing the nested for loop with a vectorised alternative, etc).

These are the plots produced by the simulation. They are not the 'perfect' S shaped curves I was thinking they might be. Do these curves look reasonable, or do they suggest there are further simulation errors? If the latter, would you have a suggestion as for how to systematically approach debugging?

Fluid enthalpy for 3 tanks:

Solid Temperature:

Mass holdup for all tanks (note mole unit):

Fluid temperature (not reason for red plot here is I used two different approaches to calculate this, just checking that they are equivalent):

The main thing I am concerned about in these plots, is that the plot shape of the first tank is not the same as the other tanks, suggesting that I am possibly doing something wrong in how I am setting up the initial/boundary conditions

 

[Chestermiller]

Hi Chet, I have removed dm/dt from the dependent variables (required some small code adjustments to the after processing etc also). The simulation takes close to 10 minutes to run 1000s of simulation time, so this evening I am looking into optimising this code (replacing the nested for loop with a vectorised alternative, etc).

These are the plots produced by the simulation. They are not the 'perfect' S shaped curves I was thinking they might be. Do these curves look reasonable, or do they suggest there are further simulation errors? If the latter, would you have a suggestion as for how to systematically approach debugging?

Fluid enthalpy for 3 tanks:
View attachment 297676
Solid Temperature:
View attachment 297677
Mass holdup for all tanks (note mole unit):
View attachment 297678
Fluid temperature (not reason for red plot here is I used two different approaches to calculate this, just checking that they are equivalent):
View attachment 297679
The main thing I am concerned about in these plots, is that the plot shape of the first tank is not the same as the other tanks, suggesting that I am possibly doing something wrong in how I am setting up the initial/boundary conditions

Check the 1st tank response against the 1 tank model with 1/n the volume.

 

[casualguitar]

Check the 1st tank response against the 1 tank model with 1/n the volume.

Yes the plot for the 1 tank model is in the same ballpark as the full model (semi surprised that they line up this well). This is using an averaged value for the T dependent parameters like heat capacity. The constant temperature zone is more pronounced in the main model, however at all times the temperature discrepancy between the two models is fairly reasonable, which suggests that the main model does not set up the boundary/initial conditions incorrectly, suggesting that the model may be bug-free?

There does remain the issue of the simulation taking a very long time to run. I am currently timing it but it is definitely 10+ minutes (EDIT: it took 15 mins to run to completion). It does look like there is something going on in the saturation zone where thermo takes a while to converge. If this can be fixed then simulation time will come down by a lot because the gas and liquid zones are relatively quick to simulate

 

[casualguitar]

Hi Chet, some short updates -

It is the thermo functions that are computationally expensive. I am looking into how these can be made less expensive to use for now (mostly to speed up model development)

We previously discussed tuning the model to experimental data. I have checked with the lab and the experimental data will be in excel format (straightforward to deal with from a programming point of view). Python has some nice curve fitting libraries that take x-y data, the model, and the tuning parameter(s), and return 'optimised' parameters using say least squares for example.

We discussed tuning parameters earlier on in the model (screenshot below is from post #99). This dispersivity parameter no longer seems to be accessible as it is not on the 'top layer' of the model anymore. Should this still be used as the tuning parameter, or are we now using another parameter that is at the top level of the model?

And lastly, I have no model functionality left on my list that seems to be very important to implement (pressure gradient using relative permeability, internal particle temperature gradients). Do you see any functionality that should be implemented, before moving to the next section of this model (adding CO2 to the air mixture)?
Edit: Before moving forward, I will do a write up of everything we have discussed on the final version of the model tomorrow (late here now)

 

[casualguitar]

Hi Chet, short write up of the discussion we've had to date (on the most complex version of this model only so far) :https://docs.google.com/document/d/...ouid=111822275555236809516&rtpof=true&sd=true

The link is to a word doc of a summarised version of this discussion, plus output from the model and the actual Python script I made. I'll leave this link active indefinitely in case this is useful to any future passers by on this forum. I will also be updating it regularly until the packed bed modelling is finished so I'll probably change the location to an open google drive file

If it is ok with you I would like to discuss adding CO2 to the existing model. The idea here is that an ambient temperature air-CO2 mixture is passed through a cold packed bed, liquefying the air and freezing the CO2. The solid CO2 will move through the bed as a 'plug' which accumulates in mass as it passes through the bed, and will be extracted from the line at the end of the packed bed therefore 'cleaning' the air stream of CO2. If this is ok with you, I would like to attempt to lay out the mass and energy balances for this model (similar to what you did initially for the air model). I will definitely have questions on the mass balance for CO2, but I think it would be useful anyway to make an attempt so we have an initial discussion point!

 

[casualguitar]

Hi Chet, just posting the updated plots here. These plots include the heat transfer coefficient correlation and the $d\rho/dH$ correlation, and now runs a lot faster than it previously did (about 5 times faster):

Currently working towards allowing non-constant enthalpy profiles to be supplied to the model

 

[casualguitar]

Hi Chet, updating here as discussed. So here is the position vs gas temperature plot for a molar flow of 0.001mol/s, and then below for 0.008mol/s.

The lower molar flow rate allows for a sharper temperature front that I think you mentioned was desirable

0.001mol/s:

0.008mol/s:

Note the 0.001mol/s flow also had n=10, I just zoomed in on the relevant section. But yes as you can see the temperature front is much more pronounced for the lower flow

 

[Chestermiller]

Hi Chet, updating here as discussed. So here is the position vs gas temperature plot for a molar flow of 0.001mol/s, and then below for 0.008mol/s.

The lower molar flow rate allows for a sharper temperature front that I think you mentioned was desirable

0.001mol/s:

View attachment 314228

0.008mol/s:

View attachment 314229

Note the 0.001mol/s flow also had n=10, I just zoomed in on the relevant section. But yes as you can see the temperature front is much more pronounced for the lower flow

This is a very difficult comparison because the advance of the fronts is different because of the different flow rates (IT'S NOT ON A COMMON BASIS). Compare the cases on a common basis in which the injected volumes are equal,, and see how they compare. So run the 0.001 for 8x the times of the 0.008.

 

[casualguitar]

Compare the cases on a common basis in which the injected volumes are equal

Ah yes, I didn't think of that (not a common basis)

Running for m = 0.001mol/s and t = 400s, and also m = 0.008mol/s and t=50s gives these plots which match almost exactly:

In addition, these are the final arrays, in the format [h0,Ts0,h1,Ts1,...hn,Tsn] where h is the fluid enthalpy and Ts is the solid temperature. They are almost identical:

So returning to the initial question about which parameters most affect axial dispersion, it seems we can conclude from the above that molar flowrate does not affect axial dispersion significantly?

 

[casualguitar]

Hey Chet, just in addition, I've got a working 'stopping condition' now. Meaning that the simulation stops once the outlet temperature of the fluid reaches the inlet temperature of the fluid

The reason this is useful is because I can now sort of quantify the beds performance degradation over time (unless I'm wrong) in this way:
1) assume the bed is initially at ambient temperature
2) pass cold air through the bed at 0.008 mol/s until the bed reaches steady state
3) save the time taken to reach this stopping condition
4) pass warm air through the bed at 0.008 mol/s (liquefying it), until the stopping condition is met, meaning that liquefaction can no longer occur
5) check the difference in time between the hot and cold 'passes', and multiply this by the molar flowrate. Then we can divide this value by the total amount of air liquefied and multiply by 100 to quantify the percentage of 'liquefying power' lost per pass

Does this approach sound reasonable?

 

[Chestermiller]

Ah yes, I didn't think of that (not a common basis)

Running for m = 0.001mol/s and t = 400s, and also m = 0.008mol/s and t=50s gives these plots which match almost exactly:
View attachment 314374
View attachment 314375

In addition, these are the final arrays, in the format [h0,Ts0,h1,Ts1,...hn,Tsn] where h is the fluid enthalpy and Ts is the solid temperature. They are almost identical:
View attachment 314372

So returning to the initial question about which parameters most affect axial dispersion, it seems we can conclude from the above that molar flowrate does not affect axial dispersion significantly?

It affects dispersion to about the same extent for a given total volume of flow.

 

[Chestermiller]

Hey Chet, just in addition, I've got a working 'stopping condition' now. Meaning that the simulation stops once the outlet temperature of the fluid reaches the inlet temperature of the fluid

The reason this is useful is because I can now sort of quantify the beds performance degradation over time (unless I'm wrong) in this way:
1) assume the bed is initially at ambient temperature
2) pass cold air through the bed at 0.008 mol/s until the bed reaches steady state
3) save the time taken to reach this stopping condition
4) pass warm air through the bed at 0.008 mol/s (liquefying it), until the stopping condition is met, meaning that liquefaction can no longer occur
5) check the difference in time between the hot and cold 'passes', and multiply this by the molar flowrate. Then we can divide this value by the total amount of air liquefied and multiply by 100 to quantify the percentage of 'liquefying power' lost per pass

Does this approach sound reasonable?

I'm not able to follow this.

 

[casualguitar]

I'm not able to follow this.

Yes I explained it badly really

So we are modelling the charging (cooling) and discharging (heating) of the bed in a cyclical operation i.e. cold air is passed through to cool the bed, then hot air is passed through with the goal of liquefying it. Liquid air being the product we want here. Then the produced liquid air is passed back through the bed to retool it. The idea is to quantify how the liquefaction performance of the bed degrades over the charge/discharge cycles due to losses etc. i.e. does the amount of liquid we get after each charge/discharge cycle stay the same, or does it gradually decrease?

The question I'm asking is - how can we quantify the degradation in liquefaction performance of the packed bed over time?

One thought I had was to simply check how much cold air is required to charge the bed, then check how much liquid air this cold bed can produce, then pass this produced liquid air back through the bed to cool the bed, then pass ambient air though to liquefy it, and repeat this cycle. We might see that each pass produces slightly less liquid air, due to losses in the system

Is there any aspect of this explanation that is unclear or can I continue with how I was thinking of calculating these losses?

 

[casualguitar]

I'm not able to follow this.

Hi Chet,

I've got the cooling/heating of the packed bed fully working now. And here is a plot of a number of charge/discharge (cooling/heating) cycles:

The simulation details:
Bed initial temperature = 80K
Fluid inlet temperature = 293K

The discharge (bed heating) is run until the last position in the bed (position 4 in this case) heats up to a temperature of 1K above 80K i.e. 81K. This signals that the bed can no longer cool the gas to 80K and liquefy it.

Then, from the outlet (position 4) liquid gas at 80K is passed through for the same amount of time, and the bed recools. It won't recoil all the way to 80K though because of some losses.

Then again hot air is passed through until the last position increases by 1K etc this repeats for a few cycles and the final temperature profiles after each cooling phase are plotted above. If this cycle was repeated a large number of times then eventually we would no longer be able to liquefy air (and the bed would undergo a full cooling phase to fully reset the temperature to 80K)

My question is how can I quantify how the performance of the bed degrades over time?

I was thinking of using the time required to increase the last position by 1K. If I store this value for each cycle, we should see that this time decreases (as the bed performance degrades). We can do $moles-of-liquid-lost = \Delta t * \dot{m}$ to calculate the amount of moles we lose per cycle, compared to the theoretical max. Does this sound reasonable?

 

[Chestermiller]

Yes I explained it badly really

So we are modelling the charging (cooling) and discharging (heating) of the bed in a cyclical operation i.e. cold air is passed through to cool the bed, then hot air is passed through with the goal of liquefying it. Liquid air being the product we want here. Then the produced liquid air is passed back through the bed to retool it. The idea is to quantify how the liquefaction performance of the bed degrades over the charge/discharge cycles due to losses etc. i.e. does the amount of liquid we get after each charge/discharge cycle stay the same, or does it gradually decrease?

The question I'm asking is - how can we quantify the degradation in liquefaction performance of the packed bed over time?

One thought I had was to simply check how much cold air is required to charge the bed, then check how much liquid air this cold bed can produce, then pass this produced liquid air back through the bed to cool the bed, then pass ambient air though to liquefy it, and repeat this cycle. We might see that each pass produces slightly less liquid air, due to losses in the system

Is there any aspect of this explanation that is unclear or can I continue with how I was thinking of calculating these losses?

I have no idea what you are driving at here. Please show graphs of what you have in mind both for the heating and cooling operations.

 

[Chestermiller]

Hi Chet,

I've got the cooling/heating of the packed bed fully working now. And here is a plot of a number of charge/discharge (cooling/heating) cycles:
View attachment 314754
The simulation details:
Bed initial temperature = 80K
Fluid inlet temperature = 293K

The discharge (bed heating) is run until the last position in the bed (position 4 in this case) heats up to a temperature of 1K above 80K i.e. 81K. This signals that the bed can no longer cool the gas to 80K and liquefy it.

Then, from the outlet (position 4) liquid gas at 80K is passed through for the same amount of time, and the bed recools. It won't recoil all the way to 80K though because of some losses.

Then again hot air is passed through until the last position increases by 1K etc this repeats for a few cycles and the final temperature profiles after each cooling phase are plotted above. If this cycle was repeated a large number of times then eventually we would no longer be able to liquefy air (and the bed would undergo a full cooling phase to fully reset the temperature to 80K)

My question is how can I quantify how the performance of the bed degrades over time?

I was thinking of using the time required to increase the last position by 1K. If I store this value for each cycle, we should see that this time decreases (as the bed performance degrades). We can do $moles-of-liquid-lost = \Delta t * \dot{m}$ to calculate the amount of moles we lose per cycle, compared to the theoretical max. Does this sound reasonable?

I am really confused about what is going on here. Please show graphs for the heating and then the subsequent cooling.

 

[casualguitar]

I am really confused about what is going on here. Please show graphs for the heating and then the subsequent cooling.

Hi Chet, my apologies for the delay. I was running a few simulations for various values of n. Anything greater than n=10 takes considerable time now. n=30 took 24 hours actually.

Anyway, I wanted to do a 'grid independence' study of sorts. If you run hot air through the bed for x amount of time, we get a thermocline in the bed. I did this for n=5,10,20, and 30 and here are the results. Clearly increasing n results in a 'tighter' thermocline or a sharper drop from high temperature to low:
n=5:

n=10:

n=20:

n=30:

I'm running higher n now (n=100) however this will take a few days. The goal of this analysis is to 1) find the value of n where the output plots stop changing i.e. grid independence and 2) to quantify the losses in the bed across a charge/discharge cycle i.e. if we discharge for 100s and then charge for 100s, how close to we get to the original temperature profile?

No questions on this, just updating