Modelling of two phase flow in packed bed (continued)

[casualguitar]

We are looking at the overall mass balance of CO2 in the column, not the instantaneous amounts entering and leaving.

Yes I was also looking at the overall mass balance. However from the plot earlier (where the total CO2 out was greater than the total CO2 in) it seems that it isn't working. Therefore I moved to the instantaneous balance to search for bugs. There is definitely a bug in the instantaneous mass balance and this may be causing the bug in the overall mass balance

In the post above I was looking for the bug in the $\dot{m}$ equation (i.e. why is the outlet flow greater than the inlet flow) in the hope that fixing this would fix the overall mass balance bug. Is this a reasonable approach?

 

[Chestermiller]

Yes I was also looking at the overall mass balance. However from the plot earlier (where the total CO2 out was greater than the total CO2 in) it seems that it isn't working. Therefore I moved to the instantaneous balance to search for bugs. There is definitely a bug in the instantaneous mass balance and this may be causing the bug in the overall mass balance

In the post above I was looking for the bug in the $\dot{m}$ equation (i.e. why is the outlet flow greater than the inlet flow) in the hope that fixing this would fix the overall mass balance bug. Is this a reasonable approach?

Don't forget to take into account the number of moles of N2 that was in the column to begin with.

 

[casualguitar]

Don't forget to take into account the number of moles of N2 that was in the column to begin with.

Ah so to account for initial N2 I've calculated the molar holdup of gas in the column initially ($m_j$), is there somewhere else I should account for the initial N2?

What I also did this morning for this model was switch to average $C_p$ in the gas phase heat balance rather than just using $C_{p,N2}$

 

[Chestermiller]

Ah so to account for initial N2 I've calculated the molar holdup of gas in the column initially ($m_j$), is there somewhere else I should account for the initial N2?

This only has to be done if you are looking at the cumulative in- and out flows, not the instantaneous values.

What I also did this morning for this model was switch to average $C_p$ in the gas phase heat balance rather than just using $C_{p,N2}$

This should not have a major effect.

 

[casualguitar]

This only has to be done if you are looking at the cumulative in- and out flows, not the instantaneous values.

This should not have a major effect.

I think I've narrowed it down to an error with the MDR/MLR functions (molar desublimation/liquefaction rates) i.e. this equation: $$M_i'' = \frac{k_i(Py_i - P(T_I))}{RT}$$

Reason being is the follows - when I set the mass transfer rate equal to zero ($k_i$ = 0 for CO2 and h2o), everything works:
- the mole balance works (see below moles in > moles out):

The CO2 mole fraction equation works (see below it levels out at 0.2, whereas with mass transfer turned on it levels out at 0.35 which is impossible since the inlet mole fraction is 0.2):

Lastly the temperature profiles work (see below the temperature levelling out at the inlet gas temperature, rather than some temperature below this):

No other update besides this. I'll investigate the MDR/MLR functions and see what's going on

 

[casualguitar]

No other update besides this. I'll investigate the MDR/MLR functions and see what's going on

In regards to calculating the sublimation pressure - I've got a correlation from Perrys working currently. One question on this - what happens to the sublimation pressure above the triple point? Currently I have something like this:

if temp < triple point temp:
-> use the correlation to calculate sublimation pressure
else:
-> return a very large number for the sublimation pressure

Do we have a sublimation rate associated with temperatures above the triple point temperature or is it 'infinitely fast'?

 

[Chestermiller]

I think I've narrowed it down to an error with the MDR/MLR functions (molar desublimation/liquefaction rates) i.e. this equation: $$M_i'' = \frac{k_i(Py_i - P(T_I))}{RT}$$

Reason being is the follows - when I set the mass transfer rate equal to zero ($k_i$ = 0 for CO2 and h2o), everything works:
- the mole balance works (see below moles in > moles out):
View attachment 312602

At long times, why doesn't the CO2 gas holdup plus the total CO2 out equal the total CO2 in?

 

[Chestermiller]

In regards to calculating the sublimation pressure - I've got a correlation from Perrys working currently. One question on this - what happens to the sublimation pressure above the triple point? Currently I have something like this:

if temp < triple point temp:
-> use the correlation to calculate sublimation pressure
else:
-> return a very large number for the sublimation pressure

Do we have a sublimation rate associated with temperatures above the triple point temperature or is it 'infinitely fast'?

Above the triple point, you use the equilibrium vapor pressure of the liquid. This is very close to the equilibrium vapor pressure for sublimation, so you can use that too. We are not allowing for the CO2 deposited to become liquid.

 

