[preis]

Hello dear Community,

this is my first thread so first of all I want to thank everybody contributing to this forum! You guys helped me a lot!

So this is going to be a longer one but I still hope someone is willing to sacrifice some time. Now to my problem:

I am trying to simulate a pipe flow through a expanding and afterwards narrowing geometry, a diffuser-pipe-nozzle configuration so to say. Both times with a 45deg angle and a inlet Re-number of 50000 (calculated with bulk velocity) and developed inlet conditions. I used different turbulence models which all showed similar wrong results so I will refer to the kOmegaSST version below.
It is an incompressible case (water-nu=1e-6) and I solved for a steady state solution with simpleFoam. The simulation was 2D as a wedge. CheckMesh was satisfied with my mesh quality and the average yPlus value has been approx. 0.5.

In my first configuration I chose a rather high area Ratio of the connected pipes of 16 (diam. 0.025m to 0.1m). The solution converged but some results were questionable. First of all I expected a static pressure rise after the diffuser along the axis which was hardly achieved (k-epselon even showed a drop) but only a huge pressure drop in the nozzle (see plots attached).
The second issue was the eddy viscosity ratio nut/nu which seems to be unreasonable high - around 1200 (see attached picture). I expected to see some value well below 500.

To verify my case setup I simulated the ERCOFTAC Case 75 (sudden Geometry expansion). Which is basically a diffuser with a 90deg half angle. The Re number was with 15600 lower but the flow type seemed to be similar. The results were very good compared to experimental data. And the eddy viscosity was approx. 180 (as expected by the cfd-online turbulence tool). In that magnitude i expected nut/nu to be for my initial case
After that I simulated my case again with a lower Re number but the results didn't get better.

My last try was to reduce the area-ratio of my case and adapted it to the one from the ERCOFTAC-case (much lower with 2.56). Inlet Re number was still 50000. And I got very much better Results. There was a higher pressure rise along the Diffuser (although much smaller area ratio) and the linear pressure drop in the big pipe indicates developed pipe flow (pressure was nearly constant along the big pipe in the first configuration). The eddy viscosity has a maximum of about 400 with what I could live but the contour plot may also suggest that there is much more convection of k and omega involved whereas in the first place there seems to be a strange eddyViscosity bubble. I attached a pressure and eddyViscosity plot as well.

So this is basically why I assume that the first case is wrong whereas the lower area ratio leads to reasonable results with the same setup. I would like to find out why?! I would appreciate any help and suggestions! I can hardly believe that the area ratio is the factor that leads to such wrong results.

I also tried out a lot different solver configurations and schemes, even one simulation with an unstructured mesh. But nothing gave me better results. I will upload my setup files and the complete case as well. I used second order linearUpwind schemes.

Best wishes
Patrick

...ahh yes and OpenFOAM 3.0.1 is my version.

Initial-files:

k:

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

internalField   uniform 0.015;

boundaryField
{
    inlet
    {
	type		fixedValue;
        value           nonuniform List<scalar>



45
(
0.00933555
0.00983949
0.0107606
0.0119197
0.0131983
0.0145152
0.0158171
0.0170712
0.0182588
0.0193703
0.0204024
0.0213552
0.0222311
0.0230335
0.0237665
0.0244342
0.0250407
0.0255896
0.0260845
0.0265282
0.0269234
0.0272721
0.0275758
0.0278355
0.0280514
0.0282229
0.0283486
0.0284258
0.0284504
0.0284162
0.0283142
0.0281319
0.027851
0.0274454
0.0268765
0.026085
0.0249641
0.0233483
0.0209806
0.0175078
0.0125525
0.00641524
0.0018269
0.000262463
0.000213873
);
    }
    outlet
    {
	type		inletOutlet;
        value           $internalField;
        inletValue      uniform 0;
    }
    walls
    {
	type		fixedValue;
	value		uniform 1e-12;
    }
    "(front|back)"
    {
        type            wedge;
    }
    axis
    {
        type            empty;
    }
}

U:

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

