
[Sponsors] 
January 22, 2006, 10:41 
Hello!
I'm writing an appli

#1 
Member
Masashi IMANO
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 34
Rep Power: 10 
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, 14:06 
Hi,
This is about right, bu

#2 
Senior Member
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,810
Rep Power: 25 
Hi,
This is about right, but not easily extensible to support the turbulence modelling runtime selection. In reality, we need a decomposition of turbulence>divR into the implicit and explicit part, which is modeldependent. 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 modelspecific), wall functions specify the source term and you will need to dig it out.  for epsilon, the nearwall 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 modelspecific :) 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, 07:28 
Thanks for the reply, Dr. Hrvo

#3 
Member
Masashi IMANO
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 34
Rep Power: 10 
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, 08:14 
Hello and sorry about this lon

#4 
Member
Masashi IMANO
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 34
Rep Power: 10 
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/OpenFOAM1.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 adhoc hack, but anyway it supports the turbulence modelling runtime selection on my application code... 

February 16, 2006, 09:00 
Hello,
Actually, you did al

#5 
Senior Member
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,810
Rep Power: 25 
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/OpenFOAM1.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 20, 2006, 00:12 
Hi,
Thanks for the kind rep

#6 
Member
Masashi IMANO
Join Date: Mar 2009
Location: Tokyo, Japan
Posts: 34
Rep Power: 10 
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: 10 
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] ] + [(((resistnu)*gh)*ddt(rho))[0 2 2 0 0 0 0] ] From function checkMethod(const fvMatrix<type>&, const GeometricField<type,>&) in file /home/liu/OpenFOAM/OpenFOAM1.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: 10 
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] ] + [(((resistnu)*gh)*ddt(rho))[0 2 2 0 0 0 0] ] From function checkMethod(const fvMatrix<type>&, const GeometricField<type,>&) in file /home/liu/OpenFOAM/OpenFOAM1.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: 125
Rep Power: 10 
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. 

Thread Tools  
Display Modes  


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 