
[Sponsors] 
November 3, 2010, 16:15 

#41 
Senior Member
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,768
Rep Power: 21 
I am getting a bit sick of this discussion: WHAT ON EARTH ARE YOU TALKING ABOUT???
So, I took a 1D convectiondiffusion equation and solved it in 1D, with U = 1 and DT = 0.01. The boundary conditions are T = 1 at inlet (left) and T = 2 at outlet (right). Aaaaand: the solution is attached, not a single problem, no single issue, nothing: it just works! So: what are you talking about? and Why do I have to waste my time with this? Hrv
__________________
Hrvoje Jasak Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk 

November 3, 2010, 16:32 

#42 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 420
Rep Power: 15 
Hrv, this is we were talking about:
Test case explanation and exact solution is in #1 (FOAM case attached), examples of output in #2, correct formulation by Prof. Versteeg in #37. BTW, What are your mesh step (in order to calculate Pe number) and div schemes? Regards. P.S. Another previous independent post on the same issue Convectiondiffusion in 1D : wrong solution for a large Delta x
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar 

November 3, 2010, 16:34 

#43 
New Member
Duncan Roy van der Heul
Join Date: May 2010
Posts: 15
Rep Power: 7 
Dear Hrvoje,
I am really happy that you do get the right solution to the problem and was really surprised to initially read here that it was not, as I told you this afternoon (really enjoyed the course by the way). So apparantly, the original posters input (and also my 'total novice OpenFOAM user' understanding of it) does not correspond to your input. Would you care to share the files? Kind regards, Duncan 

November 3, 2010, 17:26 

#44 
Senior Member
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,768
Rep Power: 21 
The length is 1, number of cells is 100, U = 1, DT = 0.01, deltaT = 0.001. It is 1D. I have increased the timestep to 0.1 (10 steps, maxCo = 10) and it still works perfectly fine.
If you have a VERY LARGE diffusivity, and a very large timestep, you will lose the diagonal dominance in the last cell, and the iterative solver will blow up (eventually). In fact, I am using ILU_0 preconditioning, so not even that will happen. If you are having trouble, either use a direct solver (tridiagonal matrix in 1D) or reduce the timestep to recover diagonal dominance. If you have very small diffusivity (0.0001), you need to be careful to preserve the matrix because with zero diffusivity the correct answer is a "shock" at the outlet  in fact this is the first solution reported as "wrong" and in fact it is right. Hrv P.S. The case is deadtrivial to set up and it is attached. I really do not understand how and why can someone get this wrong. P.P.S. Where did you all go to school?
__________________
Hrvoje Jasak Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk 

November 3, 2010, 17:51 

#45 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
The original had a relatively low diffusivity, which probably made Santiago think to an instability. Then the discussion deviated on something else, admittedly not correct.
P.S. Please let us know when you will stop being so arrogant. You have been asked by email about this problem before, and the answer you gave (posted above) was not exactly "it works", so it is possible to confuse things. 

November 3, 2010, 18:59 

