2D isothermal cylinder not converging

UPengineer · March 9, 2014, 03:08

Hello Everyone,

I am using OpenFOAM version 2.2.2 with HelyxOS 2.0 on Ubuntu 13.10. I am looking at a 2D natural convective flow from an isothermally heated cylinder.

The mesh I have may be seen below:

There is a 1.5 m by 1.5 m box surrounding the cylinder with 4X refinement.

The boundary conditions are:

Code:

Wall -->Fixed, no slip - zero gradient temperature for: front face, back face, left face and right face.

Wall --> 350 K isothermal temperature for the cylinder

Patch --> Velocity: Pressure Inlet velocity
               Pressure: Total pressure = 0 m^2/s^2
               For the top face

Patch --> Velocity: Pressure Inlet Outlet
                   Pressure: Zero gradient
                  For the Bottom face

The Solution Modeling is:

Code:

SIMPLE FOAM algorithm
2nd Order Bounded Linear Upwind for velocity and temperatureSteady
Incompressible
Laminar
With thermal energy on and buoyancy on (gravity set to -9.81 m/s^2)

The residuals are:

As you can see, the residuals are anything but convergent which makes me think the boundary conditions may be a little messed up.

An example of the temperature post processing of the case is:

which appears to be correct, however because of the residuals I cannot trust this data.

Any advice on how to get my case to converge? Any help is appreciated.

jherb · March 9, 2014, 11:31

What solver are you using? With or without Boussinesq approximation? Have you tried the same setup with a transient simulation (pimple)?

Perhaps also the area below the cylinder is to short.

UPengineer · March 9, 2014, 15:30

Quote:

Originally Posted by jherb

What solver are you using? With or without Boussinesq approximation? Have you tried the same setup with a transient simulation (pimple)?

Perhaps also the area below the cylinder is to short.

Hello jherb,

Thank you for your help and quick reply.

I did as you suggested by extending the area below the cylinder.
Unfortunately, the residuals did not converge any better.

The mesh is now:

Which output residuals seen:

For your other questions/suggestions:

I am using buoyantBoussinesqSimpleFoam as the solver, so yes I am using the Boussinesq approximation.

Would transient simulation be too overkill since I do not have any time dependent variables? There is nothing cooling, heating, or filling a space. It is just the heat from an isothermally heated cylinder.

I hope this may help pinpoint what my problem is. If there is any other questions for my case setup that may help OpenFOAM run this case I will do my best to explain everything I know.

jherb · March 10, 2014, 05:14

Unfortunately, all pictures in this thread are gone.

Your problem might be a common problem with the SIMPLE algorithm, RANS and long wavelength/low frequency oscillations. These are not handled correctly by SIMPLE because the eddies are not correctly treated by the RANS averaging. So you get some kind of vortex shedding in stationary simulations. See also this thread:
http://www.cfd-online.com/Forums/ope...m-results.html

With a pimple simulation you can check, if your boundary conditions works. Perhaps you can also use a converged solution of (one time step of) a pimple simulation as initial guess for your simple simulation.

Quote:

Originally Posted by UPengineer

Hello jherb,

Thank you for your help and quick reply.

I did as you suggested by extending the area below the cylinder.
Unfortunately, the residuals did not converge any better.

The mesh is now:

Which output residuals seen:

For your other questions/suggestions:

I am using buoyantBoussinesqSimpleFoam as the solver, so yes I am using the Boussinesq approximation.

Would transient simulation be too overkill since I do not have any time dependent variables? There is nothing cooling, heating, or filling a space. It is just the heat from an isothermally heated cylinder.

I hope this may help pinpoint what my problem is. If there is any other questions for my case setup that may help OpenFOAM run this case I will do my best to explain everything I know.

UPengineer · March 11, 2014, 13:53

Thank you for your suggestion jherb, though I am not really sure how to implement it.

After continuously changing boundary conditions I saw that the hole where the cylinder is causes a lot of disruption for the incoming velocity of the bottom face. This made me try again with a new fully meshed model.

However after troubleshooting it became clear that I have a problem with retaining heat even with the patch on bottom and top since it seems that the heat circulates as if it were a cyclic boundary.

Is there a way to make the top face just expel all heat and velocity? To make either OpenFOAM think it is an infinitely tall wall or just allow everything to pass through? An open boundary condition?

