Modelling of two phase flow in packed bed using conservation equations

[casualguitar]

For the bed

I don’t know which model has the error (or errors). This is part of doing the validation of the model.

Got it my bad working on this now.

One approach to do this could be to reduce to 3 tanks and do some manual calculations to see which model values don't match the manually calculated versions. Trying this now

 

[Chestermiller]

What were the units of that approximate time scale calculation?

 

[casualguitar]

What were the units of that approximate time scale calculation?

$\frac{kg*\frac{kJ}{kgK}}{\frac{W}{m2K}*m2}$ = $\frac{kJ}{W}$ = $1000*s$ == kilo seconds I suppose?

Changing the heat capacity units to $\frac{J}{kgK}$ gives a value of 11,000s for the time scale of the response

 

[Chestermiller]

$\frac{kg*\frac{kJ}{kgK}}{\frac{W}{m2K}*m2}$ = $\frac{kJ}{W}$ = $1000*s$ == kilo seconds I suppose?

Changing the heat capacity units to $\frac{J}{kgK}$ gives a value of 11,000s for the time scale of the response

You got to be kidding. Changing the units of heat capacity doesn’t change the value of this quantity. The result is still 1100 sec.

 

[casualguitar]

You got to be kidding. Changing the units of heat capacity doesn’t change the value of this quantity. The result is still 1100 sec.

Ah no I agree it doesn't change the value, just the unit of time i.e. from kilo seconds to seconds. My original value of 11 was in kilo seconds (according to the above). Multiplying by 1000 to convert to seconds. Is this incorrect?

Apologies for the confusion here. I think it might be best for me to spend some time looking at the output of the spatial variation model (and doing some manual calculations to check the model output). Otherwise I'm wasting your time and mine

 

[Chestermiller]

1100 is correct

 

[casualguitar]

1100 is correct

So 1100s is the approx time it would take for the solid temperature to reach steady state i.e. the timescale. The solid temperature will converge exponentially towards 300 but will never reach 300. Where would the cut-off point be here? Say within 0.1K or so? I can check the exact time required for the cutoff to be reached.

 

[Chestermiller]

So 1100s is the approx time it would take for the solid temperature to reach steady state i.e. the timescale. The solid temperature will converge exponentially towards 300 but will never reach 300. Where would the cut-off point be here? Say within 0.1K or so? I can check the exact time required for the cutoff to be reached.

We have got to take a step back, and make sure that the single tank model of the entire bed is getting the right answer. This can be checked by running a calculation in which the model can be solved analytically. Such a case would be if the stream entering the tank is a liquid at a temperature below the saturation temperature, and the bed is initially at a colder temperature. Then we can compare the analytic solution to the time-dependent numerical solution to make sure that the time scale we get numerically is correct.

After that, I have developed an even simpler version of the multi-tank model that is easier to implement, and less prone to coding error. I have shown that the multi-tank model is equivalent to n tanks in series identical in every respect to that being used in the single tank model of the entire bed, but with a mass flow to the first tank in the series equal to n times the actual mass flow rate to the actual bed. I will document the new development once we have validated the single tank model.

 

[casualguitar]

We have got to take a step back, and make sure that the single tank model of the entire bed is getting the right answer. This can be checked by running a calculation in which the model can be solved analytically. Such a case would be if the stream entering the tank is a liquid at a temperature below the saturation temperature, and the bed is initially at a colder temperature. Then we can compare the analytic solution to the time-dependent numerical solution to make sure that the time scale we get numerically is correct.

Great! I follow the first paragraph idea. So focusing on the first paragraph - we will have $h<=0$ at all times, meaning we can use the liquid correlations for temperature and mass holdup i.e. $T = T_{SAT} + \frac{h}{C_{PL}}$ and $m = \rho_{L}*V$.

Letting Tsat = 105, Tin = 90, and T0 = 80 (bed/fluid initial temperature), for the single tank model, this is the plot for the liquid phase only model:

If we are expecting a timescale of 1100s here this seems to line up relatively ok.

I did the same for gas phase only - where Tsat = 105, Tin = 300 and T0 = 200 (bed/fluid initial temperature), for the single tank model. Gas phase plot looks like this:

Seems to look off, due to that fluid temperature jump at the start. Will do some diagnostic testing on the gas phase section

 

[Chestermiller]

Great! I follow the first paragraph idea. So focusing on the first paragraph - we will have $h<=0$ at all times, meaning we can use the liquid correlations for temperature and mass holdup i.e. $T = T_{SAT} + \frac{h}{C_{PL}}$ and $m = \rho_{L}*V$.

