#111 Bug in T1 sequence analysis with sideslip, affecting mainly CL and XCP

[lct]

In PanelAnalysis.cpp

Function void PanelAnalysis::computePlane(double Alpha, double QInf, int qrhs)

cosa, sina and wind vectors WindNormal, WinDirection, WindSide are calculated using Alpha before it is updated.

As such, the wind vector is always corresponding to the first AoA point in the sequence of the analysis.

Fix: Move the Alpha update in the begining of the function.

BEFORE:

void PanelAnalysis::computePlane(double Alpha, double QInf, int qrhs)
{

int pos;
double *Mu, *Sigma;
double cosa, sina;
double Lift, IDrag, VDrag ,XCP, YCP, ZCP, WingVDrag;
Vector3d WindNormal, WindDirection, WindSide;
Vector3d Force;
QString OutString;
//   Define wind (stability) axis
cosa = cos(Alpha*PI/180.0);
sina = sin(Alpha*PI/180.0);

WindNormal.set(  -sina, 0.0, cosa);
WindDirection.set(cosa, 0.0, sina);
WindSide = WindNormal * WindDirection;

Mu     = m_Mu    + qrhs*m_MatSize;
Sigma  = m_Sigma + qrhs*m_MatSize;

m_QInf      = QInf;

if(m_pWPolar->bTilted() || m_pWPolar->isBetaPolar() || fabs(m_pWPolar->Beta())>PRECISION)
    Alpha = m_OpAlpha;
else
    m_OpAlpha = Alpha;

   ....

AFTER:

void PanelAnalysis::computePlane(double Alpha, double QInf, int qrhs)
{
int pos;
double Mu, Sigma;
double cosa, sina;
double Lift, IDrag, VDrag ,XCP, YCP, ZCP, WingVDrag;
Vector3d WindNormal, WindDirection, WindSide;
Vector3d Force;
QString OutString;

/* ------------------- * /

if (m_pWPolar->bTilted() || m_pWPolar->isBetaPolar() || fabs(m_pWPolar->Beta())>PRECISION)
Alpha = m_OpAlpha;
else
m_OpAlpha = Alpha;

/* ------------------- * /

//   Define wind (stability) axis
cosa = cos(Alpha*PI/180.0);
sina = sin(Alpha*PI/180.0);

...

[André]

I see that you have had a deep look at the code.
Note that on line 1452 in the computeAeroCoefs() method, computePlane() is called with V0+q*Velta for first argument. In the case of a type 1 or 2 polar, this is the current aoa. You can double check that by outputing alpha in the computeAeroCoefs() method:

    //   Define wind (stability) axis
    qDebug()<<Alpha;
    cosa = cos(Alpha*PI/180.0);
    sina = sin(Alpha*PI/180.0);

I agree that this is complex and confusing. It caused indirectly and error which was corrected with commit 987.
I'll double check everything again before I release the next version.

[lct]

Hej Andre,

I am currently hacking away at a 971 revision, where I made some changes to be able to multiple planes times multiple analysis from a batch script, based on the existing XML features. A lot of the code was already in place, thank you for that ! And also to analyse a fin alone (without a main wing).

Maybe I should update to the newer code, but I have not done so yet.

This is the stack frame of the bug - T1 with sideslip analysis - as I saw it:

 PanelAnalysis::loop()
 **********************
 627: if(m_pWPolar->polarType()<XFLR5::FIXEDAOAPOLAR)
    {
        if(m_pWPolar->bTilted() || fabs(m_pWPolar->Beta())>PRECISION) 
           ->  return unitLoop();

 PanelAnalysis::unitLoop()   // because of sideslip is not alphaloop()
 *************************
2368: case XFLR5::FIXEDSPEEDPOLAR:
            case XFLR5::FIXEDLIFTPOLAR:
       ->    computeAeroCoefs(m_vMin, m_vDelta, 1);   // e.g (-5,1,1)

 PanelAnalysis::computeAeroCoefs( double V0, double VDelta, int nrhs )
 ********************************
 ..
 1448: for (q=0; q<nrhs ; q++) 
 // nrhs is hardcoded 1 from previous call, so q=0 always
         ...
      ->   computePlane(V0+q*VDelta, m_3DQInf[q], q); 
                                    // V0 is always m_vMin from previous call

   PanelAnalysis::computePlane(double Alpha, double QInf, int qrhs)
   ********************************
   ...
    cosa = cos(Alpha*PI/180.0); // Alpha is always m_vMin at this point !
    sina = sin(Alpha*PI/180.0);

    // And then only 

    if (m_pWPolar->bTilted() || m_pWPolar->isBetaPolar() || fabs(m_pWPolar->Beta())>PRECISION)
        Alpha = m_OpAlpha;  // m_opAlpha is actually now correct iterating on the sequence, but the wind has been calculated with m_vMin 

best regards
Luminita

[lct]

So I quess claim is that, refering to:

called with V0+q*Velta for first argument. In the case of a type 1 or 2 polar, this is the current aoa.

For Type 1 Polar with sideslip V0+q VDelta is not the currect AoA, is always the first AoA in the sequence as per unitLoop() code.

[André]

Not a bug