I am currently trying to implement Source Panel method as described in Katz and Plotkin in Low Speed Aerodynamics. I have successfully implemented two previous methods. However, I am fully blocked on the source panel method in implementing the SOR2DC function as is called in the book in section 11.2.1. enter image description here

Theta 1 and Theta 2 Varaible are the same as in the image the same can be said about R1 and R2. P is Point_Collocation.

def Constant_Source_Induced_Vel(Panel_Start_End,Point_Collocation,Source_Strength=1):
    return np.array([up,wp])

This is my current implementation I am not using any of the change of reference system transformations to make the troubleshooting easier.Furthermore, it states in the text that this can be done. My issue currently is that the diagonal of the coefficient matrix should be 0.5 but it is not.This is due to the imposition of a boundary condition which make the surface not be crossed by the flow.

enter image description here

My question is is there any way to make this work?(Make the diagonal 0.5)

I have tried changing the way of obtaining the angles using the arctan instead of the dot product. Furthermore, I have tried changing the reference frame but this seems to cause more issues.



I have edited the code to use the local frame.

def Constant_Source_Induced_Vel(Panel_Start_End, Point_Collocation, Source_Strength=1):
    Panel_Start = Panel_Start_End[0]
    Panel_End = Panel_Start_End[1]

    a=np.arctan2(Panel_End[1]-Panel_Start[1],Panel_End[0]-Panel_Start[0])#Panel Angle

    #Local Frame
    P_Coll_LxT=Point_Collocation[0]-Panel_Start[0] #Translate x cord Coll Point
    P_Coll_LyT=Point_Collocation[1]-Panel_Start[1] #Translate y cord Coll Point
    P_Coll_Lx=P_Coll_LxT*np.cos(a) + P_Coll_LyT*np.sin(a)#Rotate to Local Frame
    P_Coll_Ly=-P_Coll_LxT*np.sin(a) + P_Coll_LyT*np.cos(a)#Rotate to Local Frame
    P_End_LxT=Panel_End[0]-Panel_Start[0] #Translate x cord Panel End Point
    P_End_LyT=Panel_End[1]-Panel_Start[1] #Translate y cord Panel End Point
    P_End_Lx=P_End_LxT*np.cos(a) +P_End_LyT*np.sin(a)#Rotate to Local Frame
    #Calculate R1,R2,TH1,TH2
    R1=(P_Coll_Lx**2 + P_Coll_Ly**2)**0.5
    R2=((P_Coll_Lx-P_End_Lx)**2 + P_Coll_Ly**2)**0.5

    if Point_Collocation[0] == (Panel_Start[0]+Panel_End[0])/2 and Point_Collocation[1] == (Panel_Start[1]+Panel_End[1])/2:
        up = 0
        wp = 0.5
    u = up *np.cos(a) - wp*np.sin(a)
    w = up *np.sin(a) + wp*np.cos(a)
    print(u,w,'Final Velocity')
    return np.array([u, w])

It results in the following matrix.

enter image description here

Although it seems that the matrix is now ok the issues persists as with an increase in the number of panels the sum of source strenghts does not get closer to zero. enter image description here

The panels in this graph start from rear bottom and go in a clockwise direction towards the top surface.The pressure on top and bottom should be the same since this is a symmetric airfoil.

  • 1
    $\begingroup$ I'm not sure what you hope we can provide in terms of help. Your question is not specific enough to actually contain a question (i.e., there is no question mark), and this isn't the kind of site where you can provide a link to several hundred or thousand lines of Python and hope that someone will go through them to point out your bug. $\endgroup$ – Wolfgang Bangerth Mar 10 at 11:42
  • $\begingroup$ I think that the panel method is the name given to the boundary element method in potential flow in. That's not my area, though. So, take it with a grain of salt. $\endgroup$ – nicoguaro Mar 10 at 18:29
  • 1
    $\begingroup$ Start with figuring out why some of the entries are NaNs. That should be easy enough to track down. $\endgroup$ – Wolfgang Bangerth Mar 11 at 9:47
  • $\begingroup$ I have added another method which should also work. $\endgroup$ – serviola_29 Mar 11 at 13:18
  • $\begingroup$ This contains the code which I am trying to replicated in the fortan. github.com/alwinw/KatzPlotkinPy/blob/master/tests/f77/SOR2DC.f $\endgroup$ – serviola_29 Mar 11 at 16:56

After the last edit there are several updates. Using Lorena Barbas work as reference it was confirmed that the source time its strength is supposed to be zero(She has excellent jupiter notebooks on source panel method using an integration approach). Therefore, the last bit of code was correct. Furthermore, it can also be observed that the pressure coefficients are not correct in the bottom surface. This is due to the fact in the calculation of the pressure coefficient the velocity due to was calculated in this surface was used with a dot product. The sign on this was incorrect and this gave the coefficients of pressure on the bottom surface a large distortion. The streamlines also look pretty good.

enter image description here enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.