Modelling of two phase flow in packed bed using conservation equations

[casualguitar]

Run some diagnostic calculations where you print out things like h and m vs t to see if they are behaving the way you expect them to behave.

I ran some tests to see if this is expected behaviour, and it seems to be. These plots are of $T_{fluid}$, $T_{solid}$, $H_{fluid}$ and mass holdup versus time. I plotted it like this to see what other variables changed significantly around the time temperature does the weird jump:

Zoomed in version around the saturation zone:

The 'resolution' of my plot was too low. I was calculating one data point for each time value, nothing in between. So the temperature seemed to jump up randomly in this case. However when I increased the number of points (10 points between every second), the jump is really a curve and it happens exactly when the enthalpy goes above the heat of vaporisation (h>200):

This seems to answer the question about the sudden temperature jump. Looking into the time constant difference now

 

[Chestermiller]

I ran some tests to see if this is expected behaviour, and it seems to be. These plots are of $T_{fluid}$, $T_{solid}$, $H_{fluid}$ and mass holdup versus time. I plotted it like this to see what other variables changed significantly around the time temperature does the weird jump:
View attachment 292572
Zoomed in version around the saturation zone:
View attachment 292573
The 'resolution' of my plot was too low. I was calculating one data point for each time value, nothing in between. So the temperature seemed to jump up randomly in this case. However when I increased the number of points (10 points between every second), the jump is really a curve and it happens exactly when the enthalpy goes above the heat of vaporisation (h>200):

View attachment 292576

This seems to answer the question about the sudden temperature jump. Looking into the time constant difference now

Very nice. Excellent

For the mass holdup, it would more instructive to examine it on a logarithmic scale.

 

[casualguitar]

For the mass holdup, it would more instructive to examine it on a logarithmic scale.

Understood yes the values are tiny in comparison to Temperature/Enthalpy. I'll change the scale for the mass holdup. Although it seems to do what is expected i.e. start at 8kg (V*rhoL = 0.01*800), gradually reduce across the phase change zone and then level out in the vapour zone

For the time constant, are we assuming this equation produces an accurate time constant?:
$$\tau=\frac{m_SC_{PS}}{\dot{m}C_{PV}}$$

And that the model should produce a value close to this if its correct? i.e. close to 200

I'll check the slope now with the higher resolution model

 

[casualguitar]

I get a temperature difference change by a factor of 10 every 700 seconds. I guess I made an error in specifying the slope I wanted. What I really wanted was Δln⁡(Tin−T)/Δt. This gives a time constant of about 300 seconds.

So using this equation:
$$\tau=\frac{m_SC_{PS}}{\dot{m}C_{PV}}$$

And the updated $C_s$ value of 1.1 (I was actually still using the old value up until recently as you spotted):
$$\tau=22s$$
And with these equations:
$$ slope = \Delta \ln{(T_{in}-T)}/\Delta t$$
$$\tau = 1/slope$$
I get:
$$\tau=18s$$
So there is relatively good similarity now it seems

Heres the code snippet:

The questions I would have here are:
1) Where do the two equations for tau come from?
2) Is it ok that I've only calculated the slope in the vapour region? i.e. I haven't found an 'averaged' slope across all regions?

 

[Chestermiller]

So using this equation:
$$\tau=\frac{m_SC_{PS}}{\dot{m}C_{PV}}$$

And the updated $C_s$ value of 1.1 (I was actually still using the old value up until recently as you spotted):
$$\tau=22s$$
And with these equations:
$$ slope = \Delta \ln{(T_{in}-T)}/\Delta t$$
$$\tau = 1/slope$$
I get:
$$\tau=18s$$
So there is relatively good similarity now it seems

Heres the code snippet:
View attachment 292577

The questions I would have here are:
1) Where do the two equations for tau come from?

See post #71. The following is an approximation for the vapor region only, based on the approximation in post #71:
$$m_SC_{PS}\frac{dT_S}{dt}=\dot{m}C_{PV}(T_{in}-T_S)$$The solution to this equation is $$\frac{T_{in}-T_S}{T_{in}-T_{sat}}=e^{-t/\tau}$$where t is measured from the time that we have all vapor, with $$\tau=\frac{m_SC_{PS}}{\dot{m}C_{PV}}$$

In this approximation, the vapor temperature T is approximately equal to TS.

2) Is it ok that I've only calculated the slope in the vapour region? i.e. I haven't found an 'averaged' slope across all regions?

Yes. The approximation applies only to the pure vapor.

 

[casualguitar]

