[Raza Javed]

Hello Everyone,


I have few questions related to relaxation factor


1. What is the relaxation factor?
2. How it helps in the simulation of our case?
3. Based on what, we modify the relaxation factor for each parameter?


I searched for these doubts, but couldn't find much about it.



I shall be thankful, if someone can guide me through this.


Thank you

[hjasak]

Pick up the old book by Patankar - he explains it well

[ingraban]

A second sub-dictionary of fvSolution that is often used in OpenFOAM is relaxationFactors which controls under-relaxation, a technique used for improving stability of a computation, particularly in solving steady-state problems. Under-relaxation works by limiting the amount which a variable changes from one iteration to the next, either by modifying the solution matrix and source prior to solving for a field or by modifying the field directly. An under-relaxation factor α,0 < α ≤ 1 specifies the amount of under-relaxation.

• No specified α: no under-relaxation.
• α = 1: guaranteed matrix diagonal equality/dominance.
• α decreases, under-relaxation increases.
• α = 0: solution does not change with successive iterations.

[yambanshee]

To follow on from the above the relaxation factor blends the previous iteration's solution and the current iteration's solution. If the factor is 1, then the new solution is entirely used, if the relaxation factor is 0, the old solution is entirely used. A simple example of a single cell could be:

iteration 1: p = 10
iteration 2: p is found to be 30, but relaxation factor is 0.6. So P = 0.6(30)+0.4(10) = 22

It's use is to ensure that the solution does not change too wildly between iterations, which could cause instability else wise.

Theoretically, as the solution converges, the difference in the answer between two iterations are minimal, and so the effect the relaxation factor has also becomes minimal. Therefor it shouldn't interfere with the answer of the final solution. It may, however, cause the solution to take more iterations to converge, as you are telling the cfd code to change less between each iteration. It is therefor ideal to find a relaxation factor that is as large as possible, whilst causing a converging and stable solution

[Raza Javed]

Originally Posted by yambanshee

To follow on from the above the relaxation factor blends the previous iteration's solution and the current iteration's solution. If the factor is 1, then the new solution is entirely used, if the relaxation factor is 0, the old solution is entirely used. A simple example of a single cell could be:

iteration 1: p = 10
iteration 2: p is found to be 30, but relaxation factor is 0.6. So P = 0.6(30)+0.4(10) = 22

It's use is to ensure that the solution does not change too wildly between iterations, which could cause instability else wise.

Theoretically, as the solution converges, the difference in the answer between two iterations are minimal, and so the effect the relaxation factor has also becomes minimal. Therefor it shouldn't interfere with the answer of the final solution. It may, however, cause the solution to take more iterations to converge, as you are telling the cfd code to change less between each iteration. It is therefor ideal to find a relaxation factor that is as large as possible, whilst causing a converging and stable solution

Hi,


Indeed, it is an excellent explanation of relaxation factor. Thank you so much.


Regards,
Raza

[Raza Javed]

Originally Posted by yambanshee

To follow on from the above the relaxation factor blends the previous iteration's solution and the current iteration's solution. If the factor is 1, then the new solution is entirely used, if the relaxation factor is 0, the old solution is entirely used. A simple example of a single cell could be:

iteration 1: p = 10
iteration 2: p is found to be 30, but relaxation factor is 0.6. So P = 0.6(30)+0.4(10) = 22

It's use is to ensure that the solution does not change too wildly between iterations, which could cause instability else wise.

Theoretically, as the solution converges, the difference in the answer between two iterations are minimal, and so the effect the relaxation factor has also becomes minimal. Therefor it shouldn't interfere with the answer of the final solution. It may, however, cause the solution to take more iterations to converge, as you are telling the cfd code to change less between each iteration. It is therefor ideal to find a relaxation factor that is as large as possible, whilst causing a converging and stable solution

Hi,


I have one more question:


Can we change the relaxation factor while running the simulation?


Thank you

[jherb]

Yes, you can change it during runtime and it will be applied in the next time step (if you have set runTimeModifiable true; in system/controlDict)

[Raza Javed]

Originally Posted by jherb

Yes, you can change it during runtime and it will be applied in the next time step (if you have set runTimeModifiable true; in system/controlDict)

Hi,


Thank you so much for your reply.



Actually I am having a problem of "maximum number of iterations exceeded" error in my simulation. When I reduce the relaxation factor from 0.7 to 1e-3 or 1e-4, then the solver RUNS for longer time but at the end it again gives the same error and stops.


I try to explain my case, may be you can help.


I am using chtMultiRegionSimpleFoam with OpenFoam 4.1.


my geomtetry has multiple regions (a box containing some heaters and one fluid region).


There is turbulent flow in the fluid region and its compressible also. I am using kEpsilon model.


