[WhiteW]

Hi to everybody.
I'm running some simulation on an airplane using both Fluent and OpenFoam 2.3.0 in order to do a comparison. I'm using the compressible solver rhoSimplecFoam and I have found a resolution strategy for OF in order to obtain a convergence:
- 100 iterations with energy equation deactivated (through fvOptions)
- first order schema until the residual are 1e-4
- second-like order schema
This strategy grants me a good accordance in the global pressure and velocity fields. The pressure field on the airplane surfaces also seems to be correct.
However the force coefficients, the drag and lift forces are very different in OpenFoam and I don't understand the reason. I attached an image of the force report for Fluent and OpenFoam analysis(I converted the OF results in fluent format using foamMeshToFluent and foamDataToFluent commands).
I report the setting used in controlDict, fvSolution and fvSchemes.
The boundary I'm using are:
p - inlet: totalPressure - outlet: fixedValue
U - inlet_ pressureInletVelocity - outlet: pressureInletOutletVelocity
T - inlet: totalTemperature - outlet: inletOutletTotalTemperature
Turbulence model - the k-omega SST.
Has someone some advice in order to fit the forces output?
Thanks in advance,
[WhiteW

ControlDict

application     rhoSimpleFoam;
startFrom       latestTime;
startTime       0;
stopAt          endTime;
endTime         6000;
deltaT          1;
writeControl    timeStep;
writeInterval   200;
purgeWrite      0;
writeFormat     ascii;
writePrecision  9;
writeCompression off;   
timeFormat      general;
timePrecision   9;
graphFormat     raw;
runTimeModifiable true;
libs
(
 "libforces.so"
 );
functions
{
	forces
	{
		type forces;
		functionObjectLibs ( "libforces.so" );
		outputControl timeStep;
		outputInterval 1;
		patches
			(
			 wall_fuselage
			 wall_nose
			 wall_sponson
			 wall_junction
			 wall_tail
			 wall_wing_nacelle
			 wall_nacelle
			 wall_pianetto_orizzontale1
			 wall_raccordo_pianetti
			 wall_pianetto_verticale
			);
		pName p;
		UName U;
		rhoName rhoInf;
		log true;
		rhoInf 1.16;
		CofR (0.98460 0 0.35195);            
		pitchAxis (0 -1 0);
	}
}

fvOptions

source1
{
    type            temperatureLimitsConstraint;
    selectionMode   all;
    active          true;

        temperatureLimitsConstraintCoeffs
        {
            Tmin     294.18;
            Tmax     294.19;
        }
}

fvSchemes 1st order:

ddtSchemes
{
	default             steadyState;
}

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

divSchemes
{
	default             bounded Gauss upwind;
	div(phi,U)          bounded Gauss upwind;
	div((muEff*dev2(T(grad(U)))))      Gauss linear;
}

laplacianSchemes
{
	default Gauss linear corrected;
}

interpolationSchemes
{
	default         linear;
}

snGradSchemes
{
	default         corrected;
}

fluxRequired
{
	default         no;
	p;
}

fvSchemes 1st-2nd order

ddtSchemes
{
	default             steadyState;
}

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

divSchemes
{
	default             bounded Gauss linearUpwind cellLimited Gauss linear 1;
	div(phi,U)          bounded Gauss linearUpwindV cellMDLimited Gauss linear 1;
	div((muEff*dev2(T(grad(U)))))      Gauss linear;
}

laplacianSchemes
{
	default Gauss linear corrected;
}

interpolationSchemes
{
	default         linear;
}

snGradSchemes
{
	default         corrected;
}

fluxRequired
{
	default         no;
	p;
}

fvSolution

solvers
{
p
    {
        solver          GAMG;
        tolerance       1e-09;
        relTol          0.1;
        smoother        GaussSeidel;
        cacheAgglomeration off;
        nCellsInCoarsestLevel 800;
        agglomerator    faceAreaPair;
        mergeLevels     1;
        maxIter         40;
    }
rho
    {
        solver          GAMG;
        tolerance       1e-09;
        relTol          0.1;
        smoother        GaussSeidel;
        cacheAgglomeration false;
        nCellsInCoarsestLevel 800;
        agglomerator    faceAreaPair;
        mergeLevels     1;
        maxIter         20;
    }
U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps 	2;
        tolerance       1e-09;
        relTol          0.1;
        maxIter         40;
    }

h
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-09;
        relTol          0.1;
        maxIter         20;
    }