jherb · March 11, 2014, 19:15

Could you post your fvSolution and fvSchemes files? And also the files in 0 time directory? Have you tried zeroGradient for the outlet (top) as boundary condition for U?

UPengineer · March 11, 2014, 20:38

Quote:

Originally Posted by jherb

Could you post your fvSolution and fvSchemes files? And also the files in 0 time directory? Have you tried zeroGradient for the outlet (top) as boundary condition for U?

Thank you very much for your continued help and time with my project. The case that matches most close to reality is slightly different than the original case.

I found that the P_(rgh) that openFOAM uses is in units m^2/s^2 so it does not actually include any density. So bottom condition is a total pressure of (35m)(9.81m/s^2) = 343.35 m^2/s^2

The case is as follows:

Mesh:

FvSolutions:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
|       o          |                                                          |
|    o     o       | HELYX-OS                                                  |
|   o   O   o      | Version: v2.0.0                                           |
|    o     o       | Web:     http://www.engys.com                            |
|       o          |                                                          |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location system;
    object fvSolution;
}

    SIMPLE
    {
        nNonOrthogonalCorrectors 0;
        pressureImplicitPorousity false;
        pRefCell 0;
        pRefValue 0;
    }

    solvers
    {
        p
        {
            solver GAMG;
            agglomerator faceAreaPair;
            mergeLevels 1;
            cacheAgglomeration true;
            nCellsInCoarsestLevel 200;
            tolerance 1e-7;
            relTol 0.01;
            smoother GaussSeidel;
            nPreSweeps 0;
            nPostSweeps 2;
            nFinestSweeps 2;
            minIter 1;
        }

        p_rgh
        {
            solver GAMG;
            agglomerator faceAreaPair;
            mergeLevels 1;
            cacheAgglomeration true;
            nCellsInCoarsestLevel 200;
            tolerance 1e-7;
            relTol 0.01;
            smoother GaussSeidel;
            nPreSweeps 0;
            nPostSweeps 2;
            nFinestSweeps 2;
            minIter 1;
        }

        f
        {
            solver GAMG;
            agglomerator faceAreaPair;
            mergeLevels 1;
            cacheAgglomeration true;
            nCellsInCoarsestLevel 200;
            tolerance 1e-7;
            relTol 0.01;
            smoother GaussSeidel;
            nPreSweeps 0;
            nPostSweeps 2;
            nFinestSweeps 2;
            minIter 1;
        }

        U
        {
            solver smoothSolver;
            smoother GaussSeidel;
            tolerance 1e-6;
            relTol 0.1;
            minIter 1;
        }

        k
        {
            solver smoothSolver;
            smoother GaussSeidel;
            tolerance 1e-6;
            relTol 0.1;
            minIter 1;
        }

        kl
        {
            solver smoothSolver;
            smoother GaussSeidel;
            tolerance 1e-6;
            relTol 0.1;
            minIter 1;
        }

        kt
        {
            solver smoothSolver;
            smoother GaussSeidel;
            tolerance 1e-6;
            relTol 0.1;
            minIter 1;
        }

        q
        {
            solver smoothSolver;
            smoother GaussSeidel;
            tolerance 1e-6;
            relTol 0.1;
            minIter 1;
        }

        zeta
        {
            solver smoothSolver;
            smoother GaussSeidel;
            tolerance 1e-6;
            relTol 0.1;
            minIter 1;
        }

        epsilon
        {
            solver smoothSolver;
            smoother GaussSeidel;
            tolerance 1e-6;
            relTol 0.1;
            minIter 1;
        }

        R
        {
            solver smoothSolver;
            smoother GaussSeidel;
            tolerance 1e-6;
            relTol 0.1;
            minIter 1;
        }

        nuTilda
        {
            solver smoothSolver;
            smoother GaussSeidel;
            tolerance 1e-6;
            relTol 0.1;
            minIter 1;
        }

        omega
        {
            solver smoothSolver;
            smoother GaussSeidel;
            tolerance 1e-6;
            relTol 0.1;
            minIter 1;
        }

        h
        {
            solver smoothSolver;
            smoother GaussSeidel;
            tolerance 1e-6;
            relTol 0.1;
            minIter 1;
        }

        T
        {
            solver smoothSolver;
            smoother GaussSeidel;
            tolerance 1e-6;
            relTol 0.1;
            minIter 1;
        }

        v2
        {
            solver smoothSolver;
            smoother GaussSeidel;
            tolerance 1e-6;
            relTol 0.1;
            minIter 1;
        }

        rho
        {
            solver PCG;
            preconditioner DIC;
            tolerance 0;
            relTol 0;
            minIter 1;
        }

        rhoFinal
        {
            solver PCG;
            preconditioner DIC;
            tolerance 0;
            relTol 0;
            minIter 1;
        }

        e
        {
            solver PBiCG;
            preconditioner DILU;
            tolerance 1e-06;
            relTol 0.1;
            minIter 1;
        }

    }

    relaxationFactors
    {
        fields
        {
            p_rgh 0.5;
            p 0.3;
            rho 0.05;
        }

        equations
        {
            U 0.7;
            h 0.5;
            k 0.7;
            kl 0.7;
            kt 0.7;
            q 0.7;
            zeta 0.7;
            epsilon 0.7;
            R 0.7;
            nuTilda 0.7;
            omega 0.7;
            T 0.7;
            v2 0.7;
            f 0.7;
        }

    }