Letting Tsat = 105, Tin = 90, and T0 = 80 (bed/fluid initial temperature), for the single tank model, this is the plot for the liquid phase only model:View attachment 293726
If we are expecting a timescale of 1100s here this seems to line up relatively ok.

I did the same for gas phase only - where Tsat = 105, Tin = 300 and T0 = 200 (bed/fluid initial temperature), for the single tank model. Gas phase plot looks like this:
View attachment 293728
Seems to look off, due to that fluid temperature jump at the start. Will do some diagnostic testing on the gas phase section

Do you know how to solve this analytically (two linear coupled ODEs in two unknowns)?

 

[casualguitar]

Do you know how to solve this analytically (two linear coupled ODEs in two unknowns)?

I don't. I can figure this out though. So to confirm, the two equations are the solid and energy balance here:

Where $h$ and $T_s$ are the two unknowns, and where the $h<=0$ correlations will be used for $T$ and $m$?

 

[Chestermiller]

I don't. I can figure this out though. So to confirm, the two equations are the solid and energy balance here:
View attachment 293731

Where $h$ and $T_s$ are the two unknowns, and where the $h<=0$ correlations will be used for $T$ and $m$?

Yes, the liquid region, and, for the analytic solution, you can use the constant heat capacity version of the equations. You are taking V as the pore volume, right? What values of the parameters are you using:

Mass of liquid in bed
Mass of solid bed
Mass flow rate of liquid to bed
U, A,
Heat capacity of solid
Heat capacity of liquid

I've done part of the analysis (but I want you to also do it), and I've got equations for the time constant(s) in terms of the above parameters. This should tell us much more about the long time response.

 

[casualguitar]

Yes, the liquid region, and, for the analytic solution, you can use the constant heat capacity version of the equations. You are taking V as the pore volume, right? What values of the parameters are you using:

Mass of liquid in bed
Mass of solid bed
Mass flow rate of liquid to bed
U, A,
Heat capacity of solid
Heat capacity of liquid

I've done part of the analysis (but I want you to also do it), and I've got equations for the time constant(s) in terms of the above parameters. This should tell us much more about the long time response.

Understood, so these equations, for the liquid phase i.e. T<Tsat:

Yes V is the pore volume for this model

Parameter Values:
Mass of liquid in bed = 8kg
Mass of solid in bed = 1kg
Mass flow rate of liquid to bed = 0.05kg/s
U = 1 W/m2K
A = 0.1m2
Heat capacity of solid = 1.1 kJ/kg.K
Heat capacity of liquid = 2 kJ/kg.K

In regards to the analytical solution, what makes up an analytical solution in this case?

Are we talking about solving for the eigenvalues, or do we mean solving the linear system, i.e. to get h(t) and Ts(t)?

The general approach I was planning to take here is:
1) Differentiate one of the energy balances w.r.t. time to get the second derivative
2) Substitute the other equation in
3) Solve the first equation for h (or Ts)
4) Substitute h (or Ts) back into the second derivative equation
5) We'll get the characteristic equation from this as we'll have the second and first derivative of h, and h
6) then assuming $\lambda_1$ and $\lambda_2$ are the roots

$$h(t) = c_{1}e^{\lambda_1t} + c_{2}e^{\lambda_2t}$$
$$Ts(t) = c_{3}e^{\lambda_1t} + c_{4}e^{\lambda_2t}$$

Where c1-c4 come from h(0), Ts(0), h'(0) and Ts'(0). If this is the route we're going then one question I would have is what would the values of h'(0) and Ts'(0) be?

 

[Chestermiller]

I’m having a medical issue presently, and am not going to be able to participate as much as I’d like to for a while. It is very frustrating because I have lots that I’d like to cover with you.

For the analytic solution, I would work with T and TS, and work with a matrix formulation, starting with subtracting out the final steady state temperatures equal to Tin.

 

[casualguitar]

I’m having a medical issue presently, and am not going to be able to participate as much as I’d like to for a while. It is very frustrating because I have lots that I’d like to cover with you.

Absolutely, take the time needed to recover from this, please! Sorry to hear this. Forget about me I'll be fine!

I can and will do this analytic solution using a matrix formulation. I guess the next step would be to check if it matches the model and if not figure out why it does not match. Plenty to work on!

Thanks Chet for everything so far. It is hugely appreciated. And take your time with the recovery!

 

[Chestermiller]

Understood, so these equations, for the liquid phase i.e. T<Tsat:
View attachment 293748
Yes V is the pore volume for this model