See post #71. The following is an approximation for the vapor region only, based on the approximation in post #71:
$$m_SC_{PS}\frac{dT_S}{dt}=\dot{m}C_{PV}(T_{in}-T_S)$$The solution to this equation is $$\frac{T_{in}-T_S}{T_{in}-T_{sat}}=e^{-t/\tau}$$where t is measured from the time that we have all vapor, with $$\tau=\frac{m_SC_{PS}}{\dot{m}C_{PV}}$$

In this approximation, the vapor temperature T is approximately equal to TS.Yes. The approximation applies only to the pure vapor.

Ah I understand, and slope here is $\frac{-1}{\tau}$?

I'll clean up the code so far and split it into the various models done so far. Are there other things I should add before we spatial variation?

Actually I'm also going to map this timeline of model development into a powerpoint

 

[Chestermiller]

Ah I understand, and slope here is $\frac{-1}{\tau}$?

I'll clean up the code so far and split it into the various models done so far. Are there other things I should add before we spatial variation?

Actually I'm also going to map this timeline of model development into a powerpoint

I think the next step is to move on to spatial variation. I'll write something up tomorrow.

 

[casualguitar]

I think the next step is to move on to spatial variation. I'll write something up tomorrow.

Looking forward to it! Again thank you Chet for walking through the modelling process with me, this is hugely appreciated.

 

[Chestermiller]

FORMULATION

Considering the fact that the fluid density is varying substantially with time and spatial position within this system, in my judgment, there is asignificant advantage to formulating and solving the heat balance equation for the fluid using the so-called divergence representation of the equation (see below); this also provides the added advantage of automatically converting mass not only in the PDE form of the equation, but also when it is expressed in finite difference form. The divergence form of the mass and energy balance for the fluid are thus written as: $$\frac{\partial \rho}{\partial t}+\frac{\partial (\rho u)}{\partial x}=0$$and$$\frac{\partial (\rho h)}{\partial t}+\frac{\partial (\rho uh)}{\partial x}=\frac{\partial}{\partial x}\left(\rho D\frac{\partial h}{\partial x}\right)+\frac{Ua}{\epsilon}(T_S-T)$$where D is the axial thermal dispersion coefficient (assumed to be much larger than the thermal diffusivity of the fluid at all typical operating conditions of the bed).

In these equations, if we substitute $\phi=\rho u$ for the axial mass flux, and express the thermal dispersion coefficient in the typical approximate form D = ul (where l is the characteristic thermal dispersion length parameter), these equations become:
$$\frac{\partial \rho}{\partial t}+\frac{\partial \phi}{\partial x}=0$$and$$\frac{\partial (\rho h)}{\partial t}+\frac{\partial (\phi h)}{\partial x}=\frac{\partial}{\partial x}\left(\phi l\frac{\partial h}{\partial x}\right)+\frac{Ua}{\epsilon}(T_S-T)$$
In lieu of any definite information on the effects of fluid density on the dispersivity parameter l, it will probably be necessary to assume that l is a constant depending only on packing geometry and to treat it as an adjustable 'tuning" parameter.

Thoughts so far?

 

[casualguitar]

in my judgment, there is asignificant advantage to formulating and solving the heat balance equation for the fluid using the so-called divergence representation of the equation (see below); this also provides the added advantage of automatically converting mass not only in the PDE form of the equation, but also when it is expressed in finite difference form.

Got it

∂ρ∂t+∂ϕ∂x=0and∂(ρh)∂t+∂(ϕh)∂x=∂∂x(ϕl∂h∂x)+Uaϵ(TS−T)

To describe the above equations I would say the following, where every term is per unit volume:
1) The rate of increase of mass plus the rate at which mass leaves in the $x$ direction $= 0$
2) The rate of increase of energy plus the rate at which energy leaves in the $x$ direction = the rate of axial dispersion in the x direction plus the heat transfer between the fluid and solid

(I notice we have arrived at the same equation I started the post with (!), just using enthalpy instead of temperature. Nice)

In these equations, if we substitute ϕ=ρu for the axial mass flux, and express the thermal dispersion coefficient in the typical approximate form D = ul (where l is the characteristic thermal dispersion length parameter), these equations become

Understood

In lieu of any definite information on the effects of fluid density on the dispersivity parameter l, it will probably be necessary to assume that l is a constant depending only on packing geometry and to treat it as an adjustable 'tuning" parameter.

Understood. I think I have seen some correlations for $D$ the axial dispersion coefficient but I guess this is of secondary importance for now.

