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)