[shereez234]

Hello Every one;

I have been trying to compare some CFD results (OpenFoam vs. Fluent)
for a flat plate case which I am guessing most of you are familiar with.

Right! So the case is incompressible and run with Mach = 0.1 ( V =34.6 m/s)).
I get the same Pressure and Velocity Field for both Sofwares. Also, I get the same results for Kinetic Energy, Turbulent Viscocity and Omega.

However, my lift and drag coefficients are off. I am expecting a CD = 0.01 and Cl = -0.0094872 (fluent results) and in OpenFoam I get CD =0.000534504 and CL -0.000464995.

With a lot of repetition I came to notice the pressure forces is correct and equal to fluent's results. However the viscous force is 3 times lower. I checked the wallGradU.exe utility and the gradient at the wall matches Fluent's results. However when I checked the Wall shear stress I have negative values in the same data range [ for example - 100 to -3] and in Fluent I have [+100 to +3].

Solver is SimpleFoam or PisoFoam, Turbulence is on, Model is SSTvorticity(modified SST).

fvScheme

ddtSchemes
{
    default         steadyState;
}

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

}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwind grad(U);
    div(phi,k)      bounded Gauss limitedLinear 1;
    div(phi,omega)  bounded Gauss limitedLinear 1;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
    snGrad(U) corrected;
}

fluxRequired
{
    default         no;
    p               ;
    Phi		    ;
}

fvSolution

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-08;
        relTol          0.001;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
	maxIter	        1000;
    }

    pFinal
    {
        solver 		GAMG;
	tolerance       1e-10;
	relTol		0.001;
        smoother        GaussSeidel;
    }

    Phi
    {
        $p;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         1;
        tolerance       1e-10;
        relTol          0.1;
	maxIter	         1000;
    }
    k
    {
        solver           smoothSolver;
        smoother         GaussSeidel;
        tolerance        1e-8;
        relTol           0.1;
        nSweeps          1;
	maxIter	         1000;
    }

    omega
    {
        solver           smoothSolver;
        smoother         GaussSeidel;
        tolerance        1e-10;
        relTol           0.1;
        nSweeps          1;
	maxIter	         1000;
    }
}

SIMPLE
{
    nCorrectors     0;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;

    residualControl
    {
        p               1e-6;
        U               1e-6;
    }
}
potentialFlow
{
    nNonOrthogonalCorrectors 10;
    pRefCell        0;
    pRefValue       0;
    PhiRefCell      0;
    PhiRefPoint     0;
    PhiRefValue     0;
    Phi		    0;
}
relaxationFactors
{
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.5;
	omega		0.5;
	k		0.5;
    }
}

forces in ControlDict

functions
{
 forces
  {
	type forceCoeffs;
	functionObjectLibs ( "libforces.so" );
	outputControl timeStep;
	outputInterval 10;
	patches
	(
	CRVS
	);
	pName p;
	UName U;
	rhoName rhoInf;
	rhoInf 1.2;
	log true;
	CofR (0 0 0 );
	liftDir (0 1 0);
	dragDir (1 0 0);
	pitchAxis (0 0 1);
	magUInf 34.7;
	lRef 2;
	Aref 2;
  }

Residuals at the End of Simulation

smoothSolver:  Solving for Ux, Initial residual = 1.82304e-008, Final residual = 5.40526e-010, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 6.87114e-008, Final residual = 3.0058e-009, No Iterations 1
GAMG:  Solving for p, Initial residual = 2.69686e-006, Final residual = 8.28443e-009, No Iterations 29
time step continuity errors : sum local = 4.57401e-009, global = -7.21234e-010, cumulative = -6.15044e-006
smoothSolver:  Solving for omega, Initial residual = 4.97458e-006, Final residual = 3.4698e-007, No Iterations 1
smoothSolver:  Solving for k, Initial residual = 1.53917e-005, Final residual = 9.17097e-007, No Iterations 2
bounding k, min: 1e-018 max: 0.240867 average: 0.0297818
ExecutionTime = 116.615 s  ClockTime = 117 s

These were the ultimate settings that helped me improve the pressure field.

Since, this is a flat plat plate there max skewness of the mesh is = 0. So mesh quality issues are eliminated from this discussion I guess.

So is there any thing I am doing wrong? How can I improve my drag and Lift forces for this simple case?

Many thanks

Shereez

[shereez234]

Dear Fellow Users of OpenFoam;


I quote Allesenadro:

http://www.cfd-online.com/Forums/ope...-cylinder.html


If the solver outputs forces then your drag will be:

Cd = (Force per unit span)/(1/2 * rho V^2 * chord )=(Drag/(z-depth))/(1/2 * rho * v^2 * S).

I mean your S, which I assume is the diameter of your cylinder.

Alternatively, you can compute CD as:

CD = (Force)/(1/2 * rho V^2 * wettedArea)=(Drag)/(1/2 * rho * v^2 * (z-depth)*pi*S).

By S I mean your S.

I hope this helps,

Alessandro

More Proof:
Forces.C

if Turbulence Model is on:

        scalar pRef = pRef_/rho(p);

        forAllConstIter(labelHashSet, patchSet_, iter)
        {
            label patchI = iter.key();

            vectorField Md
            (
                mesh.C().boundaryField()[patchI] - coordSys_.origin()
            );
            // pressure force
            vectorField fN
            (
                rho(p)*Sfb[patchI]*(p.boundaryField()[patchI] - pRef)
            );
             //viscous force
                //inner product of Surface Normal of boundary cells and devReff even though here it is devRhoReff
            vectorField fT(Sfb[patchI] & devRhoReffb[patchI]);

            vectorField fP(Md.size(), vector::zero);

            applyBins(Md, fN, fT, fP, mesh.C().boundaryField()[patchI]);
        }
    }

I changed in my polymesh points file Z direction depth from 0.02 to 1 and -0.02 to -1

and here is my new Drag Coefficient in OpenFoam:

Cd = 0.019851


Many thanks to every one. Case closed.

Best regards

Shereez