The mass balance is straightforward, however I do have questions on the energy balance to clear up.

Points of confusion:
1) Are you using the term 'heat balance' and 'energy balance' interchangeably?
2) Is this an advection term: $\frac{\partial (\phi h)}{\partial x}$? I'm not familiar with advection, other than its definition as the transfer of heat through fluid flow, so I guess this is it
3) So the axial dispersion coefficient is a 'catch all' term for a number of energy transport processes as far as I know (diffusion and conduction within the fluid, convection by the fluid in the axial direction, axial and transverse mixing of the fluid, etc, as far as I know). What terms are included in this term in our case?
4) The final term seems to be the convection between the fluid and the solid term. Why does this term not have a time or space derivative, because it will surely depend on space and time? Is it because 'time/space dependence' is built in, in the sense that $T$ and $T_S$ will change with time and space, meaning that the convection term will also change with time and space indirectly?

In lieu of any definite information on the effects of fluid density on the dispersivity parameter l, it will probably be necessary to assume that l is a constant depending only on packing geometry and to treat it as an adjustable 'tuning" parameter.

Got it

So now instead of solving a fluid and solid energy balance, we will be solving the mass balance for the fluid, and the energy balance for the fluid/solid instead?

 

[Chestermiller]

Got it

To describe the above equations I would say the following, where every term is per unit volume:
1) The rate of increase of mass plus the rate at which mass leaves in the $x$ direction $= 0$

I would call it the net rate at which mass leaves

2) The rate of increase of energy plus the rate at which energy leaves in the $x$ direction = the rate of axial dispersion in the x direction plus the heat transfer between the fluid and solid

I would call it "rate of increase of enthalpy per unit volume" and "the net rate at which enthalpy per unit volume leaves the control volume in the x direction"

(I notice we have arrived at the same equation I started the post with (!), just using enthalpy instead of temperature. Nice)Understood

Understood. I think I have seen some correlations for $D$ the axial dispersion coefficient but I guess this is of secondary importance for now.

The mass balance is straightforward, however I do have questions on the energy balance to clear up.

Points of confusion:
1) Are you using the term 'heat balance' and 'energy balance' interchangeably?

Yes.

2) Is this an advection term: $\frac{\partial (\phi h)}{\partial x}$? I'm not familiar with advection, other than its definition as the transfer of heat through fluid flow, so I guess this is it

As we said before, it represents the net advection of enthalpy per unit volume out of the control volume.

3) So the axial dispersion coefficient is a 'catch all' term for a number of energy transport processes as far as I know (diffusion and conduction within the fluid, convection by the fluid in the axial direction, axial and transverse mixing of the fluid, etc, as far as I know). What terms are included in this term in our case?

Like you said, it is a catch all for the net of all these effects. It's probably not going to be worthwhile trying to quantify it to a more detailed extent, and, in my judgment, it will be preferable to treat the dispersion length l as an adjustable parameter that the experimental results will be used to calibrate.

4) The final term seems to be the convection between the fluid and the solid term. Why does this term not have a time or space derivative, because it will surely depend on space and time? Is it because 'time/space dependence' is built in, in the sense that $T$ and $T_S$ will change with time and space, meaning that the convection term will also change with time and space indirectly?

No. The heat transfer coefficient will depend on the fluid parameters, velocities, and bed geometry, so, in reality, it can be a function of location and time. I too will be difficult to quantify, and I recommend, at least initially, treating it as a constant and, like the dispersibity, calibrating it to the experimental results. (This is how we treated it in the lumped model) Later, if necessary, we can "sharpen our pencils."

Got it

So now instead of solving a fluid and solid energy balance, we will be solving the mass balance for the fluid, and the energy balance for the fluid/solid instead?

No. The solid energy balance is still used. We have just had the add the mass balance because the density is a function of x and t (which needs to be accounted for, since the density varies by a factor of 1000 from liquid to vapor).

 

[casualguitar]

I would call it the net rate at which mass leaves

I would call it "rate of increase of enthalpy per unit volume" and "the net rate at which enthalpy per unit volume leaves the control volume in the x direction"

It is really the microscopic energy balance expressed in terms of enthalpy, assuming negligible viscous heating. See Eqn. S in Table 11.4-1 of Transport Phenomena by Bird, et al.

As we said before, it represents the net advection of enthalpy per unit volume out of the control volume.

Like you said, it is a catch all for the net of all these effects. It's probably not going to be worthwhile trying to quantify it to a more detailed extent, and, in my judgment, it will be preferable to treat the dispersion length l as an adjustable parameter that the experimental results will be used to calibrate.

