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

Convection-diffusion in 1D : wrong solution for a large Delta x

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

Reply
 
LinkBack Thread Tools Display Modes
Old   June 4, 2010, 09:09
Default Convection-diffusion in 1D : wrong solution for a large Delta x
  #1
Senior Member
 
Emanuele
Join Date: Mar 2009
Posts: 110
Rep Power: 8
nuovodna is on a distinguished road
Hi, i wrote a solver to solve a convection-diffusion equation on a 1D mesh

fvScalarMatrix noUEqn
(
fvm::div(v,noU)
== fvm::laplacian(d1,noU)
);

where v is a constant value surfaceScalarField (value is 4), d1 is a constant value surfaceScalarField (value is 0.1) and noU is a volScalarField

Mesh: 1 m x 1 m x 1m [10 x 1 x 1] -> Delta x = 0.1
Boundary condition : inlet: noU=0 outlet: noU=1

Using UPWIND scheme on div(v,noU) i obtain a wrong noU field : all values are negative (the only positive value remains at the outlet like boundary condition says)!! If i reduce Delta x to 0.01 i obtain a positive noU field according to analitical solution.
I printed noUEqn.H() and noUEqn.A() and they are like the values manually computed with a simple spreadsheet.
How can i fix it?? Why at lower resolution it gives me negative values?

Thanks in advance
Regards

Emanuele

Last edited by nuovodna; June 7, 2010 at 09:42.
nuovodna is offline   Reply With Quote

Old   June 7, 2010, 09:45
Default
  #2
Senior Member
 
Emanuele
Join Date: Mar 2009
Posts: 110
Rep Power: 8
nuovodna is on a distinguished road
If i use UPWIND scheme on div, I have this noU cell values:

Code:
 -3.41333e-07
 -2.38933e-06
 -1.26293e-05
 -6.38293e-05
 -0.000319829
 -0.00159983
 -0.00799983
 -0.0399998
 -0.2
 -1
These values are totally different from analytical solution: they should be positive
nuovodna is offline   Reply With Quote

Old   June 18, 2010, 05:47
Default
  #3
Senior Member
 
Emanuele
Join Date: Mar 2009
Posts: 110
Rep Power: 8
nuovodna is on a distinguished road
I invert the matrix manually and it returns the same negative value. Perhaps, the effect of boundary is dominant
nuovodna is offline   Reply With Quote

Old   July 1, 2010, 13:25
Default
  #4
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Emanuele: are you inverting the noUEqn matrix

Code:
fvScalarMatrix noUEqn
        (                    
            fvm::div(v,noU)            
            == fvm::laplacian(d1,noU)
        );
inmediately after it is assembled or once solve is called?. Before calling solve noUEqn doesn't include the BC contribution.

Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   July 2, 2010, 05:45
Default
  #5
Senior Member
 
Emanuele
Join Date: Mar 2009
Posts: 110
Rep Power: 8
nuovodna is on a distinguished road
I printed the noU values in this way

fvScalarMatrix noUEqn
(
fvm::div(v,noU)
== fvm::laplacian(d1,noU)
);

solve(noUEqn);

OFstream outCellValues("outCellValues.dat");
forAll(mesh.cells(), celli)
{
outCellValues << mesh.C()[celli].component(0) << " " << noU[celli] << endl;
}
nuovodna is offline   Reply With Quote

Old   July 2, 2010, 11:48
Default Typical behaviour I guess
  #6
Senior Member
 
Antonio Martins
Join Date: Mar 2009
Location: Porto, Porto, Portugal
Posts: 112
Rep Power: 8
titio is on a distinguished road
Send a message via MSN to titio Send a message via Skype™ to titio
I believe you are getting the typical behaviour. If DeltaX is too high and using upwind you will get to much numerical dispersion, that disapear if you define a lower value of deltaX. Consult Versteeg for more information.

Regards,

Antonio Martins
titio is offline   Reply With Quote

Old   October 2, 2010, 18:47
Default
  #7
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Emanuele, I had forgotten this thread until a few days ago when I started to have the same issues with a similar problem, I reported them into the FOAM bug tracker and in the bug section of this forum:

Wrong fvm::div assembling

problems seems to be related to the way of upwind scheme is implemented in FOAM.

Could post some info related to your problem, including books and papers where it is analyzed?

Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   October 2, 2010, 19:56
Default
  #8
Member
 
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 8
ovie is on a distinguished road
Hi,

I know this isnt related exactly to the problem in this thread, but bears some similarities nonetheless.

My question is: has anyone tried to implement InterFOAM in 1-D. I am trying to simulate Stefan problems using interFoam and apparently it looks like the solution in 1-D makes no sense while 2-D formulation appears to influence the results somewhat due to the additional boundary conditions.

If anyone has any ideas it would be nice.

Thanks.
ovie is offline   Reply With Quote

