Numerically solving a transport equation

In summary, the conversation discussed the use of a numerical model with a downwind approximation for the spatial derivative and the usual approximation for the time derivative. The stencil for the model was shown, along with the boundary conditions at the endpoints. The conversation also touched on the importance of using the correct counters for the different directions and the possibility of a sign error in the equation. The suggestion was made to consider using a first order scheme for a more robust solution.
  • #1
hunt_mat
Homework Helper
1,798
33
TL;DR Summary
I have a simple PDE that I want to solve:
[tex]\frac{\partial v}{\partial t}=-\frac{\partial v}{\partial x}[/tex] along with boundary conditions:
[tex]\frac{\partial v}{\partial x}\Bigg|_{x=0,1}=0[/tex]
I know this has a simple analytical solution, but that's not the point.
I'm using a ``downwind'' approximation for the spatial derivative:
[tex]\frac{\partial v}{\partial x}\approx -\frac{3}{2h}v_{j}+\frac{2}{h}v_{j-1}-\frac{1}{2h}v_{j-2}[/tex]
I'm using the usual approximation for the time derivative, I get the following for a stencil:
[tex]v_{i+1,j}=\left(1+\frac{3\alpha}{2}\right)v_{i,j}-2\alpha v_{i,j-1}+\frac{\alpha}{2}v_{i,j-1}[/tex]
where [itex]\alpha=\delta t/\delta x[/itex]. To deal with the endpoints I evaluate the governing equation at the boundary points to get
[tex]\frac{\partial v}{\partial t}\Bigg|_{x=0,1}=0[/tex], and so this simply makes:
[tex]v_{i+1,1}=v(i,1),\quad v_{i+1,N}=v_{i,N}[/tex]
To deal with the point at [itex]j=2[/itex], I write the BC at [itex]x=0[/itex] as
[tex]\frac{v_{i,1}-v_{i,0}}{\delta x}=0[/tex]
This means that [itex]v_{i,0}=v_{i,1}[/itex]. This can be used in the stencil at [itex]j=2[/itex] to yield
[tex]v_{i+1,2}=\left(1+\frac{3\alpha}{2}\right)v_{i,2}-\frac{3\alpha}{2}v_{i,1}[/tex]
This completes the numerical model. It looks okay to me but it's unstable, and I don't see how. Can anyone suggest anything?
 
Physics news on Phys.org
  • #2
What is the Courant Number you use? (##C = \frac{u\Delta t}{\Delta x}##) Since this is an explicit scheme, it should be lower than one. Also, you want to be 'upwinding' not 'downwinding', i.e. the point you calculate should be dependent on the points upwind from it, since that is what will be influencing that point.
 
  • #3
My Courant number is 0.1. I should use:
[tex]\frac{\partial v}{\partial x}\Bigg|_{x=x_{j}}\approx-\frac{3}{2\delta x}v_{j}+\frac{2}{\delta x}v_{j+1}-\frac{1}{2\delta x}v_{j+1}[/tex]
for my approximation?
 
  • #4
Just as a side remark:
normally the counter ##i## would be used for x-direction, ##j## for the y-direction, ##k## for the z-direction and ##n## for time. Same for the ##\delta## sign. If you mean a small finite distance then use ##\Delta## instead.
It would be very convenient if you use the same counters since that makes everything you write down much more readable. I will definitely use it this way since otherwise I confuse myself.

I based my remark on your own remark in that you use the 'downwind' approximation. But it really depends on the sign of ##a## in ##\frac{\partial v}{\partial t} + a \frac{\partial v}{\partial x} = 0##. In your case ##a=1## and thus the upwind side is at the 'left' side of the query point (i.e. where x is lower). So if indeed ##x_{i-1} < x_{i}## then you are already using upwind, although you called it downwind.

Actually, it seems that you have a sign error in the equation below:
hunt_mat said:
I'm using a ``downwind'' approximation for the spatial derivative:
[tex]\frac{\partial v}{\partial x}\approx -\frac{3}{2h}v_{j}+\frac{2}{h}v_{j-1}-\frac{1}{2h}v_{j-2}[/tex]

Since for a second order scheme it would look like this (which is the negative of what you got):
[tex]
\frac{\partial v}{\partial x}\Bigg|_{x=x_{i}}\approx \frac{3v_{i} -4v_{i-1} +v_{i-2}}{2\Delta x}
[/tex]Lastly, you now use a second order spatial derivative. If all fails, go back to a first order scheme since that is usually more robust.
 
Back
Top