Modelling of two phase flow in packed bed using conservation equations

casualguitar said:

Hi Chet, my apologies for the delay. I was running a few simulations for various values of n. Anything greater than n=10 takes considerable time now. n=30 took 24 hours actually.

Anyway, I wanted to do a 'grid independence' study of sorts. If you run hot air through the bed for x amount of time, we get a thermocline in the bed. I did this for n=5,10,20, and 30 and here are the results. Clearly increasing n results in a 'tighter' thermocline or a sharper drop from high temperature to low:
n=5:
View attachment 315477
n=10:
View attachment 315478
n=20:
View attachment 315479
n=30:
View attachment 315480
I'm running higher n now (n=100) however this will take a few days. The goal of this analysis is to 1) find the value of n where the output plots stop changing i.e. grid independence and 2) to quantify the losses in the bed across a charge/discharge cycle i.e. if we discharge for 100s and then charge for 100s, how close to we get to the original temperature profile?

No questions on this, just updating

I think they will continue to change as you increase n. Our idea was to calibrate the value of n to the actual data on column operation to match the actual amount of dispersion occurring experimentally.

Chestermiller said:

I think they will continue to change as you increase n

Are you saying that you don't think there will be a point where increasing n no longer changes the output i.e. grid independence?

Chestermiller said:

Our idea was to calibrate the value of n to the actual data on column operation to match the actual amount of dispersion occurring experimentally.

I understand this yes. I see how we could use n as a calibration parameter

I also understand that the grid independence value of n is likely to be very different from the calibrated value of n

Edit: Just reading back over previous posts. Didn't you say that one big advantage of this model is that there is actually no numerical dispersion (I might be referencing you incorrectly)? However I don't think all upwind schemes have zero numerical dispersion, so why does ours?

Why would the output change as it does above, if there is no numerical dispersion present in the system? Apologies for my misunderstanding

casualguitar said:

Are you saying that you don't think there will be a point where increasing n no longer changes the output i.e. grid independence?

Yes. Or, more precisely, the temperature profile will approach a step function at large n.

casualguitar said:

I understand this yes. I see how we could use n as a calibration parameter

I also understand that the grid independence value of n is likely to be very different from the calibrated value of n

Edit: Just reading back over previous posts. Didn't you say that one big advantage of this model is that there is actually no numerical dispersion (I might be referencing you incorrectly)? However I don't think all upwind schemes have zero numerical dispersion, so why does ours?

All upwind schemes feature numerical dispersion. This is a key feature of our scheme. We are using the numerical dispersion of the unwinding scheme to simulate the actual dispersion. They will match when ##\Delta x=l/2## (or wherever we determined the 2 belongs; I con't remember), where l is the dispersion length parameter.

casualguitar said:

Why would the output change as it does above, if there is no numerical dispersion present in the system? Apologies for my misunderstanding

There is no numerical dispersion if unwinding is not used. We have included dispersion explicitly in the model, but have found that the mathematical representation of the dispersion is such that, when the value of ##\Delta x## is chosen in a certain specific relation to the dispersion length l, the equations become exactly the same as with numerical dispersion of the advection terms, but without unwinding included in the equations.

Chestermiller said:

Yes. Or, more precisely, the temperature profile will approach a step function at large n.

Yes I'm running n=100 currently (will take a few days and it seems to be about half way through). So hopefully this value of n will be high enough to see something resembling this

Chestermiller said:

