Modelling of two phase flow in packed bed (continued)

[casualguitar]

Hi Chet, just to be clear on the $k_i$ units question, what I mean is that all of the correlations say that the units of $k_i$ are $mol/m^2.s$, but, the desublimation rate equation says the units of $k_i$ are $m/s$.

Equations that say $k_i$ unit is mol/m2.s:
$Sh_{loc,i} = \frac{k_id_p}{\rho_mD_{ab}}$
$Sh_{loc,i}=\frac{k_id_p}{(1-\epsilon_g)D_i\psi}$

Equation that says $k_i$ unit is m/s:
$\dot{M}_i^"=k_i\left(\frac{Py_i-p_i(T_I)}{RT_I}\right)$

 

[casualguitar]

Hi Chet, by any chance would you have a source for the molar deposition rate equation (its possibly a variation of the Hertz-Knudsen equation)? I want to look into why the units don't seem to match up for ki in the Sherwood number correlations and ki in the molar deposition rate equation

Edit: any variation of this equation I have found has had a square root and pi in it which seems quite different to what we're using currently M˙i"=kiPyi−pi(TI)RTI

Edit: The closest equation I could find to our one is this:

v is in kg/s rather than mol/m2.s here, I haven't managed to find the equation you used just yet (still looking)

 

[casualguitar]

My apologies for the multiple posts. Is it possible that this equation $\dot{M}_i^"=k_i\left(\frac{Py_i-p_i(T_I)}{RT_I}\right)$ is missing a molar volume term on the right hand side (units of m3/mol). If so, then it essentially becomes $\dot{M}_i^"=k_i\left(\frac{Py_i-p_i(T_I)}{RT_I}\right)*\rho_m$ which is $k_i*\frac{PV}{nRT}$. The brackets are dimensionless and this would mean that $k_i$ would have the same units as $M_i''$, and this would line up with the units of $k_i$ in the sherwood number correlations

 

[Chestermiller]

Lastly, regarding the units of ki, any equation with ki in it shows that the unit of ki is $mol/m^2.s$. However, the molar desublimation rate formula $\dot{M}_i^"=k_i\left(\frac{Py_i-p_i(T_I)}{RT_I}\right)\tag{1}$ suggests that the units are $m/s$. If this is correct, is there an equivalent molar desublimation rate equation that would use units of mol/m2.s for ki?

For reference, my units are:
$k_i$ = $mol/m^2.s$
$P$ = $Pa$
$y_i$ = $mol/mol$
$R$ = $J/mol.K$
$T$ = $K$

If we write the equation $\dot{M}_i^"=k_i\left(\frac{Py_i-p_i(T_I)}{RT_I}\right)\tag{1}$ as $$\dot{M}_i^"=k_ic_I(y_i-p_i(T_I)/P)$$where $c=\frac{P}{RT_I}$, then we can write $$\dot{M}_i^"=K_i(y_i-p_i(T_I)/P)$$

By analogy to Eqn. 14.5-6 of BSL, $$Sh=\frac{K_iD_p}{cD_i(1-\epsilon)\psi}$$Combining these gives $$\dot{M}_i^"=Sh \frac{c_ID_i(1-\epsilon)\psi}{D_p}(y_i-p_i(T_I)/P)$$

 

[casualguitar]

If we write the equation $\dot{M}_i^"=k_i\left(\frac{Py_i-p_i(T_I)}{RT_I}\right)\tag{1}$ as $$\dot{M}_i^"=k_ic_I(y_i-p_i(T_I)/P)$$where $c=\frac{P}{RT_I}$, then we can write $$\dot{M}_i^"=K_i(y_i-p_i(T_I)/P)$$

By analogy to Eqn. 14.5-6 of BSL, $$Sh=\frac{K_iD_p}{cD_i(1-\epsilon)\psi}$$Combining these gives $$\dot{M}_i^"=Sh \frac{c_ID_i(1-\epsilon)\psi}{D_p}(y_i-p_i(T_I)/P)$$

