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:
Code:
/*--------------------------------*- 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:
Code:
liftDir (0.606265086952 0 5.76822699314);
dragDir (5.76822699314 0 -0.606265086952);
Or, when omitting the magnitude of velocity:
Code:
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)):
Code:
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)):
Code:
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):
Code:
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
