[nusivares]

Hello all,

I am fairly new with OpenFOAM and as my first study I am trying to compare OpenFOAM results to the ones obtained to FLUENT.

Software versions: OpenFOAM: 2.4.0, FLUENT: 16.1.1

The case I am attempting is a sinusoidal channel - the one which is commonly used in heat exchangers. To eliminate the entrance effects - the flow should be modelled as periodic.

The case I ran was at Re=100. Within fluent I set up a Inlet/Outlet as interfaces and connected them using periodic boundary conditions. In terms of periodic conditions I set them us as default in Fluent - as mass flow rate.
fluent_velocity.png fluent_pressure.png


As of OpenFOAM I imported the fluent mesh, created cyclic boundary conditions and used the fvOptions file to set the momentum source terms using mean velocity to make the things flow. The result I got this was:

Openfoam_velocity.png Openfoam_pressure.png

As you can see there is a clear difference between the cases. Would anyone have any suggestions why it is happening? Do you have any solutions how to solve this issue? is there a way to put the mass flow rate as a source term to make the correlation easier?

The case files are below(The file was too big, thus, had to stick it in dropbox):

https://www.dropbox.com/s/8cfncge72e...yclic.zip?dl=0

Thanks for your help in advance! Documentation on such flows is not the best.


The code for calculating both mass flow and velocity from Re (matlab):

clear
clc
format long

W=2.23e-3; % Width of the channel/fin, m
H=1; % Height of the channel/fin, m for fluent it is 1 metre for 2D cases
L=50e-3; % Lenght of the channel/fin, m

rho = 1.225; % Density, kg/m^3
mu=1.7894e-5; % Dyncamic Viscosity kg/ms

nu=mu/rho; % Kinematic viscosity

Re=100; % Reynolds number as usual in the papers


U_calc=(Re*nu)/(2*W) % If input is velocity they make dh=2W
m_dot=rho*U_calc*W*H % Mass flow

[RodriguezFatz]

For iteration 3000 the yplus for the wall is
Patch 1 named WALL y+ : min: 0.28997018732504709071 max: 0.43480356182426332934 average: 0.34463384097139043716
So it is a low-Re mesh. Your boundary conditions are for a high-Re mesh. This might explain the difference in velocity plots.
Pressure of the incompressible solver is real pressure / density. So you need to multiply the openFoam pressure values by your operating density to compare it to fluent.

[nusivares]

Hi,

Thanks for reply. I was thinking for a while what you meant by that sentence.. Because I had kEpsilon set as a RAS model with turbulence switched off.

Then I actually tried sticking in the laminar RAS model and it made the whole lot of difference.

I was mildly surprised by that because I though that when I turned the turbulence off the solver came back to laminar equations, but apparently it did not.

The result is:

openfoam_vel_fixed.png Openfoam_p_fixed.png



As you can see the results are fairly identical - the small differences are cause by discretisation schemes, which are a bit lower accuracy currently set on OpenFOAM. I tried doing that for fluent and it gave me identical result.


Thanks for giving me a clue with the pressure! helped a lot.

[RodriguezFatz]

I didn't even recognize that you had turbulence switched off by that setting. I never used that. But anyway now it works
Normally what you did should work. I don't understand why it wasn't switched off.

[RodriguezFatz]

Well, maybe with your setting it just doesn't solve for k and epsilon / omega equation but uses the given (\0 directory)-entries as fixed settings to calculate the turbulent viscosity...
Edit: I guess this is true. If you set internal k to zero in your \0 directory and restart you will get the correct result.

[nusivares]

Thanks a lot for your explanation!


In terms of the laminar solver, I managed to get to around Re=400(u=1.310075958634575) where it decided not to converge anymore.

Do you have any ideas how to make it more robust? Fluent managed to cope it by using coupled instead of simple.

[RodriguezFatz]

Hi, what do you mean by "more robust"?

[nusivares]

Hi,

by that I meant make it to work for broader range of Re. When I increased Re number for the case to Re = 400 (Ubar ~ 1.31). It caused fluctuations in the solution, whilst the FLUENT with coupled algorithm managed to converge.

Do you know any tips how to avoid it?

[RodriguezFatz]

"fluctuations" in the residuals or where?

[nusivares]

Hi,

I will explain it better. The residuals:

Re400-laminar.png

Screenshots of contours every thousand of iterations:

8000:

RE_400_iter_8000.png

9000:

RE_400_iter_9000.png


10000:

RE_400_iter_10000.png

[RodriguezFatz]

Well that doesn't converge.
Can you post your fvSolution, fvSchemes and the last iterations of the log file?