All upwind schemes feature numerical dispersion. This is a key feature of our scheme. We are using the numerical dispersion of the unwinding scheme to simulate the actual dispersion. They will match when Δx=l/2 (or wherever we determined the 2 belongs; I con't remember), where l is the dispersion length parameter.

At low values of n the plots produced are obviously 'low resolution' in that they are somewhat jagged (not smooth). Is it not true that this will be unavoidable if if it turns out that a low value of n matches numerical dispersion with actual dispersion?

Chestermiller said:

but have found that the mathematical representation of the dispersion is such that, when the value of Δx is chosen in a certain specific relation to the dispersion length l, the equations become exactly the same as with numerical dispersion of the advection terms, but without unwinding included in the equations.

Understood

casualguitar said:

At low values of n the plots produced are obviously 'low resolution' in that they are somewhat jagged (not smooth). Is it not true that this will be unavoidable if if it turns out that a low value of n matches numerical dispersion with actual dispersion?

I don't think. that will be the case.o. Typically, the dispersion length will be on the order of a few particle diameters.

Chestermiller said:

I don't think. that will be the case.o. Typically, the dispersion length will be on the order of a few particle diameters.

Ah interesting so for a bed length of 0.1m, and a particle diameter of 0.001m, and assuming a value of 5 particle diameters for the dispersion length, then the value of n (where n = L/dz) would be n=20? Or have I miscalculated here?

Edit: Just one other note - it seems a bit unreasonable that n=100 will take a few days to run, so I'm currently looking into replacing the thermo library with algebraic correlations where possible

casualguitar said:

Ah interesting so for a bed length of 0.1m, and a particle diameter of 0.001m, and assuming a value of 5 particle diameters for the dispersion length, then the value of n (where n = L/dz) would be n=20? Or have I miscalculated here?

yes.

casualguitar said:

Edit: Just one other note - it seems a bit unreasonable that n=100 will take a few days to run, so I'm currently looking into replacing the thermo library with algebraic correlations where possible

Good idea. This seems to be the right direction to go.

Chestermiller said:

yes.

Good idea. This seems to be the right direction to go.

Is that a yes to n=20, or a yes to it being a miscalculation?

Great I’ll work on this. There is actually a possibility that the thermo library will allow us to trade accuracy for time. I’ll check this first.

Anyway, we’re using thermo to get temperature as a function of enthalpy and pressure, and also density as a function of enthalpy and pressure, for the liquid, gas, and two phase regions

casualguitar said:

Is that a yes to n=20, or a yes to it being a miscalculation?

n=20

Chestermiller said:

n=20

Definitely noticed differences going from n=20 to n=30. n=100 is actually still running believe it or not. As discussed I'm working on speeding up this script now

Chestermiller said:

n=20

Hi Chet,
I'm about to start on that "modelling axial heat dispersion as numerical dispersion" analysis. But I did spot one thing in the model derivation -

If we're saying a given ##\Delta x## introduces a given amount of axial heat dispersion into the system, and we have tuned our model to one set of process conditions, can we really say that if we change the length/flow/pressure then we will still see approximately the same amount of axial heat dispersion? The model I tuned to uses much larger flows and a much larger diamter than we have in our model

casualguitar said:

Hi Chet,
I'm about to start on that "modelling axial heat dispersion as numerical dispersion" analysis. But I did spot one thing in the model derivation -

If we're saying a given ##\Delta x## introduces a given amount of axial heat dispersion into the system, and we have tuned our model to one set of process conditions, can we really say that if we change the length/flow/pressure then we will still see approximately the same amount of axial heat dispersion? The model I tuned to uses much larger flows and a much larger diamter than we have in our model

According to our derivation, yes. The dispersion length l is some constant times the particle size only. So, if we have the same particle size, the l and ##\Delta x## should be the same.

Chestermiller said:

According to our derivation, yes. The dispersion length l is some constant times the particle size only. So, if we have the same particle size, the l and ##\Delta x## should be the same.

Hi Chet, back on this. In terms of looking further into this "modelling numerical dispersion as physical dispersion" observation, there seems to be a few options for progression.
1. Check if this happens for other equations in other domains in an equivalent sense
2. Consider the observation for this specific case in more depth.

I think option 2 is a clearer path but if you like the idea of option 1 I'm happy to work on that.

Option 2:
I think we have two main possibilities here -
1. Performance Analysis - comparing the performance of the central differencing scheme to that of the upwind scheme for a range of conditions, by comparing computation speed, accuracy, stability, and maybe memory usage
2. We could look into the characteristics of this PDE to see if there is any underlying mechanism that allows this to occur.

I like the performance analysis approach, it would be straightforward to get a sense of how much faster the upwind scheme is (if at all), and how computationally efficient it is in comparison to the 2nd order scheme

casualguitar said:

Hi Chet, back on this. In terms of looking further into this "modelling numerical dispersion as physical dispersion" observation, there seems to be a few options for progression.
1. Check if this happens for other equations in other domains in an equivalent sense

I don't know what this means.

casualguitar said:

2. Consider the observation for this specific case in more depth.

I think option 2 is a clearer path but if you like the idea of option 1 I'm happy to work on that.

Option 2:
I think we have two main possibilities here -
1. Performance Analysis - comparing the performance of the central differencing scheme to that of the upwind scheme for a range of conditions, by comparing computation speed, accuracy, stability, and maybe memory usage

We are already using central differencing, with a specific value of the finite difference length. For this specific value of the finite difference length, the equations reduce to the upwind differencing equations without explicit dispersion included in the PDEs. Are you saying you want to use a specific value of l, and vary ##\Delta x## to test the accuracy of the central difference finite difference scheme and the accuracy of using our relationship with ##Delta x=2l##?

Chestermiller said:

We are already using central differencing, with a specific value of the finite difference length.

Sorry yes my mistake I suppose there is no "comparing upwind and central differencing schemes", because only have central differencing.

Chestermiller said:

Are you saying you want to use a specific value of l, and vary Δx to test the accuracy of the central difference finite difference scheme and the accuracy of using our relationship with Deltax=2l?

I am saying this, but we may have already done this unless I'm misunderstanding.

We calculated ##l## with that Ruthven et al correlation, where ##l## is a function of particle diameter. This allowed us to calculate a ##\Delta x## value (I think it was about ##0.027m##). Then separately we matched our model to experimental data by varying the number of tanks, n, until our temperautre profiles matched experimental data, and then calculated another ##\Delta x## value from that using ##\Delta x## = L/n (the value came to be about ##0.03m##. So the two values calculated by different methods are almost equal. Is this equivalent to what you're saying above?

casualguitar said:

Sorry yes my mistake I suppose there is no "comparing upwind and central differencing schemes", because only have central differencing.I am saying this, but we may have already done this unless I'm misunderstanding.

We calculated ##l## with that Ruthven et al correlation, where ##l## is a function of particle diameter. This allowed us to calculate a ##\Delta x## value (I think it was about ##0.027m##). Then separately we matched our model to experimental data by varying the number of tanks, n, until our temperautre profiles matched experimental data, and then calculated another ##\Delta x## value from that using ##\Delta x## = L/n (the value came to be about ##0.03m##. So the two values calculated by different methods are almost equal. Is this equivalent to what you're saying above?

##\Delta x## is always 2l, where l = kDp; n has to be chosen such ##n=\frac{L}{\Delta x}=\frac{L}{2l}=\frac{L}{2kD_p}##

The real issues are if l is a constant, or whether k is a function of Re and Pr or Sc. We are taking it to be constant based on the groundwater porous media references. Other researchers are taking k to be a function of Re and Pr or Sc based on packed bed literature. There really should be no difference except for the flow in porous media being laminar and the flow in packed beds sometimes being turbulent (because of the larger tortuous flow channel dimensions). That's why we were checking to see if, over the range if Re in the beds of interest or our bed, Re varies enough that their correlations to provide substantial variation in k (and whether, in our model, a constant value of k is adequate to match the observational data). If k is a substantial function of Re and Pr or Sc, we can't use ##\Delta x=2l=2kD_p## along with our upwind differencing approximation (to our central differencing equations) to model the dispersion.

Another issue is numerical: whether ##\Delta x## is small enough using ##\Delta x=2l=2kD_p## for the finite difference equations to be sufficiently accurate numerically. This can be checked numerically by fixing l, varying ##\Delta x## independently, and not reducing the central difference equations to the forward difference version obtained if ##l=\Delta x/2##. (or wherever the 2 goes; I can never remember). But this would require some additional coding with additional chance of making numerical errors; it would be a real pain in the butt.

Chestermiller said:

If k is a substantial function of Re and Pr or Sc, we can't use Δx=2l=2kDp along with our upwind differencing approximation (to our central differencing equations) to model the dispersion.

To check this assumption for the CO2 model, we did something like the plot below:


If suitable, I will redo this calculation for this system and run for a range of Re values. It seems like we will at least get an Re 'range of validity' for the ##\Delta x = 2l## assumption, and won't have to completely throw it out. Is this similar to what you had in mind?

casualguitar said:

To check this assumption for the CO2 model, we did something like the plot below:


If suitable, I will redo this calculation for this system and run for a range of Re values. It seems like we will at least get an Re 'range of validity' for the ##\Delta x = 2l## assumption, and won't have to completely throw it out. Is this similar to what you had in mind?

A direct comparison with their numerical results and yours would be more powerful.

casualguitar said:

Some advice. When you are making a technical argument (like the "l" value being close to constant as a function of the Reynolds number), and you are confident in your contention, it is best to present your argument using visuals that most strongly emphasize your contention. In the case of this graph, even though you are using a log scale for the ordinate, I think it would be much better if you used a full decade on the ordinate scale. This would better illustrate the small percentage variation in the function, compared to the variation shown which only goes from 1 to 2.4, but spans the entire vertical domain of your graph. How you present things visually is very important to making your case.

Chestermiller said:

Some advice. When you are making a technical argument (like the "l" value being close to constant as a function of the Reynolds number), and you are confident in your contention, it is best to present your argument using visuals that most strongly emphasize your contention. In the case of this graph, even though you are using a log scale for the ordinate, I think it would be much better if you used a full decade on the ordinate scale. This would better illustrate the small percentage variation in the function, compared to the variation shown which only goes from 1 to 2.4, but spans the entire vertical domain of your graph. How you present things visually is very important to making your case.

Thanks Chet. Something like this (going from 0.1 to 10)?

Chestermiller said:

A direct comparison with their numerical results and yours would be more powerful.

If you mean calculate an ##\frac{l_h}{d_p}## value for their data, then I calculated this as being equal to about 1.33 which is inside our range of 1.2 - 2.4. I can post this calculation if necessary.

So from the above we can possibly conclude that:
1) ##l## is not a significant function of Re inside our Re range of interest (0 to about 6000)
2) The range of ##l## values we get is similar to that of gas phase experimental data for air