Parameter Values:
Mass of liquid in bed = 8kg
Mass of solid in bed = 1kg
Mass flow rate of liquid to bed = 0.05kg/s

How come the mass of solid packing is less than the mass of liquid in the bed? Isn't the density of the solid higher, and isn't there more volume fraction packing?

And, if I recall correctly, previously, wasn't the mass flow rate of fluid.

 

[casualguitar]

How come the mass of solid packing is less than the mass of liquid in the bed? Isn't the density of the solid higher, and isn't there more volume fraction packing?

And, if I recall correctly, previously, wasn't the mass flow rate of fluid.

Hi Chet, hope your medical issue is going well. As I say absolutely zero pressure to respond here! I can work away

The solid mass is less than the liquid mass because I suppose I only needed to define the pore volume in the code, not the overall volume of the tank. The total tank volume is not defined. I can do this now though to get a more reasonable value for solid mass

Taking a pore fraction of say 0.3, the tank total volume would be $0.01*10/3 = 0.033m^3$, where $0.01m^3$ is the pore volume.

Liquid mass = $0.033*0.3*800 = 8kg$, where $800kg/m^3$ is the liquid density
Solid mass = $0.033*0.7*2500 = 57.75kg$, where $2500kg/m^2$ is the solid density

I'll sub in the new solid mass value now

And, if I recall correctly, previously, wasn't the mass flow rate of fluid.

I assume this is a typo, and you're asking about previous mass flow rate values. I did use $0.01kg/s$ in the original $dT/dt$ model. This value was changed to $0.05kg/s$ for no particular reason (I was likely changing the flow to get a nicer plot), and is now $0.05kg/s$ in all models

 

[Chestermiller]

