Modelling of two phase flow in packed bed (continued)

[Chestermiller]

k/30:
View attachment 305128
View attachment 305130

View attachment 305129
Here is a plot of the rate of desublimation of CO2 for the same selection of times:
View attachment 305150

The MCO2 plot seems to suggest that there are two plugs of CO2 moving through the bed rather than one. However the dM_CO2_dt plot doesn't suggest this

So what is your interpretation of what is happening here?

Please provide a documentation of the finite difference equations for heat and mass transfer that are currently being used in the model.

 

[casualguitar]

So what is your interpretation of what is happening here?

I had turned both the MCO2 value and the dM_CO2_dt values to zero in the case where the solid buildup goes below zero, but really its just the MCO2 value that needs to be zero I think

Changing this (so now in the case of negative solid buildup the value is set to zero), here is the output for k/30:

Just checking some more ki/x values now to try the solid deposition to peak at a similar value to Tuinier et al

 

[casualguitar]

Just checking some more ki/x values now to try the solid deposition to peak at a similar value to Tuinier et al

After turning off the if statement that forced the rate of sublimation to be zero the plots look fairly close to Tuinier et al with some differences (k/100):

Adjusting the heat transfer coefficients:
If Ug and Ub are both = 100000 we get the plot below for Tb (I think very high Ub and Ug approximates the Tg=Tb assumption made in the Tuinier model?):

And here's Ub = 100000 and Ug = between 20 and 50 (calculated for each tank):

For the MCO2 plot, is their 'fudge factor' as you previously called it coming into play here? The max solid deposited seems to be in the ballpark of the Tuinier model, however our rate of sublimation of this solid seems to be a lot lower (their peaks are sharper).

Lastly, when Ug is not a very large number the constant temperature zone seems to appear for Tb, but not for Tg just yet

 

[Chestermiller]

After turning off the if statement that forced the rate of sublimation to be zero the plots look fairly close to Tuinier et al with some differences (k/100):
View attachment 305166
View attachment 305167
View attachment 305168

Adjusting the heat transfer coefficients:
If Ug and Ub are both = 100000 we get the plot below for Tb (I think very high Ub and Ug approximates the Tg=Tb assumption made in the Tuinier model?):
View attachment 305169

And here's Ub = 100000 and Ug = between 20 and 50 (calculated for each tank):

View attachment 305171
View attachment 305173

For the MCO2 plot, is their 'fudge factor' as you previously called it coming into play here? The max solid deposited seems to be in the ballpark of the Tuinier model, however our rate of sublimation of this solid seems to be a lot lower (their peaks are sharper).

Lastly, when Ug is not a very large number the constant temperature zone seems to appear for Tb, but not for Tg just yet

Refresh my memory. What do you use for Ub in the standard situation, and show its calculation.

 

[casualguitar]

Refresh my memory. What do you use for Ub in the standard situation, and show its calculation.

Apologies for the delay (public holiday here) -

$U_b$ is calculated by:
$$U_b = \frac{k_p}{\frac{d_p}{\beta}}$$
where $k_p$ = 0.8W/m., $d_p$ is 0.00404m and $\beta$ = 10

Using the above values, $U_b$ = 1980 w/m2.K

I checked the $k_p$ value (ranged it from 0.1 to 50) and changing the value of $k_p$ does not seem to affect the plots visually

 

[Chestermiller]

Apologies for the delay (public holiday here) -

$U_b$ is calculated by:
$$U_b = \frac{k_p}{\frac{d_p}{\beta}}$$
where $k_p$ = 0.8W/m., $d_p$ is 0.00404m and $\beta$ = 10

Using the above values, $U_b$ = 1980 w/m2.K

I checked the $k_p$ value (ranged it from 0.1 to 50) and changing the value of $k_p$ does not seem to affect the plots visually

Looks good.

 

[casualguitar]

Looks good.

Here is the output for $U_b$ calculated as above, and $U_g$ calculated at each interval