[casualguitar]

Above the triple point, you use the equilibrium vapor pressure of the liquid. This is very close to the equilibrium vapor pressure for sublimation, so you can use that too. We are not allowing for the CO2 deposited to become liquid.

The vapour pressure curve for CO2 is included now. If the mass transfer is turned way down (k_i = 8x10^-10), then the plots look fine
The temperature levels out at the expected temperature:

The molar desublimation rate of CO2 looks fine:

The rise and fall of the solid CO2 deposited on the bed looks good (the sublimation is rapid and not like previous plots where the sublimation was much slower than the desublimation):

Lastly, the mole fraction of CO2 levels out at the correct value of 0.2:

However when the mass transfer is turned up (k_i = 8*10^-4) these trends no longer work:
Temperature levels out at a low value:

The sublimation is much slower than the desublimation:

And lastly the CO2 mole fraction levels out at an impossible value:

It seems odd to me that there is a 'cut off' k_i value where the trends are no longer right, because k_i doesn't change the sign of MDR/MLR (molar desublimation/liquefaction rates), just the magnitude.

There is a fault in my MDR/MLR logic which potentially causes large molar flow rates at the inlet. Close to t=0, the y_CO2 value is so low that the MDR is actually negative. Could I add in some correction to say that if the y_CO2 value is extremely low then assume MDR = 0, or similar?

 

[casualguitar]

So initially, the MDR value is very slightly below zero (doesn't look like it below but it is if I zoom in). This negative value actually causes $\dot{m}_j$ to be much higher than $\dot{m}_{j-1}$, and this molar flow then propagates through the bed

 

[Chestermiller]

I think I've narrowed it down to an error with the MDR/MLR functions (molar desublimation/liquefaction rates) i.e. this equation: $$M_i'' = \frac{k_i(Py_i - P(T_I))}{RT}$$

Reason being is the follows - when I set the mass transfer rate equal to zero ($k_i$ = 0 for CO2 and h2o), everything works:
- the mole balance works (see below moles in > moles out):
View attachment 312602
The CO2 mole fraction equation works (see below it levels out at 0.2, whereas with mass transfer turned on it levels out at 0.35 which is impossible since the inlet mole fraction is 0.2):

I calculate a final total CO2 gas holdup of 0.0006 moles, which is roughly half the value on this graph.

 

[casualguitar]

I calculate a final total CO2 gas holdup of 0.0006 moles, which is roughly half the value on this graph.

Exactly yes the reason its half the value is because the temperature profile (with mass transfer turned on) is the following (below), where the temperature levels out at a value approximately half of the correct value (shown further below):

However to get the correct profile above I have to turn off mass transfer, meaning that the issue is definitely related to the mass transfer.

Looking at the desublimation rate at short time, I see that it is actually negative initially (indicating sublimation). This 1) doesn't seem to make sense as the bed is cold and 2) causes $\dot{m}_{out}$ > $\dot{m}_{in}$ which isn't possible.

So why is there apparently sublimation occurring at short times? Well if (in the equation below) $Py_i$ is less than $P(T_I)$ then this will occur. The issue is that at short times $y_i$ is so small that this equation says there is sublimation happening not desublimation. Is this correct or should I add a further constraint to stop this?

The equation for reference:
$$M_i'' = \frac{k_i(Py_i - P(T_I))}{RT}$$

 

[Chestermiller]

Exactly yes the reason its half the value is because the temperature profile (with mass transfer turned on) is the following (below), where the temperature levels out at a value approximately half of the correct value (shown further below):
View attachment 312750
View attachment 312749

However to get the correct profile above I have to turn off mass transfer, meaning that the issue is definitely related to the mass transfer.

Looking at the desublimation rate at short time, I see that it is actually negative initially (indicating sublimation). This 1) doesn't seem to make sense as the bed is cold and 2) causes $\dot{m}_{out}$ > $\dot{m}_{in}$ which isn't possible.

So why is there apparently sublimation occurring at short times? Well if (in the equation below) $Py_i$ is less than $P(T_I)$ then this will occur. The issue is that at short times $y_i$ is so small that this equation says there is sublimation happening not desublimation. Is this correct or should I add a further constraint to stop this?

The equation for reference:
$$M_i'' = \frac{k_i(Py_i - P(T_I))}{RT}$$

Let's. try the trick factor that Tuinier et al use, appropriately corrected for our units: $\frac{m_i}{m_i+0.1}$

 

[casualguitar]

Let's. try the trick factor that Tuinier et al use, appropriately corrected for our units: $\frac{m_i}{m_i+0.1}$

