Results comparision

[rafaellehmkuhl]

First of all, congratulations and thank you for the project. My aerodesign team is currently using it in it's low fidelity aircraft optimizer.

One problem we are having is that the VLM results for the inclination of the CL_alpha curve is being overestimated (by 10% or more), for both rectangular or trapezoidal wings, when compared to XFLR5 analysis.

Have you ever made a comparision for results of this code?

[aqreed]

Hello Rafael, thank you for your interest.

Yes indeed you are right, the code is not ready (please note that it has not been improved since more than one year ago). I have not yet compared it to XFLR5 or AVL, as the thing that bugged me the most when I was working on it was the no-converge of CD. I have just uploaded a notebook to illustrate this issue.

Now that you have pointed out this new issue, I might take another look at it. Please feel free to propose modifications.

[rafaellehmkuhl]

About the CD problem, a friend of mine maybe have addressed it on your code. If I remember correctly there was a gamma missing somewhere. I will check with him.

[RogerWillemsen]

Hi, thank you for this project. I'm using the latest version of the master branch. I was comparing the results to AVL and XFLR5. I found that the drag calculation was wrong and I made some changes. Below you can find the code I changed, in the PyVLM class, vlm method. I used the solved circulation to calculate the induced velocities, using a Aerodynamic Influence Coefficient matrix for total induced velocities. This can be used to calculate drag, which can be found in the same way as the lift. The results I produced with this can be found below. I hope this is helpful to you.

def vlm(self, alpha, print_output=False):
        """
        For a given set of panels applies the VLM theory:

            1. Calculates the induced velocity produced by all the
               associated horseshoe vortices of strength=1 on each panel,
               calculated on its control point where the boundary condition
               will be imposed.
            2. Computes the circulation by solving the linear equation.
            3. Calculates the aerodynamic forces.

        Parameters
        ----------
        alpha : float
            Angle of attack of the wing(degrees)
        print_output : boolean
            Prints the calculation output
        """

        Panel = self.Panels
        rho = self.rho
        alpha = np.deg2rad(alpha)
        V = 1.0

        q_inf = (1 / 2) * rho * (V**2)

        # 1. BOUNDARY CONDITION
        # To impose the boundary condition we must calculate the normal
        # components of (a) induced velocity "Wn" by horshoe vortices of
        # strength=1 and (b) upstream normal velocity "Vinf_n"

        #   (a) INDUCED VELOCITIES
        #     - "Wn", normal component of the total induced velocity by
        #       the horshoe vortices, stored in the matrix "AIC" where the
        #       element Aij is the velocity induced by the horshoe vortex
        #       in panel j on panel i
        #     - also the induced velocity by *only* the trailing vortices
        #        "Wi" on panel i is calculated and stored in the Panel object
        #       attribute "accul_trail_induced_vel"

        N = len(Panel)

        if (type(self.AIC) == int):
            # In case it is the first time the method is called, it proceeds
            # to compute the AIC matrix

            AIC = np.zeros((N, N))  # Aerodynamic Influence Coefficient matrix
            AIC_i = np.zeros((N, N)) # Aerodynamic Influence Coefficient matrix
                                     #for total induced velocities

            for i in range(N):
                panel_pivot = Panel[i]
                CP = panel_pivot.CP

                for j in range(N):
                    panel = Panel[j]
                    Wn, Wi = panel.induced_velocity(CP)
                    AIC[i, j] = Wn  # induced normal velocity by horshoe vortices
                    AIC_i[i, j] = Wi # total induced velocity by horshoe vortices

            self.AIC = AIC
            self.AIC_i = AIC_i

        #   (b) UPSTREAM NORMAL VELOCITY
        #     It will depend on the angle of attack -"alpha"- and the camber
        #     gradient at each panel' position within the local chord

        Vinf_n = np.zeros(N)  # upstream (normal) velocity

        airfoil = NACA4()
        for i in range(N):
            position = Panel[i].chordwise_position
            Panel[i].Vinf_n = -V * (alpha - airfoil.camber_gradient(position))
            Vinf_n[i] = Panel[i].Vinf_n

        # 2. CIRCULATION (Γ or gamma)
        # by solving the linear equation (AX = Y) where X = gamma
        # and Y = Vinf_n

        gamma = np.linalg.solve(self.AIC, Vinf_n)

        # 3. AERODYNAMIC FORCES
        L = 0
        D = 0
        S = 0

        #Calculate induced normal velocities
        V_i = np.matmul(self.AIC_i, gamma)

        for i in range(N):
            Panel[i].gamma = gamma[i]
            Panel[i].l = V * rho * Panel[i].gamma * Panel[i].span
            Panel[i].cl = Panel[i].l / (q_inf * Panel[i].area)

            Panel[i].accul_trail_ind_vel = V_i[i]
            Panel[i].alpha_ind = np.arctan(V_i[i]/V)  # induced AoA(rad)

