How to determine face velocity for upwind scheme?
 

Hi,

I'm new to upwind schemes. In calculating the face value using the upwind scheme, we need to determine which direction is this face value, before we can choose which values to be used for the interpolation.

But the direction and face value is currently unknown. So how do we go about it?

Do we use the direction of the face value at the previous timestep? Or do we do a simple interpolation using P and E to get the value at "e" and hence the direction?

Btw, I'm trying to implement it to a finite volume solver. If I need to calculate the eastern convective flux term c_flux = face_vel*u_e*area, are both the face_vel and u_e calculated according to the upwind scheme, or do I only use upwind scheme for the face_vel and central interpolation for u_e?

Thanks!

are you talking of first-order upwind? Then for the face "e", between nodes E and P, you will get:

u_e = 0.5*(u_P + u_E)
u_minus = 0.5*(u_e-abs(u_e))
u_plus = 0.5*(u_e+abs(u_e))
flux_e= u_minus*phi_E + u_plus*phi_P

and it automatically works

ok thanks FMDenaro, I'm starting to understand better. So phi_E/P can be u,v, is that so? However, shouldn't it be

if u_e > 0

flux_e = 0.5*(phi_E + phi_P) * u_P

else

flux_e = 0.5*(phi_E + phi_P) * u_E ?

I thought the convective velocity (u_E or u_P) should be the value which changes depending on the flow direction.

Also, is there a similar formula for 2nd-order upwind case? It makes programming much easier.

The concept of upwind involves the transported variable, thus phi must be up-winded, you cannot interpolate with a down-wind value!

The scheme can be generalized to higher order, for example you can see in literature QUICK and QUICKEST schemes by Leonard

Ok I got it now. Thanks for the clarifications!

Btw, for the west face, the formulas above must change to :

u_w = 0.5*(u_P + u_W)
u_minus = 0.5*(u_w-abs(u_w))
u_plus = 0.5*(u_w+abs(u_w))
flux_w= u_minus*phi_P + u_plus*phi_W

Is that so?

Yes, I suggest to compute first the fluxes on all the faces, then you can update the solution at time n+1 while exploting the sum of the fluxes

There's something I don't understand abt the formula.

u_e = 0.5*(u_P + u_E)
u_minus = 0.5*(u_e-abs(u_e))
u_plus = 0.5*(u_e+abs(u_e))
flux_e= u_minus*phi_E + u_plus*phi_P

if u_e<0, u_minus = -u_e, u_plus = 0.

flux_e = -u_e*phi_E

However, shouldn't flux_e = u_e*phi_E? If the negative sign is included, then the flux will now point in the wrong direction, right?

I have to consider that the sign of the flux is provided by the dot product between the normal unit to the surface and the velocity vector...Therefore, in the 1D case at section e:

flux_e = (n.v * phi )_e
being n=+i (unit vector along x) you get:
flux_e = u_e*phi_e is computed as u_minus*phi_E + u_plus*phi_P. If u_e<0 the flux enters from the right to the left and involves the value phi_E.

Note that at section w you have n=-i

I suggest defining a vector flux_x(i), defined in staggered way at location i-1/2 with respect to Phi(i). In such a way you compute first all fluxes flux_x(1:n), then

Phi_new(i) = Phi(i) - dt*(flux_x(i+1)-flux(i))/h

and the conservative property is automatically ensured.

You are wrong my friend ;)
If u_e<0, then |u_e| = -u_e then u_minus = 0.5(u_e- (-u_e)) =u_e , u_plus = 0.
thus flux_e = +u_e*phi_E

the formulation of Filippo is correct, only the surface S_e is missing in the flux which is here dy. If you are thinking in 1D it is correct.

Ops ... tks for pointing out my error. Guess those +/- signs is making me confused.

Also, for the west flux, assuming the direction is included, it should be

u_w = 0.5*(u_P + u_W)
u_minus = 0.5*(u_w-abs(u_w))
u_plus = 0.5*(u_w+abs(u_w))
flux_w= -(u_minus*phi_P + u_plus*phi_W)

The minus sign is due to n = -i.

I do not suggest to replicate the computation for sections e and w in each volume... work first sweeping faces and store the flux values, the use the unicity of the flux on each face

Yes it is correct but how Filippo mentionned it , computing all fluxes on a given control volume is not the right and efficient way to proceed. I mean do not compute flux_e AND flux_w as he stated it.
Just compute all flux_e then get all flux_w saying that
flux_w(P) = -flux_e(W)
The advantage is that you compute stuff only once, and you assure conservation up to round off error.

In addition if you do not want to struggle with the sign +/- use the normal vector and dot product in your formulation.

Start to compute mass flow rate m_e = ro_e*u_e.n_e*S_e with u_e = 0.5*(u_P + u_E)
m_w = ro_w*u_w.n_w*S_w with u_w = 0.5*(u_P + u_W)

here u_e.n_e is the dot product between velocity vector and normal vector located on _e.
n_e=+i , n_w=-i

after the flux_e = max(m_e,0)Phi_P + min(m_e, 0)Phi_E
flux_w = max(m_w,0)Phi_P + min(m_w, 0)Phi_W check the formulae is the same for 'e' and 'w' without changing anything with sign stuffs !!! ;)

that's all.

But in this case use also the sweeping procedure and compute only once the 'e' fluxes on all vertical faces
( and you compute 'n' fluxes on all horizontal faces for 2D problems
and 't' fluxes on the third direction for 3D problems )

Thanks leflix and FMDenaro! I'm also only computing flux in the forward direction since the west flux is just the negative of the east. However, i just want to be sure my concept is correct.

I tested my code with the upwind instead of central scheme for a 2D unsteady flow past a stationary cylinder at Re = 185.

I found that the ans at the same resolution as the central scheme is not correct. I try to reduce to half the resolution and the ans get much better and closer, although some difference still exist.

However, when I tried to simulate a similar case but with a vertically moving cylinder, the ans is different, even after using half the resolution. Is this expected from a 1st order upwind scheme?