[nusivares]

Sure, hopefully it does not too messy for you:

end of the log:

Time = 9295

smoothSolver:  Solving for Ux, Initial residual = 0.0030142261772487153292, Final residual = 0.00018070334041139326825, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0053867800424603354617, Final residual = 0.00033189147557002824729, No Iterations 2
Pressure gradient source: uncorrected Ubar = 1.3100603632674117094, pressure gradient = 545.82799332159220285
GAMG:  Solving for p, Initial residual = 0.06511037145398994308, Final residual = 0.0041452058269561622983, No Iterations 2
time step continuity errors : sum local = 0.51530691506081338638, global = 9.6585870312325598997e-06, cumulative = 0.089760748431006662229
Pressure gradient source: uncorrected Ubar = 1.3100815117520456798, pressure gradient = 542.91507642249086985
ExecutionTime = 120.35999999999999943 s  ClockTime = 120 s

Time = 9296

smoothSolver:  Solving for Ux, Initial residual = 0.0029965368759614871043, Final residual = 0.0001785926573439436415, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0054422825121599408688, Final residual = 0.00033430076309385209918, No Iterations 2
Pressure gradient source: uncorrected Ubar = 1.3100597716058490771, pressure gradient = 545.14540749593311375
GAMG:  Solving for p, Initial residual = 0.06451587043113798936, Final residual = 0.003530805916910612248, No Iterations 2
time step continuity errors : sum local = 0.4415437208139421088, global = 9.6587834832260265847e-06, cumulative = 0.089770407214489894576
Pressure gradient source: uncorrected Ubar = 1.3100822717510278714, pressure gradient = 542.04522315703286495
ExecutionTime = 120.37999999999999545 s  ClockTime = 120 s

Time = 9297

smoothSolver:  Solving for Ux, Initial residual = 0.003021439132626162552, Final residual = 0.00017980854255233262042, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0053910275692351083271, Final residual = 0.00033061411398559788655, No Iterations 2
Pressure gradient source: uncorrected Ubar = 1.3100620085433221096, pressure gradient = 543.96790143860755506
GAMG:  Solving for p, Initial residual = 0.06527820975546982929, Final residual = 0.0042271249612424203687, No Iterations 2
time step continuity errors : sum local = 0.55618282819879694134, global = 9.6590992031169466126e-06, cumulative = 0.089780066313693007718
Pressure gradient source: uncorrected Ubar = 1.3100792798772626924, pressure gradient = 541.58747122782585848
ExecutionTime = 120.39000000000000057 s  ClockTime = 120 s

Time = 9298

smoothSolver:  Solving for Ux, Initial residual = 0.0029947009240669690902, Final residual = 0.00017771144110895249909, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0054590914567207446742, Final residual = 0.00033386662296646844135, No Iterations 2
Pressure gradient source: uncorrected Ubar = 1.3100607748466328584, pressure gradient = 543.68067534449471623
GAMG:  Solving for p, Initial residual = 0.063088913975031510328, Final residual = 0.0039161000397144908602, No Iterations 2
time step continuity errors : sum local = 0.48745883414103380327, global = 9.6591822411424853908e-06, cumulative = 0.089789725495934155641
Pressure gradient source: uncorrected Ubar = 1.3100830731378683147, pressure gradient = 540.60668121425521804
ExecutionTime = 120.40000000000000568 s  ClockTime = 120 s

Time = 9299

smoothSolver:  Solving for Ux, Initial residual = 0.0030144707107860281223, Final residual = 0.00017813898387061542533, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0054223025449059549311, Final residual = 0.00033097957930049420189, No Iterations 2
Pressure gradient source: uncorrected Ubar = 1.3100626629309388704, pressure gradient = 542.4399464015061767
GAMG:  Solving for p, Initial residual = 0.064360624329235463503, Final residual = 0.003574292542555871071, No Iterations 2
time step continuity errors : sum local = 0.47073056129697354866, global = 9.6593609382747772503e-06, cumulative = 0.089799384856872430105
Pressure gradient source: uncorrected Ubar = 1.3100810319135225424, pressure gradient = 539.90715707213769292
ExecutionTime = 120.42000000000000171 s  ClockTime = 120 s

Time = 9300

