
[Sponsors] 
November 5, 2010, 18:24 

#61 
New Member
Duncan Roy van der Heul
Join Date: May 2010
Posts: 19
Rep Power: 8 
Can you tell me where exactly you lose the diagonal dominance?
The stencil for first order upwind scheme is (u, epsilon)>0 Internal points: [u epsilon/h u+2*epsilon/h epsilon/h] First unknown: [0 u+3*epsilon/h epsilon/h] using phi_0+phi_1=2*a Last unknown [uepsilon/h u+3*epsilon/h 0] using phi_J+phi_(J+1)=2*b Only using a central scheme, the scheme loses positivity when Pe_h>2 again, I am a OF novice, but I did get the right solution using Prof. Hjasak files modified for upwind, but maybe my input is wrong. Last edited by Duncan_vdH; November 5, 2010 at 18:44. 

November 5, 2010, 18:59 

#62 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27 
In his case DT is large, that's why you get the right answer.
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

November 5, 2010, 19:04 

#63 
New Member
Duncan Roy van der Heul
Join Date: May 2010
Posts: 19
Rep Power: 8 
I am sorry but no, I also tried the smaller values 0.001. In the upwind scheme you always get an M matrix, irrespective of the value of DT.


November 5, 2010, 19:28 

#64 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27 
Do you use underrelaxation or use an unsteady solver (did you set steadyState in fvScheme?)?
The results I posted above are obtained exactly using Hrv case, with upwind scheme, obtaining a steady solution. The first picture is what I obtain without relaxing... P.S. If you have a case with DT = 0.001, steady, without relaxation, I would like to see it, since I cannot really make it work in OF 1.7.x without relaxing it extremely heavily.
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. Last edited by alberto; November 6, 2010 at 02:46. Reason: Added P.S. 

November 6, 2010, 00:19 