We do have condensing/evaporating flow though, can the above be extended to the liquid phase? I haven't managed to find much clear experimental data on two phase flow. The Brouwers et al reference was the most useful thing I had found: https://josbrouwers.bwk.tue.nl/publications/Journal21.pdf

casualguitar said:

Thanks Chet. Something like this (going from 0.1 to 10)?


If you mean calculate an ##\frac{l_h}{d_p}## value for their data, then I calculated this as being equal to about 1.33 which is inside our range of 1.2 - 2.4. I can post this calculation if necessary.

So from the above we can possibly conclude that:
1) ##l## is not a significant function of Re inside our Re range of interest (0 to about 6000)
2) The range of ##l## values we get is similar to that of gas phase experimental data for air

We do have condensing/evaporating flow though, can the above be extended to the liquid phase? I haven't managed to find much clear experimental data on two phase flow. The Brouwers et al reference was the most useful thing I had found: https://josbrouwers.bwk.tue.nl/publications/Journal21.pdf

Check Collins' book for a chapter on two phase flow in porous media.

Chestermiller said:

Check Collins' book for a chapter on two phase flow in porous media.

Hi Chet, I think I found something relevant in section 8.30 "Miscible Displacement in Porous Media", Figure 8.2 linked here:848448


Would it be reasonable to make the case that the dispersion coefficient is not a strong function of phase but velocity?

