solver.py

import numpy as np
from .setup import Parameters, Setup
from .geometry import WingGeometry

class Results():
    def __init__(self, ground_effect, n_panels, X, Xc, ds , T, N, gamma_i, lambda_i, vt, P, cp, rho, wing_chord, Vinf):
        self.ground_effect = ground_effect
        self.rho = rho
        self.Vinf = Vinf
        self.chord = wing_chord

        self.n_panels = n_panels
        self.n_points = n_panels + np.ones_like(n_panels)

        self.Xc = Xc
        self.X = X

        self.ds = ds
        self.T = T
        self.N = N

        self.vortex = gamma_i
        self.source = lambda_i
        
        self.vt = vt
        self.P = P
        self.cp = cp

        self.forces()

    def __len__(self):
        return len(self.n_panels)
    
    def show_results(self):
        print("===============================")
        print(f"   Drag   = {round(self.Fx,2)}")
        print(f"Downforce = {round(self.Fy,2)}")
        print(f"    Cl    = {round(self.Cl,4)}")
        print(f"    Cd    = {round(self.Cd,4)}")
        print(f"  cl/cd   = {round(abs(self.Fy/self.Fx), 2)}")
        print("===============================")

    def panels(self):
        """
        The function takes the number of panels and returns the total number of panels
        :return: The number of panels in the array.
        """
        return np.sum(self.n_panels, dtype='int32')

    def points(self):
        """
        It returns the total number of points 
        :return: The sum of the points in the array.
        """
        return np.sum(self.n_points, dtype='int32')

    def forces(self):
        """
        The function calculates the force on the airfoil due to the pressure distribution
        """        
        if self.ground_effect: 
            l = np.cumsum(self.n_panels[:int(len(self.n_panels)/2)])[-1]
            P = self.P[:l,:]
            Px = P[:,0]*self.N[:l, 0]
            Py = P[:,0]*self.N[:l, 1]
            fx = Px*self.ds[:l]
            fy = Py*self.ds[:l]
        else:
            P = self.P
            Px = P[:,0]*self.N[:, 0]
            Py = P[:,0]*self.N[:, 1]
            fx = Px*self.ds
            fy = Py*self.ds
    
        self.Fx = np.sum(fx)
        self.Fy = np.sum(fy)
        self.Cl = self.Fy/(0.5*self.rho*self.chord*(self.Vinf**2))
        self.Cd = self.Fx/(0.5*self.rho*self.chord*(self.Vinf**2)) 

class MultipleResults():
    def __init__(self) -> None:
        self.wings = []
        self.results = []
    
    def add_solution(self, wing:WingGeometry, sol: Results):
        self.wings.append(wing)
        self.results.append(sol)

