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?

 

[casualguitar]

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

I've simplified to n=1. Some plots look as expected. The solid deposition of CO2 follows an n shaped curve for example. I've been debugging for effectively two weeks now (crippling but a good learning experience), and I almost have plots that look as expected. I have one question currently - what should the plots look like at t=0? i.e. should there be a jump initially up to a given value, or should this jump exist at all? If it shouldn't be there then this likely indicates that I'm setting up the boundary conditions incorrectly

For example the plot for qGI (gas interface heat flux) jumps from zero to 350 initially and then follows the expected trend. Should this just start at the max value?

Also the gas temperature has a little jump at the start for example, but the bed temperature does not exhibit this same jump

So my question is just that - is it expected that there will be a sharp jump at t=0 for gas related variables, or do these jumps suggest I've set up the boundary conditions incorrectly maybe?

 

[Chestermiller]

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

#\dot{m}_0# is the molar flow rate of the stream into the column

#y_{i,0}# is the mole fraction the mole fraction of the species entering the column

#T_0# is the temperature of the stream entering the column

These quantities do not change with time

2) The boundary conditions are never used to evaluate any value (molar deposition/liquefaction rate, $q_{GI}$, $q_{IB}$, sublimation/liquefaction pressure)

No, the entering conditions are not used to evaluate deposition rates and heat transfer rates in any of the tanks. However, they are present in the differential equations for the first tank.

3) The initial condition values (initial bed temperature, initial gas temperature, etc) are never used to evaluate any property value either

Of course not.

4) It is always 'the actual value' i.e. the solution of the relevant derivative at that point, that is used to evaluate anything

Sure.

5) The initial temperature is the same as the actual temperature at n=0 and t=0

What does n = 0 signify?

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?

Of course, assuming we can agree on what n = 0 represents.

View attachment 302799

 

[casualguitar]

Of course, assuming we can agree on what n = 0 represents.

n=0 represents the first position in the bed, and not the boundary 'position'

No, the entering conditions are not used to evaluate deposition rates and heat transfer rates in any of the tanks. However, they are present in the differential equations for the first tank.

Yes this clears it up and this is what I'm currently doing (or I think so at least).

Lastly, post 161 had a question about the initial 'jump' of certain values from zero up to a maximum value, rather than them starting at this value. Is this initial jump expected (such as in the $Q_{GI}$ plot), or does it indicate I've set something up incorrectly?

 

[Chestermiller]

n=0 represents the first position in the bed, and not the boundary 'position'

I still have no idea what this means.

With regard to post #161, are you sure about that number. I don't see any jumps in the plots of this post.

With regard to the Q's, you should substitute the expressions for these into the differential equations.

 

[casualguitar]

I still have no idea what this means.

With regard to post #161, are you sure about that number. I don't see any jumps in the plots of this post.

With regard to the Q's, you should substitute the expressions for these into the differential equations.

Sorry post #181! Bad eyesight. The $Q_{GI}$ plot

Of course, assuming we can agree on what n = 0 represents.

I might be missing some potential ambiguity over what n=0 means? Not the first position in the bed that the ODEs are evaluated at?

 

[Chestermiller]

Sorry post #181! Bad eyesight. The $Q_{GI}$ plotI might be missing some potential ambiguity over what n=0 means? Not the first position in the bed that the ODEs are evaluated at?

I'm still confused as to what you are saying n represents physically. I always though it was the total number of tanks being used in the calculation (equal to the length of the bed divided by the grid spacing).

In your graph of Q, which tank does this apply to?

 

[casualguitar]

I'm still confused as to what you are saying n represents physically. I always though it was the total number of tanks being used in the calculation (equal to the length of the bed divided by the grid spacing).

I agree yes it is this. Ah I think I see the confusion now. When I wrote n=0 earlier, I was referencing the zeroth tank i.e. the first tank in which the ODEs are being evaluated. I did not mean that there were zero tanks. Does this solve the confusion?

For the Q graph, this applies to the first tank (zeroth index). Is this expected?

 

[Chestermiller]

I agree yes it is this. Ah I think I see the confusion now. When I wrote n=0 earlier, I was referencing the zeroth tank i.e. the first tank in which the ODEs are being evaluated. I did not mean that there were zero tanks. Does this solve the confusion?

So it is tank 1?

For the Q graph, this applies to the first tank (zeroth index). Is this expected?

I don’t think so; initially the bed and gas are at the same temperature, right?

 

[casualguitar]

So it is tank 1?

Yes

I don’t think so; initially the bed and gas are at the same temperature, right?