If so, I could:
1) take the range of velocities we see in our simulations
2) calculate dispersion coefficients using our previously calculated dispersion length, ##l## of about ##0.03m##
3) plot our dispersion coefficients versus the data above and see if we get a reasonable match

casualguitar said:

Hi Chet, I think I found something relevant in section 8.30 "Miscible Displacement in Porous Media", Figure 8.2 linked here:848448


Would it be reasonable to make the case that the dispersion coefficient is not a strong function of phase but velocity?

If so, I could:
1) take the range of velocities we see in our simulations
2) calculate dispersion coefficients using our previously calculated dispersion length, ##l## of about ##0.03m##
3) plot our dispersion coefficients versus the data above and see if we get a reasonable match

4) in the case that we do get a match (with the Carberry et al liquid and gas phase dispersion coefficients above), we can make the claim that our dispersion length of ##0.03m## is suitable for both liquid and gas phase in the range of velocities simulated

casualguitar said:

4) in the case that we do get a match (with the Carberry et al liquid and gas phase dispersion coefficients above), we can make the claim that our dispersion length of ##0.03m## is suitable for both liquid and gas phase in the range of velocities simulated

Hi Chet, in your view is the above approach reasonable? This would mean that I am not tuning to any liquid phase temperature profiles, but instead tuning this two phase model to gas phase temperature profiles only (of which there is lots of data), and also showing that the dispersion coefficients present in the gas phase are applicable in the liquid phase, and therefore we can expect that tuning to gas phase profiles would give us reasonable liquid phase accuracy also

