Modelling of two phase flow in packed bed using conservation equations

[Chestermiller]

Ahh I see my apologies, yes what I did above is right except I forgot to replace ##T## with the derivative ##\frac{dT}{dt}## on the right. Should it not be (adding in m below the line. I guess you're leaving it out for me to fill it in):
$$C^*=\frac{[C_{PL}(T_{sat}-T_1)+\Delta h_{vap}+C_{PV}(T_2-T_{sat})]}{m(T_2-T_1)}$$

Not really. Between T1 and T2, we have $$m\frac{dh}{dt}=mC^*\frac{dT}{dt}=RHS$$where C* is as I wrote it in post #34. h is the enthalpy per unit mass.

So m is constant outside the phase change zone. Inside the 'phase change' zone as you said its:
$$m=\frac{V}{v_lX+\frac{RT_{sat}}{PM}(1-X)}$$

I think we want the bottom line (average specific volume) to linearly increase from ##v_l## to ##v_g## over the interval ##T1 \leq T \leq{T2}##, so we could just replace ##1-X## with ##\frac{T-T1}{T2-T1}##, to give:

$$m=\frac{V}{v_l(1-\frac{T-T1}{T2-T1})+\frac{RT}{PM}(\frac{T-T1}{T2-T1})}$$

For between T1 and T2, I get $$m=\frac{V}{\left[v_l+\left(\frac{RT_2}{PM}-v_l\right)\left(\frac{T-T_1}{T_2-T_1}\right)\right]}$$

 

[casualguitar]

Not really. Between T1 and T2, we have ##mdhdt=mC∗dTdt=RHSwhere## C* is as I wrote it in post #34. h is the enthalpy per unit mass.

Got it. I rewrote the derivative as above and it makes sense now

For between T1 and T2, I get: $$m=\frac{V}{\left[v_l+\left(\frac{RT_2}{PM}-v_l\right)\left(\frac{T-T_1}{T_2-T_1}\right)\right]}$$

Hmm interesting. Ok rearranging my equation, I get your one, except instead of ##\frac{RT_2}{PM}## I have ##\frac{RT}{PM}##. Why do you use ##T2## instead of ##T## for the gas specific volume?

My guess is that you're saying the gas specific volume is a constant (doesn't vary with temperature). Also you're treating the liquid specific volume as a constant also. Could we not use T instead of T2 here if we wanted to treat the gas specific volume as temperature dependent? We would also need to treat ##V_L## as constant in that case I would guess

I'm going to start now on the code for this using your equations

 

[casualguitar]

Not really. Between T1 and T2, we have $$m\frac{dh}{dt}=mC^*\frac{dT}{dt}=RHS$$where C* is as I wrote it in post #34. h is the enthalpy per unit mass.

For between T1 and T2, I get $$m=\frac{V}{\left[v_l+\left(\frac{RT_2}{PM}-v_l\right)\left(\frac{T-T_1}{T_2-T_1}\right)\right]}$$

So summing up the last number of posts, this seems to be the current phase change model.

The general heat balance to the fluid is:
$$mC_p\frac{dT}{dt}=\dot{m}_{in}C_p(T_{in}-T)-UA(T-T_S)$$

Where ##C_P## is:

##C_{PL}## for ##T \leq T1 ##
##C_{PV}## for ##T \geq T2 ##
##C*## for ##T1 \leq T \leq{T2}## where ##(C^*=\frac{[C_{PL}(T_{sat}-T_1)+\Delta h_{vap}+C_{PV}(T_2-T_{sat})]}{(T_2-T_1)}##

and where ##m## is:

##\frac{V}{V_L}## for ##T \leq T1 ##
##\frac{PM}{RT_2}V## for ##T \geq T2 ##
##\frac{V}{\left[v_l+\left(\frac{RT_2}{PM}-v_l\right)\left(\frac{T-T_1}{T_2-T_1}\right)\right]}## for ##T1 \leq T \leq{T2}##

The general heat balance to the solid remains the same for each temperature range:
$$MC_S\frac{dT_S}{dt}=UA(T-T_S)$$

So essentially this model is two coupled differential equations that are solved for a range of times, where the values of ##C_P## and ##m## are chosen based on the temperature.

