[ajcav2]

Hello,

I'm working on solving a pipe flow problem in OpenFOAM, and I'm having trouble getting the results to agree with the solution I obtain using other CFD programs such as ANSYS and AcuSim. The problem is a 3D internal pipe flow problem through a 90 degree bend, and the flow is turbulent and incompressible. A picture of the geometry is below.



I am using the kwSST turbulence model after some research suggested it would be the most accurate. I am not using wall functions here, and the mesh has a y+ value of ~1 at the wall.

To verify this case, I am trying to measure total pressure drop through the tube. AcuSim and ANSYS agree that the total pressure drop should be ~1050 Pa. OpenFOAM gives a total pressure drop of only 335 Pa which seems extremely low.

Here are the fvSchemes I am using:

ddtSchemes
{
    default    steadyState; 
}


gradSchemes
{
     default        Gauss linear;
     grad(p)        Gauss linear;
     grad(U)        Gauss linear;

}

divSchemes
{
     default         none;
     div(phi,U)         bounded Gauss linearUpwindV grad(U); 
     div(phi,k)          bounded Gauss upwind;
     div(phi,omega)         bounded Gauss upwind;
     div(phi,nuTilda)          bounded Gauss upwind;

     div((nuEff*dev(T(grad(U))))) Gauss linear;
     div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
    laplacian((1|A(U)),p)         Gauss linear corrected;
    laplacian(nuEff,U)         Gauss linear corrected;
    laplacian(DkEff,k)         Gauss linear corrected;
    laplacian(DepsilonEff,epsilon)         Gauss linear corrected;
    laplacian(DomegaEff,omega)         Gauss linear corrected;
    laplacian(DREff,R)         Gauss linear corrected;
    laplacian(DnuTildaEff,nuTilda)         Gauss linear corrected;
    laplacian(1,p)  Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
    interpolate(U)         linear;

}

snGradSchemes
{
    default         corrected; // corrected

}

wallDist
{
    method meshWave;
}

fluxRequired
{
    default         no;
    p         yes;

}

And here is my fvSolution file:

solvers
{
    p
    {
         solver   GAMG;
         smoother       GaussSeidel;
         nPreSweeps       0;
         nPostSweeps       2;
         cacheAgglomeration       true;
         agglomerator       faceAreaPair;
         nCellsInCoarsestLevel       10;
         mergeLevels       1;
         tolerance       1e-05;
         relTol          0.01;

    }
    
    Phi
    {
        $p;
    }

    U
    {
         solver   smoothSolver;
         smoother       GaussSeidel;
         nSweeps       2;
         tolerance       1e-08;
         relTol          0.1;
    }

    k
    {
         solver   smoothSolver;
         smoother       GaussSeidel;
         nSweeps       2;
         tolerance       1e-08;
         relTol          0.1;
    }


    omega
    {
         solver   smoothSolver;
         smoother       GaussSeidel;
         nSweeps       2;
         tolerance       1e-08;
         relTol          0.1;
    }


}
SIMPLE
{
    nNonOrthogonalCorrectors             2;
    nCorrectors             0;
    pRefCell        0;
    pRefValue       0;
    residualControl
    {
        p        1.0e-5;
        U        1.0e-6;
        k        1.0e-6;
        epsilon  1.0e-3;
        omega    1.0e-6;
        nuTilda  1.0e-3;
        T        1.0e-5;
    }

}

potentialFlow
{
    nNonOrthogonalCorrectors  10;
}

relaxationFactors
{
    p        0.3; // 0.3
    U        0.5; // 0.7
    k        0.5; // 0.7
    epsilon  0.7;
    omega    0.5; // 0.7
    nuTilda  0.7;
    T        0.7;
}

cache
{
    grad(U);
}

Finally, here is a patch summary:

wall    : r_tube_air_tube_tube_air_2
    scalar              nut             zeroGradient
    scalar              k               fixedValue
    scalar              p               zeroGradient
    scalar              omega           generic
    vector              U               fixedValue

patch   : r_inlet_flow_2d_tube_air_3
    scalar              nut             calculated
    scalar              k               turbulentIntensityKineticEnergyInlet
    scalar              p               zeroGradient
    scalar              omega           generic
    vector              U               flowRateInletVelocity