Yes, hmm well if $q_{g,I}$ is this:

and $q*$ is this:

then at t=0, the second term of $q_{g,I}$ will be zero, and the first term of $q_{g,I}$ will be at a maximum. The first term is strictly decreasing as the temperature of the gas increases. And the second term increases slightly. So I suppose you'd expect $q_{g,I}$ would follow a kind of a/x curve shape?

Meaning that the initial jump is incorrect and the $q_{g,I}$ curve should intersect the y-axis at some non zero value ($q_{g,I}$ not equal to zero)? Sorry the reason I'm slightly fixated not this small thing is because it might have a knock on effect on other parameters

 

[casualguitar]

My apologies for the 'thinking out loud'. So zooming right into what happens to $q_{g,I}$ close to t=0, we see that it shoots up to the max value very quickly (then later it takes on roughly a 1/x shape). The equation for $q_{g,I}$ has two terms:

It is actually the second term that causes the huge increase, meaning that Tg increases faster than Tb at the beginning. This works out to be about a 3 degree delta T between Tg and Tb at the early stages of the simulation. Then gradually as the system reaches steady state at position zero, the bed and gas temperature differences approach zero (resulting in the 1/x shape). So I suppose this graph shape seems to make sense

 

[casualguitar]

Things seem to look a bit better now. However I'm completely puzzled as to why the gas temperature levels off at a temperature much lower than the inlet stream temperature. Could I post graphs of the intermediate variables, and we can comment on which graphs have the correct/incorrect shapes?

 

[Chestermiller]

Things seem to look a bit better now. However I'm completely puzzled as to why the gas temperature levels off at a temperature much lower than the inlet stream temperature. Could I post graphs of the intermediate variables, and we can comment on which graphs have the correct/incorrect shapes?

Sure. Also, why don't you try a test where the mass transfer is totally shut off (say by setting the mass transfer coefficient equal to zero for each of the two species). This should result in just a temperature change advancing through the bed.

 

[casualguitar]

Sure. Also, why don't you try a test where the mass transfer is totally shut off (say by setting the mass transfer coefficient equal to zero for each of the two species). This should result in just a temperature change advancing through the bed.

Ah ok interesting so I ran that model for the single tank and multi tank (3 tank) scripts. Here's the gas temperature output for both.

Single tank:

3 tank:

Notes:
1) The script will run for any time span. With the non zero mass transfer coefficients it 'broke' at about 600s
2) The temperature trends look generally as expected. The single tank script almost exactly matches the first tank in the multi tank script (not exactly but close)

The gas phase mole fractions for H2O and CO2:

Notes:
1) 0.1 and 0.01 are the inlet mole fractions of CO2 and h2o respectively, so this plot shape makes sense (with non constant heat transfer coefficients the CO2 fraction weirdly went up to 0.15)
2) The rise to the inlet conditions is almost immediate. I'm not sure if this is expected, or if this would take time. I guess an immediate jump seems reasonable though since we're effectively flushing the bed?

Do these plots seem reasonable to you?

This does narrow it down to the mass transfer coefficient related variables that are causing the script to break. And more specifically it might be the liquefaction/sublimation pressure equations as an initial guess

Also, do you know of another simplification I can make (besides mass transfer coefficient = 0) that would possibly uncover some useful information?

Edit: No, just thinking about the mole fraction graph some more, I don't think this graph makes sense. Hmm so upstream of the frozen plug we will see the inlet conditions in the gas phase (10% co2, 1% h2o), downstream of the plug we will see 0% co2/h2o. At the exit of the bed the trend will be that there will be 0% co2/h2o present in the gas phase, except for when the plugs reach the exit, there will be a temporary spike in the gas phase mole fractions (above 10% and 1%), then it will reduce to the inlet conditions again and remain here infinitely

Oh one other question - the summation term (the last term in the gas phase mole balances), is this term the sum of the molar desublimation/liquefaction rates of CO2 and H2O? Or say in the CO2 mole balance is it just the H2O liquefaction rate that goes in here i.e. we don't account for CO2 desublimating twice in the same equation?

 

[Chestermiller]

Ah ok interesting so I ran that model for the single tank and multi tank (3 tank) scripts. Here's the gas temperature output for both.

Single tank:
View attachment 303117

3 tank:
View attachment 303118
Notes:
1) The script will run for any time span. With the non zero mass transfer coefficients it 'broke' at about 600s
2) The temperature trends look generally as expected. The single tank script almost exactly matches the first tank in the multi tank script (not exactly but close)

The gas phase mole fractions for H2O and CO2:

View attachment 303120

View attachment 303121
Notes:
1) 0.1 and 0.01 are the inlet mole fractions of CO2 and h2o respectively, so this plot shape makes sense (with non constant heat transfer coefficients the CO2 fraction weirdly went up to 0.15)
2) The rise to the inlet conditions is almost immediate. I'm not sure if this is expected, or if this would take time. I guess an immediate jump seems reasonable though since we're effectively flushing the bed?

What do the mass fraction balance equations reduce to for the 1st tank with no mass transfer between phases? What is the approximate analytic solution to these equations, and does your numerical solution match these?

Oh one other question - the summation term (the last term in the gas phase mole balances), is this term the sum of the molar desublimation/liquefaction rates of CO2 and H2O?

Yes. Of course, with no mass transfer, these terms are not present.

Or say in the CO2 mole balance is it just the H2O liquefaction rate that goes in here i.e. we don't account for CO2 desublimating twice in the same equation?

No. Include both terms.

 

[casualguitar]

What do the mass fraction balance equations reduce to for the 1st tank with no mass transfer between phases? What is the approximate analytic solution to these equations, and does your numerical solution match these?

I think the approx analytic solution for the mole balance is:
$$y_{CO2} = y_{in}(1-e^{(\dot{m}t/m)})$$
where $y_{in}$ is the $CO_2$ inlet mole fraction which is $0.1$, $\dot{m}$ is the inlet molar flow which is $0.5mol/s$ and $m$ is the molar holdup

I made a comparison plot above of the CO2 mole fraction from the simulation versus the analytic CO2 mole fraction. They are the same (the blue line is right over the red one). So that confirms the fast rise in the first tankAnd great yes I was including both terms

So this reduces the error in the 'full' model (which is constant mass transfer coefficients for now) to the parameters which use the mass transfer coefficient

Is there a way I can slightly increase the model complexity in regards to the mass transfer coefficients, without fully implementing them? In an attempt to locate what isn't working

 

[Chestermiller]

I think the approx analytic solution for the mole balance is:
$$y_{CO2} = y_{in}(1-e^{(\dot{m}t/m)})$$
where $y_{in}$ is the $CO_2$ inlet mole fraction which is $0.1$, $\dot{m}$ is the inlet molar flow which is $0.5mol/s$ and $m$ is the molar holdup
View attachment 303177I made a comparison plot above of the CO2 mole fraction from the simulation versus the analytic CO2 mole fraction. They are the same (the blue line is right over the red one). So that confirms the fast rise in the first tankAnd great yes I was including both terms

So this reduces the error in the 'full' model (which is constant mass transfer coefficients for now) to the parameters which use the mass transfer coefficient

Is there a way I can slightly increase the model complexity in regards to the mass transfer coefficients, without fully implementing them? In an attempt to locate what isn't working

I have several ideas.

1. With the mass transfer turned off, show the temperature profiles of the gas and the bed as a function of position through the bed (tank number) at a selection of times, for several values of the heat transfer coefficient of the gas ( including very high values).

2. With the heat transfer completely suppressed, show the temperature profile of the gas as a function of position in the bed (tank number) at a selection of times, and and as a function of time for tank #1. Determine how the latter compares with the analytic solution. Do the latter comparison for times shorter than 1000 sec.

3. With the mass transfer turned on, use mass transfer coefficients 10X and 100X lower than you have been using.

4. Try to run the calculation for the conditions in the literature paper, and see how the results compare. Feel free to scale your values of heat transfer coefficient and mass transfer coefficient to come close to matching. Don't forget, they use an infinite heat transfer coefficient.

 

[casualguitar]

Is a value of n=3 suitable for these tests?

1. With the mass transfer turned off, show the temperature profiles of the gas and the bed as a function of position through the bed (tank number) at a selection of times, for several values of the heat transfer coefficient of the gas ( including very high values).

Can do. So here we'd be expecting to see that the 'sharp' temperature gradient happens earlier with with increased values of gas heat transfer coefficient (htc)

2. With the heat transfer completely suppressed, show the temperature profile of the gas as a function of position in the bed (tank number) at a selection of times, and and as a function of time for tank #1. Determine how the latter compares with the analytic solution. Do the latter comparison for times shorter than 1000 sec.

So we would have $U_g$ and $U_b$ = 0 here. Are we assuming the mass transfer coefficient is also turned off at this point? Also, that $q_{g,I}$ = 0. If so, the analytical solution will have the same format as the one from earlier for gas phase mole fraction

3. With the mass transfer turned on, use mass transfer coefficients 10X and 100X lower than you have been using.