How does this look?

 

[Chestermiller]

Got it. I rewrote the derivative as above and it makes sense now

Hmm interesting. Ok rearranging my equation, I get your one, except instead of ##\frac{RT_2}{PM}## I have ##\frac{RT}{PM}##. Why do you use ##T2## instead of ##T## for the gas specific volume?

In the interval between T1 and T2, I'm expressing the specific volume as a linear function of T, varying from ##v_l## at T1 to the gas specific volume at T2, namely ##\frac{RT_2}{pM}##. The linear dependence on T is controlled by ##\frac{T-T_1}{T_2-T_1}##

 

[casualguitar]

In the interval between T1 and T2, I'm expressing the specific volume as a linear function of T, varying from ##v_l## at T1 to the gas specific volume at T2, namely ##\frac{RT_2}{pM}##. The linear dependence on T is controlled by ##\frac{T-T_1}{T_2-T_1}##

Just letting you know I was sidetracked with some other project work for the last few days. I'm working on the post #38 now

 

[casualguitar]

In the interval between T1 and T2, I'm expressing the specific volume as a linear function of T, varying from ##v_l## at T1 to the gas specific volume at T2, namely ##\frac{RT_2}{pM}##. The linear dependence on T is controlled by ##\frac{T-T_1}{T_2-T_1}##

Hi Chet, phase change model code below with a graph of the fluid and solid temperatures over time. I could tweak the inputs until the graph looks nicer but in the interest of saving time I wont. Basically the code is exactly the same as the last model, except now we just take m and cp as variable rather than constant (as shown in the mass and cp functions).

So looking at the graph, initially the fluid temperature and the solid temperature jump towards one another, and then gradually both reduce to the gas inlet temperature. This happens slowly though, meaning that the mass flow in is small relative to the volume holdup. We also have the linear increase/decrease in temperature between T1 and T2 (shown by the straight line sections on the blue curve)

How does this look to you? What are your thoughts on adding in temperature variation along the direction of flow? Is this a reasonable development?

 

[Chestermiller]

Hi Chet, phase change model code below with a graph of the fluid and solid temperatures over time. I could tweak the inputs until the graph looks nicer but in the interest of saving time I wont. Basically the code is exactly the same as the last model, except now we just take m and cp as variable rather than constant (as shown in the mass and cp functions).

View attachment 292166
View attachment 292167
View attachment 292169
So looking at the graph, initially the fluid temperature and the solid temperature jump towards one another, and then gradually both reduce to the gas inlet temperature. This happens slowly though, meaning that the mass flow in is small relative to the volume holdup. We also have the linear increase/decrease in temperature between T1 and T2 (shown by the straight line sections on the blue curve)

How does this look to you? What are your thoughts on adding in temperature variation along the direction of flow? Is this a reasonable development?

This doesn't look anything like what I expected. I was expecting to see something like the previous graph. Wasn't the initial temperature of the bed supposed to be < 80 K, and the gas temperature supposed to be about 300 K? What did you use for T1 and T2, for the inlet gas temperature, and for the bed-filled-with-liquid temperature? Also, on the graph, please show times from 0 to 1000 s only.

 

[casualguitar]

This doesn't look anything like what I expected. I was expecting to see something like the previous graph. Wasn't the initial temperature of the bed supposed to be < 80 K, and the gas temperature supposed to be about 300 K? What did you use for T1 and T2, for the inlet gas temperature, and for the bed-filled-with-liquid temperature? Also, on the graph, please show times from 0 to 1000 s only.

Ah my apologies, the ##80K/300K## temperatures refer to the liquefaction of air model. I just chose water/steam for the model above to test this phase change model. I've switched it now so the initial temperature of the bed and fluid in the bed are both ##80K##, and the gas entering is ##300K##. I used ##T1 = 100K## and ##T2 = 110K## to approximate a boiling range. Range adjusted to ##1000s## also. I've also switched the input values (cp, vL, etc to reflect average air properties).

My apologies for the confusion there. What are your thoughts on adding in temperature variation along the direction of flow? Is this a reasonable development?

 

[Chestermiller]