Now, during the flow from inlet to outlet, the fluid should get heat up due to the heaters and the temperature of fluid at the outlet should be equal to the heat transferred from heaters to fluid.


After some iterations, the temperature of fluid becomes unrealistic and the solver stops.


The strange thing is that, for the same case, I have one another geometry, the only difference is that the other geometry the fluid region is straight, and in the first geometry, the fluid region has multiple curves in it.



And for the straight pipe geometry, everything is working fine.



The only different thing in both these cases, is the geometry, everything else is same (boundary conditions, initial conditions, thermophysical properties etc.)


Can you please help me, that what could be the possible reasons for this error.


The error is also attached below:

--> FOAM FATAL ERROR: 
[3] Maximum number of iterations exceeded
[3] 
[3]     From function Foam::scalar Foam::species::thermo<Thermo, Type>::T(Foam::scalar, Foam::scalar, Foam::scalar, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar) const) const [with Thermo = Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie> >; Type = Foam::sensibleEnthalpy; Foam::scalar = double; Foam::species::thermo<Thermo, Type> = Foam::species::thermo<Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie> >, Foam::sensibleEnthalpy>]
[3]     in file /home/ubuntu/OpenFOAM/OpenFOAM-4.1/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 66.
[3] 
FOAM parallel run aborting
[3] 
[3] #0  Foam::error::printStack(Foam::Ostream&) at ??:?
[3] #1  Foam::error::abort() at ??:?
[3] #2  Foam::heRhoThermo<Foam::rhoThermo, Foam::pureMixture<Foam::polynomialTransport<Foam::species::thermo<Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie, 8>, 8>, Foam::sensibleEnthalpy>, 8> > >::calculate() at ??:?
[3] #3  Foam::heRhoThermo<Foam::rhoThermo, Foam::pureMixture<Foam::polynomialTransport<Foam::species::thermo<Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie, 8>, 8>, Foam::sensibleEnthalpy>, 8> > >::correct() at ??:?
[3] #4  ? at ??:?
[3] #5  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
[3] #6  ? at ??:?
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 3 in communicator MPI_COMM_WORLD 
with errorcode 1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------

Thank you

[jherb]

Relaxation factors in the order of magnitude of 1e-3 sound very wrong. There are probably other problems with your simulation.


Does checkMesh complain about something?



Are you sure, your boundary conditions are physically valid?

[Raza Javed]

Thank you so much for your reply.

Originally Posted by jherb

Does checkMesh complain about something?

I am attaching my checkMesh below and it seems OK, but can you please have a look.

Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           1875576
    faces:            21909869
    internal faces:   21895815
    cells:            10951421
    faces per cell:   4
    boundary patches: 3
    point zones:      0
    face zones:       0
    cell zones:       10

Overall number of cells of each type:
    hexahedra:     0
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    10951421
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
    Patch               Faces    Points   Surface topology                  
    inlet               264      164      ok (non-closed singly connected)  
    outlet              234      149      ok (non-closed singly connected)  
    defaultFaces        13556    6840     ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-0.140981 -0.096774 0) (0.132 0.096774 0.03)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (2.35742e-18 3.2066e-18 -8.15718e-17) OK.
    Max cell openness = 2.85893e-16 OK.
    Max aspect ratio = 10.056 OK.
    Minimum face area = 5.60659e-08. Maximum face area = 3.6484e-05.  Face area magnitudes OK.
    Min volume = 7.32683e-12. Max volume = 4.43234e-08.  Total volume = 0.00158505.  Cell volumes OK.
    Mesh non-orthogonality Max: 64.003 average: 15.2411
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.885953 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End

Originally Posted by jherb

Are you sure, your boundary conditions are physically valid?

I tried putting boundary conditions according to my requirement, but I am not sure either they are correct or NOT.


I want to have the velocity of 1.33 at inlet, and also as an initial value.
For outlet, I don't have any specification for velocity, so I put zeroGradient.


Temperature:


at inlet I defined 300K, and at outlet I need to check the temperature(should be higher compared to inlet due to the heaters outside the fluid region).


But for alphat, k, epsilon and nut I have put boundary conditions from tutorials, but I am not sure about them.


I am attaching my boundary conditions below:

boundary
{
    inlet
    {
        type            patch;
    }
    outlet
    {
        type            patch;
    }
}

T
{
    internalField   uniform 300;

    boundaryField
    {
        inlet
        {
            type            fixedValue;
            value           $internalField;
            
        }

        outlet
        {
            type            zeroGradient;

        }

        "fluid_to_box"
        {

            type            compressible::turbulentTemperatureCoupledBaffleMixed;
            Tnbr            T;
            kappaMethod     fluidThermo;
            value           uniform 300;
        }
    }
}

