CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

How to determine face velocity for upwind scheme?

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

Reply
 
LinkBack Thread Tools Display Modes
Old   August 28, 2012, 13:44
Default How to determine face velocity for upwind scheme?
  #1
Senior Member
 
TWB
Join Date: Mar 2009
Posts: 105
Rep Power: 7
quarkz is on a distinguished road
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 14:34.
quarkz is offline   Reply With Quote

Old   August 28, 2012, 15:02
Default
  #2
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,275
Rep Power: 17
FMDenaro will become famous soon enough
Quote:
Originally Posted by quarkz View Post
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
FMDenaro is offline   Reply With Quote

Old   August 28, 2012, 15:38
Default
  #3
Senior Member
 
TWB
Join Date: Mar 2009
Posts: 105
Rep Power: 7
quarkz is on a distinguished road
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.
quarkz is offline   Reply With Quote

Old   August 28, 2012, 15:44
Default
  #4
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,275
Rep Power: 17
FMDenaro will become famous soon enough
Quote:
Originally Posted by quarkz View Post
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
FMDenaro is offline   Reply With Quote

Old   August 28, 2012, 16:16
Default
  #5
Senior Member
 
TWB
Join Date: Mar 2009
Posts: 105
Rep Power: 7
quarkz is on a distinguished road
Ok I got it now. Thanks for the clarifications!
quarkz is offline   Reply With Quote

Old   August 28, 2012, 16:39
Default
  #6
Senior Member
 
TWB
Join Date: Mar 2009
Posts: 105
Rep Power: 7
quarkz is on a distinguished road
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?
quarkz is offline   Reply With Quote

Old   August 28, 2012, 17:14
Default
  #7
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,275
Rep Power: 17
FMDenaro will become famous soon enough
Quote:
Originally Posted by quarkz View Post
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
FMDenaro is offline   Reply With Quote

Old   September 1, 2012, 10:34
Default
  #8
Senior Member
 
TWB
Join Date: Mar 2009
Posts: 105
Rep Power: 7
quarkz is on a distinguished road
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?
quarkz is offline   Reply With Quote

Old   September 1, 2012, 12:00
Default
  #9
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,275
Rep Power: 17
FMDenaro will become famous soon enough
Quote:
Originally Posted by quarkz View Post
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.
FMDenaro is offline   Reply With Quote

Old   September 2, 2012, 10:21
Default
  #10
Senior Member
 
Join Date: Aug 2011
Posts: 251
Rep Power: 6
leflix is on a distinguished road
Quote:
Originally Posted by quarkz View Post
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.
leflix is offline   Reply With Quote

Old   September 3, 2012, 10:28
Default
  #11
Senior Member
 
TWB
Join Date: Mar 2009
Posts: 105
Rep Power: 7
quarkz is on a distinguished road
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.
quarkz is offline   Reply With Quote

Old   September 3, 2012, 10:58
Default
  #12
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,275
Rep Power: 17
FMDenaro will become famous soon enough
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
FMDenaro is offline   Reply With Quote

Old   September 3, 2012, 11:29
Default
  #13
Senior Member
 
Join Date: Aug 2011
Posts: 251
Rep Power: 6
leflix is on a distinguished road
Quote:
Originally Posted by quarkz View Post
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 )
leflix is offline   Reply With Quote

Old   September 6, 2012, 06:28
Default
  #14
Senior Member
 
TWB
Join Date: Mar 2009
Posts: 105
Rep Power: 7
quarkz is on a distinguished road
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?
quarkz is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


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 11: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 11:05
regarding velocity macro for interior face ashokme FLUENT 0 July 31, 2009 17:17
determine face position Ralf Schmidt FLUENT 1 April 24, 2009 10:37


All times are GMT -4. The time now is 16:52.