[CerpinTaxt]

Hey guys!

I have programmed the first order upwind scheme on a cartesian, equidistant grid to solve the stationary continuity equation
nabla(v*z(x,y))=0
in 2D for a given velocity field.

To test my implementation, I used the velocity field v=(vx,vy)=(y,-x). The analytical solution is therefore z(x,y)=1/2*(x^2+y^2). I wanted to solve it in the domain [-1,1]x[-1,1] and created the inflow boundary condition from the analytical solution.
My problem is, that the upwind scheme seems to break down in the presence of a rotational symmetry in the veloctiy field - the solution gets constant in a circular area around the center of the vortex and the solution does not converge with a finer grid.
Even if I shift the domain for example to [-0.5, 1.5]x [-1.1] or the velocity field, the solution breaks down around the vortex center.

The only possibility is to break the symmetry by solving the equation for example on [0,1]x[-1,1]. Than the solution is okay and converges with a finer grid.
My code has no problems at all with velocity fields without rotational symmetry like (vx,vy)=(y,x), which I've tested with the corresponding analytical solution.

Unfortionately I want to use the upwind scheme for more complex velocity fields which can have several vortices at random positions.
I've read in a few CFD books, which all say, that the upwind scheme suffers from numerical diffusion, but this doesn't seems to be the problem, since the solution doesn't converge with a finer grid.

Since it is difficult to find information about my problem in the literature, I wanted to ask you guys, if this problem is maybe inherent to the first order upwind scheme and if there are maybe better discretizaton schemes for stationary veloctiy field with vortices.

Thank you for your help!
CerpinTaxt

[FMDenaro]

Originally Posted by CerpinTaxt

Hey guys!

I have programmed the first order upwind scheme on a cartesian, equidistant grid to solve the stationary continuity equation
nabla(v*z(x,y))=0
in 2D for a given velocity field.

To test my implementation, I used the velocity field v=(vx,vy)=(y,-x). The analytical solution is therefore z(x,y)=1/2*(x^2+y^2). I wanted to solve it in the domain [-1,1]x[-1,1] and created the inflow boundary condition from the analytical solution.
My problem is, that the upwind scheme seems to break down in the presence of a rotational symmetry in the veloctiy field - the solution gets constant in a circular area around the center of the vortex and the solution does not converge with a finer grid.
Even if I shift the domain for example to [-0.5, 1.5]x [-1.1] or the velocity field, the solution breaks down around the vortex center.

The only possibility is to break the symmetry by solving the equation for example on [0,1]x[-1,1]. Than the solution is okay and converges with a finer grid.
My code has no problems at all with velocity fields without rotational symmetry like (vx,vy)=(y,x), which I've tested with the corresponding analytical solution.

Unfortionately I want to use the upwind scheme for more complex velocity fields which can have several vortices at random positions.
I've read in a few CFD books, which all say, that the upwind scheme suffers from numerical diffusion, but this doesn't seems to be the problem, since the solution doesn't converge with a finer grid.

Since it is difficult to find information about my problem in the literature, I wanted to ask you guys, if this problem is maybe inherent to the first order upwind scheme and if there are maybe better discretizaton schemes for stationary veloctiy field with vortices.

Thank you for your help!
CerpinTaxt

I don't understand your "continuity" equation that should write simply Div v = 0.

I suggest a simple modification in your test. You can use the same velocity field (a rigid rotation at angular velocity = -1) in the domain [-1,1]x[-1,1].
However, solve for the time-dependent convective equation

d z/dt + v .grad z = 0


while prescribing an initial condizion z(x,y,0) given by a cosine hill having a center at (0,-0.5) and radius of 0.2. This way you do not need inflow condition to be specified and you can check only for the upwind discretization.
Then plot the solution at T = pi, 2*pi and post the imagines

[CerpinTaxt]

Hi!

Thanks for your answer!

If I interpret your answer right, you think, that my implementation is wrong and that it is not a general problem of the first order upwind scheme?!?

Unfortunately i am not able to upload pictures at the moment, so maybe it helps if I show some of the intermediate results of my example I described above, here for a 6x6 grid, to see for you if at least the structure is correct.

Description:

Ap, An, Ae, As, Aw correspond to the usual definitions of the matrix coefficients of the equation system at each grid point like you can find it in the book of Patankar. So Ap is written into the main diagonal, Aw, An, Ae and As in the other four diagonals.
Sp contains the boundary values, is written into a vector and corresponds to the RHS of the equation system. And z is the solution of the equation system.


The data:

An =

-0.2778 -0.1667 -0.0556 0 0 0
-0.2778 -0.1667 -0.0556 0 0 0
-0.2778 -0.1667 -0.0556 0 0 0
-0.2778 -0.1667 -0.0556 0 0 0
-0.2778 -0.1667 -0.0556 0 0 0
0 0 0 0 0 0


Ae =

