CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Manual limiter of velocity doesn't work

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree9Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   March 13, 2013, 08:40
Default Manual limiter of velocity doesn't work
  #1
Member
 
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 4
batta31 is on a distinguished road
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.
batta31 is offline   Reply With Quote

Old   March 14, 2013, 04:22
Default
  #2
Member
 
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 4
batta31 is on a distinguished road
Please guy, is there someone that can help?
batta31 is offline   Reply With Quote

Old   March 14, 2013, 10:29
Default
  #3
Senior Member
 
Join Date: Nov 2009
Location: Michigan
Posts: 135
Rep Power: 7
doubtsincfd is on a distinguished road
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
doubtsincfd is offline   Reply With Quote

Old   March 15, 2013, 08:09
Default
  #4
Member
 
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 4
batta31 is on a distinguished road
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.
batta31 is offline   Reply With Quote

Old   March 15, 2013, 08:46
Default
  #5
Senior Member
 
immortality's Avatar
 
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,173
Rep Power: 16
immortality is on a distinguished road
hi Simone
I think limiting velocity can also resolve my problem.
Could you send me your modified solver?
Thanks.
immortality is offline   Reply With Quote

Old   March 16, 2013, 19:29
Default
  #6
Senior Member
 
immortality's Avatar
 
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,173
Rep Power: 16
immortality is on a distinguished road
where should i add them exactly?
immortality is offline   Reply With Quote

Old   March 17, 2013, 21:24
Default
  #7
Senior Member
 
immortality's Avatar
 
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,173
Rep Power: 16
immortality is on a distinguished road
can use velocity limiters only on a patch not entire the domain?
immortality is offline   Reply With Quote

Old   March 17, 2013, 21:26
Default
  #8
Senior Member
 
immortality's Avatar
 
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,173
Rep Power: 16
immortality is on a distinguished road
hi Omkar
Could you send me the code with added expressions for rhoPimpleFoam?
immortality is offline   Reply With Quote

Old   March 17, 2013, 21:37
Default
  #9
Senior Member
 
Join Date: Nov 2009
Location: Michigan
Posts: 135
Rep Power: 7
doubtsincfd is on a distinguished road
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
doubtsincfd is offline   Reply With Quote

Old   March 18, 2013, 05:04
Default
  #10
Member
 
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 4
batta31 is on a distinguished road
Hi immortality!
With respect to the adjointShapeOptimizationFoam solver, you should put the " two" limiter just above the lines

Quote:
// Explicitly relax pressure for adjoint momentum corrector
pa.relax();

// Adjoint momentum corrector
Ua -= rAUa*fvc::grad(pa);
Ua.correctBoundaryConditions();
into the predictor-corrector loop. At least, I did in this way!

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:
word patchName = "NAME_OF_THE_PATCH";
label patchID = mesh.boundary().findPatchID(patchName2);

forAll(U.boundaryField()[patchID],faceI)
{
U.boundaryField()[patchID][faceI].component(0)=some value;
}
Hope this works

Simone
immortality likes this.
batta31 is offline   Reply With Quote

Old   March 19, 2013, 14:57
Default
  #11
Senior Member
 
immortality's Avatar
 
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,173
Rep Power: 16
immortality is on a distinguished road
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"

        // --- Pressure-velocity 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;
}


// ************************************************************************* //
where UEqn.H is:
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);
}
hEqn.H:
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();
}
pEqn.H:
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);
immortality is offline   Reply With Quote

Old   March 21, 2013, 08:06
Post
  #12
Member
 
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 4
batta31 is on a distinguished road
In my opinion after

Code:
if (pimple.turbCorr())
{
    turbulence->correct();
}
should be fine.
batta31 is offline   Reply With Quote

Old   March 26, 2013, 09:57
Default
  #13
Senior Member
 
immortality's Avatar
 
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,173
Rep Power: 16
immortality is on a distinguished road
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);
           }
immortality is offline   Reply With Quote

Old   March 26, 2013, 12:49
Default
  #14
Senior Member
 
immortality's Avatar
 
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,173
Rep Power: 16
immortality is on a distinguished road
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);
}
immortality is offline   Reply With Quote

Old   March 27, 2013, 11:59
Default
  #15
Member
 
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 4
batta31 is on a distinguished road
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.
immortality likes this.
batta31 is offline   Reply With Quote

Old   March 27, 2013, 12:20
Default
  #16
Senior Member
 
immortality's Avatar
 
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,173
Rep Power: 16
immortality is on a distinguished road
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);
           }
another question! could we do both loops into single one by for instance I counter?
thanks.
immortality is offline   Reply With Quote

Old   March 27, 2013, 12:29
Default
  #17
Senior Member
 
immortality's Avatar
 
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,173
Rep Power: 16
immortality is on a distinguished road
when compiled it this error occurred:
Code:
ehsan@Ehsan-com:~/Desktop/rhoPimpleFoamLimited$ wmakeMaking dependency list for source file rhoPimpleFoamLimited.C
SOURCE=rhoPimpleFoamLimited.C ;  g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3  -DNoRepository -ftemplate-depth-100 -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
whats error in T?
immortality is offline   Reply With Quote

Old   March 27, 2013, 12:30
Default
  #18
Member
 
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 4
batta31 is on a distinguished road
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);
since you are looping inside the domain for the component of U you can't assign its value with respect to the value of the temperature on the boundary!

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);
You could try with the line above if I understand well your intention, i.e. you want to threshold the value of the velocity with respect to the local value of the Mach number.

Hope this help.
immortality likes this.
batta31 is offline   Reply With Quote

Old   March 27, 2013, 12:37
Default
  #19
Senior Member
 
immortality's Avatar
 
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,173
Rep Power: 16
immortality is on a distinguished road
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);
           }
but this error is shown:
Code:
ehsan@Ehsan-com:~/Desktop/rhoPimpleFoamLimited$ wmakeMaking dependency list for source file rhoPimpleFoamLimited.C
SOURCE=rhoPimpleFoamLimited.C ;  g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3  -DNoRepository -ftemplate-depth-100 -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
immortality is offline   Reply With Quote

Old   March 27, 2013, 12:40
Default
  #20
Member
 
Simone
Join Date: Sep 2012
Posts: 95
Rep Power: 4
batta31 is on a distinguished road
it's

Code:
T[cellI]
not

Code:
T.[cellI]
batta31 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


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


All times are GMT -4. The time now is 17:10.