Hello Folks, Last week I asked if anyone could advise me how best to solve the linear advection equation on a stretched cartesian grid. I recieved some really helpful comments that helped me improve my current method. Unfortunately I still have not managed to solve this problem to a level I'm happy with.

I have tried several finite difference schemes to solve u_x + u_t = 0;

with u(0,t)=-sin(t) and u(x,0)=sin(x)

Basically my problem is that the grid starts off near the origin with a very small spacing (dx=0.005) and is stretched linearly { dx(I)=1.1*dx(I-1) } until approximately x=500. Unfortunately, because of the way that this problem is coupled with another one I'm solving I cannot change the grid much.

I find that my solutions match well with the analytial solution near the origin but rapidly become innacurate near x=500. The numerical solution in this region is both out of phase and the wrong amplitude when compared to the analytical solution.

I am now thinking that my only option is to solve the advection equation on a fine mesh using a simple scheme and then to intepolate those solutions onto the coarser mesh.

This feels a bit like I'm cheating ) but I'm pretty stuck trying to get decent results using my current work. I wondered what you guys thought of this idea. Cheers, F

Hi Frank,

the natural way of solving your problem would be to transform to a uniform grid on the unit, let's call each grid point there r_i. Assuming N grid points in total we have dr(i) = 1/(N-1).

On this new grid your equation u_t + u_x = 0 becomes u_t + dr(i)/dx(i)*u_r = 0 at each grid point, where dr(i)/dx(i) = (1/(N-1))/(dx(0)*s^i) (with dx(0) = 0.005 and s = 1.1)

So you would solve the standard wave equation on a regular grid but with varying velocity.

Is this wath you're doing?

/Stefan

Hi Stefan, Thanks for your reply. I did try to this method to solve my problem but at the downstream end of the domain I still get a very large discrepancy when comparing to the analytical solution. I wonder if I'm doing something wrong?

How would you suggest dealing with the varying velocity term? I might be doing something really silly in my formulation that would explain the strange results.

Many thanks, F

I would submit that the reason you may be experiencing magnitude & phase discrepancies lies partly in the shape of your mesh & mainly in your wave-timing.

To capture wave phenomena with low-level elements, or stencils, the mesh needs to be even. Mesh distortion actually plays havoc with the wave-nature of the momentum-wave you're trying to simulate. You can see for yourself how a mesh distortion affects the original equation under study, by setting dy=a.dx & plugging back into the original pde.

To capture the momentum-wave decently use the following:

1. Even mesh;

2. Wave timing: dt = (1/2)*(dx/u,in)

Have fun.

desA

I'm afraid that the mesh and time spacing are fixed by another part of my problem. The main problem (unsteady transonic flow around an airfoil) should be the complicated part and the wave equation is a downstream wake condition which SHOULD have been easy to solve )

Do you think, then, that I should define a second mesh and solve the wave equation on the new, finer mesh? I could then interpolate the results onto my coarse mesh. This would mean that I'd get the benefit of being able to calculate the wave equation well without having to use a fine mesh for my complicated transonic problem (which would take ages to calculate). Linear interpolation would be accurate enough for what I'm doing, especially given the wide mesh spacing downstream. I don't think that the interpolation would cost too much computationally - especially when compared to some other areas in my solver.

Many Thanks, F

As your grid stretches, the local dispersion is going to make it impossible to maintain the phase and amplitude of your wave. If those quantities are critical to your project, then your suggestion of using a fine grid and interpolating the solution is probably the best approach.

For advection equation, you need to look at the number of points per wavelength,required to resolve your waves.If you use a second order scheme, the number of points needed to properly resolve the waves is around 20-25 ( rule of thumb).Hence if you do not have enough points,your results will dissipate and disperse.A simple wavenumber analysis of the numerical scheme can reveal a lot of details.Try also to study the group velocity of your numerical scheme,using a simple sin wave with gauss modulation.

Harish wrote:

A simple wavenumber analysis of the numerical scheme can reveal a lot of details. Try also to study the group velocity of your numerical scheme,using a simple sin wave with gauss modulation.

desA writes:

Hi Harish, Can you please provide a little more information on the wavenumber & group-velocity analysis techniques? Would you also perhaps have a few references?

This would be really helpful in the research I'm currently working on. Thanks so much.

desA

In the technique you first take a fourier transform of your numerical schem.Consider a simple example

df/dx=1/dx*(a_j f(j+dx) ) where j can vary from some -M to N ( M=N gives central difference)and take a fourier transform of both sides.It would be

k_a dx =-i(a_j e^(i j k_n dx ))

If d(k_a)/d(k_n) =1 , group velocity of numerical scheme is same as the system and waves are propagated with correct speed.If > 1 waves propagated faster and < 1 waves propagate slower.

For determining the wavenumbers resolved by the scheme what you do is assume a solution of the form e^(ikx) and obtain the RHS of the numerical discretization and plot the rhs for different values of k. The resolved wavenumbers will fall on the straight line. References :

My spelling of the names can be wrong.

1. Trefthen , " Group velocity of finite difference schemes "

2. Whitham " Non linear waves "

3. Tam and Webb , DRP Scheme JCP 1992

4. Lele," Compact scheme with spectral like resolution ", JCP 1992

Thanks very much Harish. Much appreciated.

desA

If you use an optimized DRP scheme, you will get great answers around 7 ppw, and even lower.

This would let you accurately get derivatives for more of the domain. After that, depending on what type of scheme you are using, your wave will start acting very nonphysical. If you are using an unbiased scheme (such as 2nd order central), then theoretically, you won't get any damping. As said before, look at the numerical wavenumber to see why.

But, you wave will disperse quite horribly, meaning that it will propagate at the incorrect speed. Since you are solving the linear Burger's equation, it will globally travel at the incorrect speed since the flux terms is linear (this is rather than the nonlinear case, where the wave error is dependent on the flow value, which causes all sorts of problems).

You could also try transforming the grid to a uniform grid before solving, as said earlier.

But anyways, for more info on DRP schemes, I believe the paper was Tam and Webb from way back in the day.

I messed up in my post. Actually, for this problem, one absolutely needs to transform the grid to accurately solve the problem.

It can be shown that for a nonuniform grid, the wavenumbers change dramatically based on the differing grid spacing. For the second order central differencing scheme, the numerical wavenumber is: i/[dx1 + dx2] * [ (dx2/dx1)sin(k dx1) + (dx1/dx2)sin(k dx2) + i { (dx2/dx1)cos(k dx1) - 1 - (dx1/dx2)cos(k dx2) - 1 } ]e^ikx

Basically, if dx2 is greater than dx1, or your grid is stretching, then the dissipation term (the real part of the wavenumber) will be greater than 1. Or, this means that the error is growing exponentially. You will not be able to get an answer using normal schemes.

The grid will need to be transformed to accurately solve the problem.