So we let desublimation as it is, and add in this trick factor for sublimation?

To convert their units (kg/m3) to ours (they use mol/m3) I guess we would just need to multiply each $m_i$ term by $\frac{1000}{mW}$?

That is interesting though. I had thought that desublimation was dependent on temperature and pressure only. I did not realize that a sufficiently low $y_i$ (and therefore a sufficiently low partial pressure) would mean that sublimation actually occurs

Also a note - I'm using units of pascals for the $M_i''$ equation calculation

 

[casualguitar]

So we let desublimation as it is, and add in this trick factor for sublimation?

To convert their units (kg/m3) to ours (they use mol/m3) I guess we would just need to multiply each $m_i$ term by $\frac{1000}{mW}$?

That is interesting though. I had thought that desublimation was dependent on temperature and pressure only. I did not realize that a sufficiently low $y_i$ (and therefore a sufficiently low partial pressure) would mean that sublimation actually occurs

Also a note - I'm using units of pascals for the $M_i''$ equation calculation

Adding in the fudge factor seems to fix the temperature profile issue, and the other plots look a bit better also. Also the final molar holdup is now 0.006 also:

Edit: I'm not sure what the benefit of the fudge factor is though. What is this for? To ensure that sublimation is slightly more inhibited than desublimation?

 

[Chestermiller]

So we let desublimation as it is, and add in this trick factor for sublimation?

yes

To convert their units (kg/m3) to ours (they use mol/m3) I guess we would just need to multiply each $m_i$ term by $\frac{1000}{mW}$?

I think our units are moles in each tank. We would multiply both numerator and denominator of their factor (including the 0.1 term) by $\frac{(1000)A_C\Delta z}{molecular\ weight}$ to get ours.

That is interesting though. I had thought that desublimation was dependent on temperature and pressure only. I did not realize that a sufficiently low $y_i$ (and therefore a sufficiently low partial pressure) would mean that sublimation actually occurs

It doesn't mean this. All this does is make sure that the mass deposited never gets less than -0.1 kg/m^3 (which is not significant). If the deposited mass tries to go negative, we add mass. Their purpose of this factor was just to provide smoothing of the behavior (not to model actual behavior) so that the deposited mass can't stay negative by more than 0.1.

Also a note - I'm using units of pascals for the $M_i''$ equation calculation

The P/RT is supposed to be moles/m^3 of gas.

 

[Chestermiller]

Adding in the fudge factor seems to fix the temperature profile issue, and the other plots look a bit better also. Also the final molar holdup is now 0.006 also:
View attachment 312767
View attachment 312768View attachment 312769View attachment 312770

Edit: I'm not sure what the benefit of the fudge factor is though. What is this for? To ensure that sublimation is slightly more inhibited than desublimation?

No. To avoid in a smooth manner the deposition going negative by more the 0.1 kg/m^3 (insignificant amount). it's a numerical trick.

 

[casualguitar]

I think our units are moles in each tank. We would multiply both numerator and denominator of their factor (including the 0.1 term) by (1000)ACΔzmolecular weight to get ours.

I guess we would just be multiplying the $m_i$ terms by the above factor to convert from kg/m3 to mol/tank, or else the factor would cancel out?

It doesn't mean this.

That is interesting though. I had thought that desublimation was dependent on temperature and pressure only. I did not realize that a sufficiently low yi (and therefore a sufficiently low partial pressure) would mean that sublimation actually occurs

When I said this what I meant was that this equation $M_i'' = \frac{k_i(Py_i - P(T_I))}{RT}$ seems to say that sublimation is a function of mole fraction (in addition to P and T), meaning that even with a temperature and pressure below the sublimation temperature, sublimation may still occur if the mole fraction of the species, $y_i$ is sufficiently low. I think anyway.

I think our units are moles in each tank. We would multiply both numerator and denominator of their factor (including the 0.1 term) by (1000)ACΔzmolecular weight to get ours.

Anyway yes in regards to the factor above - we're solving for $M_{CO2}$ (mass of deposited CO2) at each space/time increment so can I simply add this factor to the sublimation equation?: $$trick factor = \frac{M_{CO2}}{M_{CO2} + 0.1}$$

If so, then here's the output with that factor included:

The molar holdup levels out at a value of 0.006 which is what you previously calculated

And the positional plots have similar peaks to the Tuinier data (the sublimation is still slower than it is in Tuinier et al though):

 

[Chestermiller]

