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.
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): Panel_Start=Panel_Start_End Panel_End=Panel_Start_End print(Panel_Start) print(Panel_End) V_Panel=Panel_End-Panel_Start V_R1=Point_Collocation-Panel_Start V_R2=Point_Collocation-Panel_End R1=np.linalg.norm(V_R1) R2=np.linalg.norm(V_R2) print(R1,R2,'Distance') TH1=np.arccos(np.dot(V_R1,V_Panel)/(np.linalg.norm(V_R1)*np.linalg.norm(V_Panel))) TH2=np.arccos(np.dot(V_R2,V_Panel)/(np.linalg.norm(V_R2)*np.linalg.norm(V_Panel))) print(TH1,TH2,'Angles\n') up=(Source_Strength/(4*np.pi))*np.log(R1**2/R2**2) wp=(Source_Strength/(2*np.pi))*(TH2-TH1) 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.
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 Panel_End = Panel_Start_End a=np.arctan2(Panel_End-Panel_Start,Panel_End-Panel_Start)#Panel Angle #Local Frame P_Coll_LxT=Point_Collocation-Panel_Start #Translate x cord Coll Point P_Coll_LyT=Point_Collocation-Panel_Start #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-Panel_Start #Translate x cord Panel End Point P_End_LyT=Panel_End-Panel_Start #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 P_End_Ly=0 #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 TH1=np.arctan2(P_Coll_Ly,P_Coll_Lx) TH2=np.arctan2(P_Coll_Ly,P_Coll_Lx-P_End_Lx) if Point_Collocation == (Panel_Start+Panel_End)/2 and Point_Collocation == (Panel_Start+Panel_End)/2: up = 0 wp = 0.5 print('Diagonal') print(TH1,TH2,'Angles') print(R1,R2,'Distances') else: up=(Source_Strength/(2*np.pi))*np.log(R1/R2) wp=(Source_Strength/(2*np.pi))*(TH2-TH1) 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.
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.