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! 
Quote:
are you talking of firstorder 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_eabs(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 2ndorder upwind case? It makes programming much easier. 
Quote:
The concept of upwind involves the transported variable, thus phi must be upwinded, you cannot interpolate with a downwind 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_wabs(u_w)) u_plus = 0.5*(u_w+abs(u_w)) flux_w= u_minus*phi_P + u_plus*phi_W Is that so? 
Quote:

There's something I don't understand abt the formula.
u_e = 0.5*(u_P + u_E) u_minus = 0.5*(u_eabs(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? 
Quote:
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 i1/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. 
Quote:
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_wabs(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

Quote:
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? 
All times are GMT 4. The time now is 20:50. 