Ah my apologies, the ##80K/300K## temperatures refer to the liquefaction of air model. I just chose water/steam for the model above to test this phase change model. I've switched it now so the initial temperature of the bed and fluid in the bed are both ##80K##, and the gas entering is ##300K##. I used ##T1 = 100K## and ##T2 = 110K## to approximate a boiling range. Range adjusted to ##1000s## also. I've also switched the input values (cp, vL, etc to reflect average air properties).
View attachment 292183
My apologies for the confusion there. What are your thoughts on adding in temperature variation along the direction of flow? Is this a reasonable development?

I assume you've run a range of test cases to put the model through its paces and get a feel for the quantitative effects of some of the parameters such as the heat transfer coefficient.

I'm hoping that you have developed an appreciation for the value of solving simple crude models of a system first before going to the more complicated versions. In this case, you have been able to develop a strategy for including the thermodynamics of the phase change in the calculations. Imagine having tried to do this with the full version of the model. You also have some results under your belt, so you already have a feel for the time scale of the behavior; the time for the system to be fully purged is going to be on the order of about 15 minutes. In my judgment, this has been of tremendous value. Please try to incorporate this kind of approach, in building from simple models to more complicated ones, as part of your routine strategy for the future.

We are going to need to have some discussions on expanding this to include the axial variation. I've been working on a finite difference scheme that I would be comfortable with for axial discretization; I'm close to finishing, and will document it for your consideration when I'm done. We also need to address dealing with numerical issues associated with advection/diffusion situations like this. Are you familiar with the concept of axial dispersion (over and above conduction ) in porous media and packed beds? Are you familiar with the "method of lines" for solving partial differential equations involving time and spatial position?

 

[casualguitar]

I assume you've run a range of test cases to put the model through its paces and get a feel for the quantitative effects of some of the parameters such as the heat transfer coefficient.

I'm hoping that you have developed an appreciation for the value of solving simple crude models of a system first before going to the more complicated versions. In this case, you have been able to develop a strategy for including the thermodynamics of the phase change in the calculations. Imagine having tried to do this with the full version of the model. You also have some results under your belt, so you already have a feel for the time scale of the behavior; the time for the system to be fully purged is going to be on the order of about 15 minutes. In my judgment, this has been of tremendous value. Please try to incorporate this kind of approach, in building from simple models to more complicated ones, as part of your routine strategy for the future.

Yes I've played with the input parameters to get an idea of how much of an effect it has on the output. I have got a more intuitive idea now as for how the system reacts to changes in input. Yes I would not have been able to do this with a full scale model, so really I suppose it can be faster to build simple models up gradually rather than go straight for the full version. I will incorporate this approach for future models!

We are going to need to have some discussions on expanding this to include the axial variation.

Great. I'm ready for this discussion now also

I've been working on a finite difference scheme that I would be comfortable with for axial discretization; I'm close to finishing, and will document it for your consideration when I'm done.

Thank you for this again, it is hugely appreciated. I won't ask too many questions yet but I thought we were avoiding manually implementing these FD schemes and instead using existing libraries to solve PDEs? I assume I'm misunderstanding what you mean by this

We also need to address dealing with numerical issues associated with advection/diffusion situations like this. Are you familiar with the concept of axial dispersion (over and above conduction ) in porous media and packed beds? Are you familiar with the "method of lines" for solving partial differential equations involving time and spatial position?

Not really. Axial dispersion refers to generally undesirable heat transfer both along the direction of flow and against the direction of flow I believe? For example over time the temperature gradient in a packed bed would reduce due to the temperature differences between the high and low temperature zones. I don't know why this would be different to conduction though. I assumed it would be a combination of conduction and convection.

I am not familiar with numerical issues associated with advection/diffusion, and no I'm not familiar with the method of lines either.

So in the meantime I will read up on the method of lines and axial dispersion and respond here with a new understanding

 

[Chestermiller]

To use the package codes to solve PDEs, you need to first discretize the equations spatially, but retain the time derivatives, which now become ordinary time derivatives. So you will be solving a set of coupled first order ordinary differential equations in time for the temperatures as the spatial grid points. This is called the method of lines. The package codes automatically solve these, and certain options guarantee stability of the time integration.

 

[casualguitar]