In the present version of the multi-tank model, the balance equations read: $$\dot{m}_{i}=\dot{m}_{i-1}-\frac{V}{n}\frac{d\rho_i}{dt}$$$$\frac{V}{n}\rho_i\frac{dh_i}{dt}=\dot{m}_{i-1}(h_{i-1}-h_i)+U\frac{A}{n}(T_{S,i}-T_i)$$and$$\frac{m_s}{n}\frac{dT_{S,i}}{dt}=-U\frac{A}{n}(T_{S,i}-T_i)$$where V is the total pore volume in the bed and A is the total heat transfer area in the bed.
If we multiply all three of these equations by the number of tanks n, we obtain:$$\dot{M}_{i}=\dot{M}_{i-1}-V\frac{d\rho_i}{dt}$$$$(V\rho_i)\frac{dh_i}{dt}=\dot{M}_{i-1}(h_{i-1}-h_i)+UA(T_{S,i}-T_i)$$and$$m_s\frac{dT_{S,i}}{dt}=-UA(T_{S,i}-T_i)$$where $$\dot{M}_i=n\dot{m}_i$$These equations tell us that the multi-tank model is equivalent mathematically to a tanks-in-series-model in which each tank in the series is identical to the single tank model of the entire packed bed, but with the mass flow rate to the first tank in the series equal to n times the actual mass flow rate to the bed. This form of the model is much easier to code than the present version of the model because much of the code can be extracted from the single tank model (plus we don't have to divide everything by n), less error prone, and it is much easier to conceptualize. What do you think?

 

[casualguitar]

These equations tell us that the multi-tank model is equivalent mathematically to a tanks-in-series-model in which each tank in the series is identical to the single tank model of the entire packed bed, but with the mass flow rate to the first tank in the series equal to n times the actual mass flow rate to the bed. This form of the model is much easier to code than the present version of the model because much of the code can be extracted from the single tank model (plus we don't have to divide everything by n), less error prone, and it is much easier to conceptualize. What do you think?

I follow the reasoning almost completely except the mass flow point here is slightly confusing. So you say the multi-tank model is equivalent mathematically to a tanks in series model. I follow this. It means we can model each individual tank in the multi-tank model the same way the single tank is modeled. So essentially we'll have n single tank models in series, and this is equivalent to the multi-tank model.

but with the mass flow rate to the first tank in the series equal to n times the actual mass flow rate to the bed

This bit though - so we'll have the mass flow to the first tank actually equal to say 0.05kg/s, however we want to multiply this value by n. So capital $M$ is not really the mass flow, its a kind of total flow through the system value. So we will be solving for capital $M$ at all spatial points in the system, and to find the mass flow array we would just need to divide by n? This would make sense to me. Does this M term have physical significance?

I can code this model up. For this model as we're extracting code from the single tank model, I think it makes sense to finish the analytical solution for the liquid phase single tank model first. I'm actually working on the analytical solution currently. Then once its clear the numerical solution to the single tank checks out, I can use this code for the multi-tank modeled you've detailed above. How does this approach sound?

 

[Chestermiller]

For the trick model, the ratio of the mass flow rate to the the total mass of packing or of fluid in the sequence of tanks must be preserved.

 

[casualguitar]

For the trick model, the ratio of the mass flow rate to the the total mass of packing or of fluid in the sequence of tanks must be preserved.

Understood. I fully follow. If there is an error in the single tank code this will probably happen in the multi-tank model also, so do you think its worth finishing out the analytical solution model/code first?

 

[Chestermiller]

Understood. I fully follow. If there is an error in the single tank code this will probably happen in the multi-tank model also, so do you think its worth finishing out the analytical solution model/code first?

Definitely.

 

[casualguitar]

Definitely.

On it. Thanks Chet

In regards to output from the analytical solution, I guess we are we looking for:
- A graph of fluid enthalpy vs time
- A graph of solid temperature versus time

To check if the timescale for the analytical solution matches that of the numerical one, for the liquid phase

 

[Chestermiller]

On it. Thanks Chet

In regards to output from the analytical solution, I guess we are we looking for:
- A graph of fluid enthalpy vs time
- A graph of solid temperature versus time

To check if the timescale for the analytical solution matches that of the numerical one, for the liquid phase

For outside presentation of the validation checks, I think it would be better to compare temperatures vs time.

 

[casualguitar]

For outside presentation of the validation checks, I think it would be better to compare temperatures vs time.

Hi Chet, update on the analytical solution:

I haven't been able to actually derive the analytical solution myself yet (seems more difficult than I first anticipated), however Mathematica does this and has solved the equations both numerically and analytically. I'll transfer the analytical solution functions (if correct) to Python. Here is the output for fluid temperature versus time (I have the solid plot also, its basically the same):

You can see the constant values I used, plus the T(t) and Ts(t) analytical solution functions output by Mathematica (the functions are blurry so I've typed them out at the bottom), and also the analytical solution graph from 0 to 3000s. Note it seems to curve off unexpectedly just before 3000s.

Here is the numerical solution, also done by Mathematica. The plot starts out similar to the analytical one, and continues to level off towards the inlet temperature indefinitely, which seems to make more physical sense. There is no unusual decrease in temperature behaviour happening in the numerical solution

Its late here so I'm going to check these plots against my own code first thing tomorrow.

Edit: The screenshooting doesn't seem to want to cooperate with the Mathematica page so I'm copying the analytical solution below. It doesn't seem to want to wrap at the end of the line (its two equations in the same line), so I will edit this tomorrow

$${{T(t) = 90 -4.37526 E^{(-0.0133365 t)}-5.62474 E^{(-0.000737726 t)}-4.44089*10^{-16} E^{(0.0125987 t)}]
and
Ts(t) = 90 +0.585555 E^{(-0.0133365 t)}
+8.88178*10{^-16} E^{(-0.0125987 t)}-10.5856 E^{(-0.000737726 t)}
-8.88178*10^{-16} E^{(0.0125987 t)}]}}$$

 

[Chestermiller]

Hi Chet, update on the analytical solution:

I haven't been able to actually derive the analytical solution myself yet (seems more difficult than I first anticipated), however Mathematica does this and has solved the equations both numerically and analytically. I'll transfer the analytical solution functions (if correct) to Python. Here is the output for fluid temperature versus time (I have the solid plot also, its basically the same):View attachment 293819

You can see the constant values I used, plus the T(t) and Ts(t) analytical solution functions output by Mathematica (the functions are blurry so I've typed them out at the bottom), and also the analytical solution graph from 0 to 3000s. Note it seems to curve off unexpectedly just before 3000s.

Here is the numerical solution, also done by Mathematica. The plot starts out similar to the analytical one, and continues to level off towards the inlet temperature indefinitely, which seems to make more physical sense. There is no unusual decrease in temperature behaviour happening in the numerical solution
View attachment 293820

Its late here so I'm going to check these plots against my own code first thing tomorrow.

Edit: The screenshooting doesn't seem to want to cooperate with the Mathematica page so I'm copying the analytical solution below. It doesn't seem to want to wrap at the end of the line (its two equations in the same line), so I will edit this tomorrow

$${{T(t) = 90 -4.37526 E^{(-0.0133365 t)}-5.62474 E^{(-0.000737726 t)}-4.44089*10^{-16} E^{(0.0125987 t)}]
and
Ts(t) = 90 +0.585555 E^{(-0.0133365 t)}
+8.88178*10{^-16} E^{(-0.0125987 t)}-10.5856 E^{(-0.000737726 t)}
-8.88178*10^{-16} E^{(0.0125987 t)}]}}$$

