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

Estimating error for the incompressible turbulent flow

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

Reply
 
LinkBack Thread Tools Display Modes
Old   January 22, 2006, 10:41
Default Hello! I'm writing an appli
  #1
Member
 
Masashi IMANO
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 34
Rep Power: 8
imano is on a distinguished road
Hello!

I'm writing an application simpleErroeEstimate for estimating error
for the incompressible turbulent CFD application simpleFoam using
Dr. Hrvoje's residual error estimation method.

So I will need to add the Reynolds stress term (i.e. turbulence->divR) to
the following lines of icoErrorEstimate.C, but I am not sure of doing that.

applications/utilities/errorEstimation/icoErrorEstimate/icoErrorEstimate.C:

errorEstimate<vector> ee
(
resError::div(phi, U)
- resError::laplacian(nu, U)
==
-fvc::grad(p)
);


(a) If turbulence->divR returns

- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*dev(fvc::grad(U)().T()))

like KEpsilon class, can I write the "ee" as follow?

errorEstimate<vector> ee
(
resError::div(phi, U)
- resError::laplacian(turbulence->nuEff(), U)
- fvc::div(turbulence->nuEff()*dev(fvc::grad(U)().T()))
==
-fvc::grad(p)
);

(b) Is there more elegant way which can handle the divR of any
turbulent model in simpleFoam?

(c) Need a special treatment for the wall function in the residual error
estimation for turbulent flow?

Thanks,
Masashi
imano is offline   Reply With Quote

Old   January 22, 2006, 14:06
Default Hi, This is about right, bu
  #2
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,761
Rep Power: 21
hjasak will become famous soon enough
Hi,

This is about right, but not easily extensible to support the turbulence modelling run-time selection. In reality, we need a decomposition of turbulence->divR into the implicit and explicit part, which is model-dependent.

In most cases, you will have something like:

fvm::laplacian(nuEff(), U) for the implicit part (fine) and the most general way of handling the explicit part would be something like:

R() - fvc::laplacian(nuEff(), U)

which is not too bad (I would like to check that the nuEff() is OK for all the models, just to make sure).

As for c), here's how things stand:
- for the momentum equation you don't need any wall function corrections: it is all built into the wall patch values of nuEff()
- for k (now we're getting model-specific), wall functions specify the source term and you will need to dig it out.
- for epsilon, the near-wall cell value is fixed and there's no error estimation possible

However, each model is specific in its own way and equations it is solving, so please handle with care. Also, beware of boundary conditions on the turbulence equations - these are VERY model-specific :-)

I'd be very interested in your results, please keep me posted. You can find more details on error estimation in turbulent flows in a PhD Thesis of my former student dr. Franjo Juretic.

Good luck,

Hrv
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   January 23, 2006, 07:28
Default Thanks for the reply, Dr. Hrvo
  #3
Member
 
Masashi IMANO
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 34
Rep Power: 8
imano is on a distinguished road
Thanks for the reply, Dr. Hrvoje.

I think the explicit part of divR would be expressed as:

fvc::div(turbulence->R() -((2.0/3.0)*I)*turbulence->k() -nu*2*symm(fvc::grad(U)))
+ fvc::laplacian(turbulence->nuEff(), U)

but the result of error estimation (i.e. resErrorU and mag(resErrorU) )
using this expression slightly differs from the old one.

Can I consider that this is due to discretisation error?


And thanks again for introducing Dr. Franjo's PhD thesis.
To tell the truth I have read it before my post, but I could not find
treatment of Reynolds stress for the residual error estimation in it.
Anyway I will check it again!

Regards,
Masashi
imano is offline   Reply With Quote

Old   February 16, 2006, 08:14
Default Hello and sorry about this lon
  #4
Member
 
Masashi IMANO
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 34
Rep Power: 8
imano is on a distinguished road
Hello and sorry about this long message!

I have started writing an application for estimating error of k using
the residual method.
As the estimator needs the effective diffusivity and the generation
term for k equation, I added two member functions to the
turbulenceModel class and it's derived classes like these:

turbulenceModel.H:
---
//- Return the effective diffusivity for k
virtual tmp<volscalarfield> DkEff() const = 0;

//- Return the generation term for k equation
virtual tmp<volscalarfield> Gk() const = 0;
---

And I changed the core of stimateScalarError.C like that:
---
volScalarField k = turbulence->k();
volScalarField epsilon = turbulence->epsilon();
volScalarField DkEff = turbulence->DkEff();
volScalarField Gk = turbulence->Gk();

errorEstimate<scalar> ee
( resError::div(phi, k)- resError::laplacian(DkEff, k) == Gk - resError::Sp(epsilon / k, k) );
---

I succeeded a compilation of this code, but on a run time I recieved a
FATAL ERROR as follows:
----
--> FOAM FATAL ERROR : incompatible dimensions for operation
[k[0 -1 -3 0 0 0 0] ] - [((((Cmu*sqr(k))|(epsilon+epsilonSmall))*2)*magSqr( symm(grad(U))))[0 2 -3 0 0 0 0] ]

From function checkMethod(const errorEstimate<type>&, const GeometricField<type,>&)
in file ~/OpenFOAM/OpenFOAM-1.2/src/errorEstimation/lnInclude/errorEstimate.C at line 436.
----

In order to check dimensions of implicit source term and other terms, I
wrote a test code like this:
----
errorEstimate<scalar> e1 (resError::div(phi, k));
Info << "e1: " << e1.dimensions() << e1.residual()().dimensions() << endl;

errorEstimate<scalar> e2 (resError::Sp(epsilon / k, k));
Info << "e2: " << e2.dimensions() << e2.residual()().dimensions() << endl;
----
And result:
----
e1: [0 5 -3 0 0 0 0][0 2 -3 0 0 0 0]
e2: [0 2 -3 0 0 0 0][0 2 -3 0 0 0 0]
----

Now I have some questions:

a) In order to avoid this error, I changed the resError::Sp functions like this:
----
new errorEstimate<type>
(
vf,
sp.dimensions()*vf.dimensions()*dimVol, // <- sp.dimensions()*vf.dimensions(),
sp.internalField()*vf.internalField(),
scalarField(vf.internalField().size(), 0)
)
----

It seems to work well, but is this consistent with other code?

b) In Dr. Hrvoje's thesis, the coefficient of implicit source Sp is
included into the normalisation factor.
But it seems that resError::Sp functions retrurn the factor
(i.e. the fourth argument in the above code) as zero.
Why don't we use sp.internalField() for scaling?

Regards,
Masashi


P.S.
As for question (b) in my first post on this thread, instead of
decomposing the explicit part of divR, I added a new member function
divRExplicit as follows:

turbulenceModel.H:
---
//- Return the explicit part of source term for the momentum equation
virtual tmp<volvectorfield> divRExplicit(volVectorField& U) const = 0;
---

e.g. kEpsilon.C:
---
tmp<volvectorfield> kEpsilon::divRExplicit(volVectorField& U) const
{
return
(- fvc::div(nuEff()*dev(fvc::grad(U)().T())));
}
---

I realize that this is ad-hoc hack, but anyway it supports the
turbulence modelling run-time selection on my application code...
imano is offline   Reply With Quote

Old   February 16, 2006, 09:00
Default Hello, Actually, you did al
  #5
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,761
Rep Power: 21
hjasak will become famous soon enough
Hello,

Actually, you did all of the first bit right and the bug fix is fine. You will need to add anoter one in

~/OpenFOAM/OpenFOAM-1.2/src/errorEstimation/lnInclude/resErrorSup.C

on line 89 as well (again *dimVolume).

Sorry - my fault (curious that it has not been found out before).

As for the normalization, I keep changing my mind on that. The error estimate really works on the transport part of the equation. If the equation is dominated by local sources/sinks, the transport error is much less visible because the (potentially huge) Sp normalization factor masks it out. The last big systematic study I did this on has shown that the absolute value of the estimated error is better without this normalization but please feel free to experiment.