I think generally the output above is similar to the Tuinier output. The main two differences seem to be:
1) the rate of sublimation of the solid buildup (solid to gas) is lower in our model (gradually tails off rather than immediately sublimates
2) there is no constant temperature gas section showing up in our model

Checking the dTg_dt plot shows that there is no section where the value of dTg_dt reaches zero (a constant temperature section):

In the Tuinier model, what is the driving force that actually keeps the gas temperature almost constant while the phase is changing?

 

[Chestermiller]

Here is the output for $U_b$ calculated as above, and $U_g$ calculated at each interval
View attachment 305254View attachment 305255View attachment 305256View attachment 305257

I think generally the output above is similar to the Tuinier output. The main two differences seem to be:
1) the rate of sublimation of the solid buildup (solid to gas) is lower in our model (gradually tails off rather than immediately sublimates
2) there is no constant temperature gas section showing up in our model

Checking the dTg_dt plot shows that there is no section where the value of dTg_dt reaches zero (a constant temperature section):
View attachment 305259
In the Tuinier model, what is the driving force that actually keeps the gas temperature almost constant while the phase is changing?

Is it correct to say that these results were obtained using the corrected mass transfer description together with all other parameters standard?

I think the difference is that we have heat transfer resistance between the gas and bed and a different mass transfer parameterization (and different dispersion parameterization). We can make our results more like theirs by increasing Ug (I think). You can see that the bed temperature is closer to constant than the gas temperature at -90.

We can also decrease the dispersion by adding more tanks, holding the length of the column constant.

 

[casualguitar]

Is it correct to say that these results were obtained using the corrected mass transfer description together with all other parameters standard?

Yes these results were obtained using the Tuinier dimensions, flow and initial conditions as much as possible, with the corrected mass transfer coefficients of $k_{CO_2}$/100 and $k_{H_2O}$/100. The $U_b$ value is calculated according to the equation in post 305, and the$U_g$ value is calculated using the correlations from BSL

We can make our results more like theirs by increasing Ug (I think). You can see that the bed temperature is closer to constant than the gas temperature at -90.

Ok I see what you're saying. If we increase $U_g$ then the gas temperature will more closely follow the bed temperature. Checking a value of $U_g$= 100000 gives this $T_g$ profile which is visually identical to the $T_b$ profile:

As a note, if I calculate the $U_g$ values at each interval (which gives a range of about 20-50W/m2.K), then $T_b$ starts to show the almost constant temperature section, which doesn't happen at high $U_g$ values. The $T_g$ profile is unchanged from the one above:

Running the opposite (calculated $U_g$ and a high $U_b$ value gives identical output to the $T_b$ and $T_g$ profiles above)

We can also decrease the dispersion by adding more tanks, holding the length of the column constant.

Checking n=100 now for calculated $U_g$ and calculated $U_b$ values (no very high values). This will take a fair while to run so I'll post these results once it finishes