I guess we would just be multiplying the $m_i$ terms by the above factor to convert from kg/m3 to mol/tank, or else the factor would cancel out?When I said this what I meant was that this equation $M_i'' = \frac{k_i(Py_i - P(T_I))}{RT}$ seems to say that sublimation is a function of mole fraction (in addition to P and T), meaning that even with a temperature and pressure below the sublimation temperature, sublimation may still occur if the mole fraction of the species, $y_i$ is sufficiently low. I think anyway.

No. This only happens to a tiny extent when we use the fudge factor. The idea is to not let the deposition amount become significantly less than 0 as a result of slight numerical overshoot.

Anyway yes in regards to the factor above - we're solving for $M_{CO2}$ (mass of deposited CO2) at each space/time increment so can I simply add this factor to the sublimation equation?: $$trick factor = \frac{M_{CO2}}{M_{CO2} + 0.1}$$

I think you mean moles of deposited CO2 here. In any event, the fudge factor should, as I indicated earlier, be $$trick factor = \frac{M_{CO2}}{M_{CO2} + 0.1\frac{(1000)A_C\Delta z}{(molecular\ weight \ of \ CO2)}}$$A corresponding equation applies to water.

This is what it would take to make your calculations equivalent to theirs, at least with regard to the fudge factor.

If so, then here's the output with that factor included:

View attachment 312826
View attachment 312827View attachment 312828

The molar holdup levels out at a value of 0.006 which is what you previously calculated

And the positional plots have similar peaks to the Tuinier data (the sublimation is still slower than it is in Tuinier et al though):

View attachment 312829
View attachment 312830

You need to work toward (1) using the correct equivalent fudge factors, (2) increasing to 30 tanks, and (3) using our mass transfer parameterization.

 

[casualguitar]

The idea is to not let the deposition amount become significantly less than 0 as a result of slight numerical overshoot.

Got it

I think you mean moles of deposited CO2 here. In any event, the fudge factor should, as I indicated earlier, be trickfactor=MCO2MCO2+0.1(1000)ACΔz(molecular weight of CO2)A corresponding equation applies to water.

This is what it would take to make your calculations equivalent to theirs, at least with regard to the fudge factor.

Ah yes I meant moles, and now I see what you mean. Got it. So you've multiplied all terms by the $\frac{(1000)A_C\Delta z}{molecular\ weight}$ factor. I'll implement this now

So with the correct equivalent fudge factor, and n=30, the output is below

(3) using our mass transfer parameterization.

Are we not already using our mass transfer parameterisation? Or is there something I have left out?

Here is the output:

The CO2 solid moles buildup looks pretty close to the Tuinier plot. I have one or two ideas as for why the sublimation is slower. I'll look into this.

But yes, what did you mean by 'our mass transfer parameterization' if not this equation: $M_i'' = \frac{k_i(Py_i - P(T_I))}{RT}$

 

[Chestermiller]

But yes, what did you mean by 'our mass transfer parameterization' if not this equation: $M_i'' = \frac{k_i(Py_i - P(T_I))}{RT}$

No. The k correlation.

 

[casualguitar]

No. The k correlation.

Uncertain about the units for $k_i$ and $D_{ab}$ (mass diffusion coefficient) -

So the correlation I've currently got for $k_i$ is:
$$k_i = (2.19Re^{1/3} + 0.78Re^{0.619})Sc^{1/3}D_{ab}\rho_m/d_p$$

The the units of $D_{ab}$ are $m^2/s$ from the correlation above. BSL provides these $D_{ab}$ values in $m^2/s$ also. The units of $k_i$ are therefore $mol/m^2.s$ then. There is probably some unit issue here somewhere because the molar deposition rate is also mol/m2.s. Can the units of $k_i$ be the same as the units of $M_i''$?

The Schmidt number I calculated seems to be ok also when compared to BSL data. The plot produced for a range of Reynolds numbers looks like this:

For references the current constant $k_{CO2}$ value I'm using is 8*10^-4, so these seem high. Anyway I'll implement the non-constant $k_j$ functionality into the current script and let's see what happens

 

[Chestermiller]

Uncertain about the units for $k_i$ and $D_{ab}$ (mass diffusion coefficient) -

So the correlation I've currently got for $k_i$ is:
$$k_i = (2.19Re^{1/3} + 0.78Re^{0.619})Sc^{1/3}D_{ab}\rho_m/d_p$$

You keep flip-flopping between k having units of m/s and having units of moles/m^2.sec. If you are using m/s, leave out the $\rho_m\ (=\frac{P}{RT_I})$ in the above equation. Otherwise, keep it. I also think there should be a factor of $(1-\epsilon)$ in the above equation. Your equation for molar flux seems to use k is in m/s. Is that your understanding?

