Upload constant doublet code

invscidnormalpython.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def constant_doublet(clockwise_coordinates, alpha, num_panels):
+    # This program finds the pressure distribution on an arbitrary aerofoil
+    # by representing the surface as a finite number of doublet panels with
+    # constant strength (Dirichlet B.C.)
+    # Original FORTRAN program by Steven Yon, 1989; Low Speed Aerodynamics, 2nd edition, p. 565
+    # Ported to Python/Numpy by Euan French, 2020.
+
+    m = num_panels
+    n = m + 1
+    al = alpha / 57.2958
+
+    # Read in the panel end points
+    ept = clockwise_coordinates  # ept[:, 0] -> x coords, ept[:, 1] -> y coords
+
+    # Convert the panelling to clockwise
+    ep = np.flip(ept, 0)
+
+    # Establish coordinates of panel end points
+    pt1 = ep[:-1, :]
+    pt2 = np.roll(ep, -1, axis=0)[:-1, :]
+
+    # Find panel angles
+    dz = pt2[:, 1] - pt1[:, 1]
+    dx = pt2[:, 0] - pt1[:, 0]
+    th = np.arctan2(dz, dx)
+
+    # Establish collocation points
+    co = ((pt2 - pt1) / 2) + pt1
+    co = co  # Needed for broadcasting
+
+    # Establish influence coefficients
+
+    # Convert collocation coordinates to local panel coordinates
+    x2t = pt2[:, 0] - pt1[:, 0]  # Note, these are equal to dz and dx at lines 20, 21
+    z2t = pt2[:, 1] - pt1[:, 1]
+    xt = co[:, 0:1] - pt1[:, 0:1].T
+    zt = co[:, 1:2] - pt1[:, 1:2].T
+
+    x2 = x2t * np.cos(th) + z2t * np.sin(th)
+    z2 = np.zeros_like(x2)
+    x = xt * np.cos(th) + zt * np.sin(th)
+    z = -xt * np.sin(th) + zt * np.cos(th)
+
+    # Save panel lengths
+    dl = x2
+
+    # Find R and theta components
+    r1 = np.sqrt(x ** 2 + z ** 2)
+    r2 = np.sqrt((x - x2) ** 2 + z ** 2)
+    th1 = np.arctan2(z, x)
+    th2 = np.arctan2(z, x - x2)
+
+    # Compute influence coeffs. A(I, J)
+    a = -0.15916 * (th2 - th1)
+    np.fill_diagonal(a, 0.5)
+
+    # Add wake influence
+    xw = co[:, 0] - pt2[m - 1, 0]  # Have to do m-1 instead of m because of the indexing offset in Python
+    zw = co[:, 1] - pt2[m - 1, 1]
+    dthw = -np.arctan(zw / xw)
+
+    a = np.concatenate((a, np.zeros(n - 1).reshape(-1, 1)), axis=1)
+    a[:, n - 1] = -0.15916 * dthw
+
+    # Build rhs vector
+    b = co[:, 0] * np.cos(al) + co[:, 1] * np.sin(al)
+    b = np.append(b, 0)
+
+    # Add an explicit Kutta condition
+    a = np.concatenate((a, np.zeros(n).reshape(1, -1)), axis=0)
+    a[n - 1, 0] = -1
+    a[n - 1, m - 1] = 1
+    a[n - 1, n - 1] = -1
+
+    # Solve for solution vector of doublet strengths
+    g = np.linalg.solve(a, b)
+
+    # Convert doublet strengths into tangential velocities
+    r = (dl[:m-2] + np.roll(dl, -1)[:m-2]) / 2
+    vel = (np.roll(g, -1)[:m-2] - g[:m-2]) / r
+    cp = 1 - vel ** 2
+
+    return cp, co
+
+
+if __name__ == "__main__":
+    vdvaerofoil = np.loadtxt('vdv15')
+    num_panels = vdvaerofoil.shape[0] - 1
+    cp, co = constant_doublet(vdvaerofoil, 5, num_panels)
+    plt.plot(co[:num_panels-2, 0], cp)
+    plt.plot(co[:, 0], co[:, 1], marker=".", linestyle='None')
+    plt.plot(vdvaerofoil[:, 0], vdvaerofoil[:, 1])
+    plt.xlim(-0.1, 1.1)
+    plt.ylim(1, -1.8)
+    plt.show()