def solve(wing_geo: WingGeometry) -> Results:
    """
    Solves the linear system given by the panel method for a given wing
    geometry. The function computes the source (lambda) and vortex (gamma)
    strengths by solving a linear system. This code impose two conditions:

        1. The flow is tangential to the surface of the airfoil - normal
            velocity for each panel is zero.
        2. The Kutta condition at the trailing edge of each airfoil - the
            tangential velocity at the upper and lower panel at the trailing
            edge must be equal.
    
    :param wing_geo: The geometry of the wing to be solved.
    :type wing_geo: WingGeometry
    :return: The solution of the problem and information for the 
        post-processing
    :rtype: Results
    """
    params = wing_geo.params
    Vinf = np.array([params.Vinf, 0])[np.newaxis]

    n_te = len(wing_geo)                            # Number of trailing edges (flaps)    
    n_pts_list = []                                 # List of points of each airfoil
    n_pan_list = []                                 # List of panels in each airfoil
    for i in range(n_te):
        n_pts_list.append(len(wing_geo.foils_geo[i]))
        n_pan_list.append(len(wing_geo.foils_geo[i])-1)
    
    n_pts = np.sum(n_pts_list, dtype='int32')       # Total number of points
    n_pan = np.sum(n_pan_list, dtype='int32')       # Total number of panels

    n_kutta = np.zeros((n_te,2),dtype='int32')      # Kutta conditions points
    for i in range(n_te):
        if i == 0:
            n_kutta[i,0] = int(0)
            n_kutta[i,1] = int(n_pan_list[i] - 1)
        else:
            n_kutta[i,0] = int(n_kutta[i-1,0] + n_pan_list[i])
            n_kutta[i,1] = int(n_kutta[i-1,1] + n_pan_list[i])
    
    X = np.zeros((n_pts, 2))
    c = 0
    for i in range(n_te):
        for j in range(n_pts_list[i]):
            X[c+ j, 0] += wing_geo.foils_geo[i].X[j]
            X[c+ j, 1] += wing_geo.foils_geo[i].Y[j]        
        c += len(wing_geo.foils_geo[i])
    
    mask = np.ones(len(X)-1, dtype=bool)
    if n_te > 1:
        mask[np.cumsum(n_pts_list[:-1])-1] = False

    Xc = ((X[:-1] + X[1:]) / 2)[mask]
    dS = np.sqrt(np.diff(X[:, 0])**2 + np.diff(X[:, 1])**2)[mask]
    if np.any(dS == 0):
        raise ValueError("Panel with zero length detected. Check foil geometry.")
    
    T = (np.diff(X, axis=0)[mask] / dS[:, np.newaxis])
    N = np.zeros_like(T)
    N[:, 0] = -T[:, 1]
    N[:, 1] = T[:, 0]
    
    R_ij = Xc[:, np.newaxis, :] - Xc[np.newaxis, :, :]
    Xc_prime_ij = np.zeros_like(R_ij)
    Xc_prime_ij[:,:, 0] = np.einsum('ijk,jk->ij', R_ij, T)
    Xc_prime_ij[:,:, 1] = np.einsum('ijk,jk->ij', R_ij, N)
    
    dVs_prime = np.zeros_like(Xc_prime_ij)
    dVs_prime[:,:, 0] = 1 / (4*np.pi) * np.log(((Xc_prime_ij[:,:,0] + 0.5*dS[np.newaxis, :])**2 + Xc_prime_ij[:,:,1]**2) /
                                               ((Xc_prime_ij[:,:,0] - 0.5*dS[np.newaxis, :])**2 + Xc_prime_ij[:,:,1]**2))
    with np.errstate(divide='ignore', invalid='ignore'): # Ignore log(0) warnings
        dVs_prime[:,:, 1] = 1/(2*np.pi) * (np.arctan((Xc_prime_ij[:,:,0] + 0.5*dS[np.newaxis, :])/Xc_prime_ij[:,:,1]) - np.arctan((Xc_prime_ij[:,:,0] - 0.5*dS[np.newaxis, :])/Xc_prime_ij[:,:,1]))
    
    dVv_prime = np.zeros_like(dVs_prime)
    dVv_prime[:,:, 0] = dVs_prime[:,:, 1]
    dVv_prime[:,:, 1] = -dVs_prime[:,:, 0]
    
    dVs = np.einsum('jk,ij->ijk', T, dVs_prime[:,:,0]) + np.einsum('jk,ij->ijk', N, dVs_prime[:,:,1])
    dVv = np.einsum('jk,ij->ijk', T, dVv_prime[:,:,0]) + np.einsum('jk,ij->ijk', N, dVv_prime[:,:,1])
    
    b1 = -Vinf.dot(N.T).T
    b2 = -(Vinf.dot(T[n_kutta[:,0], :].T) + Vinf.dot(T[n_kutta[:,1], :].T)).T
    b = np.vstack(
        (
            b1,
            b2
        )
    )

    ids_reduce = np.r_[0, np.cumsum(n_pan_list)[:-1]]
    A11 = np.einsum('ik,ijk->ij', N, dVs)
    A12 = np.add.reduceat(np.einsum('ik,ijk->ij', N, dVv), ids_reduce, axis=1)
    A21 = np.einsum('ik,ijk->ij', T[n_kutta[:,0], :], dVs[n_kutta[:,0], :, :]) + np.einsum('ik,ijk->ij', T[n_kutta[:,1], :], dVs[n_kutta[:,1], :, :])
    A22 = np.add.reduceat(np.einsum('ik,ijk->ij', T[n_kutta[:,0], :], dVv[n_kutta[:,0], :, :]) + np.einsum('ik,ijk->ij', T[n_kutta[:,1], :], dVv[n_kutta[:,1], :, :]), ids_reduce, axis=1)
    A = np.block(
        [
            [A11, A12],
            [A21, A22]
        ]
    )
    
    sol = np.linalg.solve(A,b)
    lambda_s = sol[:n_pan,0]
    gamma_v = sol[n_pan:,0]

    ### POST PROCESSING ###
    V = Vinf + np.einsum('ijk,j->ik', dVs, lambda_s) + np.einsum('ijk,j->ik', np.add.reduceat(dVv, ids_reduce, axis=1), gamma_v)
    Vt = np.einsum('ik,ik->i',V, T)[:, np.newaxis]

    n_pan_list = np.array(n_pan_list, dtype='int32')
    
    cp = (1 - Vt**2/Vinf.dot(Vinf.T))
    P = params.Pinf + 0.5*params.rho*(Vinf.dot(Vinf.T))*cp

    sol = Results(params.ground_effect, n_pan_list, X, Xc, dS , T, N, gamma_v, lambda_s, Vt, P, cp, params.rho, wing_geo.chord, params.Vinf)

    return sol