Chestermiller said:

I think they will continue to change as you increase n. Our idea was to calibrate the value of n to the actual data on column operation to match the actual amount of dispersion occurring experimentally.

Hi Chet, I hope you're doing well. I see you're very active here lately which is great (really is a god send for postgraduate students).

I wanted to respond to the above comment from you. There now exists an experimental setup that will allow for the measurement of axial temperature profiles in a packed bed. The bed can be maintained at a set pressure value that we choose, between about 1 and 50 bar. We can also measure the flow into and out of the bed.

As you say above, to measure the dispersion coefficient we just need to tune ##\Delta x## to the experimental temperature profiles and calculate the dispersion length (and therefore the dispersion coefficient range) from this.

I was planning on doing this for gas phase only, liquid phase only, and for the condensation/evaporation case, for a range of pressures.

I was only wondering if you see any other useful experiments that could be run that might give insight into the system in practice? Or experiments that might typically be run on a system like this. If not then no worries however I had to ask

casualguitar said:

Hi Chet, I hope you're doing well. I see you're very active here lately which is great (really is a god send for postgraduate students).

I wanted to respond to the above comment from you. There now exists an experimental setup that will allow for the measurement of axial temperature profiles in a packed bed. The bed can be maintained at a set pressure value that we choose, between about 1 and 50 bar. We can also measure the flow into and out of the bed.

As you say above, to measure the dispersion coefficient we just need to tune ##\Delta x## to the experimental temperature profiles and calculate the dispersion length (and therefore the dispersion coefficient range) from this.

I was planning on doing this for gas phase only, liquid phase only, and for the condensation/evaporation case, for a range of pressures.

I was only wondering if you see any other useful experiments that could be run that might give insight into the system in practice? Or experiments that might typically be run on a system like this. If not then no worries however I had to ask

I think that, at this point in your thesis work, you should be focusing on minimizing the amount of additional work that you commit to, and, instead be focusing on wrapping things up. You need to negotiate an agreed upon end point to your work.

That being said, I would run the model to see which kind to experimental test of temperature profiles would give you the best shot on an accurate answer without excessive effort.