Is this considered 'model tuned'? I guess this is a judgement call rather than anything but it seems pretty close. If this is 'model tuned' I was thinking about taking what was learned during debugging here and applying it to the original liquefaction model (I'm guessing I can find some bugs in that now, and the model is simpler in comparison)

But besides this, I was thinking about some other things. The Tuinier paper (in my view) isn't that useful for a new reader because it doesn't actually say how good the system is (or can be), but rather just shows output for the system they used. What I mean is that they don't define performance parameters that would allow someone to compare this system with others or get a feel for how good this system actually is

To do something like this would involve (I think) answering questions like "what performance parameters are actually useful to calculate here?", "what simulations best show how these performance parameters vary with varied input" and lastly "does the model need to be developed further to be 'different enough' from the Tuinier model". Or something like these anyway

Regarding the performance parameters, I guess there are standard ones. However, possibly one meaningful parameter is the separation efficiency (how good can the system get at separating CO2 from both N2 and H2O in theory).

And lastly one other thing Tuinier didn't do is vary the ICs/BCs to arrive at an optimal solution (optimising the performance parameters). This sounds interesting to do and would also be novel as far as I can see (just trying to think of ways to further separate this model from Tuinier and make it of use to a future reader)

 

[Chestermiller]

Yes these results were obtained using the Tuinier dimensions, flow and initial conditions as much as possible, with the corrected mass transfer coefficients of $k_{CO_2}$/100 and $k_{H_2O}$/100. The $U_b$ value is calculated according to the equation in post 305, and the$U_g$ value is calculated using the correlations from BSL

So you used our k divided by 100. Please double check the finite difference equation used in our model for the mass flux from gas to solid to make sure we are not missing something (like say a factor of $\Delta z$).

Ok I see what you're saying. If we increase $U_g$ then the gas temperature will more closely follow the bed temperature. Checking a value of $U_g$= 100000 gives this $T_g$ profile which is visually identical to the $T_b$ profile:View attachment 305285

As a note, if I calculate the $U_g$ values at each interval (which gives a range of about 20-50W/m2.K), then $T_b$ starts to show the almost constant temperature section, which doesn't happen at high $U_g$ values. The $T_g$ profile is unchanged from the one above:
View attachment 305288

Running the opposite (calculated $U_g$ and a high $U_b$ value gives identical output to the $T_b$ and $T_g$ profiles above)Checking n=100 now for calculated $U_g$ and calculated $U_b$ values (no very high values). This will take a fair while to run so I'll post these results once it finishes

Is this considered 'model tuned'? I guess this is a judgement call rather than anything but it seems pretty close. If this is 'model tuned' I was thinking about taking what was learned during debugging here and applying it to the original liquefaction model (I'm guessing I can find some bugs in that now, and the model is simpler in comparison)

Model tuning is when you tune it to your own data. What we are doing here is BENCHMARKING the model against the Tunier model to make sure we are consistent. It seems to me that, at present, we are using much lower values of k than the equivalent they are using.

I suggest calculating the total amount of CO2 solid deposited on the bed so that we can compare it with the amount of CO2 that entered the bed up to any time. See what you get with the standard k and with k/100.

But besides this, I was thinking about some other things. The Tuinier paper (in my view) isn't that useful for a new reader because it doesn't actually say how good the system is (or can be), but rather just shows output for the system they used. What I mean is that they don't define performance parameters that would allow someone to compare this system with others or get a feel for how good this system actually is

To do something like this would involve (I think) answering questions like "what performance parameters are actually useful to calculate here?", "what simulations best show how these performance parameters vary with varied input" and lastly "does the model need to be developed further to be 'different enough' from the Tuinier model". Or something like these anyway

Regarding the performance parameters, I guess there are standard ones. However, possibly one meaningful parameter is the separation efficiency (how good can the system get at separating CO2 from both N2 and H2O in theory).

And lastly one other thing Tuinier didn't do is vary the ICs/BCs to arrive at an optimal solution (optimising the performance parameters). This sounds interesting to do and would also be novel as far as I can see (just trying to think of ways to further separate this model from Tuinier and make it of use to a future reader)

These are all good ideas once you are comfortable that you have a working model that matches your experimental data.

 

[casualguitar]

So you used our k divided by 100. Please double check the finite difference equation used in our model for the mass flux from gas to solid to make sure we are not missing something (like say a factor of Δz).

Yes k/100

When you say finite difference equation for the gas-solid mass flux do you mean this:
$$\dot{M}_i^"=k_i\left(\frac{Py_i-p_i(T_I)}{RT_I}\right)$$
Is this considered a 'finite difference equation'? If this is the equation then I'll check for missing terms. I could well have missed something in the $Sh_{loc,i}$ equation or something

I suggest calculating the total amount of CO2 solid deposited on the bed so that we can compare it with the amount of CO2 that entered the bed up to any time. See what you get with the standard k and with k/100.

Can do

These are all good ideas once you are comfortable that you have a working model that matches your experimental data.

Sure, I'll save this and return to it after this k/100 issue is sorted

 

[Chestermiller]

Yes k/100

When you say finite difference equation for the gas-solid mass flux do you mean this:
$$\dot{M}_i^"=k_i\left(\frac{Py_i-p_i(T_I)}{RT_I}\right)$$
Is this considered a 'finite difference equation'? If this is the equation then I'll check for missing terms. I could well have missed something in the $Sh_{loc,i}$ equation or something

No, I don't think this is the equation I was concerned with. I was concerned more with $$\frac{dM_i}{dt}=\dot{M}_i^"A_s$$where $$A_s=\frac{6}{d_p}(1-\epsilon)A_c\Delta z$$But I think we already checked this out and agreed that it is correct. Anyway, please check again.

I think we should also check the overall mass balance on the CO2 at the selected times. Please calculate the total moles of CO2 that have flowed into the bed up to time t, the total number of moles of CO2 that has exited the bed up to time t, the total molar holdup of CO2 in the gas phase in the column at time t, and the total molar holdup of CO2 on the bed at time t. In addition to these giving us a check on the CO2 mass balance, they will also tell the fraction of the CO2 that is captured by the bed up to a given time.

 

[casualguitar]

Please calculate the total moles of CO2 that have flowed into the bed up to time t

$$M_{IN,TOTAL} = \dot{m} * t\tag{1}$$

the total number of moles of CO2 that has exited the bed up to time t

Hmm. I think we could do that using the total molar holdup and the inlet flow, something like this:
$$M_{OUT,TOTAL} = \sum_{t=1}^{t}{\dot{m}_{IN} - m_{HOLDUP}}\tag{3}$$
Which I think is the sum of the outlet flows at all times in the bed

the total molar holdup of CO2 in the gas phase in the column at time t

Hmm so we have to consider each tank individually here given the temperature profile.
The molar holdup of a single tank is $(P/RT)(A_cdz\epsilon)$
So the total molar holdup would be:
$$M_{GAS} = \sum_{j=0}^{n}(P/RT_j)(A_cdz\epsilon)\tag{4}$$

and the total molar holdup of CO2 on the bed at time t

I guess this is solid phase only, so:
$$M_{SOLID} = \sum_{j=0}^{n}(M_j)\tag{5}$$

they will also tell the fraction of the CO2 that is captured by the bed up to a given time.

When you say 'fraction of CO2 captured' what does this mean? Do we define 'captured' as meaning solid phase? So the fraction would be the solid mass of CO2 at a point divided by the total mass of CO2 at that point?

I'll work away on the above calculations and if they're incorrect I can adjust

 

[casualguitar]

$$M_{IN,TOTAL} = \dot{m} * t\tag{1}$$

Hmm. I think we could do that using the total molar holdup and the inlet flow, something like this:
$$M_{OUT,TOTAL} = \sum_{t=1}^{t}{\dot{m}_{IN} - m_{HOLDUP}}\tag{3}$$
Which I think is the sum of the outlet flows at all times in the bed

Hmm so we have to consider each tank individually here given the temperature profile.
The molar holdup of a single tank is $(P/RT)(A_cdz\epsilon)$
So the total molar holdup would be:
$$M_{GAS} = \sum_{j=0}^{n}(P/RT_j)(A_cdz\epsilon)\tag{4}$$

I guess this is solid phase only, so:
$$M_{SOLID} = \sum_{j=0}^{n}(M_j)\tag{5}$$

When you say 'fraction of CO2 captured' what does this mean? Do we define 'captured' as meaning solid phase? So the fraction would be the solid mass of CO2 at a point divided by the total mass of CO2 at that point?

I'll work away on the above calculations and if they're incorrect I can adjust

I'll separate these plots out, however just for an overview these are the values mentioned above plotted versus time (note I've assumed CO2 captured is the mass of solid CO2 in the column divided by the total mass of CO2 in the column):

 

[Chestermiller]

$$M_{IN,TOTAL} = \dot{m} * t\tag{1}$$

This needs to be multiplied by $y_{in}$ for CO2

Hmm. I think we could do that using the total molar holdup and the inlet flow, something like this:
$$M_{OUT,TOTAL} = \sum_{t=1}^{t}{\dot{m}_{IN} - m_{HOLDUP}}\tag{3}$$
Which I think is the sum of the outlet flows at all times in the bed

This assumes that the mass balance is truly satisfied. The objective of this exercise is to check to see if that is the case. If $\dot{m}_i(t)$ is the total molar flow rate out of tank i at time t, then the total molar amount of CO2 leaving the column up to time t is$$\int_0^t{\dot{m}_n(t')y_{CO2}(t')dt'}$$

Hmm so we have to consider each tank individually here given the temperature profile.
The molar holdup of a single tank is $(P/RT)(A_cdz\epsilon)$
So the total molar holdup would be:
$$M_{GAS} = \sum_{j=0}^{n}(P/RT_j)(A_cdz\epsilon)\tag{4}$$

The total molar holdup of CO2 in the gas would be $$M_{CO2} = \sum_{j=1}^{n}(P/RT_j)y_{j,CO2}(A_cdz*\epsilon)\tag{4}$$
Note that the summation is from tank 1 to tank n, the mole fraction of CO2 multiplies the total number of moles, and that you multiply by the void fraction rather than dividing.

I guess this is solid phase only, so:
$$M_{SOLID} = \sum_{j=0}^{n}(M_j)\tag{5}$$

Again, the sum is from 1 to n.

When you say 'fraction of CO2 captured' what does this mean? Do we define 'captured' as meaning solid phase?

Yes

So the fraction would be the solid mass of CO2 at a point divided by the total mass of CO2 at that point?

I don't know what this means, but it doesn't sound correct.

I'll work away on the above calculations and if they're incorrect I can adjust

 

[casualguitar]

This needs to be multiplied by $y_{in}$ for CO2

Change made

This assumes that the mass balance is truly satisfied. The objective of this exercise is to check to see if that is the case. If $\dot{m}_i(t)$ is the total molar flow rate out of tank i at time t, then the total molar amount of CO2 leaving the column up to time t is$$\int_0^t{\dot{m}_n(t')y_{CO2}(t')dt'}$$

To implement the integration in code, can I do a cumulative sum here like this:
$$\int_0^t{\dot{m}_n(t')y_{CO2}(t')dt'} = \sum_{t=0}^t{\dot{m}_{out}y_{CO2,OUT}}$$
which is just calculating the outlet CO2 flow at all times and summing them together

The total molar holdup of CO2 in the gas would be $$M_{CO2} = \sum_{j=1}^{n}(P/RT_j)y_{j,CO2}(A_cdz*\epsilon)\tag{4}$$
Note that the summation is from tank 1 to tank n, the mole fraction of CO2 multiplies the total number of moles, and that you multiply by the void fraction rather than dividing.

Change made

I don't know what this means, but it doesn't sound correct.

In regards to defining the 'fraction of CO2 captured up to a given time'. So we agreed that 'captured' means 'in the solid phase', but what is this a fraction of?

I was thinking:
Fraction of CO_2 captured = $\frac{M_{CO2,SOLID}}{M_{CO2,SOLID} + M_{CO2,GAS}}$ at a given time, which is the moles of CO2 divided by the total mass of CO2 in the packed bed, at a given time

What this definition doesn't capture is the 'fraction of CO2 captured $\textbf{up to a given time}$. The above definition can only calculate the fraction of CO2 captured $\textbf{at to a given time}$, so I guess my definition above isn't correct, unless you meant at a given time rather than up to a given time?

I'll post those changed plots once the simulation finishes

 

[casualguitar]

Heres a plot of all the above variables versus time:

And zooming in on the total column CO2 gas holdup because its much smaller than the other values:

 

[Chestermiller]

Change made

To implement the integration in code, can I do a cumulative sum here like this:
$$\int_0^t{\dot{m}_n(t')y_{CO2}(t')dt'} = \sum_{t=0}^t{\dot{m}_{out}y_{CO2,OUT}}$$
which is just calculating the outlet CO2 flow at all times and summing them together

Just use the trapezoidal rule, and don't forget to multiply by $\Delta t$

Change made

In regards to defining the 'fraction of CO2 captured up to a given time'. So we agreed that 'captured' means 'in the solid phase', but what is this a fraction of?

I was thinking:
Fraction of CO_2 captured = $\frac{M_{CO2,SOLID}}{M_{CO2,SOLID} + M_{CO2,GAS}}$ at a given time, which is the moles of CO2 divided by the total mass of CO2 in the packed bed, at a given time

What this definition doesn't capture is the 'fraction of CO2 captured $\textbf{up to a given time}$. The above definition can only calculate the fraction of CO2 captured $\textbf{at to a given time}$, so I guess my definition above isn't correct, unless you meant at a given time rather than up to a given time?

I'll post those changed plots once the simulation finishes

The fraction of CO2 captured up to any time t would be the amount on the bed at time t divided by the amount entering up to time t.

 

[Chestermiller]

Heres a plot of all the above variables versus time:
View attachment 305562

And zooming in on the total column CO2 gas holdup because its much smaller than the other values:
View attachment 305563

How can the CO2 moles out exceed the CO2 moles in? You need to multiply by $\Delta t$ in getting moles out.

At short times, how can CO2 holdup on the column exceed the CO2 that entered up to time t?

 

[casualguitar]

Just use the trapezoidal rule, and don't forget to multiply by

Yes $\Delta t$ is 1 second so I think the above applies

The fraction of CO2 captured up to any time t would be the amount on the bed at time t divided by the amount entering up to time t.

Got it, change made (wont post the plot yet because of that issue with the total solid deposited being greater than the total CO2 that entered the system)

How can the CO2 moles out exceed the CO2 moles in? You need to multiply by Δt in getting moles out.

Well if the inlet flow is less dense than the density of the fluid in the tank then we would have a negative accumulation of CO2 i.e. moles out > moles in? $\Delta t$ is 1s here

At short times, how can CO2 holdup on the column exceed the CO2 that entered up to time t?

Ah I see that. Hmm doesn't make sense. Checking my calculations. Actually if I increase the mole fraction of CO2 in the inlet flow then this issue goes away. Checking this now anyway

 

[Chestermiller]

Well if the inlet flow is less dense than the density of the fluid in the tank then we would have a negative accumulation of CO2 i.e. moles out > moles in?

No freakin' way. There is no CO2 in the column to begin with. So the cumulative amount coming out can never be greater than the cumulative amount you put in.

 

[Chestermiller]

Calculation of total moles of CO2 entering column.

Mass flow in = 0.27 kg/m^2-s
Mean molecular weight = (28)(0.75)+(44)(0.2)+(18)(0.05)=30.7
Area of column cross section =$\pi(0.035)^2/4=0.000962\ m^2$
molar flow rate of CO2 in = $\frac{(270)(0.000962)(0.2)}{30.7}=0.00169\ moles/sec$

At 250 sec, moles CO2 that entered = 0.423 moles. This is consistent with your graph.

 

[casualguitar]

Calculation of total moles of CO2 entering column.

Mass flow in = 0.27 kg/m^2-s
Mean molecular weight = (28)(0.75)+(44)(0.2)+(18)(0.05)=30.7
Area of column cross section =$\pi(0.035)^2/4=0.000962\ m^2$
molar flow rate of CO2 in = $\frac{(270)(0.000962)(0.2)}{30.7}=0.00169\ moles/sec$

At 250 sec, moles CO2 that entered = 0.423 moles. This is consistent with your graph.

Meaning that the error is in the solid CO2 buildup graph. I'll manually recreate the solid CO2 buildup array and see if it matches the solution from the integrator

 

[casualguitar]

I regenerated the solid CO2 buildup array after the integrator (rather than extracting the solution from the integrator), and plotted the integrator solution (blue) versus the manually recreated solution (orange):

The manually created orange one has a better shape in that the buildup declines to zero at a rate closer to that seen in Tuinier et al.

I made the orange array with the assumption that the total solid CO2 on the bed at any time will be equal to the total solid CO2 on the bed at the previous time plus the net desublimation rate across the bed (i.e. summing the desublimation rate for each position), times the surface area available for desublimation in each bed. So something like this:
$$M_{CO2, t} = M_{CO2,t-1} + \sum_{j=1}^nMDR_t*SA$$
where MDR is the molar desublimation rate in a tank and SA is the surface area available for desublimation in the tank

 

[casualguitar]

One further note:

Checking the CO2 captured fraction up to time t using your definition above gives the following output (plotting time vs CO2 captured fraction using both the integrator CO2 buildup and the manually calculated CO2 buildup). I plotted both against the gas temperature (in tank 1) just to give a reference point

Integrator:

Manually calculated:

The captured fraction for the integrator values goes above 1 so there is something wrong there. The manually calculated version looks somewhat better. Once the temperature approximately exceeds the desublimation temperature the captured fraction starts to reduce

Looks like an error in how I've set up the integrator to calculate the mass buildup

 

[casualguitar]

Looks like an error in how I've set up the integrator to calculate the mass buildup

Just looking into this integrator integration a bit more - I'm printing out the value of both dMCO2/dt and MCO2 at each time, to check that $M_{CO2,t} = M_{CO2,t-1} + dMCO2dt$ when MCO2>0 and it doesn't seem to be:

There are two issues here. Firstly the rate of deposition is negative at the start, and secondly the above equation doesn't seem to apply anywhere.

Seeing that the sign is negative at the beginning, its possible that these are actually sublimation rates (which would be negative at the beginning), meaning that my logic to deal with MCO2<0 might be incorrect

 

[Chestermiller]

Just looking into this integrator integration a bit more - I'm printing out the value of both dMCO2/dt and MCO2 at each time, to check that $M_{CO2,t} = M_{CO2,t-1} + dMCO2dt$ when MCO2>0 and it doesn't seem to be:
View attachment 305615
There are two issues here. Firstly the rate of deposition is negative at the start, and secondly the above equation doesn't seem to apply anywhere.

Seeing that the sign is negative at the beginning, its possible that these are actually sublimation rates (which would be negative at the beginning), meaning that my logic to deal with MCO2<0 might be incorrect

keep at it. you'll get there. This is why I suggested checking the overall balances.

 

[casualguitar]

keep at it. you'll get there. This is why I suggested checking the overall balances.

Plotting the total solid buildup (somewhat corrected) for both the integrator (blue) and manually (orange), actually shows the integrator perfectly matches the manual recreation up until t=200, and then the manual one starts to decrease while the integrator increases further.

Anyway just working through this. Just said I'd post in case a plot like this is indicative of some error type in particular. It may be a sign error at this point (t=200) since they somewhat mirror one another from t=200 through y=0.13 (approx).

 

[casualguitar]

I'm not sure how this didn't register with me earlier but yes as you said before $\dot{m}_{out}$ is clearly greater than $\dot{m}_{in}$ and this is not possible

Some thinking out loud here-
To confirm, the equation I'm using to calculate the outlet flowrate is:
$$\dot{m}_{OUT} = \dot{m}_{IN} + \frac{m_j}{T_g}*\frac{\partial T_g}{\partial t} - ((MDR_{CO2} + MLR_{H2O})*SA)$$ where MDR and MLR are the desublimation and liquefaction rates of CO2 and h2o.

So in what cases will $\dot{m}_{out}$ be greater than $\dot{m}_{in}$? Looking at the terms, $m_j$ and $T_g$ will always be positive, as will \frac{\partial T_g}{\partial t} unless its zero (steady state)

Checking the last term (the $((MDR_{CO2} + MLR_{H2O})*SA)$ bit), printing out the integrator $MLR_{H2O}$ values show that it is positive in all cases in this simulation. This is ok for now because the gas temperature only gets up to 200K. The $MDR_{CO2}$ values start out positive and then go negative as the temperature goes up, as expected.

So to have $\dot{m}_{out}$ < $\dot{m}_{in}$ we need $MDR_{CO2}$ + $MLR_{H2O}$ to be positive? Hmm. So $MLR_{H2O}$ is always positive (for now), therefore all we need is for $MDR_{CO2}$ to be positive, or have a smaller absolute value than $MLR_{H2O}$

Looking at the equation for these, the general form is:
$$M_i'' = \frac{k_i(Py_i - P(T_I))}{RT}$$
So $k_i$, $P$ and $R$ are constants, and T is the same for both h2o and co2, so the only questionable terms here are the desublimation/liquefaction pressures and $y_i$.

Two questions come up at this point - is it not true that varying $y_{CO2}$ will mean that in some cases $MDR_{CO2}$ will be greater than $MLR_{H2O}$, and sometimes less than? Meaning that $\dot{m}_{out}$ < $\dot{m}_{in}$ is dependent on $y_{CO2}$ being in a certain range?

 

[Chestermiller]

I'm not sure how this didn't register with me earlier but yes as you said before $\dot{m}_{out}$ is clearly greater than $\dot{m}_{in}$ and this is not possible

Some thinking out loud here-
To confirm, the equation I'm using to calculate the outlet flowrate is:
$$\dot{m}_{OUT} = \dot{m}_{IN} + \frac{m_j}{T_g}*\frac{\partial T_g}{\partial t} - ((MDR_{CO2} + MLR_{H2O})*SA)$$ where MDR and MLR are the desublimation and liquefaction rates of CO2 and h2o.

So in what cases will $\dot{m}_{out}$ be greater than $\dot{m}_{in}$? Looking at the terms, $m_j$ and $T_g$ will always be positive, as will \frac{\partial T_g}{\partial t} unless its zero (steady state)

Checking the last term (the $((MDR_{CO2} + MLR_{H2O})*SA)$ bit), printing out the integrator $MLR_{H2O}$ values show that it is positive in all cases in this simulation. This is ok for now because the gas temperature only gets up to 200K. The $MDR_{CO2}$ values start out positive and then go negative as the temperature goes up, as expected.

So to have $\dot{m}_{out}$ < $\dot{m}_{in}$ we need $MDR_{CO2}$ + $MLR_{H2O}$ to be positive? Hmm. So $MLR_{H2O}$ is always positive (for now), therefore all we need is for $MDR_{CO2}$ to be positive, or have a smaller absolute value than $MLR_{H2O}$

Looking at the equation for these, the general form is:
$$M_i'' = \frac{k_i(Py_i - P(T_I))}{RT}$$
So $k_i$, $P$ and $R$ are constants, and T is the same for both h2o and co2, so the only questionable terms here are the desublimation/liquefaction pressures and $y_i$.

Two questions come up at this point - is it not true that varying $y_{CO2}$ will mean that in some cases $MDR_{CO2}$ will be greater than $MLR_{H2O}$, and sometimes less than? Meaning that $\dot{m}_{out}$ < $\dot{m}_{in}$ is dependent on $y_{CO2}$ being in a certain range?

I have no idea what you are doing here. We are looking at the overall mass balance of CO2 in the column, not the instantaneous amounts entering and leaving.

(Total moles of CO2 that entered column up to time t) = (Total moles of CO2 that exited column up to time t) + (Total moles of CO2 accumulated in column at present time t)

 

[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}$