Thats interesting. What you did above shows that $M_i''$ has the same units as $K_i$. Issue solved then. Great! It looks like you did the same thing I did (there was a density factor I had not considered). I have this written up in code form so we're almost there on the functions front also. I have pretty much everything plotted for ranges of dependent variables also

Regarding the gas interface heat transfer coefficient $U_g$, you mentioned I should continue to use the Nu correlation for flow in a packed bed, and add the wall viscosity correction. By this, do you mean:
$U_{g,NEW} = U_{g,OLD}*\frac{\mu_b}{\mu_0}^{0.14}$

Besides that question I think that's all issues solved for now. I just need to code up the last few heat transfer related functions (QgI, QIb, etc)

 

[Chestermiller]

Thats interesting. What you did above shows that $M_i''$ has the same units as $K_i$. Issue solved then. Great! It looks like you did the same thing I did (there was a density factor I had not considered). I have this written up in code form so we're almost there on the functions front also. I have pretty much everything plotted for ranges of dependent variables also

Regarding the gas interface heat transfer coefficient $U_g$, you mentioned I should continue to use the Nu correlation for flow in a packed bed, and add the wall viscosity correction. By this, do you mean:
$U_{g,NEW} = U_{g,OLD}*\frac{\mu_b}{\mu_0}^{0.14}$

Besides that question I think that's all issues solved for now. I just need to code up the last few heat transfer related functions (QgI, QIb, etc)

It's only an approximate thing. The viscosity term is only going to be a minor correction.

 

[casualguitar]

It's only an approximate thing. The viscosity term is only going to be a minor correction.

Perfect. Sure. This correction is in there now anyway. I'm tempted to post plots of all the main functions (Re, Sc, Pr, ki, U_g, etc) for suitable ranges of dependent variables later today to be sure they're ok before moving on

Actually one thing I noticed in writing the code for these functions is that there are 'layers' of functions. For example, viscosity is dependent on temperature, Reynolds number is dependent on viscosity, the mole transfer coefficient is dependent on the Reynolds number, and the molar deposition rate is dependent on the mole transfer coefficient.

This makes defining functions a bit messy in that I have for example the molar deposition rate function as a function of T, the sherwood number, the mass transfer coefficient, the mole fraction and the diffusion coefficient. I am tempted to look for a set of 'basis' variables that all functions are dependent on. For example, if I have every function as a function of temperature and the molar flow, would this be enough to calculate everything in any 'layer'? (assuming constants like P and $y_i$ are defined elsewhere). For now I will finish off these plots but it looks like finding this common basis would be useful for code simplicity

 

[casualguitar]

Hi Chet, issue above resolved after chatting with a postdoc. Looks like the standard/OOP approach is to not do what I was thinking of doing above. This is fine as my code is already in this format anyway

One question - did we agree that $Ug = h_{fs}$? I know in the last model we had the lumped parameter approach where we had $\frac{1}{h} = \frac{1}{h_{fs}} + \frac{dp}{k*beta}$ however here we are saying $h = h_{fs}$?

Every function is written up now, and most plots are done. I've been able to check any of the thermodynamic property functions via literature data, however there are some other plots I'd like to check with you if possible. Also, most parameters have a number of dependent variables, so I would like to ask you which dependent variable is most insightful to vary for the plots. Here are my thoughts on possible plots and dependent variables I would vary within a range to create the plots. In the format '(The function output -> the variable being varied)':

Reynolds number -> molar flow
Prandtl number -> viscosity
Sherwood number -> mass transfer coefficient
Nusselt number -> $h_{fs}$
Schmidt number -> viscosity
Mass transfer coefficient -> Reynolds number
Molar deposition rate -> mass transfer coefficient
$h_{fs}$ -> Reynolds number

Working on these plots this morning. If you think there are other useful plots at this stage, or that there are better dependent variables to vary than the ones I chose just let me know

 

[Chestermiller]

Hi Chet, issue above resolved after chatting with a postdoc. Looks like the standard/OOP approach is to not do what I was thinking of doing above. This is fine as my code is already in this format anyway