To use the package codes to solve PDEs, you need to first discretize the equations spatially, but retain the time derivatives

Understood. I did some reading on solving PDEs with the solve_ivp package I've used above and yes they generally said something similar to what you have said.

So you will be solving a set of coupled first order ordinary differential equations in time for the temperatures as the spatial grid points. This is called the method of lines.

I have seen solutions of PDEs similar to the ones I included in my original question that use methods like Crank Nicolson etc to solve for the temperature profiles at each time increment, where both the time and space domains are discretised. Why choose the method of lines approach over an approach like this?

 

[casualguitar]

To use the package codes to solve PDEs, you need to first discretize the equations spatially, but retain the time derivatives, which now become ordinary time derivatives. So you will be solving a set of coupled first order ordinary differential equations in time for the temperatures as the spatial grid points. This is called the method of lines. The package codes automatically solve these, and certain options guarantee stability of the time integration.

Hi Chet, just a quick update - I have now read up on the method of lines and have followed an example solution in MATLAB. To answer my own previous question regarding the use cases for the MOL, it seems the MOL would be chosen to take advantage of the methods and software that have been developed for numerically integrating ODEs

 

[Chestermiller]

You need to read up on axial heat dispersion in packed beds (over and above axial heat conduction). We need to discuss how to handle this.

 

[casualguitar]

You need to read up on axial heat dispersion in packed beds (over and above axial heat conduction). We need to discuss how to handle this.

Understood. I will read up on this now. Any aspects of porous media axial dispersion in particular? If not I will read up on how axial dispersion in packed beds is modeled. Thank you

 

[casualguitar]

I've spent the last few hours reading papers that include axial dispersion. I understand now why this isn't trivial. There are lots of correlations we could use to find a thermal axial dispersion coefficient (as a function of thermal conductivity, porosity etc).

This paper is quite similar to what we're doing, and it uses an axial dispersion coefficient correlation: https://www.sciencedirect.com/science/article/abs/pii/0038092X82902067

So essentially it seems to be 'lumped in' with the conduction term, by using a pseudoproperty ##K_Z## rather than the actual thermal conductivity.

Is this what you had in mind?

Also while reading I found this paper, which is sort of close to what we're doing: https://www.sciencedirect.com/science/article/pii/000925097480059X?via=ihub

I will take a close read of both of these tomorrow morning to prepare for this discussion

 

[Chestermiller]

I've spent the last few hours reading papers that include axial dispersion. I understand now why this isn't trivial. There are lots of correlations we could use to find a thermal axial dispersion coefficient (as a function of thermal conductivity, porosity etc).

This paper is quite similar to what we're doing, and it uses an axial dispersion coefficient correlation: https://www.sciencedirect.com/science/article/abs/pii/0038092X82902067

So essentially it seems to be 'lumped in' with the conduction term, by using a pseudoproperty ##K_Z## rather than the actual thermal conductivity.

Is this what you had in mind?

Yes. You may have to make this an adjustable parameter when you tune your model to the experimental data.

Also while reading I found this paper, which is sort of close to what we're doing: https://www.sciencedirect.com/science/article/pii/000925097480059X?via=ihub

I will take a close read of both of these tomorrow morning to prepare for this discussion

This paper doesn't seem to include a 2nd fluid phase.

 

[casualguitar]

Yes. You may have to make this an adjustable parameter when you tune your model to the experimental data.

Understood

This paper doesn't seem to include a 2nd fluid phase.

Yes you're absolutely right I linked the wrong one. This is the correct one: https://www.sciencedirect.com/science/article/pii/001793109090255S

This paper does assume axial dispersion is negligible so it won't be of use for that. It does model two phase flow though

 

[Chestermiller]

Understood

Yes you're absolutely right I linked the wrong one. This is the correct one: https://www.sciencedirect.com/science/article/pii/001793109090255S

This paper does assume axial dispersion is negligible so it won't be of use for that. It does model two phase flow though

I have access to only the abstract.

It looks like they are taking into account pressure variations, while we are assuming constant pressure.

In our analysis, a complexity in writing down the finite difference expressions for the spatial direction is that the density and mass rate of flow are both functions of time and spatial position; this needs to be accounted for in the mass balance equation and in the energy balance equation. This is what I've been working on.

 

