[hhuang84]

I am trying to simulate fully developed turbulent flow in a square duct with cyclic BC in simpleFoam with K-Epsilon model. The mesh was generated in pointwise, with cyclic BCs automatically specified. here attach my setup files, but it fail to get a converged solution. Any hits? Thanks in advance.
0/U &p
dimensions [0 1 -1 0 0 0 0];

internalField uniform (0.329 0 0);

boundaryField
{
Cyclic
{
type cyclic;
value uniform (0.329 0 0);
}

Wall
{
type fixedValue;
value uniform (0 0 0);
}

}
dimensions [0 2 -2 0 0 0 0];

internalField uniform 0;

boundaryField
{
Cyclic
{
type cyclic;
}

Wall
{
type zeroGradient;
}


}

and I also set Ubar in constant/transportProperties
Ubar Ubar [ 0 1 -1 0 0 0 0 ] (0.329 0 0 );
transportModel Newtonian;

nu nu [ 0 2 -1 0 0 0 0 ] 1.003e-06;
in system/fvSchemes I set
SIMPLE
{

nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
I has been trying for days, but got no good results, even if I turned off turbulence. the maximum stream velocity from simpleFoam is 0.329, so the mass is not conserved.

[lakeat]

How about this,

Cyclic 
{
type cyclic;
}

and why did you use simplefoam to simulate a fully developed turbulence???

[hhuang84]

Any other solvers that can implement the periodic boundary conditions?

[santos]

Have you looked into channelFoam? It can handle cyclics, while simpleFoam doesnt by default (I think).

Regards,
Jose Santos

[hhuang84]

However, if I just wanna run a laminar or RANS , what solver can I chose?

[lakeat]

Try pisoFoam.

[santos]

For laminar, change LESModel in channelFoam to laminar. For RANS I guess you need to make your own solver, combining eg pisoFoam (as suggested above) and channelFoam.

Regards,
Jose Santos

[ouafa]

hi everybody,

i try to run the biperiodic tutorial channelOoddles initially proposed in LES, with the solver simpleFoam for steady RANS using the K-omega SST model and taking the same parameters Ubar and nu.

so i add the term gradP like in the channelOoddles solver in the Ueqn but the result was a laminar regime even after long time simulations. k and omega are not zero. can someone help me please, here after my proposed:

fvVectorMatrix UEqn
(
fvm::div(phi, U)
+ turbulence->divDevReff(U)
==
flowDirection*gradP
);

UEqn.relax();

eqnResidual = solve
(
UEqn == -fvc::grad(p)
).initialResidual();

maxResidual = max(eqnResidual, maxResidual);

-----------------------------------------------------------------------------

p.boundaryField().updateCoeffs();

volScalarField AU =UEqn.A();
volScalarField rUA = 1.0/UEqn.A();

U = UEqn.H()/AU;

phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, p);

// Non-orthogonal pressure corrector loop
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(1.0/AU, p) == fvc::div(phi)
);

pEqn.setReference(pRefCell, pRefValue);
// retain the residual from the first iteration
if (nonOrth == 0)
{
eqnResidual = pEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
}
else
{
pEqn.solve();
}

if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}

# include "continuityErrs.H"

// Explicitly relax pressure for momentum corrector
p.relax();

// Momentum corrector
U -= fvc::grad(p)/AU;
U.correctBoundaryConditions();


// Correct driving force for a constant mass flow rate

// Extract the velocity in the flow direction
dimensionedScalar magUbarStar =
(flowDirection & U)().weightedAverage(mesh.V());

// Calculate the pressure gradient increment needed to
// adjust the average flow-rate to the correct value
dimensionedScalar gragPplus =
(magUbar - magUbarStar)/rUA.weightedAverage(mesh.V());

U += flowDirection*rUA*gragPplus;

gradP += gragPplus;

Info<< "Uncorrected Ubar = " << magUbarStar.value() << tab
<< "pressure gradient = " << gradP.value() << endl;

[chegdan]