No. The heat transfer coefficient will depend on the fluid parameters, velocities, and bed geometry, so, in reality, it can be a function of location and time. I too will be difficult to quantify, and I recommend, at least initially, treating it as a constant and, like the dispersibity, calibrating it to the experimental results. (This is how we treated it in the lumped model) Later, if necessary, we can "sharpen our pencils."

All understood. Ok so we're leaving the dispersion length and the heat transfer coefficient as calibration parameters

No. The solid energy balance is still used. We have just had the add the mass balance because the density is a function of x and t (which needs to be accounted for, since the density varies by a factor of 1000 from liquid to vapor).

Understood yes this is what I was trying to communicate also

I follow all of that. No further questions on any of the above

I do have two short questions on the assumptions we will make here:
1) Are we assuming here that saturation occurs at a specific temperature (rather than find the temperature with bubble/dew point correlations)
2) Should I prepare to start using correlations for any parameters previously kept constant (heat capacity, enthalpy, etc), or will we be using the equations used previously

 

[Chestermiller]

All understood. Ok so we're leaving the dispersion length and the heat transfer coefficient as calibration parameters

Understood yes this is what I was trying to communicate also

I follow all of that. No further questions on any of the above

I do have two short questions on the assumptions we will make here:
1) Are we assuming here that saturation occurs at a specific temperature (rather than find the temperature with bubble/dew point correlations)
2) Should I prepare to start using correlations for any parameters previously kept constant (heat capacity, enthalpy, etc), or will we be using the equations used previously

A foundation of success in doing modeling is to start with simple models and gradually build to more advanced versions incorporating additional features of scomplexity until one arrives at a final version that satisfies the modeler (and matches whatever experimental data has been developed to characterize the system and test the model). So, typically, the data development and the model development are done in tandem. So you should have in mind going to more advanced version of the model as needed.

Depending on your judgment, it may be necessary to employ more accurate correlations of physical property parameters as the model evolves to future versions.

Understand that it is a very bad idea to try to include everything in the model all at once. The concept is that, if you can't solve the more simple versions, you'll never be able to solve the most complex version. Moreover, for each version along the way, you will already have model results to show for your effort (which can be compared with the previous results to quantify the effect of the added complexity and to compare with experimental data).

Right now, as initially specified by you, the model assumes a pure component phase change. But there is no reason you can't start doing separate calculations by hand or by computer to evaluate the VLE behavior at a specified pressure of the effect of temperature in the 2-phase region. Start with one mole of an air mixture, and, using Raout's law, as a function of temperature, calculate the molar split of vapor and liquid, the specific volume of the vapor, the specific volume of the liquid, the average specific volume of the combination of vapor and liquid, the specific enthalpy of the liquid, the specific enthalpy of the vapor, the average specific enthalpy of the combination of vapor and liquid. Get an idea of what this all looks like. Maybe it can be incorporated pretty easily into the model.

One thing we haven't discussed yet is the possibility that the liquid will be traveling through the bed at a different velocity than the vapor (because of viscous drag on the liquid). You should also be thinking about the possibility of pressure variations due to drag, and how that might come into play for a two phase mixture in a packed bed (as well as for the single phase regions). (You have chosen a very complicated problem to work on). All these things should be in the back of your mind, and you should be starting to think about whether they need to be included in later versions of your model (or, at least, how to estimate whether to include them).

 

[Chestermiller]

RECOMMENDED FINITE DIFFERENCE SCHEME

The finite difference scheme I am recommending employs a finite difference grid with spacing $\Delta x$, grid cell boundaries at $(x=0,\ \Delta x,\ 2\Delta x...)$, and grid cell centers at $(x=\Delta x/2,\ 3\Delta x/2,\ 5\Delta x/2,\ ...$). The mass flux into the bed $\phi_0$ is known at the cell boundary x = 0, and we calculate the mass flux at all the other cell boundaries. On the other hand, we calculate the fluid temperatures, enthalpies, and densities at the cell centers.