[casualguitar]

I've set up a WeTransfer (allows sending of the PDF) link if you want to look at the paper: https://we.tl/t-9ebfBtGiuR

Yes correct they are. The library I am using can do temperature and pressure dependent parameters, this is an overview of the documentation if you want to take a look at the functionality: https://thermo.readthedocs.io

In our analysis, a complexity in writing down the finite difference expressions for the spatial direction is that the density and mass rate of flow are both functions of time and spatial position; this needs to be accounted for in the mass balance equation and in the energy balance equation. This is what I've been working on.

Ok so are you saying we will not just have two energy balance equations, but also a mass balance to the fluid equation, meaning our 'system' will involve solving three ODEs?

 

[Chestermiller]

I've set up a WeTransfer (allows sending of the PDF) link if you want to look at the paper: https://we.tl/t-9ebfBtGiuR

Yes correct they are. The library I am using can do temperature and pressure dependent parameters, this is an overview of the documentation if you want to take a look at the functionality: https://thermo.readthedocs.ioOk so are you saying we will not just have two energy balance equations, but also a mass balance to the fluid equation, meaning our 'system' will involve solving three ODEs?

Not exactly. I'll get to that soon.

But, before getting to the full problem, I'd like you to try something, please.

Try modifying the stirred tank version of the numerical model by solving the fluid heat balance equation ##m\frac{dh}{dt}=\dot{m}_{in}(h_{in}-h)-UA(T(h)-T_S)## (from post #12) in terms of the enthalpy h rather than the temperate T. (The heat balance for the solid can be left in terms of T and TS). The reason I'm asking this is that I think there would be an advantage to working in terms of h in the full model because the heat capacity is discontinuous, which would lead to logistical complexities in evaluating spatial derivatives numerically; such complexities could be avoided by working in terms of h (which is continuous). If this doesn't make a lot of sense, please take me at my word. Once the enthalpy is known, the temperature follows directly (in our scheme). Please make comparison plots so we can see how the results for the temperature dependence compare. Thanks.

 

[casualguitar]

Try modifying the stirred tank version of the numerical model by solving the fluid heat balance equation mdhdt=m˙in(hin−h)−UA(T(h)−TS) (from post #12) in terms of the enthalpy h rather than the temperate T. (The heat balance for the solid can be left in terms of T and TS).

Will do. Just some questions to clarify:
1) Am I to solve this for a single phase only, or for two phase flow?
2) The equation you supplied: ##m\frac{dh}{dt}=\dot{m}_{in}(h_{in}-h)-UA(T(h)-T_S)## seems to be in a useful form already. So if we are working on a single phase model then I actually only need an equation of state here I think, H(T) and T(H), and then we can solve using the ODE solver as usual

Is the above correct? i.e. your equation is in a useful form and I just need an EOS (thermodynamic property library) to solve?

 

[Chestermiller]

Will do. Just some questions to clarify:
1) Am I to solve this for a single phase only, or for two phase flow?

Two phase.

2) The equation you supplied: ##m\frac{dh}{dt}=\dot{m}_{in}(h_{in}-h)-UA(T(h)-T_S)## seems to be in a useful form already. So if we are working on a single phase model then I actually only need an equation of state here I think, H(T) and T(H), and then we can solve using the ODE solver as usual

Is the above correct? i.e. your equation is in a useful form and I just need an EOS (thermodynamic property library) to solve?

You use the property equations for M and h that we have written for the 2 phase model. It's easiest to express T and M as functions of h over the three intervals.

 

[casualguitar]

You use the property equations for M and H that we have written for the 2 phase model.

Are we assuming ##C_{PL}## and ##C_{PG}## are constant?

So the differences in terms of calculations here would be:
1) We are first calculating the inlet enthalpy using the previously developed property equations
2) The mass holdup will now be a function of enthalpy instead of temperature, meaning that we should solve for ##T_{sat}## in each of the three enthalpy property equations
3) At each time step we solve for the fluid ##H##, and then convert to ##H## to ##T## for input into the solid equation
4) Also we would need the enthalpy at ##T1## and ##T2## rather than the temperature bounds