#65 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
Hi all, please let me post again (see post # 54) the spirit which I started this thread:
The main matter is, What one must do if had a downstream BC and is doing upwind? To use the BC value or to see upstream to the cell centroid to take a value for the boundary? If you follow Jasak you have to do one thing (the first one), if you follow Versteeg do the opposite. Which is conceptually correct? Can anybody answer this fundamental question? Then let's talk about other related topics. Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

November 6, 2010, 09:06 

#66 
New Member
Duncan Roy van der Heul
Join Date: May 2010
Posts: 19
Rep Power: 8 
I think you should formulate your requirements, and then decide what to do:
The choice of how you incorporate the BC in the first order upwind scheme should:
A reasoning that 'upwind' is a scheme that should use information from 'upwind' is in my opinion not that valid. You are solving an elliptic equation, that does not have characteristics like a hyperbolic equation does. 'Upwind' is a first order scheme that has positivity irrespective of the mesh peclet number would be my definition. Just my opinion. Unfortunately I cant post my results before monday, when I am back in my office. 

November 6, 2010, 12:35 

#67  
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
Duncan, respect of:
Quote:
I think things would became much more clear if all people in the thread takes only 10 minutes, a pencil and a paper and solve the 5 cells problem proposed in Versteeg book. Then do the same using the info in Hrv. thesis. Systems obtained will be the same, except for the last equation. Assembling only the divergence term: Versteeg S_f*U*(phi_4+phi_5)=0 Hrv S_f*U*phi_4+0*phi_5=S_f*U*BC the issue is not centered in solving but in assembling. Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

November 6, 2010, 21:24 

#68 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 13 
It seems clear (as mud) to me that there is no right or wrong here, just two different ways of approximating the same thing. I am not aware of any physical systems that are even conceptually similar to a fixed value outlet as is being discussed here, so I am not even sure "right" and "wrong" have any meaning in the sense that one approach can be said to approximate reality better than the other. Are the Versteeg boundaries better than the OpenFOAM fixed value boundaries? For most practical applications, probably not. If you want to compare artificial systems, then the two implementations can give different results. But what does this mean given that the systems being compared are artificial? The argument becomes something like "OpenFOAM behaves incorrectly when I try to model nonphysical systems". It all gets a bit philosophical if you ask me.
The OpenFOAM fixedValue boundaries are what they are and do what they say on the tin. This happens not to be identical to what is detailed in some text books, but it is an approach that works just fine most of the time. If you want or need Versteeg type boundaries, they can be implemented with a modicum of effort, but they should not replace the existing boundaries. Both can coexist in the same framework. Continuing to debate the correctness of the different approaches is not going to solve any problems. For those who are interested, the simplest way to implement Versteeglike fixed value outlets, is to make a new boundary derived from fixedValue with valueInternalCoeffs and valueBoundaryCoeffs replaced by those used in zeroGradient boundaries. Of course, these boundaries will not be convection scheme aware and should only be used with pure upwind and on outlets where all flux is going out of the domain. If I implement anything more general at some point in the future, I will leave a notice here. 

November 6, 2010, 22:30 

#69 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
Eugene, I think things are not so polarized, Versteeg way to assemble BC's is the same as FOAM, except for the case of downstream BC and upwind scheme. For this particular case I've found that FOAM solution is wrong for high Pe.
I only want an explanation about why to use downwind info if I'm doing upwind. Some people argues that if you have a BC you must use it, but the fact is that this lead to an unbounded solution. There is a related topic to this one. If one would use node centered FVM, like in Patankar's book, it allows to put the unknown just on the boundary, using a "half cell". In this formulation Hrv's and Verteeg's upwind schemes performs equal and resembles a Finite Difference solution. In cell centered FVM, BC' aren't imposed in real unknown positions but in faces of volumes, so from my point of view the strong sense of a "BC" like in FDM or FEM cannot directly be applied, particularly if the way of reconstructing the face value leads to an unstable scheme. Anyway... Imagine you are teaching and you want to explain why FOAM assembles the downstream BC with upwind in the way what is done. What would you say? (It's very probable that I have to give this explanation in the next semester...) In Versteeg way I have the book and the facts supporting me, at least. Regards. P.S. You said that this issue appears in industrial cases too, didn't you?
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

November 7, 2010, 01:31 

#70  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27 
Quote:
Anyway, I did what you said (took paper, pencil... and dig out some old code, which I should have done before ), and tried the two implementations. Well, results are indeed quite interesting, since with Versteeg implementation, I can easily obtain a solution with Pe' = 100, whatever linear solver I use (I tried GaussSeidel and TDMA), while with the "OpenFOAMlike" implementation, a Pe' of 10 the correct solution is not obtained without altering the conditions in some way. Santiago: do the result I posted above (see plots / codes in one of my previous post) match what you found with OF (you have a private message btw ) ?
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

November 7, 2010, 08:32 

#71 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 13 
Santiago, if you look in the code, Versteeg's way and the OpenFOAM way are not equivalent. They might produce equivalent results most of the time, but there is no consideration of the convection scheme in any noncoupled boundaries in FOAM. This is a fundamentally different approach  the FOAM way says the boundary conditions take precedence over the discretisation scheme, Versteeg's way says it is the other way around. In FOAM, the value on the boundary is thus the same for all operators. In Versteeg, the boundary simultaneously has multiple values. You could easily argue that Versteeg is the one who is confused!
Here is what I would say in a classroom: OpenFOAM does not implement the Versteeg fixed value boundary. It has an alternative approach. In OpenFOAM you are not using upwind or downwind info, in fact you are not considering the interpolation scheme (which is by definition approximate) at all. Rather you are using the more accurate exact information, i.e. the specified fixed value on the boundary for all purposes. What is the justification in Versteeg for using upwind interpolation to obtain a less accurate value at the boundary than the one already supplied? Numerical stability? From my perspective this falls into the realms of expediency and certainly is no more "correct" than the alternative, but neither is it "wrong". I think we must try to make a distinction between 1) Not matching accepted dogma on a subject and 2) Containing a fundamental flaw  i.e. absolute correctness. FOAM does not match the generally accepted way of doing things and thus does not produce results in line with Versteeg. However, assuming that Versteeg and other established sources must be correct and everything that differs from them must be wrong is not necessarily a supportable position. Having the book on your side does not carry that much weight in my book. If someone can unequivocally prove to me that there is only one way to assemble the div operator in boundary cells, then I will gladly put my hands up, but so far all I have seen are "he said/she said" arguments. All instances where "correct" results have been obtained so far is when comparisons have been made to other numerical results using the same implementation. This is a circular argument and doesn't prove anything. Neither do I think you can argue Versteeg is correct because his solutions remain bounded  who says boundedness is an arbiter of correctness? Don't get me wrong, I think we should implement Versteeglike boundaries and that they are appropriate for a good range of problems, where they represent a more convenient and stable approximation than the alternative. Consider for example a porous jump outlet boundary with infinite resistance in the tangential direction (i.e. Utangential = (0 0 0) ). For everyday operating conditions, this outlet will not behave stably in conjunction with upwind differencing using a single value of U on the boundary (irrespective of whether it is provided by a gradient or a fixed value). So you might say "Ah, this proves Versteeg is correct because this real physical system is modelled more accurately by his approach than the FOAM approach." But it does not. Applying the equivalent of infinite resistance over zero space is in no way physical. In reality there is no such thing as a porous jump, it is only an approximation we use to make modelling some systems more convenient. I think it is the same for Versteeg  his div assembly is an approximation we should use to make modelling some systems easier. Saying that it is "correct" for this reason is like saying upwind differencing is correct and central differencing is wrong. That's just how I see it at the moment. 