Their analytic solution is very puzzling because there are only two eigenvalues to the characteristic equation. Here is the eigenvalue equation I came up with: $$\lambda^2+\left[\frac{\dot{m}_{in}}{m}+\frac{UA}{mC}+\frac{UA}{m_sC_s}\right]\lambda+\frac{\dot{m}_{in}}{m}\frac{UA}{m_sC_s}=0$$
For the parameter values you showed in your printout, this gives $$\lambda=-0.0133$$and$$\lambda=-0.000738$$which matches their values (at least on the two main terms). This implies a time constant at long times of 1/0.000738=1355 sec.

Let's see how your numerical solution compares to their analytic solution..

 

[casualguitar]

Their analytic solution is very puzzling because there are only two eigenvalues to the characteristic equation. Here is the eigenvalue equation I came up with: $$\lambda^2+\left[\frac{\dot{m}_{in}}{m}+\frac{UA}{mC}+\frac{UA}{m_sC_s}\right]\lambda+\frac{\dot{m}_{in}}{m}\frac{UA}{m_sC_s}=0$$
For the parameter values you showed in your printout, this gives $$\lambda=-0.0133$$and$$\lambda=-0.000738$$which matches their values (at least on the two main terms). This implies a time constant at long times of 1/0.000738=1355 sec.

Let's see how your numerical solution compares to their analytic solution..

The analytical and numerical solutions from Mathematica now match. The problem was in how Mathematica deals with decimal numbers. So the graphs of T and Ts vs time now match.

And here is the analytical solution:

Editing my numerical solution now to compare.

There we go, my numerical solution does seem to match the Mathematica one:

And here is the plot for the liquid phase temperature, using the spatial variation model, assuming n=2 (so the inlet conditions plus one spatial point). Ignoring the jumps, This also matches the other plots:

Note the time is x1000 (due to 1000 points calculated per time step) so dividing the x-axis values by 1000 gives the same values as the previous plots

And lastly, I've set n = 32 and plotted the final position (index 31) versus time:

This also roughly matches the timescale of the other plots. Is this all reasonable to you?

 

[Chestermiller]

The analytical and numerical solutions from Mathematica now match. The problem was in how Mathematica deals with decimal numbers. So the graphs of T and Ts vs time now match.

And here is the analytical solution:
View attachment 293840

Editing my numerical solution now to compare.

There we go, my numerical solution does seem to match the Mathematica one:
View attachment 293843

And here is the plot for the liquid phase temperature, using the spatial variation plot, assuming n=2 (so the inlet conditions plus one spatial point). Ignoring the jumps, This also matches the other plots:
View attachment 293844
Note the time is x1000 (due to 1000 points calculated per time step) so dividing the x-axis values by 1000 gives the same values as the previous plots

And lastly, I've set n = 32 and plotted the final position (index 31) versus time:
View attachment 293845
This also roughly matches the timescale of the other plots. Is this all reasonable to you?

How about running the 3 tank version and plotting temperature vs time on the same graph for all three tanks? Don’t you expect a time delay between subsequent tank time responses as the wave travels through the bed?

 

[casualguitar]

Don’t you expect a time delay between subsequent tank time responses as the wave travels through the bed?

Yes there seems to be a slight delay looking at the n=32 plot

How about running the 3 tank version and plotting temperature vs time on the same graph for all three tanks?

Doing this now

 

[casualguitar]

How about running the 3 tank version and plotting temperature vs time on the same graph for all three tanks? Don’t you expect a time delay between subsequent tank time responses as the wave travels through the bed?

Here we go, this is for n=4. I went with n=4 so we can see three inner-tank Tf vs t graphs, position 0 is just the inlet condition.

As you mentioned we can see the time delay across the tank. The plot is quite messy due to these 'numerical bumps'. The output does look similar to the analytical and numerical output from Mathematica though.

Reverting back to the two phase model (assuming liquid in the tank, with a gas 300K input). Heres the graph for n=4. This also seems to approximately matches the timescale of the others:

Am I correct in saying this means that the spatial variation model is correct? i.e. no particular need to move to the model you outlined yesterday? I will code this up anyway but it seems the original spatial variation model is ok?

 

[Chestermiller]