I will start on this solution this evening however I'm not sure why this enthalpy model would output different results to the temperature model. Why would it?

I'll do the necessary equation adjusting and comment with a clear algorithm before I start coding

 

[Chestermiller]

Are we assuming ##C_{PL}## and ##C_{PG}## are constant?

So the differences in terms of calculations here would be:
1) We are first calculating the inlet enthalpy using the previously developed property equations
2) The mass holdup will now be a function of enthalpy instead of temperature, meaning that we should solve for ##T_{sat}## in each of the three enthalpy property equations
3) At each time step we solve for the fluid ##H##, and then convert to ##H## to ##T## for input into the solid equation
4) Also we would need the enthalpy at ##T1## and ##T2## rather than the temperature bounds

I will start on this solution this evening however I'm not sure why this enthalpy model would output different results to the temperature model. Why would it?

I'll do the necessary equation adjusting and comment with a clear algorithm before I start coding

OK. Let's do the interval thing with this approach.

1. If ##h\leq h_{sat,L}## , then $$T=T_{sat}+(h-h_{sat,L})/C_{PL}$$ and $$m=\rho_L V$$

2. If ##h_{sat,L}\leq h \leq h_{sat,V}##, then $$T=T_{sat}$$and $$m=\frac{V}{\left[\frac{1}{\rho_L}+\left(\frac{RT_{sat}}{PM}-\frac{1}{\rho_L}\right)\left(\frac{(h-h_{sat,L})}{(h_{sat,V}-h_{sat,L})}\right)\right]}$$

3. If ##h_{sat,V}\leq h##, then $$T=T_{sat}+(h-h_{sat,V})/C_{PV}$$ and $$m=\frac{PM}{RT} V$$ with T given by the previous equation.

Hmmm...This method gets us out of the artificial approach of using T1 and T2, and is "exact" for the relations between enthalpy, temperature. and mass.

We expect that the results of this calculation will match the results from the temperature-version, so it should provide a good validation of the model. Plus, it will show the effect of the temperature remaining exactly constant with time at the saturation temperature while the fluid is converting from vapor to liquid.

Note finally, that, with no loss of generality, we can simplify further by taking the reference temperature for zero enthalpy of the liquid as ##T_{sat}##, such that ##h_{sat,L}=0## and ##h_{sat,V}=\Delta h_{vaporization}##.

 

[casualguitar]

Note finally, that, with no loss of generality, we can simplify further by taking the reference temperature for zero enthalpy of the liquid as Tsat, such that hsat,L=0 and hsat,V=Δhvaporization.

Got it, so just adding that in:

1. If ##h\leq h_{sat,L}##, then
$$T=T_{sat}+h/C_{PL}$$
$$m=\rho_L V$$

2. If ##h_{sat,L}\leq h \leq h_{sat,V}## then
$$T=T_{sat}$$
$$m=\frac{V}{\left[\frac{1}{\rho_L}+\left(\frac{RT_{sat}}{PM}-\frac{1}{\rho_L}\right)\left(\frac{h}{\Delta h_{vap}}\right)\right]}$$

3. If ##h_{sat,V}\leq h## then
$$T=T_{sat}+(h-\Delta h_{vap})/C_{PV}$$
$$m=\frac{PM}{RT} V$$

We expect that the results of this calculation will match the results from the temperature-version, so it should provide a good validation of the model.

The post #38 model yes?

This method gets us out of the artificial approach of using T1 and T2, and is "exact" for the relations between enthalpy, temperature. and mass.

It does, I guess dealing with a mixture means that we should get an average value for ##T_{sat}## here (N2 boils at 77K, O2 boils at 90K) so maybe a Tsat of 73.5K is suitable here? The same for any other constant I could take the average of the N2 and O2 values.

I think I follow this model fully so, I will aim to implement this tomorrow morning!

 

[Chestermiller]

Got it, so just adding that in:

More like

1. If ##h\leq 0##, then
$$T=T_{sat}+h/C_{PL}$$
$$m=\rho_L V$$

2. If ##0\leq h \leq \Delta h_{vap}## then
$$T=T_{sat}$$
$$m=\frac{V}{\left[\frac{1}{\rho_L}+\left(\frac{RT_{sat}}{PM}-\frac{1}{\rho_L}\right)\left(\frac{h}{\Delta h_{vap}}\right)\right]}$$

