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

조회 수: 36 (최근 30일)
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));
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));
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.
  댓글 수: 8
Torsten 2022년 5월 5일
편집: Torsten 2022년 5월 5일
I have the feeling that keeping velocity constant over time forbids to take into account the terms -h(du/dx) - h(dv/dy) of your equation because thickening or thinning of the rock layer due to compressional or extensional ice flow will itself influence the velocity field. But this was not your question.
So concerning your discretization: the discretization you suggest is a 2nd order upwind method, but maybe not stable since your velocity fields look quite diffuse. I suggest starting with the usual 1st-order upwind scheme. This produces numerical diffusion, but is stable. Further I suggest using (in the case u,v >= 0)
d(u*h)/dx ~ ((u*h)(i)-(u*h)(i-1))/dx instead of u(i)*(h(i)-h(i-1))/dx + h(i)*(u(i)-u(i-1))/dx
d(v*h)/dy ~ ((v*h)(j)-(v*h)(j-1))/dy instead of v(j)*(h(j)-h(j-1))/dy + h(j)*(v(j)-v(j-1))/dy
if you really want to solve
dh/dt + d(u*h)/dx +d(v*h)/dy = 0
instead of
dh/dt + u*dh/dx + v*dh/dy
as I assume to be more correct if u and v are prescribed and kept constant over time (as noted above).
For a stable 2nd order upwind scheme that can cope with discontinuities in u and v see e.g.
Furthermore, you will have to think about the boundary conditions for h. For u,v >=0, they have to be prescribed at x=xmin, ymin<=y<=ymax and xmin <=x<=xmax, y=ymin, e.g.
Yoni Verhaegen -WE-1718-
Yoni Verhaegen -WE-1718- 2022년 5월 6일
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 2022년 5월 6일
편집: Torsten 2022년 5월 6일
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.
I strongly recommend using clawpack for your task where a second-order upwind implementation for the 2-dimensional advection equation is already implemented:
  댓글 수: 5
Torsten 2022년 5월 6일
편집: Torsten 2022년 5월 6일
Then your equation must include a source term s in its formulation:
dh/dt = -d(u*h)/dx - d(v*h)/dy + s(x,y)
And did you check the validity of the global balance :
integral integral h(x,y,t) dx dy - integral integral h(x,y,0) dx dy = t* integral integral s(x,y) dx dy
Yoni Verhaegen -WE-1718-
Yoni Verhaegen -WE-1718- 2022년 5월 6일
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?

댓글을 달려면 로그인하십시오.

추가 답변 (0개)


Help CenterFile Exchange에서 Eigenvalue Problems에 대해 자세히 알아보기




Community Treasure Hunt

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

Start Hunting!

Translated by