What are the actual times on the time axis? Is 1 the same as 1000 seconds? If so, can you please show the results only up to the 1000 seconds part (so 1000 is the full horizontal range)? Also, can you please include the single tank model results for the entire bed for comparison? All this for the liquid only case. Thanks.

 

[casualguitar]

What are the actual times on the time axis? Is 1 the same as 1000 seconds?

Yes correct

If so, can you please show the results only up to the 1000 seconds part (so 1000 is the full horizontal range)?

Can do, this is both graphs above, to 1000s, for the liquid phase only (90K inlet stream, bed initial temperature of 80K).

Single tank:

Spatial variation plot: for n=4:

The single tank plot seems to best match the first position here.

Also, can you please include the single tank model results for the entire bed for comparison? All this for the liquid only case. Thanks.

To clarify by this do you mean plot both the fluid and solid temperature for the single tank, for the liquid phase?

Note: Its just the spatial variation plot that is time*1000 for the x axis

 

[Chestermiller]

Yes correct

Can do, this is both graphs above, to 1000s, for the liquid phase only (90K inlet stream, bed initial temperature of 80K).

Single tank:
View attachment 293854

Spatial variation plot: for n=4:
View attachment 293855
The single tank plot seems to best match the first position here.To clarify by this do you mean plot both the fluid and solid temperature for the single tank, for the liquid phase?

Note: Its just the spatial variation plot that is time*1000 for the x axis

Any chance of showing the single tank model results for the entire bed on the same plot as the multi-tank model results for this same case? Also show the single tank model for the entire bed with 3X the mass flow rate in?

 

[casualguitar]

Any chance of showing the single tank model results for the entire bed on the same plot as the multi-tank model results for this same case?

I think so, I'm not exactly sure how to do that yet but yes

Also show the single tank model for the entire bed with 3X the mass flow rate in?

I can do this now. When you say 'entire bed' what do you mean here? As there is only one tank in the single tank model. Or do you mean plot each position for the n=4 model? For 3x flowrate

 

[Chestermiller]

I think so, I'm not exactly sure how to do that yet but yes

I can do this now. When you say 'entire bed' what do you mean here? As there is only one tank in the single tank model. Or do you mean plot each position for the n=4 model? For 3x flowrate

The former.

 

[casualguitar]

The former.

Here we go, this is a 3x flowrate plot for the single tank model. So the flow here is 0.15kg/s. Everything else is the same:

I will look into this overlay plot now. There is almost definitely a convenient way to do this

 

[Chestermiller]

Here we go, this is a 3x flowrate plot for the single tank model. So the flow here is 0.15kg/s. Everything else is the same:

View attachment 293867
I will look into this overlay plot now. There is almost definitely a convenient way to do this

This is what the first tank in the 3 tank model should give. Does it?

 

[casualguitar]

This is what the first tank in the 3 tank model should give. Does it?

Hmm, this is the Mathematica analytical solution for 3x flowrate (zoomed to appropriate section because its fairly blurry otherwise):

Taking a reference point of 1000s, the temperature is just above 89, say 89.1K

Here is my numerical solution for the 3x flowrate single tank model (from above):

At t = 1000s, the temperature seems to be around 89K also, so its similar to the plot above. I think we can now say the Mathematica plot, and my single tank numerical solution are both correct

Multi-tank plots:
Here is the multi-tank plot for n=3, and m = 0.05kg/s:

At t = 1000s, the temperature is 88K, so its 1K lower than the single tank plots.

Am I correct to use n = 3 here?

Here is the plot for n= 4:

So t = 1000 and T = 88.3K or so. about 0.7K off.

I don't know if these are 'close enough' or not, but my guess is that they are not close enough. I'll run some diagnostics now in case it is not close enough

 

[Chestermiller]

Hmm, this is the Mathematica analytical solution for 3x flowrate (zoomed to appropriate section because its fairly blurry otherwise):
View attachment 293875
Taking a reference point of 1000s, the temperature is just above 89, say 89.1K

Here is my numerical solution for the 3x flowrate single tank model (from above):
View attachment 293881
At t = 1000s, the temperature seems to be around 89K also, so its similar to the plot above. I think we can now say the Mathematica plot, and my single tank numerical solution are both correct

Multi-tank plots:
Here is the multi-tank plot for n=3, and m = 0.05kg/s:
View attachment 293880
At t = 1000s, the temperature is 88K, so its 1K lower than the single tank plots.

Am I correct to use n = 3 here?

Here is the plot for n= 4:
View attachment 293879
So t = 1000 and T = 88.3K or so. about 0.7K off.

I don't know if these are 'close enough' or not, but my guess is that they are not close enough. I'll run some diagnostics now in case it is not close enough

