
[Sponsors] 
March 13, 2013, 07:40 
Manual limiter of velocity doesn't work

#1 
Member
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 7 
Hi to everyone guys!
I'm experiencing this problem and I really don' t know what to do: I'm running with the solver adjointShapeOptimizationFoam, but now I'm trying to modify it a little, since sometimes the equations diverge. What I want to do is to use a manual limiter, that limits the value of the velocity in each cell if it raises too much, so I've inserted these lines at the bottom of the solver: forAll(Ua,cellI) { Ua[cellI].component(0)=min(Ua[cellI].component(0), 200); Ua[cellI].component(1)=min(Ua[cellI].component(1), 200); Ua[cellI].component(2)=min(Ua[cellI].component(2), 200); } forAll(Ua,cellJ) { Ua[cellJ].component(0)=max(Ua[cellJ].component(0), 200); Ua[cellJ].component(1)=max(Ua[cellJ].component(1), 200); Ua[cellJ].component(2)=max(Ua[cellJ].component(2), 200); } I'm expecting that these lines behave like a threshold value for each component of the velocity at each iteration, but I've run a simulation and I discovered that the values of the velocity were bigger than the range [200:200]. Maybe the lines I've added are wrong? Any help is really appreciated. Thanks in advance Simone P.s. I'm running in parallel, but I hope this is not a problem for the code I've added. 

March 14, 2013, 03:22 

#2 
Member
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 7 
Please guy, is there someone that can help?


March 14, 2013, 09:29 

#3 
Senior Member
Join Date: Nov 2009
Location: Michigan
Posts: 135
Rep Power: 10 
You will need to limit the values on the patches also (all faces on patches)
forAll(U.boundaryField(),patchI) { forAll(U.boundaryField()[patchI],faceI) { U.boundaryField()[patchI][faceI].component(0)=some value; } } Also your solver might adjust U values not only after solving momentum equation. So make sure that your are limiting after each calculation step of U 

March 15, 2013, 07:09 

#4 
Member
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 7 
Thank you Omkar! The problem was exactly on the patches, since I had forgotten to loop on them. By adding your lines now it seems to work perfectly.


March 15, 2013, 07:46 

#5 
Senior Member
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,209
Rep Power: 20 
hi Simone
I think limiting velocity can also resolve my problem. Could you send me your modified solver? Thanks. 

March 16, 2013, 18:29 

#6 
Senior Member
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,209
Rep Power: 20 
where should i add them exactly?


March 17, 2013, 20:24 

#7 
Senior Member
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,209
Rep Power: 20 
can use velocity limiters only on a patch not entire the domain?


March 17, 2013, 20:26 

#8 
Senior Member
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,209
Rep Power: 20 
hi Omkar
Could you send me the code with added expressions for rhoPimpleFoam? 

March 17, 2013, 20:37 

#9 
Senior Member
Join Date: Nov 2009
Location: Michigan
Posts: 135
Rep Power: 10 
My code consists of combination of above two codes.
Like I said, it did not work for me so I erased that code. The codes in this forum clearly explain how to limit velocity in the domain and the patches 

March 18, 2013, 04:04 

#10  
Member
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 7 
Hi immortality!
With respect to the adjointShapeOptimizationFoam solver, you should put the " two" limiter just above the lines Quote:
If you want to "limit" only the patches you should insert only the lines that doubtsincfd suggested. If, instead, you asked to limit only one specific patch you can use something like this: Quote:
Simone 

March 19, 2013, 13:57 