Can do

4. Try to run the calculation for the conditions in the literature paper, and see how the results compare. Feel free to scale your values of heat transfer coefficient and mass transfer coefficient to come close to matching. Don't forget, they use an infinite heat transfer coefficient.

Understood

So in summary the tests are:
1) Tank number vs gas/bed temperature (for $k_i$ = 0 and a wide range of $U_g$ values)
2) Tank number vs gas temperature (for $k_i$,$U_g$,$U_b$ = 0)
3) Time vs gas temperature (numeric and analytical solution comparisons for t<1000)
4) Check output for $k_i$ 10 and 100 times smaller
5) Replicate Tuinier et al model

I'll start on these now

 

[casualguitar]

So there is a bug related to the heat transfer coefficients. I don't know what it is exactly yet though.

Position vs gas temperature for range of $U_g$ values:
Varying the gas and bed heat transfer coefficients from $10^{-4}$ up to $10^4$ has no effect on the position vs temperature (gas and bed) plots. The plot looks like this for all values of U:

Similarly when the heat transfer is suppressed (I did this by letting $U_g$ and $U_b$ equal to 0.000001, as there would be division by zero otherwise), the plot looks the same as above.

Time versus gas temperature for tank 1 (numeric vs analytic):
The analytic and numerical solutions match for HTCs set to zero. Also suggesting that adding in heat transfer terms has no effect (some error in the code)

Checking the output for $k_i$ 10 and 100 times smaller:
Ok slightly interestingly here I can get normal output for $10^{-8}$. The $10^{-6}$ value I was using before did not allow the simulation to run. However, I compare this output to the output where $k_i$ = 0 and it is the same! For both of these tests the heat transfer coefficients were a constant 100W/m2.KSo in summary, the mass and heat transfer coefficients currently have no effect on the gas or bed temperature (assuming the values used allow the simulation to finish)! I will look into this now however any thoughts appreciated as usual

 

[Chestermiller]

Is a value of n=3 suitable for these tests?

Please try more, say at least 10.

Can do. So here we'd be expecting to see that the 'sharp' temperature gradient happens earlier with with increased values of gas heat transfer coefficient (htc)

If you are talking about spatial position, increasing U should bring the bed and gas temperatures much closer together and make the temperature wave travel more slowly through the bed.

So we would have $U_g$ and $U_b$ = 0 here. Are we assuming the mass transfer coefficient is also turned off at this point?

Yes.

Also, that $q_{g,I}$ = 0. If so, the analytical solution will have the same format as the one from earlier for gas phase mole fraction

Yes. But please expand the scale so that we can see what is happening in more detail at short times.

Can doUnderstood

So in summary the tests are:
1) Tank number vs gas/bed temperature (for $k_i$ = 0 and a wide range of $U_g$ values)
2) Tank number vs gas temperature (for $k_i$,$U_g$,$U_b$ = 0)

ki =0

3) Time vs gas temperature (numeric and analytical solution comparisons for t<1000)

Yes.

4) Check output for $k_i$ 10 and 100 times smaller
5) Replicate Tuinier et al model

I'll start on these now

 

[Chestermiller]

So there is a bug related to the heat transfer coefficients. I don't know what it is exactly yet though.

Position vs gas temperature for range of $U_g$ values:
Varying the gas and bed heat transfer coefficients from $10^{-4}$ up to $10^4$ has no effect on the position vs temperature (gas and bed) plots. The plot looks like this for all values of U:
View attachment 303206
Similarly when the heat transfer is suppressed (I did this by letting $U_g$ and $U_b$ equal to 0.000001, as there would be division by zero otherwise), the plot looks the same as above.

Maybe Ug is not high enough. What is the bed temperature doing?

Time versus gas temperature for tank 1 (numeric vs analytic):
The analytic and numerical solutions match for HTCs set to zero. Also suggesting that adding in heat transfer terms has no effect (some error in the code)View attachment 303209

View attachment 303210

I can't see the comparison on this scale. Spread out the time axis.

Checking the output for $k_i$ 10 and 100 times smaller:
Ok slightly interestingly here I can get normal output for $10^{-8}$. The $10^{-6}$ value I was using before did not allow the simulation to run. However, I compare this output to the output where $k_i$ = 0 and it is the same! For both of these tests the heat transfer coefficients were a constant 100W/m2.KSo in summary, the mass and heat transfer coefficients currently have no effect on the gas or bed temperature (assuming the values used allow the simulation to finish)! I will look into this now however any thoughts appreciated as usual