smoothSolver:  Solving for Ux, Initial residual = 0.0030060085023729749346, Final residual = 0.0001775962916378959059, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0054377854809080318141, Final residual = 0.00033143736345043990155, No Iterations 2
Pressure gradient source: uncorrected Ubar = 1.3100633286004443878, pressure gradient = 541.64885297837497546
GAMG:  Solving for p, Initial residual = 0.062612587898651206331, Final residual = 0.0040836100267523317911, No Iterations 2
time step continuity errors : sum local = 0.50924408601294590682, global = 9.659438363008983355e-06, cumulative = 0.08980904429523543786
Pressure gradient source: uncorrected Ubar = 1.3100836275253004981, pressure gradient = 538.84960841913800778
ExecutionTime = 120.43000000000000682 s  ClockTime = 120 s

Time = 9301

smoothSolver:  Solving for Ux, Initial residual = 0.0029993729454928868237, Final residual = 0.00017639142097321101497, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0054679932764793289998, Final residual = 0.00033218215220257886738, No Iterations 2
Pressure gradient source: uncorrected Ubar = 1.3100633575023452604, pressure gradient = 540.58746802936445874
GAMG:  Solving for p, Initial residual = 0.062273976018851491532, Final residual = 0.0030614791476954047143, No Iterations 2
time step continuity errors : sum local = 0.40569798671617252062, global = 9.6595026991233798341e-06, cumulative = 0.08981870379793456538
Pressure gradient source: uncorrected Ubar = 1.3100788330487438405, pressure gradient = 538.45318941518178235
ExecutionTime = 120.45000000000000284 s  ClockTime = 120 s

Time = 9302

smoothSolver:  Solving for Ux, Initial residual = 0.0030126219283441268808, Final residual = 0.00017685873923883887755, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0054212643616981508882, Final residual = 0.0003293912110210666733, No Iterations 2
Pressure gradient source: uncorrected Ubar = 1.3100658910909372956, pressure gradient = 539.84167044233481647
GAMG:  Solving for p, Initial residual = 0.062754410783897957016, Final residual = 0.0041134422295080593987, No Iterations 2
time step continuity errors : sum local = 0.51270005829692355537, global = 9.6596181861437968927e-06, cumulative = 0.089828363416120715623
Pressure gradient source: uncorrected Ubar = 1.3100831747255314852, pressure gradient = 537.45797094441934405
ExecutionTime = 120.45999999999999375 s  ClockTime = 120 s

Time = 9303

smoothSolver:  Solving for Ux, Initial residual = 0.002990507618087509447, Final residual = 0.00017533903922666086968, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0054850112576477674103, Final residual = 0.00033195624508595321855, No Iterations 2
Pressure gradient source: uncorrected Ubar = 1.3100648542525064055, pressure gradient = 538.98941952984546333
GAMG:  Solving for p, Initial residual = 0.061633026457751864668, Final residual = 0.0029219085350524580391, No Iterations 2
time step continuity errors : sum local = 0.38715982644265090062, global = 9.6595322213938960524e-06, cumulative = 0.08983802294834210278
Pressure gradient source: uncorrected Ubar = 1.3100749450238227123, pressure gradient = 537.59776195152642231
ExecutionTime = 120.46999999999999886 s  ClockTime = 120 s

fvSchemes:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
    grad(U)         cellLimited Gauss linear 1;
}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss upwind;
//    div(phi,k)      bounded Gauss upwind;
//    div(phi,epsilon) bounded Gauss upwind;
    div(phi,R)      bounded Gauss upwind;
    div(R)          Gauss linear;
//   div(phi,nuTilda) bounded Gauss upwind;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

fluxRequired
{
    default         no;
    p;
    Phi;
}


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

fvSolution:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

   solvers
{
    Phi
    {
        $p;
}
potentialFlow
{
    nNonOrthogonalCorrectors 10;
}

cache
{
    grad(U);
}

   p
   {
        solver          GAMG;
        tolerance       1e-20;
        relTol          0.1;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        cacheAgglomeration on;
        agglomerator    faceAreaPair;
        nCellsInCoarsestLevel 10;
        mergeLevels     1;
    }

    "(U|k|epsilon|R|nuTilda)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-20;
       relTol          0.1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    pRefCell        1001;
    pRefValue       0;
    residualControl
    {
        p               1e-20;
        U               1e-20;
        "(k|epsilon|omega)" 1e-20;
    }
}

relaxationFactors
{
    fields
    {
        p               0.5;
    }
    equations
    {
        U               0.7;
        k               0.7;
        epsilon         0.7;
        R               0.7;
       nuTilda         0.7;
   }
}


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

[RodriguezFatz]

For me it looks like you already have pretty safe settings.
You are using cyclic boundary conditions? I tryed to get them running for different very similar cases, namely a straight pipe and a channel flow. Both of them did not converge properly no matter what numerical settings I used. I eventually used the "mapped" boundary condition instead of the cyclic, which does pretty much the same for this kind of set up. That worked without any problems. I suspect that there is some issue with the solver and cyclic b.c. - or I am simply incompetent, which could be taken into account.