Chestermiller said:

I think that, at this point in your thesis work, you should be focusing on minimizing the amount of additional work that you commit to, and, instead be focusing on wrapping things up. You need to negotiate an agreed upon end point to your work.

That being said, I would run the model to see which kind to experimental test of temperature profiles would give you the best shot on an accurate answer without excessive effort.

Hi Chet, my apologies yes I have effectively completed thesis writing at this point, I was more curious about future work on this than anything.

I did have one question on the modelling of physical dispersion using numerical dispersion (I was hoping to flesh this out a small bit more as it is very central to the model) - is it true to say that for this set of equations, if we use a second order accurate finite difference representation of the mass/energy balances, and let ##l## equal to ##\Delta x/2##, it will always mean that the amount of numerical dispersion we introduce is equal to the amount of physical dispersion we see? i.e. are there conditions where the physical dispersion cannot be modelled using numerical dispersion, for this set of equations? Or can it be mathematically shown that numerical dispersion must equal physical dispersion when the grid spacing = ##2l##?

The equations:
\begin{equation}
\frac{\partial h_x}{\partial t} = \frac{\phi_{x-\Delta x/2}(h_{x-\Delta x} - h_x)}{\rho_x \Delta x} + \frac{Ua}{\rho_x\epsilon}(T_{s,x}-T_x)
\end{equation}

\begin{equation}
m_j\frac{dh_j}{dt} = \dot{m}_{j-1}(h_{j-1} - h_j) + UA(T_{S,j}-T_j)
\end{equation}

casualguitar said:

Hi Chet, my apologies yes I have effectively completed thesis writing at this point, I was more curious about future work on this than anything.

I did have one question on the modelling of physical dispersion using numerical dispersion (I was hoping to flesh this out a small bit more as it is very central to the model) - is it true to say that for this set of equations, if we use a second order accurate finite difference representation of the mass/energy balances, and let ##l## equal to ##\Delta x/2##, it will always mean that the amount of numerical dispersion we introduce is equal to the amount of physical dispersion we see? i.e. are there conditions where the physical dispersion cannot be modelled using numerical dispersion, for this set of equations? Or can it be mathematically shown that numerical dispersion must equal physical dispersion when the grid spacing = ##2l##?

The equations:
\begin{equation}
\frac{\partial h_x}{\partial t} = \frac{\phi_{x-\Delta x/2}(h_{x-\Delta x} - h_x)}{\rho_x \Delta x} + \frac{Ua}{\rho_x\epsilon}(T_{s,x}-T_x)
\end{equation}

\begin{equation}
m_j\frac{dh_j}{dt} = \dot{m}_{j-1}(h_{j-1} - h_j) + UA(T_{S,j}-T_j)
\end{equation}

We are setting this up with a central difference scheme for better numerical accuracy. But the way that we set up these central differences matters too. For fluxes into a cell, we use the 2nd order central difference approximation for the 1st derivative on the left face of the cell and for fluxes leaving the cell, we use the 2nd order central difference approximation for the 1st derivative on the right face of the cell. This is based on knowing the overall flows into and out of the cell on the left and right faces. This approach automatically conserves mass, not only for the differential equation but also numerically in the finite difference scheme. In my judgment, if the finite difference scheme is not set up in the specific way that we have done it, it would be impossible to recognize the relationship we have elucidated between numerical dispersion and actual dispersion. That's what so neat about what we have done. Of course, for this to hold together, the dispersion coefficient has to be expressible as in the groundwater literature (##D=v_ol##).

It is important to recognize that the first order picture of what is happening is provided by the sharp front model. All that the diffusion does is smooth the transitions between zones, and we can see from the results that this is not a major effect. So getting the dispersion included with super accuracy is not going to be an imperative (in my judgment) from an engineering standpoint.

Chestermiller said:

We are setting this up with a central difference scheme for better numerical accuracy. But the way that we set up these central differences matters too. For fluxes into a cell, we use the 2nd order central difference approximation for the 1st derivative on the left face of the cell and for fluxes leaving the cell, we use the 2nd order central difference approximation for the 1st derivative on the right face of the cell. This is based on knowing the overall flows into and out of the cell on the left and right faces. This approach automatically conserves mass, not only for the differential equation but also numerically in the finite difference scheme. In my judgment, if the finite difference scheme is not set up in the specific way that we have done it, it would be impossible to recognize the relationship we have elucidated between numerical dispersion and actual dispersion. That's what so neat about what we have done. Of course, for this to hold together, the dispersion coefficient has to be expressible as in the groundwater literature (##D=v_ol##).

It is important to recognize that the first order picture of what is happening is provided by the sharp front model. All that the diffusion does is smooth the transitions between zones, and we can see from the results that this is not a major effect. So getting the dispersion included with super accuracy is not going to be an imperative (in my judgment) from an engineering standpoint.

Hi Chet, is the claim that the numerical dispersion from the upwind scheme is exactly the same as the physical dispersion in the central differencing scheme, or is the claim only that they are of the same order of magnitude?

We are comparing the numerical dispersion introduced by upwind differencing to the physical dispersion captured by the central differencing scheme.

Given the energy balance:
\begin{equation}
\frac{d(\rho_x h_x)}{dt} = \frac{\phi_{x-\Delta x/2}(\frac{h_{x-\Delta x} + h_x}{2})- \phi_{x+\Delta x/2}(\frac{h_{x+\Delta x} + h_x}{2})}{\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}{\epsilon}(T_s-T)
\end{equation}

Using Taylor series expansion about point x:
\begin{align}
h_{x+\Delta x} &= h_x + \Delta x \frac{\partial h}{\partial x} + \frac{\Delta x^2}{2!}\frac{\partial^2 h}{\partial x^2} + O(\Delta x^3) \\
h_{x-\Delta x} &= h_x - \Delta x \frac{\partial h}{\partial x} + \frac{\Delta x^2}{2!}\frac{\partial^2 h}{\partial x^2} - O(\Delta x^3)
\end{align}

The upwind scheme for the first derivative at x can be expressed as:
##\frac{\partial h}{\partial x} \approx \frac{h_x - h_{x-\Delta x}}{\Delta x}##

Substituting in the Taylor expansion for ##( h_{x-\Delta x})##, we get:
## \frac{\partial h}{\partial x} \approx \frac{\partial h}{\partial x} + \frac{\Delta x}{2} \frac{\partial^2 h}{\partial x^2} + O(\Delta x^2)##

The truncation error, introducing numerical dispersion, is:
##\frac{\Delta x}{2} \frac{\partial^2 h}{\partial x^2} + O(\Delta x^2)##

Matching this with the dispersion term in the 2nd order representation:
\begin{equation}
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}
\end{equation}

Showing equivalence:
\begin{equation}
\frac{\Delta x}{2} \frac{\partial^2 h}{\partial x^2} = 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}
\end{equation}

I may have done this in a roundabout way. Does the above somewhat make sense?

casualguitar said:

Hi Chet, is the claim that the numerical dispersion from the upwind scheme is exactly the same as the physical dispersion in the central differencing scheme, or is the claim only that they are of the same order of magnitude?

To terms of 2nd order accuracy, they are the same when ##\Delta x## and l are related in the way we have identified.

Chestermiller said:

To terms of 2nd order accuracy, they are the same when ##\Delta x## and l are related in the way we have identified.

Thanks Chet. Have I identified this relationship correctly in my post #450 above?

Hi Chet, my apologies please ignore the above as its not correct. I will revert later this morning on this

Hi Chet, am I correct in saying -
"If the truncation error on the advection term in the upwind scheme is the same as the dispersion term in the central differencing scheme (when l = delta x/2), then the numerical dispersion associated with the upwind scheme will be equal to the physical dispersion in the central differencing scheme"