Old   October 2, 2010, 23:07
Default
  #9
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Ovie, I've tried an 1D example of interFoam in order to study how compressive term works. It's a very interesting case because, momentum equation vanishes and only non-linear equation for phase is solved really. It allowed me to compare interFoam results with simple Matlab/Octave code. I'll be working on that again next weeks.

Maybe you can start a new thread about this topic to share some results.

Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   October 2, 2010, 23:40
Default
  #10
Member
 
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 8
ovie is on a distinguished road
Thanks Santiago,

As it is, I wasnt sure if my point had any merit to it thats why I was reluctant to start a new thread. But on your suggestion, I might just consider doing that to see if others with similar challenges can report their findings.

Thanks for your response all the same.
ovie is offline   Reply With Quote

Old   October 3, 2010, 00:59
Default
  #11
Member
 
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 8
ovie is on a distinguished road
Quote:
Originally Posted by santiagomarquezd View Post
Ovie, I've tried an 1D example of interFoam in order to study how compressive term works. It's a very interesting case because, momentum equation vanishes and only non-linear equation for phase is solved really. It allowed me to compare interFoam results with simple Matlab/Octave code. I'll be working on that again next weeks.

Maybe you can start a new thread about this topic to share some results.

Regards.
I have started a thread already.

But I am curious. How did you formulate the problem? I mean the mesh and what boundary conditions did you impose on patches? I have tried a simple approach of declaring only inlet and outlet patches for a 1-D rectangular block but when I run the computation, it is just completely useless results I get. So I would really like to see how you managed to pull this off...

Thanks
ovie is offline   Reply With Quote

Old   October 4, 2010, 05:53
Default one dimensional convection-diffusion results
  #12
Senior Member
 
Emanuele
Join Date: Mar 2009
Posts: 110
Rep Power: 8
nuovodna is on a distinguished road
Hi Santiago, these are my results.

Equation:

fvm::div(v,U) == fvm::laplacian(d,U)

with
d= 0.001
v = 1
U(0) = 0
U(1) = 1

Four schemes combination on div :
-) Gauss linear on div + Gauss linear corrected on laplacian (in figure as CDS)
-) Gauss upwind on div + Gauss linear corrected on laplacian (UPWIND)
-) Gauss linearUpwind cellMDLimited Gauss linear 1 on div + Gauss linear corrected on laplacian (LIMITERS)
-) Gauss limitedLinear 1 on div + Gauss linear corrected on laplacian (limitedLinear)


Classified according to Peclet number:

Pe = v * Delta x / d


Regards

Emanuele

Last edited by nuovodna; October 4, 2010 at 06:42. Reason: errors on schemes definition
nuovodna is offline   Reply With Quote

Old   October 4, 2010, 07:29
Default What is the value of the time increment
  #13
Senior Member
 
Antonio Martins
Join Date: Mar 2009
Location: Porto, Porto, Portugal
Posts: 112
Rep Power: 8
titio is on a distinguished road
Send a message via MSN to titio Send a message via Skype™ to titio
Hi,

What is the value of the value increment? According to the figures, it looks you have a larger courant than allowed by stability requirements. Also, try to use more precise interpolation methods, such as gamma, minmod, or SFCD. Although non limited, they are much more precise than upwind...

Also, more information on this problem can be found in the book of Malaskera.

Regards,

António Martins
titio is offline   Reply With Quote

Old   October 4, 2010, 07:58
Default
  #14
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Emanuele. your results seems as I expected, particularly for Central Difference and Upwind. In CD the oscillations are normal due large Pe, this scheme is unstable for Pe>1. With respect of upwind it should be stable in all range of Pe numbers, but in FOAM it is not the case due the form of divergence assembling (check the link I posted a few days ago). Respect of limited schemes I would have expected a bounded behavior but it seem to have the same problems.
Could you post a paper o book in which this problem is shown? I'm searching too, but the only references I've got are for FEM.

Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   October 4, 2010, 08:18
Default expected results
  #15
Senior Member
 
Emanuele
Join Date: Mar 2009
Posts: 110
Rep Power: 8
nuovodna is on a distinguished road
Hi Santiago, i read your bug report and my results correspond to your alert on fvm::div assembling. This problem is described here

Finite Volume Method for Convection-Diffusion Problems

and in famous book like Peric or Veersteg
nuovodna is offline   Reply With Quote

Old   October 20, 2010, 13:36
Default
  #16
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Emanuele, Professor Tung gives the right answer in slide 47 as you said, the explanation is taken from Versteeg.

Thx.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd 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
Solution in the Large for Nonlinear Hyperbolic CFD_student Main CFD Forum 0 March 19, 2007 13:04
Analytical Solution of 1D Convection Diffusion Eq. Suman Kumar Main CFD Forum 7 July 15, 2003 14:05
CFL Condition Matt Umbel Main CFD Forum 14 January 12, 2001 15:34
Numerical diffusion error Z.Zeng Main CFD Forum 8 October 22, 1999 09:06
Wall functions Abhijit Tilak Main CFD Forum 6 February 5, 1999 02:16


All times are GMT -4. The time now is 13:00.