#46 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 420
Rep Power: 15 
Hrv, Alberto, problem was:
V*dphi/dxk*d^2phi/dx^2=0, phi(0)=1, phi(L)=0 exact solution, phi(x)=(e^(V/k*L)+e^(V/k*x))/(1e^(V/k*L)) completely bounded between 0 and 1. Until I understand this problem has no shock solution while diffusivity is not zero, it has a "boundary layer" at outlet whichs is so thinner as the magnitude of Pe (Problem Pe not mesh Peclet). We checked the matrix assembled in FOAM and it is clear that using the downwind value (when you really are using upwinded scheme) to assemble last equation is exactly the opposite that proffesor Versteeg explains in his book (see reference in post #37). In other words in Hrv Thesis is clear that BC's are the boss even in upwind scheme, but professor Versteeg get rid off the BC in the last cell and see upstream (to the upstream value at cell centroid), then diagonal dominance is assured again. It is quite interesting because he explains the issue exactly with this problem (i knew about this reference some weeks after I started this thread). I ever maintained the topic centered in this issue, and Eugene proposed a solution. I agree with Alberto, following the answer he posted, it seemed to be something like "we kwnow about the issue but is not really important". I'll download the Hrv's case at home and will test it, I hope a really academic and fundamented explanation could be given. Regards. P.S. My academic data is below, our alma mater is Sergio R. Idelsohn.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar 

November 3, 2010, 22:50 

#47 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 420
Rep Power: 15 
Hrv, I tested you case, try this:
ddt(T) steadyState div(phi,T) Gauss Upwind DT={0.1, 0.01, 0.001}, so Pe_mesh={0.1, 1, 10} Solutions are in the figure. Due the steady state, no tricks can be done with the timestep (i.e. increase the diagonal dominance by V/dt contribution), so if div is not correctly assembled for the upwind, solution goes wrong for high Peclet. The only contribution to diagonal you have in the last cell is the one from diffusivity, which is small. If you use the upwind cell centroid value for the face in the outlet (right border), the contribution to diagonal is Sf_{right}*U. The same is explained by Versteeg in his book, he uses F_B*phi_P not F_B*phi_B (as in Eq. 3.58 of your thesis), see attached figures. Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar 

November 5, 2010, 09:03 

#48 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12 
Lets not get all excited.
Hrv, I can second Santiago's latest analysis. This is a fundamental issue, but it is only apparent if you use upwind or an upwindbiased scheme (more than 50% upwind) in conjunction with a fixed value for the transported property at the outlet. Without upwind face interpolation of the convected property there is no issue at all. I have personally had problems with this in several practical applications, so it is not just an academic concern. PS. I went to the same school as you. 

November 5, 2010, 10:25 

#49 
Senior Member
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,768
Rep Power: 21 
... and here's one with Upwind as well.
Would you care to offer an explanation why this works for me and does not for you  having in mind we went to the same school? Hrv
__________________
Hrvoje Jasak Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk 

November 5, 2010, 10:26 

#50 
Senior Member
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,768
Rep Power: 21 
To reproduce this result, take the 1D case above and change the convection discretisation to Gauss upwind.
I made no other changes and... it works!
__________________
Hrvoje Jasak Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk 

November 5, 2010, 10:37 

#51 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 420
Rep Power: 15 
Hrv, could you maintain all your settings and change
ddt(T) steadyState DT={0.01, 0.001}, so Pe_mesh={1, 10} and post your results? Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar 

November 5, 2010, 10:54 

#52 
Senior Member
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,768
Rep Power: 21 
Only if I add underrelaxation to recover diagonal dominance in the matrix that you have just thrown away with the removal of the ddtterm.
If I do not do that, I will no longer preserve the Mmatrix properties, and the ITERATIVE linear equation solver would blow up, JUST LIKE THEY TAUGHT ME AT SCHOOL!!! My other option is to use a direct solver, so even that may work with TDMA. the point is that there is nothing wrong with discretisation nor with linear solvers nor with boundary conditions and you guys have been talking about this for a month.
__________________
Hrvoje Jasak Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk 

November 5, 2010, 10:56 

#53 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
Santiago, the case with DT = 0.01 works fine.
The case with DT = 0.001 will give you problems (you have no points in the BL). Refine the grid around the outlet (I know I am changing the conditions :P), and it will still works. Best, 

November 5, 2010, 11:37 

#54 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 420
Rep Power: 15 
Hrv, I resembled the FOAM problem in octave with a mate and we obtained the same matrix as in FOAM, then solving with an exact solver issue appears too (values goes out of bounds). Diagonal dominance is lost because FOAM assembles the last equation using downstream info instead of really doing upwind.
This is quite interesting, if one use Prof. Versteeg formulation for the last cell (post #47) all the things are good and solution from his book is obtained (bounded solution). We coded this in Matlab and FOAM and things were good. The main matter is, What one must do if had a downstream BC and is doing upwind? To use the BC or to see upstream to the cell centroid to take the BC value? If you follow Jasak you have to do one thing, if you follow Versteeg do the opposite. Which is conceptually correct? Alberto, if you change the dx you are changing the mesh Peclet, then is a completely different numerical problem. Upwind ought give you a bounded solution for all finite Pe. Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar 

November 5, 2010, 12:56 

#55 
New Member
Duncan Roy van der Heul
Join Date: May 2010
Posts: 15
Rep Power: 7 
I now get the correct solutions:
With central discretisation, when the meshPeclet number <2, to guarantee the operator is of positive type, like expected With upwind discretisation, irrespective of the meshPeclet number, like expected (operator will be of positive type irrespective of the meshPeclet number). In both cases using the steadystate solver. Great, I am happy with this result! 

November 5, 2010, 13:33 

#56  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
Quote:
Reducing the step just brings back the matrix to a better situation, exactly as underrelaxing. Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image. OpenQBMM  An opensource implementation of quadraturebased moment methods Last edited by alberto; November 5, 2010 at 14:52. 

November 5, 2010, 16:06 

#57  
New Member
Duncan Roy van der Heul
Join Date: May 2010
Posts: 15
Rep Power: 7 
Quote:
Underrelaxation will make the iterative process converge....to a discrete solution that is not monotonic and therefore not physical. 

November 5, 2010, 16:16 

#58 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
I think we are mixing things (sorry if I contributed to that).
We are talking about upwind schemes, not central schemes, which we all know have a stability condition on the Peclet number. Now, if you consider the upwind scheme, you have no limitation on Pe, but you observe, as Hrv pointed out, a problem due to the numerical iterative solvers of the associated linear system (loss of diagonal dominance => the iterative solver is not guaranteed to converge or be stable), when Pe is large. If you reduce the mesh size where this problem happens, you lower Pe, so you keep the problem under control (which does not mean upwind schemes have a limitation on Pe, it is just a trick to recover in some way the diagonal dominance). That was my observation. Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image. OpenQBMM  An opensource implementation of quadraturebased moment methods Last edited by alberto; November 5, 2010 at 16:20. Reason: Added clarification 

November 5, 2010, 16:28 

#59 
New Member
Duncan Roy van der Heul
Join Date: May 2010
Posts: 15
Rep Power: 7 
I don't know exactly how it is implemented in OF, but normally the first order upwind scheme for the scalar CD equation will lead to an Mmatrix irrespective of the mesh peclet number.
The solution will always be monotonic. The thickness of the boundary layer will be overestimated if the meshpeclet number is large, but the qualitative behavior will be right (monotonic). I used Prof. Hjasak input files, put the div scheme to upwind and got the correct solution for the upwind scheme. Last edited by Duncan_vdH; November 5, 2010 at 16:29. Reason: added 'first order' to upwind 

November 5, 2010, 17:57 

#60 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
If you set DT = 0.001, with that case, without relaxing, refining the grid next to the outlet and/or using a small time stepping in the case of unsteady runs, you won't get the correct solution for the reason explained above (you lose diagonal dominance).
The pictures attached show you what I get with OF 1.7.x, using scalarTransportFoam with T0 = 1 T1 = 2 DT = 0.001 Upwind and with the same settings, but relaxing with URF 0.001 (still wrong solution), and 0.00001. The case and code are attached, since I had to modify scalarTransportFoam to include relaxation. P.S. Plots are made in paraview, so the scales are what they are... 

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 