This is an old thread, but its close to what I would like to discuss. I have made some changes to simpleFoam using channelFoam as a guide to model steady-state channel flow using RANS. I noticed this thread that was left hanging....has anyone managed to get periodic boundary conditions in simpleFoam to work WITHOUT direct mapped patches? Ubar and gradP diverge after a few hundred iterations. Thoughts? Below is my code so far.

Dan

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Application
    simpleChannelFoam

Description
    Steady-state solver for incompressible, laminar or turbulent flow with 
    periodic boundary conditions

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "IFstream.H"  //added
#include "OFstream.H"  //added
#include "RASModel.H"

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

int main(int argc, char *argv[])
{
#   include "setRootCase.H"
#   include "createTime.H"
#   include "createMesh.H"
#   include "readTransportProperties.H" //added from channelFoam
#   include "createFields.H"
#   include "initContinuityErrs.H"
#   include "createGradP.H" //added from channelFoam

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

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

#       include "readSIMPLEControls.H"
#       include "initConvergenceCheck.H"

        p.storePrevIter();

        // Pressure-velocity SIMPLE corrector
	
	turbulence->correct();//moved here

//solve the U-equation

	    tmp<fvVectorMatrix> UEqn
	    (
	        fvm::div(phi, U)
	      + turbulence->divDevReff(U)
	      ==
	        flowDirection*gradP
	    );
	
	    UEqn().relax();
	
	    eqnResidual = solve
	    (
	        UEqn() == -fvc::grad(p)
	    ).initialResidual();
	
	    maxResidual = max(eqnResidual, maxResidual);

//end solving the U-equation
//solve the p-equation with SIMPLE Algorithm
	    p.boundaryField().updateCoeffs();

	    volScalarField rUA = 1.0/UEqn().A();//added
	    U = rUA*UEqn().H();
	    UEqn.clear();
	    phi = fvc::interpolate(U) & mesh.Sf();
	    adjustPhi(phi, U, p);

	    // Non-orthogonal pressure corrector loop
	    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
	    {
	        fvScalarMatrix pEqn
        	(
	            fvm::laplacian(rUA, p) == fvc::div(phi)
        	);

        	pEqn.setReference(pRefCell, pRefValue);

        	// Retain the residual from the first iteration
	        if (nonOrth == 0)
        	{
        	    eqnResidual = pEqn.solve().initialResidual();
        	    maxResidual = max(eqnResidual, maxResidual);
        	}
        	else
        	{
        	    pEqn.solve();
        	}

        	if (nonOrth == nNonOrthCorr)
        	{
        	    phi -= pEqn.flux();
        	}
    	   }

	#   include "continuityErrs.H"

	    // Explicitly relax pressure for momentum corrector
	    p.relax();

	    // Momentum corrector
	    U -= rUA*fvc::grad(p);
	    U.correctBoundaryConditions();

//end solve the p-equation with SIMPLE Algorithm

	// Correct driving force for a constant mass flow rate

        // Extract the velocity in the flow direction
        dimensionedScalar magUbarStar =
            (flowDirection & U)().weightedAverage(mesh.V());

        // Calculate the pressure gradient increment needed to
        // adjust the average flow-rate to the correct value
        dimensionedScalar gragPplus =
            (magUbar - magUbarStar)/rUA.weightedAverage(mesh.V());

        U += flowDirection*rUA*gragPplus;

        gradP += gragPplus;

        //- Solve the turbulence equations and correct the turbulence viscosity after U has been updated
        //turbulence->correct();

        Info<< "Uncorrected Ubar = " << magUbarStar.value() << tab
            << "pressure gradient = " << gradP.value() << endl;

        runTime.write();

        #include "writeGradP.H"	//added from channelFoam

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;

#       include "convergenceCheck.H"
    }

    Info<< "End\n" << endl;

    return 0;
}


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

[tom ato]

Dear Dan,

in the thread above about „channel flow with simpleFoam aka simpleChannelFoam“ you wrote
about having problems with convergence of Ubar and gradP.

I'm having the same problem right now (not surprising since I'm using a quite similar solver)
and I'm not able to find the problem.

Have you been able to identify or even solve the problem of the solver you posted?
It would be great to get some feedback :-)