November 8, 2010, 16:16 

#72  
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
Quote:
Quote:
Quote:
I think the base problem is what I posted before: "There is a related topic to this one. If one would use node centered FVM, like in Patankar's book, it allows to put the unknown just on the boundary, using a "half cell". In this formulation Hrv's and Verteeg's upwind schemes performs equal and resembles a Finite Difference solution. In cell centered FVM, BCs aren't imposed in real unknown positions but in faces of volumes, so from my point of view the strong sense of a "BC" like in FDM or FEM cannot directly be applied, particularly if the way of reconstructing the face value leads to an unstable scheme." Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

November 9, 2010, 13:29 

#73 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 13 
Well I think we are pretty much in agreement then. FOAM imposes boundaries in the strong sense, most other approaches do so weakly to ensure stability for upwind and Pe > 1. This doesn't mean the upwind scheme is different or that it does central differencing in the last cell (although the effect might be equivalent), it just means boundary conditions are applied in the strong sense. This applies to all schemes and all boundary conditions.
I don't see this as right or wrong, just different and as long as the difference is understood, absolutely fine. The problem comes in when you assume the boundaries are applied weakly and then get numerical instabilities because they are applied in a strong sense. So perhaps some information along these lines should be put in a public place and weak alternative boundaries developed for applications that require them. For interest sake, are the results for strong and weak boundaries identical for Pe < 1? 

December 6, 2010, 12:19 

