|
[Sponsors] |
Estimating error for the incompressible turbulent flow |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
January 22, 2006, 09:41 |
Hello!
I'm writing an appli
|
#1 |
Member
Masashi IMANO
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 34
Rep Power: 17 |
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 |
|
January 22, 2006, 13:06 |
Hi,
This is about right, bu
|
#2 |
Senior Member
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,905
Rep Power: 33 |
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 |
|
January 23, 2006, 06:28 |
Thanks for the reply, Dr. Hrvo
|
#3 |
Member
Masashi IMANO
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 34
Rep Power: 17 |
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 |
|
February 16, 2006, 07:14 |
Hello and sorry about this lon
|
#4 |
Member
Masashi IMANO
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 34
Rep Power: 17 |
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... |
|
February 16, 2006, 08:00 |
Hello,
Actually, you did al
|
#5 |
Senior Member
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,905
Rep Power: 33 |
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 |
|
February 19, 2006, 23:12 |
Hi,
Thanks for the kind rep
|
#6 |
Member
Masashi IMANO
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 34
Rep Power: 17 |
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 |
|
June 20, 2006, 23:20 |
Hi friend,
I also have prob
|
#7 |
Senior Member
Guoxiang
Join Date: Mar 2009
Posts: 109
Rep Power: 17 |
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; } |
|
June 20, 2006, 23:25 |
Hi friend,
I also have prob
|
#8 |
Senior Member
Guoxiang
Join Date: Mar 2009
Posts: 109
Rep Power: 17 |
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; } |
|
October 9, 2013, 16:04 |
|
#9 |
Senior Member
M. Montero
Join Date: Mar 2009
Location: Madrid
Posts: 153
Rep Power: 17 |
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. |
|
|
|
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 |