One question - did we agree that $Ug = h_{fs}$? I know in the last model we had the lumped parameter approach where we had $\frac{1}{h} = \frac{1}{h_{fs}} + \frac{dp}{k*beta}$ however here we are saying $h = h_{fs}$?

$h_{fs}$ is the heat transfer coefficient on the gas side, right?
Also, I think in this model we are calling the heat transfer coefficient on the gas side, the heat transfer coefficient on the solid bed side, and the overall heat transfer coefficient as $U_g$, $U_b$, and U. I'm not really sure. I lost track of the nomenclature we are using. You need to go back to my posts.

Every function is written up now, and most plots are done. I've been able to check any of the thermodynamic property functions via literature data, however there are some other plots I'd like to check with you if possible. Also, most parameters have a number of dependent variables, so I would like to ask you which dependent variable is most insightful to vary for the plots. Here are my thoughts on possible plots and dependent variables I would vary within a range to create the plots. In the format '(The function output -> the variable being varied)':

Reynolds number -> molar flow
Prandtl number -> viscosity
Sherwood number -> mass transfer coefficient
Nusselt number -> $h_{fs}$
Schmidt number -> viscosity
Mass transfer coefficient -> Reynolds number
Molar deposition rate -> mass transfer coefficient
$h_{fs}$ -> Reynolds number

Working on these plots this morning. If you think there are other useful plots at this stage, or that there are better dependent variables to vary than the ones I chose just let me know

The comparisons that make sense to me are in terms of the dimensionless groups, say $Nu/Pr^{1/3}$ vs Re (say comparing the correlation for turbulent flow in a tube to that in a packed bed).

 

[casualguitar]

Also, I think in this model we are calling the heat transfer coefficient on the gas side, the heat transfer coefficient on the solid bed side, and the overall heat transfer coefficient as Ug, Ub, and U. I'm not really sure. I lost track of the nomenclature we are using. You need to go back to my posts.

Yep looking back at those posts we are doing this

The comparisons that make sense to me are in terms of the dimensionless groups, say Nu/Pr1/3 vs Re (say comparing the correlation for turbulent flow in a tube to that in a packed bed).

Almost finished with those other plots, so I'll add this one in then also. What information does a Nu/Pr^1/3 vs Re plot actually tell us? I suppose it can tell you how dominant convective heat transfer is for a range of Re values?

Edit: Will add all of these plots first thing tomorrow. If there are any other useful ones I can easily plot these two as I have functions for everything we would need to calculate now

 

[casualguitar]

Just posting some plots here:

Molar flow vs Reynolds number (using system specific dimensions):

Reynolds Number vs Mass Transfer Coefficient for CO2 and H2O:

Reynolds Number vs Nu/Pr^1/3:

If these look ok to you I think it's about time to start adding these functions into the simulation? Would you recommend any general order for adding these functions in? i.e. should I leave out any functions (and leave some constants) for now, or just go straight for adding them all in?

 

[Chestermiller]

Just posting some plots here:

Molar flow vs Reynolds number (using system specific dimensions):
View attachment 301954
Reynolds Number vs Mass Transfer Coefficient for CO2 and H2O:
View attachment 301955

Reynolds Number vs Nu/Pr^1/3:
View attachment 301956

If these look ok to you I think it's about time to start adding these functions into the simulation? Would you recommend any general order for adding these functions in? i.e. should I leave out any functions (and leave some constants) for now, or just go straight for adding them all in?

What is your judgment on this? My rule of thumb for modeling is always "start simple and build in complexity."

 

[casualguitar]

What is your judgment on this? My rule of thumb for modeling is always "start simple and build in complexity."

I do agree with this, and in addition it would be nice to start with a simple model that we can 'validate' against another model or against what we would intuitively expect to happen.

If we were to look at the mole fraction of co2/h2o in a single tank versus time, what trend would we expect?

I guess given that the ambient stream is hitting a cold bed we would see a yCO2 and yH2O spikes, and then when this bed heats up again the mole fractions would gradually decrease to zero?

If this is true, it at least offers some form of sanity checking.

As for simplifying it down initially though, we would need to implement the minimum amount of functions that allows for this trend to be seen.