#74 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 13 
Just a quick update on a practical application that benefits from the approach discussed here. I have a boundary similar to Fluent's porous jump outlet. The normal velocity is zero gradient, but the tangential velocity should be zero. I used a fixedGradient boundary with the normal gradient set to zero and the tangential gradient set such that the velocity at the outlet would be zero. Convection scheme for U is upwind. Re_inlet ~ 1000, laminar.
I ran two cases, one using upwind for the gauss convection scheme outlet face value and the other using the standard fixedGradient specification. The difference between the two is simply setting valueBoundaryCoeffs to zero to use the upwind approach. The normal approach produces large discontinuities and unphysical behaviour near the outlet, while the modified formulation produces a smooth and much more plausible solution. (I would upload some pictures, but embarrassingly, I don't know how!) For this kind of boundary the discretisation linked face formulation is thus of critical importance. 

December 7, 2010, 03:12 

#75  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27 
Quote:
Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

December 8, 2010, 19:07 

#76 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 13 
Thanks Alberto. Here are the images in question (I hope!).
The one with the high velocity streak along the bottom (outlet) and left wall is the normal convection assembly formulation. The upwind formulation is just much better. Now admittedly this situation was set up to favour the upwind formulation, but it does illustrate the point. 

December 9, 2010, 14:35 

#77  
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
Quote:
I will prepare more pics solving the problem for upwind scheme using finite difference (node based scheme) to compare. Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

December 11, 2010, 12:06 
problem with scalar when BC fixed at outlet

#78 
Member
George Pichurov
Join Date: Jul 2010
Posts: 39
Rep Power: 8 
As I am solving a SS problem with scalar transport and BC fixedValue at outlet, I can confirm that problem with divergence exists and that changing just the BC to zerogradient fixes things perfectly. The equation solved:
Code:
dimensionedScalar LMRAsource ( "LMRAsource", dimensionSet(0, 0, 1, 0, 0, 0, 0), 1.0 ); fvScalarMatrix LMRAEqn ( fvm::div(phi, LMRA)  fvm::laplacian(kappaEff, LMRA) == LMRAsource ); LMRAEqn.solve(); LMRAEqn.relax(); Code:
LMRA { solver PBiCG; preconditioner DILU; //smoother DILU; tolerance 1e05; relTol 0.1; } It is an issue indeed, because I only use stable cshemes: upwind for div and Gauss linear corrected for Laplace. And I do believe that the matrix coefficients are not composed correct for the near outlet cell(s). It is these cells that receive fantastic values in the process of itteration and flip flop the values in all other cells. This happens even for strong relaxation 0.1. But above all, just changing the BC to fixedgradient makes things sing. 

December 11, 2010, 21:54 

#79 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
George, maybe changing the outlet BC is proper for the description of your particular physical problem, but in the case of the problem of this thread, using another outlet BC is solving a complete different problem. The solution comes by changing the way of BC are applied from 'strong' BCs to 'weak' BCs.
Respect of the inspection of matrix coefficients you have to go with a debugger (gdb for example) into the solve method until reach the point where BC's are imposed and system is finally solved. In case of scalarTransportFoam this point is after line 141 of fvScalarMatrix.C. Then in order to examine the coefficients: p this>diag().v_[start]@end p this>lower().v_[start]@end p this>upper().v_[start]@end where start and end are the starting and ending indices of the piece of vector for diagonal, lower diagonal and upper diagonal that you want to inspect. There are many other ways to do this but it is simple and powerful. Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

December 20, 2010, 04:05 

#80  
Member
George Pichurov
Join Date: Jul 2010
Posts: 39
Rep Power: 8 
Quote:
I have to disagree that fixed values at outlets are to be rendered purely artificial or philosophical. I need to solve a parameter called "residual age of air" in a room, which has neatly defined theretically proven equation (standard convection, diffusion and fixed source terms), but also a very sound physical sense. It shows how much time is left for air molecules until they leave the domain. If you are familiar with the age of air parameter, this one is kind if opposite to it. And of course, due to the pure definition of this parameter, it has to be zero everywhere where it leaves the room. This is precisely the BC I need: zero value at all outlets. The problem I get with OF: unrealistic fluctuations in this and only this parameter, regardles of scheme and solution method. The fluctuations are so wild in the cells adjacent to outlet, that all domain values are ripped off for that parameter. Another example I would give for using fixed outlet BC is the pressure BC. It is (almost) always fixed at outlets. There is no problem in solving the pressure though, because of the SIMPLE and others which create very different final matrix, and that is for pressure correction, not pressure itself. I am still on a quest of composing stable matrix with zero outlet BC. Will post if anything comes around. Last edited by jorkolino; December 20, 2010 at 05:48. 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Warning: Dynamic zone with wrong CG using 6DOF  Manoj Kumar  FLUENT  1  August 11, 2012 04:03 
BuoyantBoussinesqSimpleFoam and axialsymmetric results wrong mass flow  Thomas Baumann  OpenFOAM  6  December 21, 2009 11:31 
udf error  srihari  FLUENT  0  February 9, 2009 10:00 
Pressure contour seems wrong???  Harry Qiu  FLUENT  1  June 29, 2001 05:53 
What's wrong with my UDF?  olivia  FLUENT  1  June 23, 2001 17:06 