U
{
    internalField   uniform (1.33 0 0);

    boundaryField
    {
        inlet
        {
            type            fixedValue;
            value           uniform (1.33 0 0);
        }

        outlet
        {                   
            type            zeroGradient;
            
        }
        "fluid_to_box"
        {
            type            noSlip;
        }
    }
}

epsilon
{
    internalField   uniform 0.064;

    boundaryField
    {
        inlet
        {
            type            fixedValue;
            value           uniform 0.064;
        }

        outlet
        {
            type            zeroGradient;
        }

        ".*"
        {
            type            epsilonWallFunction;
            value           uniform 0.064;
        }
    }
}

k
{
    internalField   uniform 2.45e-3;

    boundaryField
    {
        inlet
        {
            type            fixedValue;
            value           uniform 2.45e-3;
        }

        outlet
        {
            type            zeroGradient;
            //value           uniform 0;
        }

        "fluid_to_box"
        {
            type            kqRWallFunction;
            value           uniform 2.45e-3;
        }
    }
}


p_rgh
{
    internalField   uniform 0;

    boundaryField
    {
        inlet
        {
            type            zeroGradient;
            value           uniform 0;
        }

        outlet
        {
            type            fixedValue;
            value           uniform 0;
        }

        ".*"
        {
            type            fixedFluxPressure;
            value           uniform 0;
        }
    }
}

p
{
    internalField   uniform 100000;

    boundaryField
    {
        ".*"
        {
            type            calculated;
            value           $internalField;
        }
    }
}

alphat
{
    internalField   uniform 0;

    boundaryField
    {
        inlet
        {
            type            calculated;
            value           uniform 0;
        }

        outlet
        {
            type            calculated;
            value           uniform 0;
        }

        "fluid_to_box"
        {
            type            compressible::alphatWallFunction;

        }

    }
}

nut
{
    internalField   uniform 0;

    boundaryField
    {
        inlet
        {
            type            calculated;
            value           uniform 0;
        }

        outlet
        {
            type            calculated;
            value           uniform 0;
        }

        "fluid_to_box"
        {
            type            nutkWallFunction;
            value           uniform 0;
        }

    }
}

// ************************************************************************* //

I shall be very thankful, if you can guide me through this.


Thank you

[jherb]

One thing is probably not correct. You use chtMultiRegionSimpleFoam? Then you need to set a physically correct pressure, because it is a "compressible" solver, calculating the density from some equation of state. So p = 0 at the outlet cannot be correct.



Actually, what is the message of the solver, when it crashes?

[Raza Javed]

Originally Posted by jherb

Then you need to set a physically correct pressure, because it is a "compressible" solver, calculating the density from some equation of state. So p = 0 at the outlet cannot be correct.

I am putting pressure (p) equals to calculated on all boundaries, but yes p_rgh, I am putting zero at outlet. Do I need to change the boundary conditions for p_rgh?

Originally Posted by jherb

Actually, what is the message of the solver, when it crashes?

The log file of the solver is given below:

Time = 0.1


Solving for fluid region fluid
DILUPBiCG:  Solving for Ux, Initial residual = 1, Final residual = 0.00248155, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 1, Final residual = 0.00420011, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 1, Final residual = 0.00300392, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 0.999619, Final residual = 0.00174638, No Iterations 1
Min/max T:299.968 300.002
GAMG:  Solving for p_rgh, Initial residual = 1, Final residual = 0.0294703, No Iterations 10
time step continuity errors : sum local = 1833.44, global = 63.8567, cumulative = 63.8567
Min/max rho:1000 1000
DILUPBiCG:  Solving for epsilon, Initial residual = 0.0127737, Final residual = 1.03336e-05, No Iterations 1
bounding epsilon, min: -6.47425 max: 198.948 average: 0.11565
DILUPBiCG:  Solving for k, Initial residual = 1, Final residual = 0.00132192, No Iterations 1
bounding k, min: -0.0598646 max: 1.36506 average: 0.00247787

Solving for solid region box
DICPCG:  Solving for h, Initial residual = 0.999899, Final residual = 0.0537438, No Iterations 2
Min/max T:299.999 300

Solving for solid region p1
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.0154276, No Iterations 2
Min/max T:300 300.056

Solving for solid region p2
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.0150834, No Iterations 2
Min/max T:300 300.044

Solving for solid region p3
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.0139625, No Iterations 2
Min/max T:300 300.056

Solving for solid region p4
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.0142509, No Iterations 2
Min/max T:300 300.047

Solving for solid region p5
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.0159586, No Iterations 2
Min/max T:300 300.049

Solving for solid region p6
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.0153529, No Iterations 2
Min/max T:300 300.056