Kind Regards

Tom

[chegdan]

Originally Posted by tom ato

Dear Dan,

in the thread above about „channel flow with simpleFoam aka simpleChannelFoam“ you wrote
about having problems with convergence of Ubar and gradP.

I'm having the same problem right now (not surprising since I'm using a quite similar solver)
and I'm not able to find the problem.

Have you been able to identify or even solve the problem of the solver you posted?
It would be great to get some feedback :-)

Kind Regards

Tom

Dear Tom,

Not yet.

Regards, Dan

[tom ato]

Dear Dan,
Dear Foamers,

the solver posted above by Dan should be correct. A very similar solver
I used is running quite fine by considering the following aspects:

1. it is absolutely vital to use an upwind scheme for solving U (can be changed in the fvSchemes-dictionary if not set by default)

2. relaxationFactors for pressure and velocity must be reduced and should not be bigger than

p 0.15;
U 0.5;

(The relaxation factors can be changed in the fvSolution dictionary)

3. one should take at least 8000 iteration-steps ( set in the controldict-file), although in some cases 4000 iterations might be sufficient as well. I made best experiences with 10000 steps.


In particular, I posted this for Dan in the hope of giving him, and anybody else
having the same problem and reading this thread, an helpful answer, but I fear this is coming to late.

Kind Regards

Tom

[chegdan]

Tom,

Thanks. I'll give it a try and If I can add anything I'll get back here and post it.

Thanks.

Dan

[Mirage]

Originally Posted by chegdan

This is an old thread, but its close to what I would like to discuss. I have made some changes to simpleFoam using channelFoam as a guide to model steady-state channel flow using RANS. I noticed this thread that was left hanging....has anyone managed to get periodic boundary conditions in simpleFoam to work WITHOUT direct mapped patches? Ubar and gradP diverge after a few hundred iterations. Thoughts? Below is my code so far.

Dan

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Application
    simpleChannelFoam

Description
    Steady-state solver for incompressible, laminar or turbulent flow with 
    periodic boundary conditions

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "IFstream.H"  //added
#include "OFstream.H"  //added
#include "RASModel.H"

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

int main(int argc, char *argv[])
{
#   include "setRootCase.H"
#   include "createTime.H"
#   include "createMesh.H"
#   include "readTransportProperties.H" //added from channelFoam
#   include "createFields.H"
#   include "initContinuityErrs.H"
#   include "createGradP.H" //added from channelFoam

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

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

#       include "readSIMPLEControls.H"
#       include "initConvergenceCheck.H"

        p.storePrevIter();

        // Pressure-velocity SIMPLE corrector
    
    turbulence->correct();//moved here

//solve the U-equation

        tmp<fvVectorMatrix> UEqn
        (
            fvm::div(phi, U)
          + turbulence->divDevReff(U)
          ==
            flowDirection*gradP
        );
    
        UEqn().relax();
    
        eqnResidual = solve
        (
            UEqn() == -fvc::grad(p)
        ).initialResidual();
    
        maxResidual = max(eqnResidual, maxResidual);

//end solving the U-equation
//solve the p-equation with SIMPLE Algorithm
        p.boundaryField().updateCoeffs();

        volScalarField rUA = 1.0/UEqn().A();//added
        U = rUA*UEqn().H();
        UEqn.clear();
        phi = fvc::interpolate(U) & mesh.Sf();
        adjustPhi(phi, U, p);

        // Non-orthogonal pressure corrector loop
        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
        {
            fvScalarMatrix pEqn
            (
                fvm::laplacian(rUA, p) == fvc::div(phi)
            );

            pEqn.setReference(pRefCell, pRefValue);

            // Retain the residual from the first iteration
            if (nonOrth == 0)
            {
                eqnResidual = pEqn.solve().initialResidual();
                maxResidual = max(eqnResidual, maxResidual);
            }
            else
            {
                pEqn.solve();
            }

            if (nonOrth == nNonOrthCorr)
            {
                phi -= pEqn.flux();
            }
           }

    #   include "continuityErrs.H"

        // Explicitly relax pressure for momentum corrector
        p.relax();

        // Momentum corrector
        U -= rUA*fvc::grad(p);
        U.correctBoundaryConditions();

//end solve the p-equation with SIMPLE Algorithm

    // Correct driving force for a constant mass flow rate

        // Extract the velocity in the flow direction
        dimensionedScalar magUbarStar =
            (flowDirection & U)().weightedAverage(mesh.V());

        // Calculate the pressure gradient increment needed to
        // adjust the average flow-rate to the correct value
        dimensionedScalar gragPplus =
            (magUbar - magUbarStar)/rUA.weightedAverage(mesh.V());

        U += flowDirection*rUA*gragPplus;

        gradP += gragPplus;

        //- Solve the turbulence equations and correct the turbulence viscosity after U has been updated
        //turbulence->correct();

        Info<< "Uncorrected Ubar = " << magUbarStar.value() << tab
            << "pressure gradient = " << gradP.value() << endl;

        runTime.write();

        #include "writeGradP.H"    //added from channelFoam

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;

#       include "convergenceCheck.H"
    }

    Info<< "End\n" << endl;

    return 0;
}


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

