# Implementation of 2D advection equation with second order upwind scheme?

43 views (last 30 days)
Yoni Verhaegen -WE-1718- on 5 May 2022
Hi all,
I am trying to numerically discretize a 2D advection equation to model the transport of rocks with thickness (h_debris) on top of glacier ice with velocity components (velx_mod and vely_mod). The equation to be discretized is as follows:
dh/dt = -d(u*h)/dx - d(v*h)/dy = -u(dh/dx) - v(dh/dy) - h(du/dx) - h(dv/dy)
% in the model:
% u = velx_mod (x component of ice surface velocity)
% v = vely_mod (y component of ice surface velocity)
% h = h_debris (thickness of the rock layer)
I therefore did the following for a second order upwind scheme in 2D:
if velx_mod(i,j) < 0
term1_debris_x(i,j) = -velx_mod(i,j).*((-h_debris(i+2,j)+4*h_debris(i+1,j)-3*h_debris(i,j))./(2*deltax_d)) - h_debris(i,j).*((-velx_mod(i+2,j)+4*velx_mod(i+1,j)-3*velx_mod(i,j))./(2*deltax_d));
elseif velx_mod(i,j) >= 0
term1_debris_x(i,j) = -velx_mod(i,j).*((3*h_debris(i,j)-4*h_debris(i-1,j)+h_debris(i-2,j))./(2*deltax_d)) - h_debris(i,j).*((3*velx_mod(i,j)-4*velx_mod(i-1,j)+velx_mod(i-2,j))./(2*deltax_d));
end
if vely_mod(i,j) < 0
term1_debris_y(i,j) = -vely_mod(i,j).*((-h_debris(i,j+2)+4*h_debris(i,j+1)-3*h_debris(i,j))./(2*deltax_d)) - h_debris(i,j).*((-vely_mod(i,j+2)+4*vely_mod(i,j+1)-3*vely_mod(i,j))./(2*deltax_d));
elseif vely_mod(i,j) >= 0
term1_debris_y(i,j) = -vely_mod(i,j).*((3*h_debris(i,j)-4*h_debris(i,j-1)+h_debris(i,j-2))./(2*deltax_d)) - h_debris(i,j).*((3*vely_mod(i,j)-4*vely_mod(i,j-1)+vely_mod(i,j-2))./(2*deltax_d));
end
h_debris(i,j) = h_debrisini(i,j) + deltat_d*(term1_debris_x(i,j) + term1_debris_y(i,j));
I am not an expert in this field and this is the first time I implement this advection equation in 2D. Could someone check my code and assure I did the right thing here? Thanks already.
Yoni Verhaegen -WE-1718- on 6 May 2022
Thank you for your input, it really helps me a lot and I appreciate it. So if the velocity field were to change over time, mny original implementation (as stated in the code of my question here), would be correct?

Torsten on 6 May 2022
Edited: Torsten on 6 May 2022
Formally, your implementation is correct as it is (I did not check the one-sided difference formulae in detail).
The problems I see are:
1. the stability of your scheme (usually, 2nd order schemes which are not stabilized (as your scheme) tend to become instable/oscillatory). That's why you should start with first-order upwind or better: with my suggestion from below.
2. the velocity field. Usually, a velocity field is divergence free to ensure mass-conservation (i.e. the condition du/dx + dv/dy = 0 automatically holds). In this case, your terms h*du/dx + h*dv/dy automatically cancel out and you are left with the equation dh/dt + u*dh/dx + v*dh/dy = 0. You should test this before starting. If you don't get du/dx + dv/dy approximately 0, you must be very careful with the results you get from your simulation.
3. The boundary condition and where you have to set it.
4. The implementation of the scheme. Note that one cell apart from an inflow boundary, you must use a value for h beyond the boundary points (h_debris(i+2,j),h_debris(i-2,j),(-h_debris(i,j+2),h_debris(i,j-2)). So, an extrapolation is necessary.
Yoni Verhaegen -WE-1718- on 6 May 2022
Unfortunately, I am not an expert in this field, so I don't really know how I neede to check and handle these issues. I put a lot of effort in it already to get it where it is now, and for that I am really happy that you have been helping me. Would it be an idea to share the model code somehow so that we can figure it out together?

### Categories

Find more on MATLAB Compiler in Help Center and File Exchange

R2022a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!