[MBttR]

Hi all,

I'm doing some force calculations on an aircraft under angle of attack and up to this point I have always filled in U*sin(alpha) and U*cos(alpha) into the vectors for liftDir and dragDir. I have now tried a calculation with just sin(alpha) and cos(alpha) as I saw someone else do it and my forces are completely different for both simpleFoam and pisoFoam. In case it matters, I am running OpenFOAM 4.0 (but as far as I understood, I compared the code in github with 2.4 and 3.0 and it seems the same).

The data:
Flying with a velocity of 5.8m/s at an angle of attack of -6 degrees (negative). This results in the following 0/U file:

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

flowVelocity    (5.76822699314 0 -0.606265086952);

pressure        0;

turbulentKE     0.24;

turbulentOmega  1.78;

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

internalField   uniform (5.76822699314 0 -0.606265086952);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (5.76822699314 0 -0.606265086952);
    }
    outlet
    {
        type            inletOutlet;
        inletValue      uniform (0 0 0);
        value           uniform (5.76822699314 0 -0.606265086952);
    }
    lowerWall
    {
        type            inletOutlet;
        inletValue      uniform (0 0 0);
        value           uniform (5.76822699314 0 -0.606265086952);
    }
    md43000
    {
        type            noSlip;
    }
    upperWall
    {
        type            fixedValue;
        value           uniform (5.76822699314 0 -0.606265086952);
    }
    "(front|back)"
    {
        type            symmetryPlane;
    }
    "proc.*"
    {
        type            processor;
    }
}


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

As a result I pass the following vectors to calculate the force coefficients. FYI, my x-axis is pointing to the rear of the aircraft (hence against direction of travel), z-axis points up:

liftDir (0.606265086952 0 5.76822699314);
dragDir (5.76822699314 0 -0.606265086952);

Or, when omitting the magnitude of velocity:

liftDir (0.104528463267 0 0.994521895368);
dragDir (0.9945218953680 -0.104528463267);

So the size is different, but the direction is the same. I was starting to think I'm stupid after seeing the results but Wolfram Alpha assured I'm right.

Here are the last entries in the log.simpleFoam files. pisoFoam shows similarly different results but for the sake of the length of this post I'm only posting simpleFoam.

For Case #1 (with U*cos(alpha) and U*sin(alpha)):

smoothSolver:  Solving for Ux, Initial residual = 0.000159418, Final residual = 1.53285e-05, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0032299, Final residual = 0.00028852, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 0.000345581, Final residual = 3.07293e-05, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.000389341, Final residual = 1.69226e-05, No Iterations 2
time step continuity errors : sum local = 2.10409e-06, global = -1.20322e-07, cumulative = 0.00206658
smoothSolver:  Solving for nuTilda, Initial residual = 5.87586e-05, Final residual = 5.13931e-06, No Iterations 2
ExecutionTime = 760.82 s  ClockTime = 797 s

forceCoeffs forces write:
    Cm    = -0.00646708
    Cd    = 0.909421
    Cl    = -0.174108
    Cl(f) = -0.093521
    Cl(r) = -0.0805869

For Case #2 (with cos(alpha) and sin(alpha)):

smoothSolver:  Solving for Ux, Initial residual = 0.000159418, Final residual = 1.53285e-05, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0032299, Final residual = 0.00028852, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 0.000345581, Final residual = 3.07293e-05, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.000389341, Final residual = 1.69226e-05, No Iterations 2
time step continuity errors : sum local = 2.10409e-06, global = -1.20322e-07, cumulative = 0.00206658
smoothSolver:  Solving for nuTilda, Initial residual = 5.87586e-05, Final residual = 5.13931e-06, No Iterations 2
ExecutionTime = 2367 s  ClockTime = 2449 s

forceCoeffs forces write:
    Cm    = -0.00646708
    Cd    = 0.156797
    Cl    = -0.0300186
    Cl(f) = -0.0214764
    Cl(r) = -0.00854222

Note that the massive difference in runtime is probably due to me running two virtual machines at the same time while doing the latter case (I thought I was doing that as well before but can't be sure).

As you can see the residuals are exactly the same. The forces are completely different, by orders of magnitude.

I wanted to find out what caused this difference (and which approach is correct), and found the following in the code (but can't quite understand it):

coeffs[0] = (totForce & liftDir_)/(Aref_*pDyn);
coeffs[1] = (totForce & dragDir_)/(Aref_*pDyn);

scalar Cl = sum(coeffs[0]);
scalar Cd = sum(coeffs[1]);

What does (totForce & liftDir_) mean? is it some sort of multiplication? Anyway, the fact that liftDir_ is in this formula worries me and makes me think not just the direction but the magnitude is used as well, without some sort of normalization.

What is the right way? Should I just use cos(alpha) and sin(alpha)? Why "code-wise" does this matter? The direction stays the same. Is this something that should be addressed to the OF-devs?

I hope I have given enough information and someone can clear this up and help me out

Thanks,
Cheers,

Bruno

[MBttR]

Nobody?

A simple yes or no answer to whether I should only use sine and cosine would already suffice.

Cheers,
Bruno

[alexeym]

Hi,

1. For OpenFOAM vectors & is dot product.
2. In general liftDir and dragDir should be normalized (i.e. ||liftDir|| = 1 and ||dragDir|| = 1).

[MBttR]

Originally Posted by alexeym

Hi,

1. For OpenFOAM vectors & is dot product.
2. In general liftDir and dragDir should be normalized (i.e. ||liftDir|| = 1 and ||dragDir|| = 1).

Thanks a lot Alexey! Good to know & is dot product