Solving for solid region p7
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.0134462, No Iterations 2
Min/max T:300 300.05

Solving for solid region p8
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.0144757, No Iterations 2
Min/max T:300 300.053
ExecutionTime = 133.06 s  ClockTime = 134 s

Time = 0.2


Solving for fluid region fluid
DILUPBiCG:  Solving for Ux, Initial residual = 0.504383, Final residual = 0.00162629, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.300932, Final residual = 0.00101981, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.270728, Final residual = 0.000718876, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 0.316352, Final residual = 0.000658972, No Iterations 1
Min/max T:299.961 300.002
GAMG:  Solving for p_rgh, Initial residual = 0.776721, Final residual = 0.0236817, No Iterations 10
time step continuity errors : sum local = 2063.36, global = 36.928, cumulative = 100.785
Min/max rho:1000 1000
DILUPBiCG:  Solving for epsilon, Initial residual = 0.480833, Final residual = 6.92497e-05, No Iterations 1
bounding epsilon, min: -102.556 max: 5835.89 average: 0.137401
DILUPBiCG:  Solving for k, Initial residual = 0.844893, Final residual = 0.000528436, No Iterations 1
bounding k, min: -0.451503 max: 51.2525 average: 0.00274481

Solving for solid region box
DICPCG:  Solving for h, Initial residual = 0.9999, Final residual = 0.02301, No Iterations 2
Min/max T:299.993 300.015

Solving for solid region p1
DICPCG:  Solving for h, Initial residual = 0.689583, Final residual = 0.00986339, No Iterations 2
Min/max T:300.006 300.096

Solving for solid region p2
DICPCG:  Solving for h, Initial residual = 0.698901, Final residual = 0.00942908, No Iterations 2
Min/max T:300.006 300.078

Solving for solid region p3
DICPCG:  Solving for h, Initial residual = 0.679693, Final residual = 0.00839357, No Iterations 2
Min/max T:300.005 300.1

Solving for solid region p4
DICPCG:  Solving for h, Initial residual = 0.692854, Final residual = 0.00897344, No Iterations 2
Min/max T:300.006 300.085

Solving for solid region p5
DICPCG:  Solving for h, Initial residual = 0.686304, Final residual = 0.0101575, No Iterations 2
Min/max T:300.006 300.088

Solving for solid region p6
DICPCG:  Solving for h, Initial residual = 0.68812, Final residual = 0.00966565, No Iterations 2
Min/max T:300.006 300.099

Solving for solid region p7
DICPCG:  Solving for h, Initial residual = 0.696385, Final residual = 0.00854814, No Iterations 2
Min/max T:300.006 300.09

Solving for solid region p8
DICPCG:  Solving for h, Initial residual = 0.693705, Final residual = 0.00922509, No Iterations 2
Min/max T:300.007 300.095
ExecutionTime = 145.92 s  ClockTime = 147 s

Time = 0.3


Solving for fluid region fluid
DILUPBiCG:  Solving for Ux, Initial residual = 1, Final residual = 0.00696453, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 1, Final residual = 0.00537407, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 1, Final residual = 0.00773095, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 1, Final residual = 0.00227395, No Iterations 1
[3] 
[3] 
[3] --> FOAM FATAL ERROR: 
[3] Maximum number of iterations exceeded
[3] 
[3]     From function Foam::scalar Foam::species::thermo<Thermo, Type>::T(Foam::scalar, Foam::scalar, Foam::scalar, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar) const) const [with Thermo = Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie> >; Type = Foam::sensibleEnthalpy; Foam::scalar = double; Foam::species::thermo<Thermo, Type> = Foam::species::thermo<Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie> >, Foam::sensibleEnthalpy>]
[3]     in file /home/ubuntu/OpenFOAM/OpenFOAM-4.1/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 66.
[3] 
FOAM parallel run aborting
[3] 
[3] #0  Foam::error::printStack(Foam::Ostream&) at ??:?
[3] #1  Foam::error::abort() at ??:?
[3] #2  Foam::heRhoThermo<Foam::rhoThermo, Foam::pureMixture<Foam::polynomialTransport<Foam::species::thermo<Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie, 8>, 8>, Foam::sensibleEnthalpy>, 8> > >::calculate() at ??:?
[3] #3  Foam::heRhoThermo<Foam::rhoThermo, Foam::pureMixture<Foam::polynomialTransport<Foam::species::thermo<Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie, 8>, 8>, Foam::sensibleEnthalpy>, 8> > >::correct() at ??:?
[3] #4  ? at ??:?
[3] #5  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
[3] #6  ? at ??:?
--------------------------------------------------------------------------

please let me know if you need any further information.


Thank you