Based on this description, and employing the method of lines, the 2nd order accurate finite difference representations of the mass balance and heat balance equations for the fluid are expressed as:
$$\frac{d\rho_x}{dt}=\frac{(\phi_{x-\Delta x/2}-\phi_{x+\Delta x/2})}{\Delta x}$$and$$\frac{d(\rho_xh_x)}{dt}=\rho_x\frac{dh_x}{dt}+h_x\frac{d\rho_x}{dt}$$$$=\frac{\phi_{x-\Delta x/2}\left(\frac{h_{x-\Delta x}+h_x}{2}\right)-\phi_{x+\Delta x/2}\left(\frac{h_{x+\Delta x}+h_x}{2}\right)}{\Delta x}$$$$+l\frac{\phi_{x+\Delta x/2}(h_{x+\Delta x}-h_x)-\phi_{x-\Delta x/2}(h_x-h_{x-\Delta x})}{(\Delta x)^2}$$$$+\frac{Ua(T_{S,x}-T_x)}{\epsilon}$$
OK so far?

Are you familiar with the finite differencing concepts of upwind differencing and numerical dispersion associated with modeling advection in transport models?

 

[casualguitar]

A foundation of success in doing modeling is to start with simple models and gradually build to more advanced versions incorporating additional features of complexity until one arrives at a final version that satisfies the modeler

Might have to add this as my laptop wallpaper

Understand that it is a very bad idea to try to include everything in the model all at once.

Got it. From my own limited experience I don't think I've never actually been able to implement the full version of any model first try. I always end up fighting the urge start simply though, as it feels unintuitive and as though it would be slower overall than tackling the entire thing first go I suppose. I'll remember to avoid this approach from now on, and start simply

But there is no reason you can't start doing separate calculations by hand or by computer to evaluate the VLE behavior at a specified pressure of the effect of temperature in the 2-phase region. Start with one mole of an air mixture, and, using Raout's law, as a function of temperature, calculate the molar split of vapor and liquid, the specific volume of the vapor, the specific volume of the liquid, the average specific volume of the combination of vapor and liquid, the specific enthalpy of the liquid, the specific enthalpy of the vapor, the average specific enthalpy of the combination of vapor and liquid. Get an idea of what this all looks like. Maybe it can be incorporated pretty easily into the model.

Understood. I've actually got some dew/bubble point functionality coded up, and also flash functionality (TP,PH, etc) for pure fluids and mixtures, so really in including two components instead of 1, a boiling range instead of a single temperature etc, the difficult bit for me is conceptualising the model. The functionality in code is ready to be used

One thing we haven't discussed yet is the possibility that the liquid will be traveling through the bed at a different velocity than the vapor (because of viscous drag on the liquid). You should also be thinking about the possibility of pressure variations due to drag, and how that might come into play for a two phase mixture in a packed bed (as well as for the single phase regions). (You have chosen a very complicated problem to work on)

Got it. I suppose with the approach of gradually building in complexity, this is something that doesn't have to be included now but can be looked at later (once other things, like the boiling range, have been added)

Actually one note I want to add is that I haven't yet described the existing functionality I've got available in code. Effectively any property/basic functionality like property calculations for pure fluids/mixtures, temperature and/or pressure dependent properties, flash for mixtures to calculate liquid fractions, or the composition of the gas/liquid at given conditions, etc is all possible. So as I mentioned above the most difficult part is likely the conceptualisation of the model. I mention this mainly because I don't want to limit the model capability (where possible) because of lack of access to property/process functionality like the above

 

[casualguitar]

The finite difference scheme I am recommending employs a finite difference grid with spacing Δx, grid cell boundaries at (x=0, Δx, 2Δx...), and grid cell centers at (x=Δx/2, 3Δx/2, 5Δx/2, ...).

Got it. I have implemented some basic FD schemes manually before so I'm not starting from zero (just close to it)

The mass flux into the bed ϕ0 is known at the cell boundary x = 0, and we calculate the mass flux at all the other cell boundaries

On the other hand, we calculate the fluid temperatures, enthalpies, and densities at the cell centers.

So this bit I'm not sure about. I think these questions clear up the confusion though:
1) You say we're calculating mass flux at the cell boundaries, so at $(x=0,\ \Delta x,\ 2\Delta x...)$. But in eq.1 the flux terms are evaluated $\phi_{x-\Delta x/2}$ and $\phi_{x+\Delta x/2}$, which are not at the cell boundaries if I'm correct. How can we know these values, if we are evaluating only mass flux at the boundaries i.e. $\phi_{x}$ and $\phi_{2\Delta x}$ etc

2) Similarly in eq.2, I follow the algebra, however if we're evaluating enthalpy at the cell centre, then I thought the subscript would be $x+\Delta x/2$ etc rather than $(x=0,\ \Delta x,\ 2\Delta x...)$, which are the boundaries.