FvSchemes:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
|       o          |                                                          |
|    o     o       | HELYX-OS                                                  |
|   o   O   o      | Version: v2.0.0                                           |
|    o     o       | Web:     http://www.engys.com                            |
|       o          |                                                          |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location system;
    object fvSchemes;
}

    ddtSchemes
    {
        default steadyState;
    }

    gradSchemes
    {
        default Gauss linear;
        grad(nuTilda) cellLimited Gauss linear 1;
        grad(k) cellLimited Gauss linear 1;
        grad(kl) cellLimited Gauss linear 1;
        grad(omega) cellLimited Gauss linear 1;
        grad(epsilon) cellLimited Gauss linear 1;
        grad(q) cellLimited Gauss linear 1;
        grad(zeta) cellLimited Gauss linear 1;
        grad(v2) cellLimited Gauss linear 1;
        grad(f) cellLimited Gauss linear 1;
        grad(sqrt(kt)) cellLimited Gauss linear 1;
        grad(kt) cellLimited Gauss linear 1;
        grad(sqrt(kl)) cellLimited Gauss linear 1;
    }

    divSchemes
    {
        default Gauss linear;
        div(phi,U) bounded Gauss linearUpwindV grad(U);
        div(phi,k) bounded Gauss linearUpwind grad(k);
        div(phi,epsilon) bounded Gauss linearUpwind grad(epsilon);
        div(phi,zeta) bounded Gauss linearUpwind grad(zeta);
        div(phi,q) bounded Gauss linearUpwind grad(q);
        div(phi,omega) bounded Gauss linearUpwind grad(omega);
        div(phi,nuTilda) bounded Gauss linearUpwind grad(nuTilda);
        div(phi,T) bounded Gauss linearUpwind grad(T);
        div(phi,kl) Gauss limitedLinear 1;
        div(phi,kt) Gauss limitedLinear 1;
        div(phi,R) Gauss upwind;
        div(R) Gauss linear;
        div((nuEff*dev(grad(U).T()))) Gauss linear;
        div(phi,v2) bounded Gauss linearUpwind grad(v2);
        div(phi,f) bounded Gauss linearUpwind grad(f);
    }

    fluxRequired
    {
        default no;
        p_rgh ;
        p ;
    }

    interpolationSchemes
    {
        default linear;
        interpolate(HbyA) linear;
    }

    laplacianSchemes
    {
        default Gauss linear limited 0.333;
    }

    snGradSchemes
    {
        default limited 0.333;
    }

Residuals:

Outputs at time=1000s

Velocity:

Pressure:

Temperature:

This was the best I could do after spending many many hours changing boundary conditions, mesh, relaxation factor for pressure, and etc...

Some questions I still have for this simulation:

Why is the pressure negative at the top, when it was specified to have a total pressure of 0 m^2/s^2?
Why do I have such a high velocity? It makes me believe that this is somewhat forced convection instead of natural convection.
Any other advice to make it more realistic to natural convection?

Again thank you jherb for your extended time and help!