[RodriguezFatz]

Here is the thread:
http://www.cfd-online.com/Forums/ope...tml#post521820

[nusivares]

Ok, I will attempt to modify everything into the mapped conditions and I will post the results.

I did try using the mapped condiitons initially without too much luck,but again, my knowledge of OpenFOAM increases significantly every day.

[nusivares]

Ok.. So I did try it and I think cyclic is better than mapped in my case at least when everything is periodic.. I will explain why and give the bits and pieces of code.

So, for mapped case. I used Re=200 case which works well everywhere. Boundary file:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       polyBoundaryMesh;
    location    "constant/polyMesh";
    object      boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

4
(
    OUTLET
    {
        type            patch;
        nFaces          59;
        startFace       21937;
    }
    WALL
    {
        type            wall;
        inGroups        1(wall);
        nFaces          376;
        startFace       21996;
    }
    INLET
    {
	type            mappedPatch;
        inGroups        1(mappedPatch);
        nFaces          59;
        startFace       22372;
        sampleMode      nearestPatchFace;//nearestCell;
        sampleRegion    region0;
        samplePatch     OUTLET;
        offsetMode      uniform;
        offset          (0.00628 0 0);
        //type            patch;
        //nFaces          59;
        //startFace       22372;
    }
    frontAndBackPlanes
    {
        type            empty;
        inGroups        1(empty);
        nFaces          22184;
        startFace       22431;
    }
)

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

0/U:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    location    "0";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0.655037979317287 0 0);

boundaryField
{
     INLET
    {
        type                mapped;
        value               uniform (0.655037979317287 0 0);
        interpolationScheme cell;
        setAverage          true;
        average             (0.655037979317287 0 0);
	//type            fixedValue;
        //value           uniform (1.31 0 0);
    }

    OUTLET
    {
        type            fixedValue;
	value		$internalField;
    }

    WALL
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    frontAndBackPlanes
    {
        type            empty;
    }
}


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

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform 0;

boundaryField
{
    INLET
    {
	type                mapped;
        value               $internalField;
        interpolationScheme cell;
        setAverage          true;
        average             0;        
	//type            fixedValue;
	//value		$internalField;
    }

    OUTLET
    {
        type            fixedValue;
        value           $internalField;
    }

    WALL
    {
        type            zeroGradient;
    }
    frontAndBackPlanes
    {
        type            empty;
    }
}

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

And it gives a result:
RE200_residuals.png RE200_vel.png RE_200_p.png


Applying cyclic condition as before, meanwhile, gives:

RE_200_vel_cyclic.png RE_200_p_cyclic.png


Using FLUENT is on another message as I ran out of maximum attachment number.

[nusivares]

So, using Fluent gives:

RE_200_u_fluent.png RE_200_p_fluent.png


As you can see fluent agrees fairly well with the cyclic condition and for mapped condition.. not so much. There are some settings as you may know as well whilst using the mapped boundary condition complicated to make it truly cyclic.


As of Re = 400 both fluent and OpenFOAM struggle..

So, I have tried several turbulence models as well to see what is happening. I know RE number is low, but again, the geometry of the problem is such that flow would be disturbed to maximum.

So, I have got some Results which are in the following messages.

[nusivares]

So, using now FLUENT with COUPLED (Simple misbehaved and did not converge s in OpenFOAM) to predict flow at Re=400 and laminar:

fluent_laminar_u.png fluent_lam_p.png

Openfoam on laminar:

open_lam_u.png open_lam_p.png

[nusivares]

So, the next step was to use k-epsilon (because why not). Again, FLUENT:


fluent_k-e-u.png fluent_k-e_p.png


Openfoam

open_k-e-u.png open_k-e-p.png


Residuals:

open-k-e-res.png

Looks like a fair match between the two

[nusivares]

And finally k-omegaSST (because why not again). Fluent:

fluent_k-o-u.png fluent-k-o-p.png


OpenFOAM:

open_k-o-u.png open-k-o-p.png

and residuals:

k-o-res.png

Again, not horrific match, however, as can be seen the result did not converge that well.

I have noticed this before, that somehow the kOmegaSST is much more sensitive to the initial conditions (my guess) which consequently causes all of this.

It did work reasonably for the simple inlet/outlet case, suggesting that cyclic conditions might cause that as well to some extent.

To avoid this, I was planing to make a quick code to estimate everything. Unless your opinion of this is different. Or maybe you know the app for estimating quantities within OpenFOAM already exists?