So the above two questions are really the same and are just 'what does it mean to evaluate flux at the boundary, and the other properties at the cell centre?'

3) Final question on that is why not evaluate all properties at the boundary, or at the centre. I guess there is a limitation here in that we can't evaluate flux at the centre or the other properties at the boundary, but I can't see it yet. I will look for this limitation in the morning

Are you familiar with the finite differencing concepts of upwind differencing and numerical dispersion associated with modeling advection in transport models?

I'm familiar with upwind schemes. I do not know what numerical dispersion associated with modelling advection is though, and did not seem to find much online about this. What is this?

 

[Chestermiller]

Got it. I have implemented some basic FD schemes manually before so I'm not starting from zero (just close to it)So this bit I'm not sure about. I think these questions clear up the confusion though:
1) You say we're calculating mass flux at the cell boundaries, so at $(x=0,\ \Delta x,\ 2\Delta x...)$. But in eq.1 the flux terms are evaluated $\phi_{x-\Delta x/2}$ and $\phi_{x+\Delta x/2}$, which are not at the cell boundaries if I'm correct. How can we know these values, if we are evaluating only mass flux at the boundaries i.e. $\phi_{x}$ and $\phi_{2\Delta x}$ etc

2) Similarly in eq.2, I follow the algebra, however if we're evaluating enthalpy at the cell centre, then I thought the subscript would be $x+\Delta x/2$ etc rather than $(x=0,\ \Delta x,\ 2\Delta x...)$, which are the boundaries.

So the above two questions are really the same and are just 'what does it mean to evaluate flux at the boundary, and the other properties at the cell centre?

Let's consider the 3rd cell. The cell center is at $x=\frac{5}{2}\Delta x$. The left boundary of the cell is at $x-\frac{\Delta x}{2}=\frac{5}{2}\Delta x-\frac{\Delta x}{2}=2\Delta x$. The right boundary of the cell is at $x+\frac{\Delta x}{2}=\frac{5}{2}\Delta x+\frac{\Delta x}{2}=3\Delta x$.

3) Final question on that is why not evaluate all properties at the boundary, or at the centre. I guess there is a limitation here in that we can't evaluate flux at the centre or the other properties at the boundary, but I can't see it yet. I will look for this limitation in the morning

For accuracy of the finite difference scheme, it is better to evaluate the fluxes at the cell boundaries, and the non-fluxes at the cell centers. For example, in a conductive heat transfer finite difference scheme, we evaluate the q's ($q=-k\frac{\partial T}{\partial x}$) at the cell boundaries and the temperatures at the cell centers (or vice versa, which is equivalent). The point is that the fluxes are evaluated a half-grid-space away from the temperatures. In an advective transport problem like ours, another advantage of this approach is that the grid cells preserve the physically intuitive analogy to a series of mixed tanks, with fluid flowing in at the left cell boundary and exiting at the right cell boundary.

I'm familiar with upwind schemes. I do not know what numerical dispersion associated with modelling advection is though, and did not seem to find much online about this. What is this?

Are you saying that you are familiar with upwind (first order) finite difference schemes applied to advection-dispersion problems?

 

[casualguitar]

Let's consider the 3rd cell. The cell center is at x=52Δx. The left boundary of the cell is at x−Δx2=52Δx−Δx2=2Δx. The right boundary of the cell is at x+Δx2=52Δx+Δx2=3Δx.

Understood

The point is that the fluxes are evaluated a half-grid-space away from the temperatures. In an advective transport problem like ours, another advantage of this approach is that the grid cells preserve the physically intuitive analogy to a series of mixed tanks, with fluid flowing in at the left cell boundary and exiting at the right cell boundary.

Got it

Are you saying that you are familiar with upwind (first order) finite difference schemes applied to advection-dispersion problems?

I have seen (and understood) examples of upwind finite difference schemes applied to advection equations. The key idea is that we use information 'upstream' of the flow direction i.e. $x < x_{i}$ to solve at $x_{i}$. I have not seen upwind schemes applied to advection-dispersion equations. Are there key differences between advection, and advection-dispersion upwind schemes? Or will an understanding of upwind schemes applied to advection problems be enough here i.e. the addition of the dispersion term doesn't add any other difficulties

 

[Chestermiller]

Understood

Got it

I have seen (and understood) examples of upwind finite difference schemes applied to advection equations. The key idea is that we use information 'upstream' of the flow direction i.e. $x < x_{i}$ to solve at $x_{i}$. I have not seen upwind schemes applied to advection-dispersion equations. Are there key differences between advection, and advection-dispersion upwind schemes? Or will an understanding of upwind schemes applied to advection problems be enough here i.e. the addition of the dispersion term doesn't add any other difficulties

