forked from jte0419/Panel_Methods
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathCOMPUTE_IJ_SPM.py
79 lines (70 loc) · 4.7 KB
/
COMPUTE_IJ_SPM.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
# FUNCTION - COMPUTE I AND J GEOMETRIC INTEGRALS FOR SOURCE PANEL METHOD
# Written by: JoshTheEngineer
# YouTube : www.youtube.com/joshtheengineer
# Website : www.joshtheengineer.com
# Started : 02/03/19 - Transferred from MATLAB to Python
# - Works as expected
# : 04/28/20 - Fixed E value error handling
#
# PURPOSE
# - Compute the integral expression for constant strength source panels
# - Source panel strengths are constant, but can change from panel to panel
# - Geometric integral for panel-normal : I(ij)
# - Geometric integral for panel-tangential: J(ij)
#
# REFERENCES
# - [1]: Normal Geometric Integral SPM, I(ij)
# Link: https://www.youtube.com/watch?v=76vPudNET6U
# - [2]: Tangential Geometric Integral SPM, J(ij)
# Link: https://www.youtube.com/watch?v=JRHnOsueic8
#
# INPUTS
# - XC : X-coordinate of control points
# - YC : Y-coordinate of control points
# - XB : X-coordinate of boundary points
# - YB : Y-coordinate of boundary points
# - phi : Angle between positive X-axis and interior of panel
# - S : Length of panel
#
# OUTPUTS
# - I : Value of panel-normal integral (Eq. 3.163 in Anderson or Ref [1])
# - J : Value of panel-tangential integral (Eq. 3.165 in Anderson or Ref [2])
import numpy as np
import math as math
np.seterr('raise')
def COMPUTE_IJ_SPM(XC,YC,XB,YB,phi,S):
# Number of panels
numPan = len(XC) # Number of panels/control points
# Initialize arrays
I = np.zeros([numPan,numPan]) # Initialize I integral matrix
J = np.zeros([numPan,numPan]) # Initialize J integral matrix
# Compute integral
for i in range(numPan): # Loop over i panels
for j in range(numPan): # Loop over j panels
if (j != i): # If the i and j panels are not the same
# Compute intermediate values
A = -(XC[i]-XB[j])*np.cos(phi[j])-(YC[i]-YB[j])*np.sin(phi[j]) # A term
B = (XC[i]-XB[j])**2 + (YC[i]-YB[j])**2 # B term
Cn = np.sin(phi[i]-phi[j]) # C term (normal)
Dn = -(XC[i]-XB[j])*np.sin(phi[i])+(YC[i]-YB[j])*np.cos(phi[i]) # D term (normal)
Ct = -np.cos(phi[i]-phi[j]) # C term (tangential)
Dt = (XC[i]-XB[j])*np.cos(phi[i])+(YC[i]-YB[j])*np.sin(phi[i]) # D term (tangential)
E = np.sqrt(B-A**2) # E term
if (E == 0 or np.iscomplex(E) or np.isnan(E) or np.isinf(E)): # If E term is 0 or complex or a NAN or an INF
I[i,j] = 0 # Set I value equal to zero
J[i,j] = 0 # Set J value equal to zero
else:
# Compute I (needed for normal velocity), Ref [1]
term1 = 0.5*Cn*np.log((S[j]**2 + 2*A*S[j] + B)/B) # First term in I equation
term2 = ((Dn-A*Cn)/E)*(math.atan2((S[j]+A),E)-math.atan2(A,E)) # Second term in I equation
I[i,j] = term1 + term2 # Compute I integral
# Compute J (needed for tangential velocity), Ref [2]
term1 = 0.5*Ct*np.log((S[j]**2 + 2*A*S[j] + B)/B) # First term in I equation
term2 = ((Dt-A*Ct)/E)*(math.atan2((S[j]+A),E)-math.atan2(A,E)) # Second term in I equation
J[i,j] = term1 + term2 # Compute J integral
# Zero out any problem values
if (np.iscomplex(I[i,j]) or np.isnan(I[i,j]) or np.isinf(I[i,j])): # If I term is complex or a NAN or an INF
I[i,j] = 0 # Set I value equal to zero
if (np.iscomplex(J[i,j]) or np.isnan(J[i,j]) or np.isinf(J[i,j])): # If J term is complex or a NAN or an INF
J[i,j] = 0 # Set J value equal to zero
return I, J # Return both I and J matrices