#11 
Senior Member
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,209
Rep Power: 20 
thanks.where the codes should be added in rhoPimpleFoam?
rhoPimpleFoam.C is this: Code:
#include "fvCFD.H" #include "basicPsiThermo.H" #include "turbulenceModel.H" #include "bound.H" #include "pimpleControl.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { #include "setRootCase.H" #include "createTime.H" #include "createMesh.H" pimpleControl pimple(mesh); #include "createFields.H" #include "initContinuityErrs.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; while (runTime.run()) { #include "readTimeControls.H" #include "compressibleCourantNo.H" #include "setDeltaT.H" runTime++; Info<< "Time = " << runTime.timeName() << nl << endl; #include "rhoEqn.H" //  Pressurevelocity PIMPLE corrector loop while (pimple.loop()) { #include "UEqn.H" #include "hEqn.H" //  Pressure corrector loop while (pimple.correct()) { #include "pEqn.H" } if (pimple.turbCorr()) { turbulence>correct(); } } runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End\n" << endl; return 0; } // ************************************************************************* // Code:
// Solve the Momentum equation tmp<fvVectorMatrix> UEqn ( fvm::ddt(rho, U) + fvm::div(phi, U) + turbulence>divDevRhoReff(U) ); UEqn().relax(); volScalarField rAU(1.0/UEqn().A()); if (pimple.momentumPredictor()) { solve(UEqn() == fvc::grad(p)); K = 0.5*magSqr(U); } Code:
{ fvScalarMatrix hEqn ( fvm::ddt(rho, h) + fvm::div(phi, h)  fvm::laplacian(turbulence>alphaEff(), h) == dpdt  (fvc::ddt(rho, K) + fvc::div(phi, K)) ); hEqn.relax(); hEqn.solve(); thermo.correct(); } Code:
rho = thermo.rho(); rho = max(rho, rhoMin); rho = min(rho, rhoMax); rho.relax(); U = rAU*UEqn().H(); if (pimple.nCorrPISO() <= 1) { UEqn.clear(); } if (pimple.transonic()) { surfaceScalarField phid ( "phid", fvc::interpolate(psi) *( (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvm::div(phid, p)  fvm::laplacian(rho*rAU, p) ); pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); if (pimple.finalNonOrthogonalIter()) { phi == pEqn.flux(); } } } else { phi = fvc::interpolate(rho)* ( (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); while (pimple.correctNonOrthogonal()) { // Pressure corrector fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvc::div(phi)  fvm::laplacian(rho*rAU, p) ); pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); if (pimple.finalNonOrthogonalIter()) { phi += pEqn.flux(); } } } #include "rhoEqn.H" #include "compressibleContinuityErrs.H" // Explicitly relax pressure for momentum corrector p.relax(); // Recalculate density from the relaxed pressure rho = thermo.rho(); rho = max(rho, rhoMin); rho = min(rho, rhoMax); rho.relax(); Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; U = rAU*fvc::grad(p); U.correctBoundaryConditions(); K = 0.5*magSqr(U); dpdt = fvc::ddt(p); 

March 21, 2013, 07:06 

#12 
Member
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 7 
In my opinion after
Code:
if (pimple.turbCorr()) { turbulence>correct(); } 

March 26, 2013, 08:57 

#13 
Senior Member
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,209
Rep Power: 20 
thank you dear Simone for your help
so I added it.is it correct?why you have wrote patchName2 in findPatchID? I put U[cellI] instead of Ua[cellI] due to use in rhoPimpleFoam does it have any problem? I want not to let U becomes higher than sound speed(flow should be subsonic) will it be true that I write sqrt(1.4*287.14*T[cellI]) instead of 500 I have put now? I want to be certain to compile the modified solver. thanks again. Code:
word patchName = "left"; label patchID = mesh.boundary().findPatchID(patchName); forAll(U.boundaryField()[patchID],faceI) { U.boundaryField()[patchID][faceI].component(0)=min(U[cellI].component(0), 500); } 

March 26, 2013, 11:49 

#14 
Senior Member
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,209
Rep Power: 20 
I used it for the patch.but velocity in cells near the inflow boundary are supersonic.
where I add the expressions to limit velocity on top of domain? and what do you mean from cellI and cellJ in this expressions: forAll(Ua,cellI) { Ua[cellI].component(0)=min(Ua[cellI].component(0), 200); Ua[cellI].component(1)=min(Ua[cellI].component(1), 200); Ua[cellI].component(2)=min(Ua[cellI].component(2), 200); } forAll(Ua,cellJ) { Ua[cellJ].component(0)=max(Ua[cellJ].component(0), 200); Ua[cellJ].component(1)=max(Ua[cellJ].component(1), 200); Ua[cellJ].component(2)=max(Ua[cellJ].component(2), 200); } 

March 27, 2013, 10:59 

#15 
Member
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 7 
You have to put those line before the ones to limit the boundary. cellI and cellJ are two (unuseful) different counters. You can use the same for both cycles.


March 27, 2013, 11:20 

#16 
Senior Member
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,209
Rep Power: 20 
thank you dear Simone again
is this correct now?: Code:
forAll(U,cellI) { U[cellI].component(0)=min(U[cellI].component(0), sqrt(1.4*287.14*T.boundaryField()[patchID][faceI])30); } forAll(U,cellJ) { U[cellJ].component(1)=max(U[cellI].component(1),150); } word patchName = "left"; label patchID = mesh.boundary().findPatchID(patchName); forAll(U.boundaryField()[patchID],faceI) { U.boundaryField()[patchID][faceI].component(0)=min(U.boundaryField()[patchID][faceI].component(0), sqrt(1.4*287.14*T.boundaryField()[patchID][faceI])30); U.boundaryField()[patchID][faceI].component(1)=max(U.boundaryField()[patchID][faceI].component(1),150); } thanks. 

March 27, 2013, 11:29 

#17 
Senior Member
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,209
Rep Power: 20 
when compiled it this error occurred:
Code:
ehsan@Ehsancom:~/Desktop/rhoPimpleFoamLimited$ wmakeMaking dependency list for source file rhoPimpleFoamLimited.C SOURCE=rhoPimpleFoamLimited.C ; g++ m64 Dlinux64 DWM_DP Wall Wextra Wnounusedparameter Woldstylecast Wnonvirtualdtor O3 DNoRepository ftemplatedepth100 I/opt/openfoam211/src/thermophysicalModels/basic/lnInclude I/opt/openfoam211/src/turbulenceModels/compressible/turbulenceModel I/opt/openfoam211/src/finiteVolume/cfdTools I/opt/openfoam211/src/finiteVolume/lnInclude IlnInclude I. I/opt/openfoam211/src/OpenFOAM/lnInclude I/opt/openfoam211/src/OSspecific/POSIX/lnInclude fPIC c $SOURCE o Make/linux64GccDPOpt/rhoPimpleFoamLimited.o rhoPimpleFoamLimited.C: In function ‘int main(int, char**)’: rhoPimpleFoamLimited.C:90:79: error: ‘Foam::T’ does not have class type rhoPimpleFoamLimited.C:90:96: error: ‘patchID’ was not declared in this scope rhoPimpleFoamLimited.C:90:105: error: ‘faceI’ was not declared in this scope rhoPimpleFoamLimited.C:94:41: error: name lookup of ‘cellI’ changed for ISO ‘for’ scoping [fpermissive] rhoPimpleFoamLimited.C:94:41: note: (if you use ‘fpermissive’ G++ will accept your code) rhoPimpleFoamLimited.C:100:132: error: ‘Foam::T’ does not have class type make: *** [Make/linux64GccDPOpt/rhoPimpleFoamLimited.o] Error 1 

March 27, 2013, 11:30 

#18 
Member
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 7 
Actually I don't know if OF will allow you such a statement:
Code:
U[cellI].component(0)= min(U[cellI].component(0), sqrt(1.4*287.14*T.boundaryField()[patchID][faceI])30); I don't know if you can use something like this: Code:
U[cellI].component(0)= min(U[cellI].component(0), sqrt(1.4*287.14*T.[cellI])30); Hope this help. 

March 27, 2013, 11:37 

#19 
Senior Member
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,209
Rep Power: 20 
thank you.
yes Simone you understood well. I modified it as so: Code:
forAll(U,cellI) { U[cellI].component(0)=min(U[cellI].component(0), sqrt(1.4*287.14*T.[cellI])30); } forAll(U,cellJ) { U[cellJ].component(1)=max(U[cellI].component(1),150); } word patchName = "left"; label patchID = mesh.boundary().findPatchID(patchName); forAll(U.boundaryField()[patchID],faceI) { U.boundaryField()[patchID][faceI].component(0)=min(U.boundaryField()[patchID][faceI].component(0), sqrt(1.4*287.14*T.boundaryField()[patchID][faceI])30); U.boundaryField()[patchID][faceI].component(1)=max(U.boundaryField()[patchID][faceI].component(1),150); } Code:
ehsan@Ehsancom:~/Desktop/rhoPimpleFoamLimited$ wmakeMaking dependency list for source file rhoPimpleFoamLimited.C SOURCE=rhoPimpleFoamLimited.C ; g++ m64 Dlinux64 DWM_DP Wall Wextra Wnounusedparameter Woldstylecast Wnonvirtualdtor O3 DNoRepository ftemplatedepth100 I/opt/openfoam211/src/thermophysicalModels/basic/lnInclude I/opt/openfoam211/src/turbulenceModels/compressible/turbulenceModel I/opt/openfoam211/src/finiteVolume/cfdTools I/opt/openfoam211/src/finiteVolume/lnInclude IlnInclude I. I/opt/openfoam211/src/OpenFOAM/lnInclude I/opt/openfoam211/src/OSspecific/POSIX/lnInclude fPIC c $SOURCE o Make/linux64GccDPOpt/rhoPimpleFoamLimited.o rhoPimpleFoamLimited.C: In function ‘int main(int, char**)’: rhoPimpleFoamLimited.C:90:79: error: ‘Foam::T’ does not have class type rhoPimpleFoamLimited.C:90:96: error: ‘patchID’ was not declared in this scope rhoPimpleFoamLimited.C:90:105: error: ‘faceI’ was not declared in this scope rhoPimpleFoamLimited.C:94:41: error: name lookup of ‘cellI’ changed for ISO ‘for’ scoping [fpermissive] rhoPimpleFoamLimited.C:94:41: note: (if you use ‘fpermissive’ G++ will accept your code) rhoPimpleFoamLimited.C:100:132: error: ‘Foam::T’ does not have class type make: *** [Make/linux64GccDPOpt/rhoPimpleFoamLimited.o] Error 1 

March 27, 2013, 11:40 

#20 
Member
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 7 
it's
Code:
T[cellI] Code:
T.[cellI] 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
UDF error  parabolic velocity profile  3D turbine  Zaqie  Fluent UDF and Scheme Programming  9  June 25, 2016 19:08 
unable to get parabolic velocity profile with pimplefoam  houkensjtu  OpenFOAM  4  October 8, 2012 04:41 
Emergency:UDF for a time dependent parabolic velocity  zumaqiong  Fluent UDF and Scheme Programming  12  March 25, 2010 13:00 
Neumann pressure BC and velocity field  Antech  Main CFD Forum  0  April 25, 2006 02:15 
what the result is negatif pressure at inlet  chong chee nan  FLUENT  0  December 29, 2001 06:13 