Add some airfoil related functions

[rainman110]

From [email protected] on November 25, 2011 13:29:16

It would be useful to have a rapid function to calculate the maximum thickness of an airfoil. Such a function should require wing, segment and relative spanwidth coordinate as input and should compute the relative airfoil thickness of the thickest point (relative to the chordlength) and the corresponding ksi-coordinate. I am currently doing this iteratively using point functions (see attached code snippet), but this is quite time consuming und I could imagine, there is a faster way when having direct access to the geometry.

If this function is included, it is also necessary to have a "getChordlength" function available. This one should use the total difference (x, y and z) between leading and trailing edge (and not only the x-coordinate).

Best regards,

Carsten

Attachment: getMaximumThicknessKsi.cpp

Original issue: http://code.google.com/p/tigl/issues/detail?id=52

[rainman110]

From [email protected] on November 29, 2011 02:30:21

Status: Accepted
Owner: [email protected]
Labels: -Type-Defect Type-Enhancement Milestone-Releasex.x

[rainman110]

Anforderungen:

• Funktion für die Ausgabe der maximalen relativen Dicke eines Profils (Dicke)
• Funktion für die Ausgabe der Ksi-Koordinate der maximalen relativen Dicke eines Profils (Dickenrücklage)
• Parameter der beiden Funktionen: Flügel, Segment, Segement-Eta
• Ausgabe: Parameter „Dicke“ und „Dickenrücklage“ im Verhältnis zur örtlichen Flügeltiefe (Sehnenlänge in iso-eta)
• In dem Zusammenhang wäre auch eine Funktion zur Ausgabe der benutzten Sehnenlänge sinnvoll (wurzel(dx²+dy²+dz²) mit den Differenzen zwischen VK-punkt und HK-Punkt (Auf ISO_Eta, Mitte zwischen HK oben und HK unten)

[rainman110]

Possible implementation from @CLiersch:

/**
**************************************************************************************************************************************************************
* This const function returns the relative chordlength position of the maximum thickness on the given wing segment and relative spanwidth position.
* The function takes three parameters.
* \param  Wing                      an unsigned int    containing the number of the wing
* \param  Segment                   an unsigned int    containing the number of the segment
* \param  RelativeSpanwidthPosition a           double containing the eta position on the segment
* \return A double containing the relative chordlength position
**************************************************************************************************************************************************************
**/
double cCPACSInputBaseClass::getMaximumThicknessKsi (unsigned int Wing,
                                                     unsigned int Segment,
                                                     double       RelativeSpanwidthPosition) const
{

    bool LowerKsiFound = false;
    bool UpperKsiFound = false;
    bool ExtremumFound = false;

    double LowerKsi                = 0.;
    double LowerGradient           = 0.;
    double UpperKsi                = 1.;
    double UpperGradient           = 0.;
    double NewKsi                  = 0.;
    double NewKsiLowThickness      = 0.;
    double NewKsiHighThickness     = 0.;
    double NewGradient             = 0.;

    double MaximumThicknessKsi     = 0.;

    double ChordLength             = getProjectedChordLength (Wing, Segment, RelativeSpanwidthPosition);

    const double GRADIENT_DISTANCE = 0.00001;                                           // Distance for gradient calculation (relative to Chordlength)


    // Assure relative spanwidth between 0. and 1.
    // (Deviations greater than EPS will already have been found in getProjectedChordLength ())
    // ----------------------------------------------------------------------------------------
    RelativeSpanwidthPosition = limitToInterval (RelativeSpanwidthPosition, 0., 1.);


    // Check, if Chordlength is > 0
    if (!isGeoZero (ChordLength))
    {
        // Find two borders around the maximum thickness (or an extremum itself)
        do
        {
            // Select way to find next NewKsi
            if (!(LowerKsiFound && UpperKsiFound))
            {
                // Compute new ksi in the middle between upper and lower border
                NewKsi = (LowerKsi + UpperKsi) / 2.;
            }
            else
            {
                // Interpolate new ksi (linear "0"-search between lower and upper ksi)
                NewKsi = (LowerKsi - UpperKsi) * LowerGradient / (UpperGradient - LowerGradient) + LowerKsi;
            }

            // Get gradient around the new ksi
            NewKsiLowThickness  = getThickness (Wing, Segment, RelativeSpanwidthPosition, NewKsi - GRADIENT_DISTANCE / 2.);
            NewKsiHighThickness = getThickness (Wing, Segment, RelativeSpanwidthPosition, NewKsi + GRADIENT_DISTANCE / 2.);
            NewGradient         = (NewKsiHighThickness - NewKsiLowThickness) / (GRADIENT_DISTANCE * ChordLength);

            // Check, if it is already an extremum, or a new lower or upper border
            if (isZero (NewGradient))
            {
                // Extremum
                ExtremumFound = true;
            }
            else if (NewGradient > 0.)
            {
                // If step size becomes smaller than GRADIENT_DISTANCE, the result is considered to be close enough to the maximum
                if ((NewKsi - LowerKsi) < GRADIENT_DISTANCE)
                    ExtremumFound = true;

                // New lower ksi
                LowerKsi               = NewKsi;
                LowerGradient          = NewGradient;
                LowerKsiFound          = true;
            }
            else
            {
                // If step size becomes smaller than GRADIENT_DISTANCE, the result is considered to be close enough to the maximum
                if ((UpperKsi - NewKsi) < GRADIENT_DISTANCE)
                    ExtremumFound = true;

                // New upper ksi
                UpperKsi               = NewKsi;
                UpperGradient          = NewGradient;
                UpperKsiFound          = true;
            }

        } while (!ExtremumFound);

        MaximumThicknessKsi = NewKsi;

/*
        // Print Thicknesses
        cout << RelativeSpanwidthPosition << endl;
        for (int i = 0;i<=100;++i)
            cout << i/100. << "     " << getRelativeThickness (Wing, Segment, RelativeSpanwidthPosition,i/100.) << endl;
*/
    }

    // Return the final relative thickness ksi coordinate
    // --------------------------------------------------
    return MaximumThicknessKsi;
}

Und für die Dicke:

    // Upper point
    TIGLReturnCode = tiglWingGetUpperPoint (TIGLHandle,
                                            static_cast<int> (Wing),
                                            CPACSSegment,
                                            CPACSRelativeSpanwidthPosition,
                                            RelativeChordlengthPosition,
                                            &XUpper,
                                            &YUpper,
                                            &ZUpper);
    if (TIGLReturnCode != TIGL_SUCCESS)
    {
        ErrorMessage << TIGLHandle                     << ", "
                     << Wing                           << ", "
                     << CPACSSegment                   << ", "
                     << CPACSRelativeSpanwidthPosition << ", "
                     << RelativeChordlengthPosition    << ", "
                     << &XUpper                        << ", "
                     << &YUpper                        << ", "
                     << &ZUpper;
        throw xTIGLError (TIGLReturnCode,
                          "tiglWingGetUpperPoint",
                          ErrorMessage.str ());
    }

    // Lower point
    TIGLReturnCode = tiglWingGetLowerPoint (TIGLHandle,
                                            static_cast<int> (Wing),
                                            CPACSSegment,
                                            CPACSRelativeSpanwidthPosition,
                                            RelativeChordlengthPosition,
                                            &XLower,
                                            &YLower,
                                            &ZLower);
    if (TIGLReturnCode != TIGL_SUCCESS)
    {
        ErrorMessage << TIGLHandle                     << ", "
                     << Wing                           << ", "
                     << CPACSSegment                   << ", "
                     << CPACSRelativeSpanwidthPosition << ", "
                     << RelativeChordlengthPosition    << ", "
                     << &XLower                        << ", "
                     << &YLower                        << ", "
                     << &ZLower;
        throw xTIGLError (TIGLReturnCode,
                          "tiglWingGetLowerPoint",
                          ErrorMessage.str ());
    }

    // Calculate absolute thickness
    // ----------------------------
    RelativeThickness = sqrt (  pow (XUpper - XLower, 2.)
                              + pow (YUpper - YLower, 2.)
                              + pow (ZUpper - ZLower, 2.));