- **OpenFOAM Running, Solving & CFD**
(*http://www.cfd-online.com/Forums/openfoam-solving/*)

- - **Estimating error for the incompressible turbulent flow**
(*http://www.cfd-online.com/Forums/openfoam-solving/60154-estimating-error-incompressible-turbulent-flow.html*)

Hello!
I'm writing an appliHello!
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 |

Hi,
This is about right, buHi,
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 |

Thanks for the reply, Dr. HrvoThanks 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 |

Hello and sorry about this lonHello 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... |

Hello,
Actually, you did alHello,
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 |

Hi,
Thanks for the kind repHi,
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 |

Hi friend,
I also have probHi 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; } |

Hi friend,
I also have probHi 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; } |

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 |

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