k
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps 	2;
        tolerance       1e-09;
        relTol          0.1;
        maxIter         20;
    }
omega
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps 	2;
        tolerance       1e-09;
        relTol          0.1;
        maxIter         20;
    }
e  
    {
        solver          GAMG;
        tolerance       1e-08;
        relTol          0.1;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 800;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 4;
    pMin            pMin [1 -1 -2 0 0 0 0] 90000;
    pMax            pMax [1 -1 -2 0 0 0 0] 105000;
    rhoMin          rhoMin [1 -3 0 0 0 0 0] 1.12; //then changed to 1.08
    rhoMax          rhoMax [1 -3 0 0 0 0 0] 1.20;  //then changed to 1.24

    residualControl //tutorial
    {
        p               1e-4;
        U               1e-4;
        e               1e-4;
        "(k|epsilon|omega)" 1e-3;
    }
}

relaxationFactors
{
        p               0.3;
        rho             0.1;
        e               0.7;
        U               0.7;
        h               0.7;
        k               0.7;
        omega           0.7;
}

[WhiteW]

I noticed the output values of forces strongly depends from the initialization of turbulence parameters.
I used k=0 and omega=2, for the internalField.
Now, using different values (k=0.027 and omega=2.05) the force output in the first iterations differ of about an order of magnitude(but then the simulation diverges).

I also run a checkMesh and this is the log:

 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           4949919
    faces:            36290908
    internal faces:   35830274
    cells:            16487628
    faces per cell:   4.37426063
    boundary patches: 20
    point zones:      0
    face zones:       1
    cell zones:       3

Overall number of cells of each type:
    hexahedra:     0
    prisms:        6170670
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    10316958
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
    Patch               Faces    Points   Surface topology
    wall_fuselage       18346    9448     ok (non-closed singly connected)
    wall_pilone_lato_aereo2626     1390     ok (non-closed singly connected)
    wall_nose           21680    11037    ok (non-closed singly connected)
    wall_sponson        21869    11040    ok (non-closed singly connected)
    wall_junction       29688    15182    ok (non-closed singly connected)
    wall_tail           34684    17834    ok (non-closed singly connected)
    wall_wing_nacelle   82898    41788    ok (non-closed singly connected)
    wall_nacelle        43241    21880    ok (non-closed singly connected)
    wall_pianetto_orizzontale140740    20449    ok (non-closed singly connected)
    wall_raccordo_pianetti7184     3805     ok (non-closed singly connected)
    wall_pianetto_verticale33317    17231    ok (non-closed singly connected)
    wall_pil_fisso      6542     3432     ok (non-closed singly connected)
    inflow              146      90       ok (non-closed singly connected)
    outflow             148      91       ok (non-closed singly connected)
    symm_parete         1571     851      ok (non-closed singly connected)
    symm_floor          615      364      ok (non-closed singly connected)
    symm_ceiling_interno1550     1155     ok (non-closed singly connected)
    symm_ceiling_esterno2254     1224     ok (non-closed singly connected)
    symm_esterno        7978     4224     ok (non-closed singly connected)
    symm_interno        103557   77227    ok (non-closed singly connected)

Checking geometry...
    Overall domain bounding box (-6.95 -3.5 -2.273955) (14.65 3.33069446e-06 2.726045)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (-2.05138739e-18 -7.24646001e-16 1.31380308e-16) OK.
 ***High aspect ratio cells found, Max aspect ratio: 4925.81694, number of cells 1586
  <<Writing 1586 cells with high aspect ratio to set highAspectRatioCells
    Minimum face area = 4.70543557e-10. Maximum face area = 0.32125922.  Face area magnitudes OK.
    Min volume = 8.94286652e-14. Max volume = 0.0550260094.  Total volume = 345.427077.  Cell volumes OK.
    Mesh non-orthogonality Max: 78.9019986 average: 12.3286271
   *Number of severely non-orthogonal (> 70 degrees) faces: 69.
    Non-orthogonality check OK.
  <<Writing 69 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
 ***Max skewness = 14.8103056, 4889 highly skew faces detected which may impair the quality of the results
  <<Writing 4889 skew faces to set skewFaces
    Coupled point location match (average 0) OK.

Failed 2 mesh checks.

End

Any idea about the causes of this error in the forces estimation?
Thanks,

WhiteW