It seems that the molar deposition rate is essential here (can't have this as a constant). Other possible requirements are $Q_{IB}$ and $Q_{GI}$

What about leaving everything (all property functions and heat related functions) as constants, and just implementing variable $M_i''$ for a single tank?

Edit: Also in addition why the interest in a Nu/Pr^1/3 plot? What information does this convey?

Edit: Lastly just posting the actual equations here to bring them forward:

 

[Chestermiller]

I do agree with this, and in addition it would be nice to start with a simple model that we can 'validate' against another model or against what we would intuitively expect to happen.

If we were to look at the mole fraction of co2/h2o in a single tank versus time, what trend would we expect?

I guess given that the ambient stream is hitting a cold bed we would see a yCO2 and yH2O spikes, and then when this bed heats up again the mole fractions would gradually decrease to zero?

If this is true, it at least offers some form of sanity checking.

As for simplifying it down initially though, we would need to implement the minimum amount of functions that allows for this trend to be seen.

It seems that the molar deposition rate is essential here (can't have this as a constant). Other possible requirements are $Q_{IB}$ and $Q_{GI}$

No, but you can impose constant heat transfer- and mass transfer coefficients.

What about leaving everything (all property functions and heat related functions) as constants, and just implementing variable $M_i''$ for a single tank?

Try it and see how it plays out.

Edit: Also in addition why the interest in a Nu/Pr^1/3 plot? What information does this convey?

It would be interesting to compare this with the relationship with the corresponding variation for turbulent flow in a pipe.

Edit: Lastly just posting the actual equations here to bring them forward:
View attachment 301961

I think these are OK. Recheck the algebra to be sure.

 

[casualguitar]

No, but you can impose constant heat transfer- and mass transfer coefficients.

Ah nice, yes that sounds good. It means that temperature, heat transfer rate an molar deposition rate are all that will vary in this basic model. I guess this will be enough to see the yCO2 and yH2O trends

Try it and see how it plays out.

Will do this that morning then

It would be interesting to compare this with the relationship with the corresponding variation for turbulent flow in a pipe.

Can do. Why would this be interesting? It would quantify the effect of the packing on heat transfer coefficient I suppose

 

[casualguitar]

Ah nice, yes that sounds good. It means that temperature, heat transfer rate an molar deposition rate are all that will vary in this basic model. I guess this will be enough to see the yCO2 and yH2O trends

Will do this that morning then

Can do. Why would this be interesting? It would quantify the effect of the packing on heat transfer coefficient I suppose

Hi Chet, have all functions etc for the current 'constant heat and mass transfer coefficient' model in the model. The model is giving some errors etc so I'm currently debugging these. Will post a condensed version of the initial conditions and output if I can't debug it

 

[casualguitar]

Hi Chet, have all functions etc for the current 'constant heat and mass transfer coefficient' model in the model. The model is giving some errors etc so I'm currently debugging these. Will post a condensed version of the initial conditions and output if I can't debug it

I haven't solved this yet, but I did find something odd. Solve_ivp offers a number of integrators (LSODA, Radau, RK45, RK23, DOP853, BDF). I tried a few and they all seem to return very different results. Actually all of them show the gas temperature reaching temperatures that are not present in the system at all which indicates I'm doing something wrong anyway, but regardless I didn't expect the output from each integrator do vary this much

 

[casualguitar]

I haven't solved this yet, but I did find something odd. Solve_ivp offers a number of integrators (LSODA, Radau, RK45, RK23, DOP853, BDF). I tried a few and they all seem to return very different results. Actually all of them show the gas temperature reaching temperatures that are not present in the system at all which indicates I'm doing something wrong anyway, but regardless I didn't expect the output from each integrator do vary this much

For reference, these are the initial, boundary, and constant values I'm using. Anything stand out here to you as being a value that is outside its typical range?

n = 5 #number of tanks
rho0CO2 = 1.1 #initial density of CO2 in bed (mol/m3)
rho0H2O = 1.1 #initial density of h2o in bed (mol/m3)
y0CO2 = 0 #initial gas phase mole fraction of CO2 in bed (mol/mol)
y0H2O = 0 #initial gas phase mole fraction of h2o in bed (mol/mol)
M0CO2 = 0.0 #initial solid CO2 moles deposited on bed (mol/m2)
M0H2O = 0.0 #initial solid h2o moles deposited on bed (mol/m2)
Tg0 = 150 #initial gas temperature (K)
Tb0 = 150 #initial bed temperature (K)
A_C = 0.005 #cross sectional area (m2)
A_s = 267 #specific surface area of solid (m2/m3)
h_vap_h2o = 40650 (J/mol)
v_desublimation_co2 = 26000 (J/mol)
k_s = 18 #solid heat capacity (W/m.K)
dz = 0.01
U_b = 100 #bed heat transfer coefficient (W/m2.K)
U_g = 200 #gas phase heat transfer coefficient (W/m2.K)
cp_CO2 = 45 #co2 heat capacity (J/mol.K)
cp_H2O = 45 #h2o heat capacity (J/mol.K)
ki_co2 = 8 #co2 mass transfer coefficient (mol/m2.s)
ki_h2o = 16 #h2o Mass Transfer coefficient (mol/m2.s)
m_co2 = rho0CO2*dz*A_C (value of 5*10^-5)
m_h2o = rho0H2O*dz*A_C (value of 5*10^-5)
M_al = 72.2 (mol/tank) (about 7kg/tank)

#Boundary conditions
mol_in = 0.5 #mol/s
y_co2_in = 0.1 #mol/mol
y_h2o_in = 0.01 #mol/mol
T_in = 220 #K

 

[Chestermiller]

For reference, these are the initial, boundary, and constant values I'm using. Anything stand out here to you as being a value that is outside its typical range?

n = 5 #number of tanks

rho0CO2 = 1.1 #initial density of CO2 in bed (mol/m3)
rho0H2O = 1.1 #initial density of h2o in bed (mol/m3)

What are these supposed to represent? I thought there is no water or CO2 deposited on the bed initially, and there is no water or CO2 in the gas phase initially.

y0CO2 = 0 #initial gas phase mole fraction of CO2 in bed (mol/mol)
y0H2O = 0 #initial gas phase mole fraction of h2o in bed (mol/mol)
M0CO2 = 0.0 #initial solid CO2 moles deposited on bed (mol/m2)
M0H2O = 0.0 #initial solid h2o moles deposited on bed (mol/m2)
Tg0 = 150 #initial gas temperature (K)
Tb0 = 150 #initial bed temperature (K)
A_C = 0.005 #cross sectional area (m2)

Bed diameter = 8 cm?

A_s = 267 #specific surface area of solid (m2/m3)

Please provide the calculation which led to this.

h_vap_h2o = 40650 (J/mol)
v_desublimation_co2 = 26000 (J/mol)
k_s = 18 #solid heat capacity (W/m.K)

thermal conductivity, not heat capacity

dz = 0.01

Units are meters? So total bed length is 0.05 m = 5 cm?

U_b = 100 #bed heat transfer coefficient (W/m2.K)
U_g = 200 #gas phase heat transfer coefficient (W/m2.K)
cp_CO2 = 45 #co2 heat capacity (J/mol.K)
cp_H2O = 45 #h2o heat capacity (J/mol.K)

cp of N2?

ki_co2 = 8 #co2 mass transfer coefficient (mol/m2.s)
ki_h2o = 16 #h2o Mass Transfer coefficient (mol/m2.s)

m_co2 = rho0CO2*dz*A_C (value of 5*10^-5)
m_h2o = rho0H2O*dz*A_C (value of 5*10^-5)

what are these?

M_al = 72.2 (mol/tank) (about 7kg/tank)

what is this?

#Boundary conditions
mol_in = 0.5 #mol/s
y_co2_in = 0.1 #mol/mol
y_h2o_in = 0.01 #mol/mol
T_in = 220 #K

 

[casualguitar]

Bed diameter = 8 cm?

Yes

thermal conductivity, not heat capacity

Typo

Units are meters? So total bed length is 0.05 m = 5 cm?

No the bed length is 2m. I've changed dz to 0.4m

cp of N2?

Yes I haven't included this yet I just put in an approximate constant average value of the cp of N2,CO2,H2O

So to the problematic ones:

What are these supposed to represent? I thought there is no water or CO2 deposited on the bed initially, and there is no water or CO2 in the gas phase initially.

Yes there is no deposited co2/water initially, and there is no co2/h2o present in the gas phase initially. I have set the deposited values to zero, but we're dividing by the molar holdup of CO2 and the molar holdup of h2o in 4/6 ODEs so I don't think I can let the molar holdup equal to zero. Therefore I've let it equal to the density of the vapour*cross sectional area*dz. This is actually wrong though also. The mass holdup is density of vapour * tank volume. Is it correct to say I can't have a molar holdup of zero ever because of this division issue?

Please provide the calculation which led to this.

I think I had (and possibly still do) have an error here also. I just changed the calculation in the last few minutes. I still get an error in the simulation but the results do look better.

Specific surface area:
$S_A = V_{particle}/V_{bed} = \frac{2*d_p^3}{3*d_b^2h}$
where $d_p$ is the particle diameter, $d_b$ is the bed diameter and $h$ is the length of the bed

what is this?

$M_{al}$ is the total moles of solid alumina per tank

EDIT: So just to confirm, since that specific area fix, the gas temperature in the bed trends towards the inlet temperature, the gas fractions of co2/h2o trend towards the inlet mole fractions now

 

[casualguitar]

A positive update-
Finally getting normal output for every variable (except bed temperature weirdly but I'll fix that). The time scale is obviously too short to reflect reality but here's an example plot of time vs gas temperature for each position (n=3):

Working on that gas temperature bug now (the bed temperature seems to stay at the initial temperature)

 

[Chestermiller]

Yes there is no deposited co2/water initially, and there is no co2/h2o present in the gas phase initially. I have set the deposited values to zero, but we're dividing by the molar holdup of CO2 and the molar holdup of h2o in 4/6 ODEs so I don't think I can let the molar holdup equal to zero. Therefore I've let it equal to the density of the vapour*cross sectional area*dz. This is actually wrong though also. The mass holdup is density of vapour * tank volume. Is it correct to say I can't have a molar holdup of zero ever because of this division issue?

I don't understand this at all. Where do you divide by the molar holdup in your equations? And why does this have to include water and co2?

I think I had (and possibly still do) have an error here also. I just changed the calculation in the last few minutes. I still get an error in the simulation but the results do look better.

Specific surface area:
$S_A = V_{particle}/V_{bed} = \frac{2*d_p^3}{3*d_b^2h}$
where $d_p$ is the particle diameter, $d_b$ is the bed diameter and $h$ is the length of the bed

I don't understand this at all. The amount of particle surface area per unit volume of bed is $\frac{6}{d_p}(1-\epsilon)$ (assuming spherical particles).

$M_{al}$ is the total moles of solid alumina per tank

EDIT: So just to confirm, since that specific area fix, the gas temperature in the bed trends towards the inlet temperature, the gas fractions of co2/h2o trend towards the inlet mole fractions now

I don't understand what this is saying.

 

[casualguitar]

I don't understand this at all. Where do you divide by the molar holdup in your equations? And why does this have to include water and co2?

The first term in the gas phase mole and heat balances $m_{m,j}$ is the molar holdup which is initially zero. When I divide across by this term (I do this when I'm setting up the ODEs in the solver), I'm left with a molar holdup in the denominator which is zero initially, unless I specify otherwise. No?

I don't understand this at all. The amount of particle surface area per unit volume of bed is 6dp(1−ϵ) (assuming spherical particles).

Yes sorry this is right I just quickly derived it by hand and got it wrong. Yes I agree its $\frac{6}{d_p}(1-\epsilon)$

The gas temperature should be high at the inlet and low downstream.

Yep that's also right (completely missed that). Interesting. Weird. I'll check this out

 

[casualguitar]

Yep that's also right (completely missed that). Interesting. Weird. I'll check this out

Actually no this is fine. It shows the temperature profile vs time for the three tanks. We'd expect them all to be at T0 at time=0, and the first tank will be the first to rise up to the inlet temperature, etc

However even though the trend is fine, for this simulation my inlet gas temperature is actually 220 so It shouldn't be levelling out at about 174. Bug somewhere

 

[Chestermiller]

The first term in the gas phase mole and heat balances $m_{m,j}$ is the molar holdup which is initially zero. When I divide across by this term (I do this when I'm setting up the ODEs in the solver), I'm left with a molar holdup in the denominator which is zero initially, unless I specify otherwise. No?
View attachment 302244

No. $$m_{m,j}=\frac{P}{RT_j}A\Delta z \epsilon$$This is the total number of moles of gas contained in tank j.

 

[casualguitar]

No. $$m_{m,j}=\frac{P}{RT_j}A\Delta z \epsilon$$This is the total number of moles of gas contained in tank j.

Ahh ok so this term isn't species specific? i.e. it is an overall term for the tank (meaning that there is only one $m_{m,j}$ value that applies to both gas phase mole balances and the gas phase heat balance)

 

[casualguitar]

Ahh ok so this term isn't species specific? i.e. it is an overall term for the tank (meaning that there is only one $m_{m,j}$ value that applies to both gas phase mole balances and the gas phase heat balance)

Debugging currently. The old specific surface area value I was using was 0.0008 m2/m3. The new/correct one from $\frac{6}{d_p}(1-\epsilon)$ is 840 m2/m3, assuming epsilon = 0.3 and dp = 0.005m.

The old value let's the integrator run for the specified rime range, and produces the expected plot shapes like:

which shows the CO2 'plug' moving down the bed, and:

which is the temperature distribution for given positions.

However when I use the actual value calculated from $\frac{6}{d_p}(1-\epsilon)$ the integrator does not run for the expected time range and produces plots like this:

I guess because the specific surface area value is now much bigger, there is some other term trending to 0 or inf somewhere. However there is no specific error returned. I'm debugging this now anyway, just updating

 

[casualguitar]

Just one further check to confirm the units $q*$? Its:

So mol/m2.s * J/mol = J/m2.s?

Working back through the computational flow to find the bug(s). It seems like the star term is quite high, which results in an absurdly high dTg/dT value. For context the qstar values are 1,000,000 initially and go down to 7000 or so

This suggests the $M_i''$ calculation isn't working. For context the $M_{co2}''$ values start at about 40 mol/m2.s and go down to about -8mol/m2.s (this trend seems right but the values may be very high)

Also note: the reason I know the dTg/dt term is way too high is because the mdot (mass flow out of a tank) is 10^5 or so after 1 iteration

Just another idea: it could also be the mass transfer coefficients. Currently 8 and 16 mol/m2.s (I took an average value from the function plot of re vs ki I did to leave these as constants for now). They seem very high

 

[casualguitar]

Hi Chet, as I mentioned I think I've narrowed the last bugs to how I'm setting up the initial/boundary conditions. I think my confusion is in where to use each of the initial, boundary, and actual values in setting up the boundary and non boundary equations.

Is it correct to say:
1) The boundary (inlet flow, inlet temperature, inlet composition) values are only used for the j-1 position values at n=0 ($m_dot_{j-1}$, $y_{j-1}$, $T_{j-1}$) and nowhere else
2) The boundary conditions are never used to evaluate any value (molar deposition/liquefaction rate, $q_{GI}$, $q_{IB}$, sublimation/liquefaction pressure)

3) The initial condition values (initial bed temperature, initial gas temperature, etc) are never used to evaluate any property value either
4) It is always 'the actual value' i.e. the solution of the relevant derivative at that point, that is used to evaluate anything
5) The initial temperature is the same as the actual temperature at n=0 and t=0

In short, the question is, for the equations below, when setting them up at the boundary (n=0), what values of $y_{j-1}$, $y_j$, temperature for the molar deposition term, temperatures for the heat flux terms should be used? Is it always the 'current' temperature, and never the initial/boundary temperature?

 

[casualguitar]

Is there a case for simplifying further to just have a single element? i.e. n=1?