
[Sponsors] 
August 28, 2012, 12:44 
How to determine face velocity for upwind scheme?

#1 
Senior Member
TWB
Join Date: Mar 2009
Posts: 118
Rep Power: 9 
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! Last edited by quarkz; August 28, 2012 at 13:34. 

August 28, 2012, 14:02 

#2  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,513
Rep Power: 31 
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 

August 28, 2012, 14:38 

#3 
Senior Member
TWB
Join Date: Mar 2009
Posts: 118
Rep Power: 9 
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. 

August 28, 2012, 14:44 

#4  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,513
Rep Power: 31 
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 

August 28, 2012, 15:16 

#5 
Senior Member
TWB
Join Date: Mar 2009
Posts: 118
Rep Power: 9 
Ok I got it now. Thanks for the clarifications!


August 28, 2012, 15:39 

#6 
Senior Member
TWB
Join Date: Mar 2009
Posts: 118
Rep Power: 9 
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? 

August 28, 2012, 16:14 

#7 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,513
Rep Power: 31 
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


September 1, 2012, 09:34 

#8 
Senior Member
TWB
Join Date: Mar 2009
Posts: 118
Rep Power: 9 
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? 

September 1, 2012, 11:00 

#9  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,513
Rep Power: 31 
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. 

September 2, 2012, 09:21 

#10  
Senior Member
Join Date: Aug 2011
Posts: 251
Rep Power: 8 
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. 

September 3, 2012, 09:28 

#11 
Senior Member
TWB
Join Date: Mar 2009
Posts: 118
Rep Power: 9 
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. 

September 3, 2012, 09:58 

#12 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,513
Rep Power: 31 
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


September 3, 2012, 10:29 

#13  
Senior Member
Join Date: Aug 2011
Posts: 251
Rep Power: 8 
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 ) 

September 6, 2012, 05:28 

#14 
Senior Member
TWB
Join Date: Mar 2009
Posts: 118
Rep Power: 9 
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? 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
error message with modeling a cube with a hold at the center  hsingtzu  OpenFOAM Native Meshers: blockMesh  2  March 14, 2012 10:56 
mesh airfoil NACA0012  anand_30  OpenFOAM Meshing & Mesh Conversion  12  December 12, 2011 05:16 
how to determine friction velocity?  nickvinn  Main CFD Forum  0  May 2, 2011 10:05 
regarding velocity macro for interior face  ashokme  FLUENT  0  July 31, 2009 16:17 
determine face position  Ralf Schmidt  FLUENT  1  April 24, 2009 09:37 