If you understand that upwind difference schemes are first order accurate and introduce unwanted numerical dispersion into the calculation (over and above actual physical diffusion/dispersion), while central difference schemes (e.g., the one we are using), while sometime problematic with respect to advection (by sometimes resulting in unwanted "zig-zags" in the axial profiles) do not introduce numerical dispersion, that will be sufficient.

 

[casualguitar]

(over and above actual physical diffusion/dispersion)

Ah I didn't understand what you meant by this previously. I do now

that will be sufficient.

All understood. Ready to go

Is there anything I should do before we continue?

 

[Chestermiller]

Ah I didn't understand what you meant by this previously. I do now

All understood. Ready to go

Is there anything I should do before we continue?

Let's continue.

Here's a math problem for you: What would the finite difference version of our heat balance equation reduce to if the dispersion length $l$ in the equation just happened to be equal to half our selected grid spacing $\Delta x$? (That is $l=\frac{\Delta x}{2}$)

 

[casualguitar]

Let's continue.

Here's a math problem for you: What would the finite difference version of our heat balance equation reduce to if the dispersion length $l$ in the equation just happened to be equal to half our selected grid spacing $\Delta x$? (That is $l=\frac{\Delta x}{2}$)

We'll get:
$=\frac{\phi_{x-\Delta x/2}h_{x-\Delta x} -\phi_{x+\Delta x/2} h_{x}}{\Delta x}$ $+\frac{Ua(T_{S,x}-T_x)}{\epsilon}$

Working:

 

[Chestermiller]

We'll get:
$=\frac{\phi_{x-\Delta x/2}h_{x-\Delta x} -\phi_{x+\Delta x/2} h_{x}}{\Delta x}$ $+\frac{Ua(T_{S,x}-T_x)}{\epsilon}$

Working:
View attachment 292769

Correct. Do you realize that this is the exact same result you would get if you used upwind differencing on the advection term and if the system did not feature any real physical dispersion? Any thoughts on how you could capitalize on this?

 

[casualguitar]

Correct. Do you realize that this is the exact same result you would get if you used upwind differencing on the advection term and if the system did not feature any real physical dispersion? Any thoughts on how you could capitalize on this?

Ok so finite difference and $=\frac{\Delta x}{2}$ gives the same result as upwind differencing and $l = 0$?

No I've thought about and I have no good suggestions so far. I guess we are somehow looking to use the advantages the upwind differencing scheme gives, while not assuming physical dispersion is zero? But no, I can't see how we can capitalize. How would we do this?

I have thought about it some more and no nothing jumps out to me. Hmm any hints?

 

[Chestermiller]

Ok so finite difference and $=\frac{\Delta x}{2}$ gives the same result as upwind differencing and $l = 0$?

No I've thought about and I have no good suggestions so far. I guess we are somehow looking to use the advantages the upwind differencing scheme gives, while not assuming physical dispersion is zero? But no, I can't see how we can capitalize! How would we do this?

Your guess is correct. We are seeing that a grid spacing of $\Delta x=2l$ gives the exact same amount of numerical dispersion for the upwind differencing scheme without actual dispersion as the more accurate 2nd order scheme with actual dispersion included. So, in this case, the upwind differencing numerical dispersion exactly matches the actual dispersion (at least in the 2nd order approximation). Suppose that, in our full 2nd order scheme with dispersion, a grid spacing of $\Delta x=2l$ were found to be sufficiently adequate to accurately solve the true differential differential equation (and thus to accurately model the real physical system, assuming a constant value of $l$ in the model could do this). Then we could use the upwind difference scheme without dispersion (with its desirable advantages) to model the system, and could employ the grid spacing $\Delta x$, rather than the dispersion parameter $l$, as a calibration parameter for matching the experimental data. When we got a match, we would know the value of $l$ that would also match the data with the full 2nd order scheme including dispersion. Your thoughts? (This is entirely new. I have never seen this simple strategy for a dispersion problem with dispersivity parameter $l$ identified before.)

 

[casualguitar]

We are seeing that a grid spacing of Δx=2l gives the exact same amount of numerical dispersion for the upwind differencing scheme without actual dispersion as the more accurate 2nd order scheme with actual dispersion included. So, in this case, the upwind differencing numerical dispersion exactly matches the actual dispersion (at least in the 2nd order approximation).

Got it