You need to compare the temperature vs time of the first tank in your multi-tank 3 tank model with the temperature vs time for the single tank model of the full bed, the latter at 3x the flow rate of the former. The two results should match exactly. I hope I’m making it clear what I’m asking for. This is a validation check on the multi-tank model.

 

[casualguitar]

You need to compare the temperature vs time of the first tank in your multi-tank 3 tank model with the temperature vs time for the single tank model of the full bed, the latter at 3x the flow rate of the former. The two results should match exactly. I hope I’m making it clear what I’m asking for. This is a validation check on the multi-tank model.

The idea is clear, and I can validation check by using the values of the other parameters (h,m,Ts).

One question - I'm not clear on the correct value of n. My thoughts are that temperature is evaluated at the 'centre' of a tank. n=0 is at the inlet and n=n is at the outlet. So to have 3 'tanks', we would need n=5?

Either way none of n=3,4,5 give exactly the same solution as the analytical one, so I need to do some testing regardless

 

[Chestermiller]

The idea is clear, and I can validation check by using the values of the other parameters (h,m,Ts).

One question - I'm not clear on the correct value of n. My thoughts are that temperature is evaluated at the 'centre' of a tank. n=0 is at the inlet and n=n is at the outlet. So to have 3 'tanks', we would need n=5?

Either way none of n=3,4,5 give exactly the same solution as the analytical one, so I need to do some testing regardless

The parameter n is the total number of tanks being considered, not a counting index. The counting index “i” is 1 for the first tank, 2 for the second tank, and 3 for the 3rd tank. The location of the center of the i’th tank is at $x=\frac{(2i-1)}{2}\Delta x$. The left side of the i'th tank is at $x=(i-1)\Delta x$. The right side of the i'th tank is at $x=i \Delta x$.

 

[casualguitar]

The parameter n is the total number of tanks being considered, not a counting index. The counting index “i” is 1 for the first tank, 2 for the second tank, and 3 for the 3rd tank. The location of the center of the i’th tank is at $x=\frac{(2i-1)}{2}\Delta x$. The left side of the i'th tank is at $x=(i-1)\Delta x$. The right side of the i'th tank is at $x=i \Delta x$.

Understood now. I read the original comments you had deriving the spatial variation model and I follow now.

Ok then sounds like n=3 and I should run some diagnostics

 

[casualguitar]

Understood now. I read the original comments you had deriving the spatial variation model and I follow now.

Ok then sounds like n=3 and I should run some diagnostics

So assuming the single tank model is correct (for now), given it matches the Mathematica model, I'm working with on getting the spatial variation model to match the single tank model. Specifically, getting tank index 1 with min = 0.05kg/s to match the single tank model, for min = 0.15kg/s
Some data for n=3, tank 1:

The mass holdup versus time:

This checks out because its just 8kg/3

dm/dt:

Zero as expected

dp/dh is also zero for all times

The differential equations all look ok to me, so I reckon the error is in how I have set up the initial and boundary conditions. Checking this now

 

[Chestermiller]

I have a feeling that the thing that is causing your numerical problems is the overall mass balance for each tank.

Are you still using an automatic integrator, and calculating the time derivative of the entire solution vector all at once for the integrator to use? If so, do you have the solution vector as the column vector ordered as $h_1,T_{S,1}, h_2, T_{S,2},h_3, T_{S,3},...$?

If you are using this scheme, you calculate the time derivative of $h_1$ for the first tank using $\dot{m}_0$ which does not change with time. Then you calculate $\dot{m}_1$ exiting tank 1 and entering tank 2 using the mass balance equation with the time derivative of $h_1$ that you just calculated for the first tank. Then you use this value of $\dot{m}_1$ entering tank 2 to calculate the time derivative of $h_2$ for the second tank. Then you calculate $\dot{m}_2$ exiting tank 2 and entering tank 3 using the mass balance equation with the time derivative of $h_2$ that you just calculated for the second tank. Then you use this value of $\dot{m}_2$ entering tank 3 to calculate the time derivative of $h_3$ for the third tank. etc. In the end, you will have all the time derivatives of the h's in all the tanks, which gets loaded into the time derivative vector. And, of course, in each of these steps you are also calculating the time derivative of the bed temperature for the specific tank, and also loading this into the derivative vector.

Is this what you are doing? If not, what is the scheme that you are applying?

 

[casualguitar]

Hi Chet, I had a separate report to finish today so I'm starting this now.