0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
-0.0556 -0.0556 -0.0556 -0.0556 -0.0556 0
-0.1667 -0.1667 -0.1667 -0.1667 -0.1667 0
-0.2778 -0.2778 -0.2778 -0.2778 -0.2778 0


Aw =

0 -0.2778 -0.2778 -0.2778 -0.2778 -0.2778
0 -0.1667 -0.1667 -0.1667 -0.1667 -0.1667
0 -0.0556 -0.0556 -0.0556 -0.0556 -0.0556
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0


Ap =

0.5556 0.4444 0.3333 0.3333 0.4444 0.5556
0.4444 0.3333 0.2222 0.2222 0.3333 0.4444
0.3333 0.2222 0.1111 0.1111 0.2222 0.3333
0.3333 0.2222 0.1111 0.1111 0.2222 0.3333
0.4444 0.3333 0.2222 0.2222 0.3333 0.4444
0.5556 0.4444 0.3333 0.3333 0.4444 0.5556


As =

0 0 0 0 0 0
0 0 0 -0.0556 -0.1667 -0.2778
0 0 0 -0.0556 -0.1667 -0.2778
0 0 0 -0.0556 -0.1667 -0.2778
0 0 0 -0.0556 -0.1667 -0.2778
0 0 0 -0.0556 -0.1667 -0.2778


sp =

0.2353 0 0 0.0285 0.1042 0.2353
0.1042 0 0 0 0 0
0.0285 0 0 0 0 0
0 0 0 0 0 0.0285
0 0 0 0 0 0.1042
0.2353 0.1042 0.0285 0 0 0.2353



z =

0.7522 0.7183 0.7091 0.6765 0.6572 0.7522
0.6572 0.6618 0.6629 0.6663 0.6618 0.7183
0.6765 0.6663 0.6663 0.6663 0.6629 0.7091
0.7091 0.6629 0.6663 0.6663 0.6663 0.6765
0.7183 0.6618 0.6663 0.6629 0.6618 0.6572
0.7522 0.6572 0.6765 0.7091 0.7183 0.7522

analytical soultion

z_ana =

0.6944 0.4722 0.3611 0.3611 0.4722 0.6944
0.4722 0.2500 0.1389 0.1389 0.2500 0.4722
0.3611 0.1389 0.0278 0.0278 0.1389 0.3611
0.3611 0.1389 0.0278 0.0278 0.1389 0.3611
0.4722 0.2500 0.1389 0.1389 0.2500 0.4722
0.6944 0.4722 0.3611 0.3611 0.4722 0.6944


Thank you for your help!

[aerosjc]

Originally Posted by FMDenaro

I don't understand your "continuity" equation that should write simply Div v = 0.

I suggest a simple modification in your test. You can use the same velocity field (a rigid rotation at angular velocity = -1) in the domain [-1,1]x[-1,1].
However, solve for the time-dependent convective equation

d z/dt + v .grad z = 0


while prescribing an initial condizion z(x,y,0) given by a cosine hill having a center at (0,-0.5) and radius of 0.2. This way you do not need inflow condition to be specified and you can check only for the upwind discretization.
Then plot the solution at T = pi, 2*pi and post the imagines

I think that he just wants to solve a boundary problem. The governing equation is \nabla ( V(x,y) \rho (x, y) ) = 0 of which V(x, y) is V = (y, -x).

[aerosjc]

What is the boundary conditions?

[CerpinTaxt]

Originally Posted by aerosjc

I think that he just wants to solve a boundary problem. The governing equation is \nabla ( V(x,y) \rho (x, y) ) = 0 of which V(x, y) is V = (y, -x).

Exactly! So z(x,y) in my case can be interpreted as a mass density.

[CerpinTaxt]

Originally Posted by aerosjc

What is the boundary conditions?

I've constructed the boundary condition from the analytical solution z(x,y)=1/2*(x^2+y^2) to test the code.
So for x=-1, I take z_B=z(-1,y), y=-1..1;
x=1: z_B=z(1,y), y=-1..1;
y=-1: z_B=z(x,-1), x=-1..1;
y=1: z_B=z(x,1), x=-1..1;

[FMDenaro]

Originally Posted by CerpinTaxt

I've constructed the boundary condition from the analytical solution z(x,y)=1/2*(x^2+y^2) to test the code.
So for x=-1, I take z_B=z(-1,y), y=-1..1;
x=1: z_B=z(1,y), y=-1..1;
y=-1: z_B=z(x,-1), x=-1..1;
y=1: z_B=z(x,1), x=-1..1;

that is mathematical well posed for elliptic equations, for hyperbolic equations you can not set BC.s on all boundaries ...

[CerpinTaxt]

Originally Posted by FMDenaro

that is mathematical well posed for elliptic equations, for hyperbolic equations you can not set BC.s on all boundaries ...

With an upwind scheme, the boundary conditions i mentioned are simple inflow boundary conditions, since the boundary values at the outflow part simply doesn't matter (so you can set them to zero, if you want - it will have no effect).
And the stationary continuity equation with inflow boundary conditions is a mathematically well posed problem!