#            d0 = (Panel[i].l * Panel[i].accul_trail_ind_vel) / V
#            d1 = -rho * abs(Panel[i].gamma) * Panel[i].span * Panel[i].accul_trail_ind_vel
#            d2 = Panel[i].l * np.sin(Panel[i].alpha_ind)
            d3 = -rho * Panel[i].gamma * Panel[i].span * Panel[i].accul_trail_ind_vel

            # print("%8.3f % 8.3f %8.3f" % (d0, d1, d2))

            Panel[i].d = d3
            Panel[i].cd = Panel[i].d / (q_inf * Panel[i].area)

            L += Panel[i].l
            D += Panel[i].d
            S += Panel[i].area

        CL = L / (q_inf * S)
        CD = D / (q_inf * S)

[aqreed]

it is useful indeed @RogerWillemsen, thank you very much. Drag results have been wrong since the beginning, but I have focused on other parts of the code until now. While I do not completely understand the need to have two AIC matrices (can you referred where have you seen this, like papers or books?), results seems fine to me. Open a pull request to include your code!

I am about to release a new version that moves the AIC calculation to the _addwing method (which would be renamed as _addsurface) and some others improvements to make the code look more pythonic. If you preferred to wait to open the PR I believe by the end ot this week the new release will be out.

[RogerWillemsen]

I based my solution on Low Speed Aerodynamics by Katz and Plotkin, Chapter 12, below you can find the part where the drag calculation is described for a lifting-line solution with horseshoe elements. The drag for each element is given by equation 12.10, where w_ind is the induced velocity due to the trailing vortices. To calculate w_ind, a new AIC needs to be constructed with coefficient b_i,j.

I believe that w_ind should be perpendicular to the freestream velocity, which means a correction should be added to account for the angle of attack and camber gradient. However, I'm not yet 100% sure about this.

Also, I found that for some cases, the difference in drag compared to AVL and XFLR5 is still (relatively) large. The error increases especially for higher angles of attack. So I think this is due to the wake representation. Right now the wake is parallel to the chordline, while it should be parallel to the freestream velocity. I might try to improve this, if it is needed for the analyses I want to run.

Anyway, I'll open a pull request once I'm about done :)

[aqreed]

ok, now I understand, and I agree with your solution. Doing:

for i in range(N):
   panel_pivot = Panel[i]
   CP = panel_pivot.CP

   Wi_ = 0
   for j in range(N):
      panel = Panel[j]
      Wn, Wi = panel.induced_velocity(CP)
      AIC[i, j] = Wn  # induced normal velocity by horshoe vortices
      Wi_ += Wi  # induced normal velocity by trailing vortices

      Panel[i].accul_trail_ind_vel = Wi_
      Panel[i].alpha_ind = np.arctan(abs(Wi_)/V)  # induced AoA(rad)

      self.AIC = AIC

is wrong, Wi (induced velocity by trailing vortices, not total induced velocity by horshoe vortices btw) must be stored in a AIC_i matrix as you have done, not summed with Wi_ += Wias I did. Later on with the np.matmul(self.AIC_i, gamma) operation we can finally have the total induced velocity by trailing vortices as required to compute the drag.

I think the first step must be to fix this drag calculation issue, later we can work on the induced AOA to improve the results.

[aqreed]

https://github.com/aqreed/PyVLM/tree/release-0.0.2

@RogerWillemsen I have created a branch for the new release, you can see the changes made in the CHANGELOG file. Please feel free to pull to it, dont forget to add your name to the AUTHORS file if you want.