Are you still using an automatic integrator, and calculating the time derivative of the entire solution vector all at once for the integrator to use? If so, do you have the solution vector as the column vector ordered as h1,TS,1,h2,TS,2,h3,TS,3,...?

Yes to all this, except the solution vector is $h1,TS1,m1,etc$ where m1 is the flowrate. The way I've set up the mass holdup and mass flow does look suspect to me. I don't think this affects the solution but I could remove this from the solution vector in case I'm doing something wrong there in the setup

If you are using this scheme, you calculate the time derivative of h1 for the first tank using m˙0 which does not change with time. Then you calculate m˙1 exiting tank 1 and entering tank 2 using the mass balance equation with the time derivative of h1 that you just calculated for the first tank. Then you use this value of m˙1 entering tank 2 to calculate the time derivative of h2 for the second tank. Then you calculate m˙2 exiting tank 2 and entering tank 3 using the mass balance equation with the time derivative of h2 that you just calculated for the second tank. Then you use this value of m˙2 entering tank 3 to calculate the time derivative of h3 for the third tank. etc. In the end, you will have all the time derivatives of the h's in all the tanks, which gets loaded into the time derivative vector.

Yes to all of this, as far as I know I'm following this exact scheme unless the error is that I've coded this scheme incorrectly

And, of course, in each of these steps you are also calculating the time derivative of the bed temperature for the specific tank, and also loading this into the derivative vector.

Instead of this I'm just taking the solid energy balance, and generating the fluid temperature vector afterwards. The Ts array is known so I generate the dTs/dt array also and get the fluid temperature array from this. This image shows the flow for that:

I'll go for removing m from the solution array first

 

[Chestermiller]

Hi Chet, I had a separate report to finish today so I'm starting this now.

Yes to all this, except the solution vector is $h1,TS1,m1,etc$ where m1 is the flowrate. The way I've set up the mass holdup and mass flow does look suspect to me. I don't think this affects the solution but I could remove this from the solution vector in case I'm doing something wrong there in the setup

Yes to all of this, as far as I know I'm following this exact scheme unless the error is that I've coded this scheme incorrectly

Instead of this I'm just taking the solid energy balance, and generating the fluid temperature vector afterwards. The Ts array is known so I generate the dTs/dt array also and get the fluid temperature array from this. This image shows the flow for that:

View attachment 293952

I'll go for removing m from the solution array first

You are definitely not going to be integrating $dm_i/dt$ with respect to time. Once $dh_i/dt$ is known, $dm_i/dt$ can be evaluated and used to get $\dot{m}_i$ from $$\dot{m}_i=\dot{m}_{i-1}-\frac{dm_i}{dt}$$I guess that's what you do.

 

[Chestermiller]

Hi Chet, I had a separate report to finish today so I'm starting this now.

Yes to all this, except the solution vector is $h1,TS1,m1,etc$ where m1 is the flowrate. The way I've set up the mass holdup and mass flow does look suspect to me. I don't think this affects the solution but I could remove this from the solution vector in case I'm doing something wrong there in the setup

Yes to all of this, as far as I know I'm following this exact scheme unless the error is that I've coded this scheme incorrectly

Instead of this I'm just taking the solid energy balance, and generating the fluid temperature vector afterwards. The Ts array is known so I generate the dTs/dt array also and get the fluid temperature array from this. This image shows the flow for that:

View attachment 293952

I'll go for removing m from the solution array first

The fluid temperature and bed temperature must be solved simultaneously, not one and then the other.

 

[casualguitar]

I have a feeling that the thing that is causing your numerical problems is the overall mass balance for each tank.

Got it. Its the way I'm regenerating the derivatives. Solved now!:

So this leaves us at the point of having a working spatial variation model. We had discussed your $M_i$ model, however I haven't implemented this. Was this #M_i## model made purely to create a stepping stone between the single and multi tank model, or is there other reason to implement it?

If so I can implement this first thing tomorrow

 

[casualguitar]

The fluid temperature and bed temperature must be solved simultaneously, not one and then the other.

To answer this also, yes I am actually solving them simultaneously, but the way I've set up solve_ivp I don't have T_f in the solution array so I can't easily extract it. Thats why I regenerate it after

 

[casualguitar]

You are definitely not going to be integrating $dm_i/dt$ with respect to time. Once $dh_i/dt$ is known, $dm_i/dt$ can be evaluated and used to get $\dot{m}_i$ from $$\dot{m}_i=\dot{m}_{i-1}-\frac{dm_i}{dt}$$I guess that's what you do.

Missed this post. Yes correct I don't integrate it its just evaluated as part of the loop