[imano]

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

[hjasak]

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

[imano]

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]

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...

[hjasak]

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

[imano]

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

[liugx212]

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]

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;
}

[be_inspired]

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