Thanks for sharing your code

I d like to simulate a cyclic flow inside a ribbed Channel and i 'd like to know, if i should modify my solver to get the cyclic behavior. I thought, that I have just to use the cyclic BC in 0 in the case folder to simulate a cyclic flow. Are U editing the code to solve a convergence problem? I don't really understand the added line code (colored in red). Could you please briefly explain the meaning of the equations Or could U post a link, where I can read more about what U did Thank you

For the other Foamers, please feel free to interact, if you could help

[tom ato]

Dear Mirage,

I'm sorry not to have answered you sooner, but my mail-adress changed and I got your answer just yesterday, again it seems that an answer in this thread is coming too late.

First of all, you have to be aware, that the entries in this post are related to OpenFOAM Versions 1.7 or lower, so with your current version you might not face the same problem. But I will try to share what is left on my mind from that time and how I remember the problem. The code posted above was about to solve the following problem during the solution of a cyclic channel:

A cyclic channel has the same velocity profile over all of its length. Anyway, the channel of infinite length has also a certain pressure drop which is constant depending on viscosity and the curvature of the wall (indeed, the pressure drop is the quantitiy driving the flow through the channel). Since you discretize the channel of infinite length with an equivalent channel of finite length and applying periodic boundary conditions, you will get different values for your pressure on the outlet than on the inlet caused by the pressure drop, the difference will be the pressure drop integrated over the length of your channel.

As I remember, what the solver did at that time by using periodic boundary conditions, was to take ALL(!) quantities from the outlet and prescribe it on the inlet and repeat the solution process until convergence (this is of course just a simplified way to describe the solution process to give you a short overview over the problem). So the pressure field in the iteration step n+1 at the inlet what the pressure field from the outlet at iteration n.

So with every single iteration step, the starting pressure level was decreased by the calculated pressure drop and after a certain number of steps the flow felt "asleep" and nothing happened anymore because pressure and driving pressure drop were zero. The code above fixed this problem by adding the pressure drop of the channel at the outlet before being described on the inlet of the next iteration step (at least this is was I remember, I didn't take the time to read the code again. Let me know, if I'm wrong here).

To your rather technical questions: I implemented this solver in my own openFoam version, but I cannot tell you, if or how this is handled in newer versions of openFoam. In the version above, you had to set periodic/cyclic boundary conditions as you described it and then you had a solution which was real fine.

About the red code lines: 1. the pressure gradient is implemented as the source term in the momentum predictor on the right hand side of the equation system (see the dissertation of Prof. Jasak or something like this for information about implicit pressure correction solvers...), the terms of the left hand side are convection and diffusion-terms of ico Navier-Stokes equations.

2. In the second red snippet the pressure is increased to enforce a constant mass flow (as mentioned in the comment on top). If you need details, let me know or refer to the userGuide, it's always worth trying to understand the code.

I hope I could help you after all this time, if not simply ignore it.

Kind regards,
Tom