internalField   uniform (2 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           nonuniform List<vector>


45
(
(2.3977 -4.38574e-06 4.10496e-17)
(2.38947 -3.19018e-05 9.86728e-17)
(2.37359 -4.41901e-05 1.87704e-16)
(2.35264 -5.38954e-05 2.69189e-16)
(2.32815 -5.72397e-05 3.28321e-16)
(2.30123 -5.73149e-05 3.65412e-16)
(2.27266 -5.50426e-05 4.26e-16)
(2.24301 -5.16693e-05 4.9447e-16)
(2.21269 -4.76964e-05 5.23782e-16)
(2.18198 -4.35582e-05 5.33742e-16)
(2.15107 -3.94613e-05 5.85707e-16)
(2.12009 -3.55532e-05 6.15477e-16)
(2.08912 -3.19225e-05 6.27487e-16)
(2.0582 -2.86069e-05 6.48472e-16)
(2.02735 -2.56122e-05 6.59231e-16)
(1.99655 -2.2926e-05 6.67904e-16)
(1.9658 -2.05268e-05 6.7409e-16)
(1.93506 -1.83888e-05 6.78534e-16)
(1.90429 -1.64855e-05 6.77416e-16)
(1.87344 -1.47913e-05 6.77412e-16)
(1.84244 -1.32824e-05 6.76851e-16)
(1.81123 -1.19369e-05 6.73991e-16)
(1.77972 -1.07356e-05 6.69179e-16)
(1.74782 -9.66107e-06 6.64312e-16)
(1.71542 -8.69819e-06 6.58646e-16)
(1.6824 -7.83365e-06 6.51936e-16)
(1.64861 -7.05578e-06 6.43942e-16)
(1.61387 -6.35445e-06 6.35235e-16)
(1.57798 -5.72102e-06 6.24993e-16)
(1.54068 -5.14755e-06 6.1312e-16)
(1.50167 -4.6266e-06 6.00062e-16)
(1.46055 -4.15171e-06 5.85673e-16)
(1.41684 -3.72524e-06 5.7007e-16)
(1.36988 -3.35115e-06 5.52832e-16)
(1.31881 -2.91031e-06 5.33624e-16)
(1.26245 -2.03916e-06 5.11128e-16)
(1.19891 -2.16759e-06 4.8554e-16)
(1.12542 -2.34127e-06 4.566e-16)
(1.0384 -2.59479e-06 4.21496e-16)
(0.932893 -1.90032e-06 3.78671e-16)
(0.802023 2.63587e-07 3.25737e-16)
(0.638514 -5.14399e-07 2.58909e-16)
(0.449021 -3.81458e-07 1.81614e-16)
(0.25914 -1.19122e-07 1.04628e-16)
(0.0837018 -7.53384e-09 3.37192e-17)
);
    }
    outlet
    {
        type            inletOutlet;
        value           $internalField;
        inletValue      uniform (0 0 0);
    }
    walls
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    "(front|back)"
    {
        type            wedge;
    }
    axis
    {
        type            empty;
    }
}

omega:

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

internalField   uniform 500;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           nonuniform List<scalar>



45
(
97.309
101.697
110.073
121.108
133.957
148.115
163.337
179.571
196.894
215.471
235.524
257.313
281.134
307.316
336.229
368.294
403.998
443.904
488.674
539.093
596.098
660.819
734.631
819.218
916.67
1029.61
1161.34
1316.12
1499.49
1718.76
1983.76
2307.95
2710.23
3217.77
3870.79
4731.42
5899.09
7532.41
9892.61
13455.8
19243.6
30173.7
58660.6
174787
1.06099e+06
);
    }
    outlet
    {
        type            zeroGradient;
    }
    walls
    {
        type            omegaWallFunction;
        value           uniform 500;
    }
    "(front|back)"
    {
        type            wedge;
    }
    axis
    {
        type            empty;
    }
}

p:

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

internalField   uniform 0; 

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }
    walls
    {
        type            zeroGradient;
    }
    "(front|back)"
    {
        type            wedge;
    }
    axis
    {
        type            empty;
    }
}

nut:

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

internalField   uniform 0;

boundaryField
{
    inlet
    {
	type	calculated;
    }
    outlet
    {
	type 	calculated;
    }
    walls
    {
	type 	calculated;
	value	uniform 0;
    }
    "(front|back)"
    {
        type            wedge;
    }
    axis
    {
        type            empty;
    }
}

fvSchemes:

ddtSchemes
{
    default         steadyState;
}

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

divSchemes
{
    default	    none;
    
    /* second order

    div(phi,U)      bounded Gauss linearUpwind grad(U);
    div(phi,k)      bounded Gauss linearUpwind default;
    div(phi,omega)  bounded Gauss linearUpwind default;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
	
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method meshWave;
}

fvSolution:

solvers
{
    p
    {
        solver                GAMG;
        tolerance             1e-06;
        relTol                1e-06;
	smoother	      GaussSeidel;
	nPreSweeps	      0;
	nPostSweeps	      2;
        cacheAgglomeration    on;
	agglomerator	      faceAreaPair;
	nCellsInCoarsestLevel 10;
	mergeLevels           1;
    }
	
    U
    {
        solver          PBiCG;
    	preconditioner	DILU;
        tolerance       1e-08;
        relTol          0.0;
    }
    
    "(epsilon|R|omega|k|nut)"
    {
	solver		smoothSolver;
	smoother	GaussSeidel;
	tolerance	1e-06;
	relTol		0.0;
	nSweeps		1;
    }	
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    //pRefValue	0;
    //pRefCell 1001;

    residualControl
    {
        p               0.5e-7;
        U               0.5e-7;
        "(k|epsilon|omega)" 1e-7;
    }
}

relaxationFactors
{
    fields
    {
    	p	0.3;
    }
    equations
    {
	"(k|omega|epsilon|U)"	0.7;
    }
}

[preis]

And here are the case files for the simulation with the "wrong"results