The the units of $D_{ab}$ are $m^2/s$ from the correlation above. BSL provides these $D_{ab}$ values in $m^2/s$ also. The units of $k_i$ are therefore $mol/m^2.s$ then. There is probably some unit issue here somewhere because the molar deposition rate is also mol/m2.s. Can the units of $k_i$ be the same as the units of $M_i''$?

You need to get the units issue resolved.

 

[casualguitar]

You keep flip-flopping between k having units of m/s and having units of moles/m^2.sec. If you are using m/s, leave out the $\rho_m\ (=\frac{P}{RT_I})$ in the above equation. Otherwise, keep it. I also think there should be a factor of $(1-\epsilon)$ in the above equation. Your equation for molar flux seems to use k is in m/s. Is that your understanding?

You need to get the units issue resolved.

That is my understanding yes. I had some issues getting the non constant $k_i$ model to converge however it does converge now. The output below is for n=5. Adding in non constant $k_i$ seems to increase the runtime of the model by quite a bit for some reason. So I'll run n=30 in a bit (currently on the way home) and see how long it takes. Anyway here's the n=5 output:

I don't like how the MCO2 plot looks (the peaks should increase in size over time), so I'll look into that. But yes first I'll run n=30

 

[casualguitar]

That is my understanding yes. I had some issues getting the non constant $k_i$ model to converge however it does converge now. The output below is for n=5. Adding in non constant $k_i$ seems to increase the runtime of the model by quite a bit for some reason. So I'll run n=30 in a bit (currently on the way home) and see how long it takes. Anyway here's the n=5 output:
View attachment 312935
View attachment 312936
I don't like how the MCO2 plot looks (the peaks should increase in size over time), so I'll look into that. But yes first I'll run n=30

n = 10 with non-constant $k_i$:

No such luck getting n=30 to converge. It seems to get stuck around 180K. I'm not sure what extra difficulty would be caused be n=30 that is not present in n=10 besides maybe some convergence issues. Looking into it. I'm leaving the simulation running just on the off chance it does actually converge after long times

 

[Chestermiller]

n = 10 with non-constant $k_i$:View attachment 312987
View attachment 312988
No such luck getting n=30 to converge. It seems to get stuck around 180K. I'm not sure what extra difficulty would be caused be n=30 that is not present in n=10 besides maybe some convergence issues. Looking into it. I'm leaving the simulation running just on the off chance it does actually converge after long times

What implicit integrator are you using? What is the structure of the Jacobian matrix, and do you take advantage of this in the matrix inversion chosen option>

 

[casualguitar]

What implicit integrator are you using? What is the structure of the Jacobian matrix, and do you take advantage of this in the matrix inversion chosen option>

The LSODA integrator (from the solve_ivp package). I have seen mentions of the Jacobian matrix in the documentation so I'll take a look at this

 

[casualguitar]

The LSODA integrator (from the solve_ivp package). I have seen mentions of the Jacobian matrix in the documentation so I'll take a look at this

Hi Chet, I looked into the Jacobian matrix question and if I understand the question the structure of the Jacobian matrix is just a rectangle of shape [n,6], where n is the number of spatial positions and 6 is the number of ODEs we have. Does this answer the question?

Also I'm not sure if you got the notification of the private message I sent. I attached the report as mentioned, and this also shows the Jacobian and solution array structures

What does it mean to take advantage of the structure of the Jacobian matrix in the matrix inversion chosen option?

 

[casualguitar]

What implicit integrator are you using? What is the structure of the Jacobian matrix, and do you take advantage of this in the matrix inversion chosen option>

Hi Chet,

As mentioned, here is the model output for the new Reynolds Number formulation:

We're now getting a clear constant temperature zone around the sublimation point of CO2, similar to the Tuinier model:

The mass buildup looks odd here in that the maximum buildup amount at a position seems to be constant, and also the value (approx 100) is about double the Tuinier value. Looking into this

 

[Chestermiller]

Hi Chet,

As mentioned, here is the model output for the new Reynolds Number formulation:

We're now getting a clear constant temperature zone around the sublimation point of CO2, similar to the Tuinier model:
View attachment 313411

The mass buildup looks odd here in that the maximum buildup amount at a position seems to be constant, and also the value (approx 100) is about double the Tuinier value. Looking into this

View attachment 313412

This looks like a major improvement. Is the only difference the equation used to calculate the Re? It doesn't seem possible.