The turbulence stuff is all fine as well - this is why we released the code in full source in the first place.

Enjoy,

Hrv


Hrv
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   February 20, 2006, 00:12
Default Hi, Thanks for the kind rep
  #6
Member
 
Masashi IMANO
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 34
Rep Power: 8
imano is on a distinguished road
Hi,

Thanks for the kind reply.
I understand the reason about the normalization and I will go on experiment around that on my test study!

Regards,
Masashi
imano is offline   Reply With Quote

Old   June 20, 2006, 23:20
Default Hi friend, I also have prob
  #7
Senior Member
 
Guoxiang
Join Date: Mar 2009
Posts: 109
Rep Power: 8
liugx212 is on a distinguished road
Hi friend,

I also have problem like you. So could you please give some advice?

my error like:
--> FOAM FATAL ERROR : incompatible dimensions for operation
[pd[0 0 -1 0 0 0 0] ] + [(((resist|nu)*gh)*ddt(rho))[0 2 -2 0 0 0 0] ]

From function checkMethod(const fvMatrix<type>&, const GeometricField<type,>&)
in file /home/liu/OpenFOAM/OpenFOAM-1.2/src/OpenFOAM/lnInclude/fvMatrix.C at line 935.

FOAM aborting

Aborted


my solver is:\*--------------------------------------------------------------------------- */

#include "fvCFD.H"
#include "basicThermo.H"
//#include "compressible/turbulenceModel/turbulenceModel.H"
#include "fixedGradientFvPatchFields.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{

# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "readEnvironmentalProperties.H"
# include "createFields.H"
//# include "initContinuityErrs.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nStarting time loop\n" << endl;

while (!runTime.end())
{
# include "readPISOControls.H"
# include "CourantNo.H"

runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;

//# include "rhoEqn.H"
solve(Porous*fvm::ddt(rho) + fvc::div(phi) == rho*fvm::Sp(rhoinj,rho)); // DDD! add rho equation

//# include "UEqn.H"
fvVectorMatrix UEqn // DDD! solving Darcy velocity
(
fvm::Sp((nu/resist),U) - fvc::grad(rho)*gh
);
solve(UEqn == - fvc::grad(pd));

// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
//# include "pEqn.H"
pd.boundaryField() ==
p.boundaryField() - rho.boundaryField()*gh.boundaryField() - pRef.value();
volScalarField rUA = 1.0/UEqn.A();
U = rUA*UEqn.H();
phi = (fvc::interpolate(rho*U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
- fvc::interpolate(rUA*rho*gh)*fvc::snGrad(rho)*mesh .magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn // DDD! solving pressure equation
(
-(resist/nu)*fvm::laplacian(pd)
+ (resist/nu)*gh*fvc::ddt(rho)
);
solve(pEqn == rho*fvm::Sp(rhoinj, rho)); // DDD! solving pressure equation
if (nonOrth == nNonOrthCorr)
{
phi += pEqn.flux();
}
}
/* p = (pd + rho*g + pRef);
dpdt = fvc::ddt(p);
*/
fvScalarMatrix CEqn // DDD! solving concentration equation
(
Porous*fvm::ddt(rho, C)
+ fvm::div(phi, C)
- D*fvm::laplacian(C)
);
solve(CEqn == rho*fvm::Sp(rhoinj, rho)); // DDD! solving concentration equation
}
runTime.write();
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;
}
Info<< "End\n" << endl;
return 0;
}
liugx212 is offline   Reply With Quote

Old   June 20, 2006, 23:25
Default Hi friend, I also have prob
  #8
Senior Member
 
Guoxiang
Join Date: Mar 2009
Posts: 109
Rep Power: 8
liugx212 is on a distinguished road
Hi friend,

I also have problem like you. So could you please give some advice?

my error like:
--> FOAM FATAL ERROR : incompatible dimensions for operation
[pd[0 0 -1 0 0 0 0] ] + [(((resist|nu)*gh)*ddt(rho))[0 2 -2 0 0 0 0] ]