jherb · March 13, 2014, 05:17

I would suggest that you start over with a very simple geometry, e.g. the hot room tutorial for buoyantBoussinesqSimpleFoam. Then change the boundary conditions of this example so you have open patches at the bottom and the top. Due to the "hot spot" at the bottom (created by the setFields command in the Allrun script) you should also get natural circulation.

(Be carefull that the T boundary condition in 0/T is overwritten by the 0/T.org file)

If this works you can start with your mesh.

P.S. Your pictures/graphs are gone again.

March 9, 2014, 03:08	2D isothermal cylinder not converging	#1
UPengineer New Member Kramer Join Date: Aug 2013 Posts: 11 Rep Power: 12	Hello Everyone, I am using OpenFOAM version 2.2.2 with HelyxOS 2.0 on Ubuntu 13.10. I am looking at a 2D natural convective flow from an isothermally heated cylinder. The mesh I have may be seen below: There is a 1.5 m by 1.5 m box surrounding the cylinder with 4X refinement. The boundary conditions are: Code: Wall -->Fixed, no slip - zero gradient temperature for: front face, back face, left face and right face. Wall --> 350 K isothermal temperature for the cylinder Patch --> Velocity: Pressure Inlet velocity Pressure: Total pressure = 0 m^2/s^2 For the top face Patch --> Velocity: Pressure Inlet Outlet Pressure: Zero gradient For the Bottom face The Solution Modeling is: Code: SIMPLE FOAM algorithm 2nd Order Bounded Linear Upwind for velocity and temperatureSteady Incompressible Laminar With thermal energy on and buoyancy on (gravity set to -9.81 m/s^2) The residuals are: As you can see, the residuals are anything but convergent which makes me think the boundary conditions may be a little messed up. An example of the temperature post processing of the case is: which appears to be correct, however because of the residuals I cannot trust this data. Any advice on how to get my case to converge? Any help is appreciated.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Incorrect Drag and Drag Coefficient for flow over a cylinder	ozzythewise	Main CFD Forum	8	June 13, 2012 06:24
convection from an isothermal horizontal cylinder	engahmed	Main CFD Forum	0	April 19, 2011 12:27
Not converging flow around cylinder with Re=1000000. Help please!	Wotild	ANSYS	5	December 7, 2010 21:07
Not converging flow around cylinder with Re=1000000. Help please!	Wotild	ANSYS	0	December 3, 2010 09:50
[blockMesh] Specifying boundary faces failes in blockMesh	blaise	OpenFOAM Meshing & Mesh Conversion	0	May 10, 2010 03:56

March 9, 2014, 11:31		#2
jherb Senior Member Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 21	What solver are you using? With or without Boussinesq approximation? Have you tried the same setup with a transient simulation (pimple)? Perhaps also the area below the cylinder is to short.

March 11, 2014, 13:53		#5
UPengineer New Member Kramer Join Date: Aug 2013 Posts: 11 Rep Power: 12	Thank you for your suggestion jherb, though I am not really sure how to implement it. After continuously changing boundary conditions I saw that the hole where the cylinder is causes a lot of disruption for the incoming velocity of the bottom face. This made me try again with a new fully meshed model. However after troubleshooting it became clear that I have a problem with retaining heat even with the patch on bottom and top since it seems that the heat circulates as if it were a cyclic boundary. Is there a way to make the top face just expel all heat and velocity? To make either OpenFOAM think it is an infinitely tall wall or just allow everything to pass through? An open boundary condition?

March 11, 2014, 19:15		#6
jherb Senior Member Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 21	Could you post your fvSolution and fvSchemes files? And also the files in 0 time directory? Have you tried zeroGradient for the outlet (top) as boundary condition for U?

March 13, 2014, 05:17		#8
jherb Senior Member Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 21	I would suggest that you start over with a very simple geometry, e.g. the hot room tutorial for buoyantBoussinesqSimpleFoam. Then change the boundary conditions of this example so you have open patches at the bottom and the top. Due to the "hot spot" at the bottom (created by the setFields command in the Allrun script) you should also get natural circulation. (Be carefull that the T boundary condition in 0/T is overwritten by the 0/T.org file) If this works you can start with your mesh. P.S. Your pictures/graphs are gone again.