3. If ##\Delta h_{vap}\leq h## then
$$T=T_{sat}+(h-\Delta h_{vap})/C_{PV}$$
$$m=\frac{PM}{RT} V$$

The post #38 model yes?

Yes.

It does, I guess dealing with a mixture means that we should get an average value for ##T_{sat}## here (N2 boils at 77K, O2 boils at 90K) so maybe a Tsat of 73.5K is suitable here? The same for any other constant I could take the average of the N2 and O2 values.

Please use your best judgment for this. I guess at some point the model can be modified to more the precise VLE behavior. I think the 73.5 is a typo, and you meant 83.5, right?

 

[Chestermiller]

I think I can solve for the case with U = 0 analytically.

 

[casualguitar]

Working on the above model now I should have results for this in an hour or so. I'll change the inputs of the model we're comparing to so that they're close to equivalent. What would a U = 0 model show us? The effect of having zero heat transfer between fluid and solid?

 

[Chestermiller]

Working on the above model now I should have results for this in an hour or so. I'll change the inputs of the model we're comparing to so that they're close to equivalent. What would a U = 0 model show us? The effect of having zero heat transfer between fluid and solid?

It would tell us the effect of purging liquid initially at a low temperature from the lumped mixed tank with a vapor at a high temperature.

 

[casualguitar]

It would tell us the effect of purging liquid initially at a low temperature from the lumped mixed tank with a vapor at a high temperature.

So I made the above model, and ran it for ranges of input parameters etc. I've printed the output and I can see the expected increase of temperature to the saturation temperature, a (short) constant temperature zone, and then an increase up to the fluid inlet temperature. I've checked that the correct enthalpy/mass equations are being used for the right enthalpy ranges. I've also made this models inputs as close as possible to the old model inputs. This one produces results that are quite a bit different, in that the temperature of the fluid and gas increase at a much slower rate. I'm attaching the output to model 1 and 2, and also the code to the new model below. Note I plotted over 10000s as this was the time required to see a levelling out of the temperature curve in the new model. I'm assuming for now there is a bug in my code somewhere so I'll hopefully find it this evening!

New model:

Old model:

Enthalpy model code:

 

[Chestermiller]

I don’t see anything like a constant temperature zone in either plot. What were the parameter values you used, including the mass flow rate in, and the solid mass and heat capacity?

 

[casualguitar]

I don’t see anything like a constant temperature zone in either plot. What were the parameter values you used, including the mass flow rate in, and the solid mass and heat capacity?

There is a constant temperature zone in the new plot, however it’s very small (I think it spans only a few seconds). I’m away from my laptop currently so I will zoom in on this in the plot later this evening. In regards to mass flow in, solid mass and heat capacity, these are minL, mS and cpS in the first code image. I’ll also change these parameters this evening to get a longer constant temperature zone

 

[Chestermiller]

There is a constant temperature zone in the new plot, however it’s very small (I think it spans only a few seconds). I’m away from my laptop currently so I will zoom in on this in the plot later this evening. In regards to mass flow in, solid mass and heat capacity, these are minL, mS and cpS in the first code image. I’ll also change these parameters this evening to get a longer constant temperature zone

You omitted multiplying by the bed volume when you calculated the mass of fluid in the bed for pure vapor. This is an indication that you really need to spend some time testing and checking your code.

What are the units of the physical parameters in your code?

Heat transfer coefficeint
Heat transfer area
Heat capacity of bed
mass of bed

 

[casualguitar]

You omitted multiplying by the bed volume when you calculated the mass of fluid in the bed for pure vapor.

Thanks!

This is an indication that you really need to spend some time testing and checking your code.

Yes I agree completely. I did spend over an hour looking, however this is really a bug I should have been able to find if I'd spent some more time. I'll take more care to find these myself from now on

What are the units of the physical parameters in your code?

Units included in the comments here:

Updated graph below. Shows the constant temperature zone and levels out at the same rate:

I'm going to make a graph of this one over the previous model graph to see how well they actually overlay. However from visual inspection of the individual graphs they look ok