From function checkMethod(const fvMatrix<type>&, const GeometricField<type,>&)
in file /home/liu/OpenFOAM/OpenFOAM-1.2/src/OpenFOAM/lnInclude/fvMatrix.C at line 935.

FOAM aborting

Aborted


my solver is:\*--------------------------------------------------------------------------- */

#include "fvCFD.H"
#include "basicThermo.H"
//#include "compressible/turbulenceModel/turbulenceModel.H"
#include "fixedGradientFvPatchFields.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{

# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "readEnvironmentalProperties.H"
# include "createFields.H"
//# include "initContinuityErrs.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nStarting time loop\n" << endl;

while (!runTime.end())
{
# include "readPISOControls.H"
# include "CourantNo.H"

runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;

//# include "rhoEqn.H"
solve(Porous*fvm::ddt(rho) + fvc::div(phi) == rho*fvm::Sp(rhoinj,rho)); // DDD! add rho equation

//# include "UEqn.H"
fvVectorMatrix UEqn // DDD! solving Darcy velocity
(
fvm::Sp((nu/resist),U) - fvc::grad(rho)*gh
);
solve(UEqn == - fvc::grad(pd));

// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
//# include "pEqn.H"
pd.boundaryField() ==
p.boundaryField() - rho.boundaryField()*gh.boundaryField() - pRef.value();
volScalarField rUA = 1.0/UEqn.A();
U = rUA*UEqn.H();
phi = (fvc::interpolate(rho*U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
- fvc::interpolate(rUA*rho*gh)*fvc::snGrad(rho)*mesh .magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn // DDD! solving pressure equation
(
-(resist/nu)*fvm::laplacian(pd)
+ (resist/nu)*gh*fvc::ddt(rho)
);
solve(pEqn == rho*fvm::Sp(rhoinj, rho)); // DDD! solving pressure equation
if (nonOrth == nNonOrthCorr)
{
phi += pEqn.flux();
}
}
/* p = (pd + rho*g + pRef);
dpdt = fvc::ddt(p);
*/
fvScalarMatrix CEqn // DDD! solving concentration equation
(
Porous*fvm::ddt(rho, C)
+ fvm::div(phi, C)
- D*fvm::laplacian(C)
);
solve(CEqn == rho*fvm::Sp(rhoinj, rho)); // DDD! solving concentration equation
}
runTime.write();
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;
}
Info<< "End\n" << endl;
return 0;
}
liugx212 is offline   Reply With Quote

Old   October 9, 2013, 16:04
Default
  #9
Senior Member
 
M. Montero
Join Date: Mar 2009
Location: Madrid
Posts: 112
Rep Power: 8
be_inspired is on a distinguished road
Hi all,

I have use the simpleFoamResidual utility implemented into OF1.6ext. The implementation is different that the one stated in this post:

uResidual:
mag(fvc::div(phi,U) + fvc::div(turbulence->R())+fvc::grad(p)

and my results are:
uResidual max:2255.13 mean:0.000733257
The turbulence model that I have used in realizableKE.

If I review the uResidual field with paraFoam, there are some cells with high uResidual but I do not know how to interpretate the figures. Is needed to normalize the uResidual field to obtain a error field in meters per second?
uResidual equals to 2255 is high? Very High? Is needed to refine the mesh in those regions?

Thank you very much

Last edited by be_inspired; October 10, 2013 at 04:40.
be_inspired 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
3D Incompressible turbulent flow validation case Yu Main CFD Forum 2 April 4, 2008 03:13
2d incompressible turbulent flow na Main CFD Forum 4 June 16, 2005 08:41
2D incompressible turbulent flow na Main CFD Forum 5 June 2, 2005 11:02
2d turbulent incompressible flow aghiba Main CFD Forum 2 April 11, 2005 08:53
Turbulent incompressible flow & moving boundaries Cristian Orozco Main CFD Forum 1 July 5, 2002 12:45


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