I'm having a bit of trouble proving that they are the same mathematically. Working on it though. Just checking that I have the right idea.

Heres what I've attempted so far:

Central Differencing Scheme (CDS) Dispersion Term:
\begin{equation*}
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}
\end{equation*}

Upwind Scheme (US) Advection Term:
\begin{equation*}
\frac{\phi_{x-\Delta x/2} h_{x-\Delta x} - \phi_{x+\Delta x/2} h_x}{\Delta x}
\end{equation*}

To determine the truncation error for the advection term in the Upwind Scheme, we can use the Taylor Series expansion around x:

Expanding (h_{x-\Delta x}) using the Taylor series:
\begin{equation*}
h_{x-\Delta x} = h_x - \Delta x \frac{\partial h}{\partial x} + \frac{\Delta x^2}{2} \frac{\partial^2 h}{\partial x^2} - \dots
\end{equation*}

Substituting this into the US advection term and simplifying I think should give the truncation error term:
\begin{equation*}
\frac{\Delta x}{2} \frac{\partial^2 h}{\partial x^2}
\end{equation*}

However I havent been able to get this yet. If I did, then I'd equate this truncation error with the dispersion term in CDS when (##l = \frac{\Delta x}{2}##).

With (##l = \frac{\Delta x}{2}##), the dispersion term in CDS becomes:
\begin{equation*}
\frac{\Delta x}{2} \frac{\partial^2 h}{\partial x^2}
\end{equation*}

This term is the same as the truncation error from the US. This indicates that when ##(l = \frac{\Delta x}{2})##, the truncation error in the upwind scheme is equivalent to the physical dispersion modeled by the CDS.

Is this the right track?

casualguitar said:

Hi Chet, am I correct in saying -
"If the truncation error on the advection term in the upwind scheme is the same as the dispersion term in the central differencing scheme (when l = delta x/2), then the numerical dispersion associated with the upwind scheme will be equal to the physical dispersion in the central differencing scheme"

I'm having a bit of trouble proving that they are the same mathematically. Working on it though. Just checking that I have the right idea.

Heres what I've attempted so far:

Central Differencing Scheme (CDS) Dispersion Term:
\begin{equation*}
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}
\end{equation*}

Upwind Scheme (US) Advection Term:
\begin{equation*}
\frac{\phi_{x-\Delta x/2} h_{x-\Delta x} - \phi_{x+\Delta x/2} h_x}{\Delta x}
\end{equation*}

To determine the truncation error for the advection term in the Upwind Scheme, we can use the Taylor Series expansion around x:

Expanding (h_{x-\Delta x}) using the Taylor series:
\begin{equation*}
h_{x-\Delta x} = h_x - \Delta x \frac{\partial h}{\partial x} + \frac{\Delta x^2}{2} \frac{\partial^2 h}{\partial x^2} - \dots
\end{equation*}

Substituting this into the US advection term and simplifying I think should give the truncation error term:
\begin{equation*}
\frac{\Delta x}{2} \frac{\partial^2 h}{\partial x^2}
\end{equation*}

However I havent been able to get this yet. If I did, then I'd equate this truncation error with the dispersion term in CDS when (##l = \frac{\Delta x}{2}##).

With (##l = \frac{\Delta x}{2}##), the dispersion term in CDS becomes:
\begin{equation*}
\frac{\Delta x}{2} \frac{\partial^2 h}{\partial x^2}
\end{equation*}

This term is the same as the truncation error from the US. This indicates that when ##(l = \frac{\Delta x}{2})##, the truncation error in the upwind scheme is equivalent to the physical dispersion modeled by the CDS.

Is this the right track?

I havent yet been able to work this out. Is this a dead end by any chance? I'm not 100% sure that it is the truncation error I should be considering. My apologies for all of the questions on this