patch   : r_outlet_pres_2d_tube_air_4
    scalar              nut             zeroGradient
    scalar              k               turbulentIntensityKineticEnergyInlet
    scalar              p               fixedValue
    scalar              omega           generic
    vector              U               inletOutlet

The inlet is defined with a mass flow rate of 0.3 kg/s, and the density is constant at 1.20069 kg/m^3. The outlet has a zero gauge pressure boundary condition. I am using the most recent version of bluecape's windows ported openfoam (4.x).

Please let me know if there is any additional information I can provide. This is my first post here, but I've read through the post guidelines and tried to provide all relevant information and be as thorough as possible.

Thank you.

[jeytsav]

Dear ajcav2

Have you managed to figure out what the problem was?

Maybe you can try posting the boundary conditions for each quantity as well. I am talking about the files into the \0 folder.

Or it would have been even better if you could upload the whole case in Dropbox or elsewhere so I can have a better look.

Another thing you can do is to initially run the case as laminar to check if everything works well and then implement the turbulence model of your preference.

Kind regards!

[ajcav2]

jeytsav,

Thanks for your response. I still haven't been able to figure out what the problem could be. I took your advice and ran the case as laminar and didn't notice anything out of the ordinary. I've also tried using different versions of openfoam to run the simulation without any luck. I've uploaded the case for OpenFOAM version 4.x to Dropbox. Please let me know if you'd like me to upload the case for a different version of OpenFOAM.

Thanks so much for your help.

[jeytsav]

Dear ajcav2

Sorry for the late reply I am a bit busy with exams etc this period. Have you got any news?

I tried running several times your case, by modifying some parameters. I couldn't get a match of ~1050 Pa total pressure drop.

Surprisingly, I ran the case exactly like you gave it to me (turbulent case) and I got a total pressure drop of ~1680 Pa.

I calculated the total pressure drop as follows

Total pressure at outlet = static + dynamic = 0 + 0.5*rho*v_outlet^2

Total pressure at inlet = static + dynamic = p_inlet*rho + 0.5*rho*v_inlet^2


openFoam calculates the pressure in [m^2/s^2] so you have to multiply times density to get the static pressure at the inlet.

I also get

p_inlet = 2293
v_inlet = 49.5 (as set)
v_outlet = 65


These results are after 4200 iteration, it appears that the pressure drop is growing as I let the simulation run.. maybe it needs more time..

I tried some different settings, but I couldn't make it work and get a ~1050 pressure drop..

Maybe you can also have a look at this pdf describing an external aerodynamics case, could make some sense. Some settings are different than yours. Maybe you can have a look.

How are you calculating the pressure drop?!

I also think you do not have to specify a velocity value at the inlet since you are already specifying a flow rate. Guess it works both ways but just in case..

Please correct me If I am wrong!

Attached find the residuals.
myResiduals.jpg

P.S. I ll also try some other settings as well, let me know if you have already solved it :P



Kind regards

[ajcav2]

jeytsav,

Sorry for the late reply. I've been trying everything to get the case to work. I've tried looking over the pdf you sent and using those settings, but I get the same incorrect answer. I tried refining my mesh, and it doesn't change anything either .

To calculate pressure drop, I am calculating total pressure on each surface, then integrate the total pressure over the entire surface, then divide by the area of the surface. This should yield the total pressure at whatever patch was being used.

I'm beginning to suspect that it's the shape of the mesh that isn't ideal for this type of case. Maybe a hexahedral mesh will yield better results. I'm going to see if I can find a T-Grid license anywhere to mesh the case with hexahedral cells. I'll update as soon as I get results.

[jherb]

Have you seen this thread and the boundary conditions for turbulence properties for the SST model for small Y+ values:
One more time: wall functions - SST

[ajcav2]

Just in case anyone sees this thread, I did eventually manage to find the error, but not where I expected. I started using EnSight to post process my data and, as it turns out, the way data is read into EnSight from OpenFOAM's format results in incorrect pressure calculations which threw off my total pressure calculations. I have a script now to calculate the total pressure drop through the tube directly from OpenFOAM data, and the values are nearly exactly as expected. Thanks for your help, everyone.