Suppose that, in our full 2nd order scheme with dispersion, a grid spacing of Δx=2l were found to be sufficiently adequate to accurately solve the true differential differential equation (and thus to accurately model the real physical system, assuming a constant value of l in the model could do this). Then we could use the upwind difference scheme without dispersion (with its desirable advantages) to model the system, and could employ the grid spacing δx, rather than the dispersion parameter l, as a calibration parameter for matching the experimental data. When we got a match, we would know the value of l that would also match the data with the full 2nd order scheme including dispersion. Your thoughts?

So if I follow this correctly, you're saying:
1) Let's use $\Delta x=2l$ to build the model because we can use the upwind scheme without dispersion
2) Once experimental results are produced, we can use $\Delta x$ to match the model with the data
3) If we can match the model to the data well enough, then we can say that using an upwind scheme without dispersion is a valid approach
4) If it is a valid approach, then we can also just divide our $\Delta x$ value by two to find the value of $l$ that would produce the same results in the model with dispersion?

If this is all correct then I follow fully so far

 

[Chestermiller]

Got it

So if I follow this correctly, you're saying:
1) Let's use $\Delta x=2l$ to build the model because we can use the upwind scheme without dispersion
2) Once experimental results are produced, we can use $\Delta x$ to match the model with the data
3) If we can match the model to the data well enough, then we can say that using an upwind scheme without dispersion is a valid approach

A valid approach to simulate central differencing with dispersion included.

4) If it is a valid approach, then we can also just divide our $\Delta x$ value by two to find the value of $l$ that would produce the same results in the model with dispersion?

If this is all correct then I follow fully so far

Yes, you are right on target. But, if this approach is used, I need you to confirm that you are comfortable with it. After all, it is really your research, and you are in the drivers seat. I'm just here to help you identify the possibilities and offer suggestions.

In terms of calibrating the model, might I also suggest doing experiments not on the actual fluid, but on say, just a gas or a liquid without phase change.

Continuing the development, we recognize that, using our unwinding approach, what we are dealing with here conceptually is a series of stirred tanks with the fluid entering each tank having the same enthalpy at that present in the previous tank. For the full heat balance equation corresponding to your result in post #112, what do you get if you substitute $\phi_{x+\Delta x/2}$ from the mass balance equation into this relationship?

 

[casualguitar]

Yes, you are right on target. But, if this approach is used, I need you to confirm that you are comfortable with it. After all, it is really your research, and you are in the drivers seat. I'm just here to help you identify the possibilities and offer suggestions.

I'm comfortable with it. One hope of mine is that work I publish will be 'useful' to others. As I learn I'm hoping that I will get a better intuition for spotting 'useful' research. For example, is this approach (above) as useful as one that uses known correlations for the dispersion coefficient, for example? I think at some stage I would like to compare this approach with a model that uses correlations for the dispersion coefficient, however as you said the idea is to gradually build in complexity so I think this is a solid approach that I'm 100% happy to work on

But yes as I go I will aim to focus on doing the widely 'useful' model rather than any other

In terms of calibrating the model, might I also suggest doing experiments not on the actual fluid, but on say, just a gas or a liquid without phase change.

Can do there are plenty options available in the lab for this

Continuing the development, we recognize that, using our unwinding approach, what we are dealing with here conceptually is a series of stirred tanks with the fluid entering each tank having the same enthalpy at that present in the previous tank. For the full unwinding heat balance equation corresponding to your result in post #112, what do you get if you substitute ϕx+Δx/2 from the mass balance equation into this relationship?

Working on this now

 

[casualguitar]

what do you get if you substitute ϕx+Δx/2 from the mass balance equation into this relationship?

Hmm I'm not sure I understand. Are you asking me to solve for \phi_{x+\Delta x/2} in the fluid mass balance equation, and to sub this into the heat balance in #112?

Ah if you are saying this I think I see why. We would have the mass flux in terms of the mass flux 'upwind' and the density? The idea here being to get things in terms of the upwind parameters? If the above is true, I am not sure how to deal with the density derivative

Got a result. Posting now

 

[Chestermiller]

Hmm I'm not sure I understand. Are you asking me to solve for \phi_{x+\Delta x/2} in the fluid mass balance equation, and to sub this into the heat balance in #112?

Ah if you are saying this I think I see why. We would have the mass flux in terms of the mass flux 'upwind' and the density? The idea here being to get things in terms of the upwind parameters? If the above is true, I am not sure how to deal with the density derivative

Got a result. Posting now

Yes. That is what I am asking for.