# How to determine face velocity for upwind scheme?

 Register Blogs Members List Search Today's Posts Mark Forums Read

 August 28, 2012, 12:44 How to determine face velocity for upwind scheme? #1 Senior Member   TWB Join Date: Mar 2009 Posts: 138 Rep Power: 10 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: 3,562
Rep Power: 40
Quote:
 Originally Posted by quarkz 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

 August 28, 2012, 14:38 #3 Senior Member   TWB Join Date: Mar 2009 Posts: 138 Rep Power: 10 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.

August 28, 2012, 14:44
#4
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 3,562
Rep Power: 40
Quote:
 Originally Posted by quarkz 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

 August 28, 2012, 15:16 #5 Senior Member   TWB Join Date: Mar 2009 Posts: 138 Rep Power: 10 Ok I got it now. Thanks for the clarifications!

 August 28, 2012, 15:39 #6 Senior Member   TWB Join Date: Mar 2009 Posts: 138 Rep Power: 10 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?

August 28, 2012, 16:14
#7
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 3,562
Rep Power: 40
Quote:
 Originally Posted by quarkz 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

 September 1, 2012, 09:34 #8 Senior Member   TWB Join Date: Mar 2009 Posts: 138 Rep Power: 10 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?

September 1, 2012, 11:00
#9
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 3,562
Rep Power: 40
Quote:
 Originally Posted by quarkz 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.

September 2, 2012, 09:21
#10
Senior Member

Join Date: Aug 2011
Posts: 271
Rep Power: 9
Quote:
 Originally Posted by quarkz There's something I don't understand abt the formula. 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?
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.

 September 3, 2012, 09:28 #11 Senior Member   TWB Join Date: Mar 2009 Posts: 138 Rep Power: 10 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.

 September 3, 2012, 09:58 #12 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 3,562 Rep Power: 40 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: 271
Rep Power: 9
Quote:
 Originally Posted by quarkz 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.

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: 138 Rep Power: 10 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 Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post hsingtzu OpenFOAM Native Meshers: blockMesh 2 March 14, 2012 10:56 anand_30 OpenFOAM Meshing & Mesh Conversion 12 December 12, 2011 05:16 nickvinn Main CFD Forum 0 May 2, 2011 10:05 ashokme FLUENT 0 July 31, 2009 16:17 Ralf Schmidt FLUENT 1 April 24, 2009 09:37

All times are GMT -4. The time now is 12:01.