Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Multiple profiles update #2

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Multiple profiles update #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
60 changes: 59 additions & 1 deletion geom_discr.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,51 @@

#Classes

class Assy_wrapper:
'''
Profile assy wrapper

'''
def __init__(self, **kwargs):
self.profiles = []
self.profile_count = 0
self.Ns = []
self.Ms = []
self.N_total = 0
self.M_total = 0

def add_profile(self, profile):
self.profiles.append(profile)
self.profile_count += 1
self.Ns.append(profile.N)
self.Ms.append(profile.M)
self.N_total += profile.N
self.M_total += profile.M

def assy_mtrxs(self):
self.betas = np.copy(self.profiles[0].betas)
self.tgs = np.copy(self.profiles[0].tgs)
self.norms = np.copy(self.profiles[0].norms)
self.x_points = np.copy(self.profiles[0].x_points)
self.y_points = np.copy(self.profiles[0].y_points)
self.dx = np.copy(self.profiles[0].dx)
self.dy = np.copy(self.profiles[0].dy)
self.x_mid = np.copy(self.profiles[0].x_mid)
self.y_mid = np.copy(self.profiles[0].y_mid)
self.dL = np.copy(self.profiles[0].dL)

for prof in self.profiles[1:]:
self.betas = np.append(self.betas, prof.betas)
self.tgs = np.append(self.tgs, prof.tgs, axis = 0)
self.norms = np.append(self.norms, prof.norms, axis = 0)
self.x_points = np.append(self.x_points, prof.x_points)
self.y_points = np.append(self.y_points, prof.y_points)
self.dx = np.append(self.dx, prof.dx)
self.dy = np.append(self.dy, prof.dy)
self.x_mid = np.append(self.x_mid, prof.x_mid)
self.y_mid = np.append(self.y_mid, prof.y_mid)
self.dL = np.append(self.dL, prof.dL)

class Prof_gen:
'''
Profile object
Expand All @@ -38,20 +83,33 @@ class Prof_gen:
- update(): Calculates some geometric data (normals, tangentials, distances, lenghts...)
'''

def __init__(self,**kwargs):
def __init__(self, **kwargs):
if 'from_coords' in kwargs:
if kwargs.get('from_coords'):
self.x_coords = np.array(kwargs.get('x_points'))
self.y_coords = np.array(kwargs.get('y_points'))
if 'traslation' in kwargs:
self.x_offset = kwargs.get('traslation')[0]
self.y_offset = kwargs.get('traslation')[1]
else:
self.x_offset = 0
self.y_offset = 0

def update(self):

self._gen_p_points()
self._transform_geom()
self._calc_thet_midpoints()
self._panel_length()
self.N = len(self.x_points)
self.M = self.N-1
self._gen_norms()


def _transform_geom(self):
self.x_points += self.x_offset
self.y_points += self.y_offset

def _gen_norms(self, **kwargs):
self.norms = []
self.tgs = []
Expand Down
12 changes: 7 additions & 5 deletions plotter.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,12 +27,14 @@ def directors_plot(ax, prof, scale = 0.05):
ax.plot([p0[0],pt[0]], [p0[1], pt[1]], 'r')
ax.plot([p0[0],pn[0]], [p0[1], pn[1]], 'b')

def plot_tester(prof):
def plot_tester(profs):
fig, ax = plt.subplots()
ax.plot(prof.x_points, prof.y_points)
# ax.plot(prof.x_coords, prof.y_coords, 'ro')
# ax.plot(prof.x_mid, prof.y_mid, 'bo')
# directors_plot(ax, prof)

for prof in profs:
ax.plot(prof.x_points, prof.y_points)
# ax.plot(prof.x_coords, prof.y_coords, 'ro')
# ax.plot(prof.x_mid, prof.y_mid, 'bo')
# directors_plot(ax, prof)
ax.grid()
ax.set_aspect('equal')
return fig,ax
Expand Down
157 changes: 103 additions & 54 deletions pot_flow.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,12 @@ def VOR2D_L(gamma_1, gamma_2, X_i, x2, i, j, mode, x1 = 0):
w1 = -gamma_1*((x2-X_i[1]*dtheta)-X_i[0]*np.log(R1/R2)+x2*np.log(R1/R2))/(2*np.pi*x2)
w2 = gamma_2*((x2-X_i[1]*dtheta)-X_i[0]*np.log(R1/R2))/(2*np.pi*x2)
else:
if R1 == 0: #CHECK THIS. OR 0?
u1 = -0.5*gamma_1*(X_i[0]-x2)/x2
u2 = 0.5*gamma_2*X_i[0]/x2
w1 = -gamma_1/(2*np.pi)
w2 = gamma_2/(2*np.pi)
else:
u1 = -gamma_1*(X_i[1]*np.log(R2/R1)+X_i[0]*dtheta - x2*dtheta)/(2*np.pi*x2)
u2 = gamma_2*(X_i[1]*np.log(R2/R1) + X_i[0]*dtheta)/(2*np.pi*x2)

Expand All @@ -54,81 +60,124 @@ def VOR2D_L(gamma_1, gamma_2, X_i, x2, i, j, mode, x1 = 0):

return(u1, u2, w1, w2)

def EVAL_FIELD(x, y, prof, gamma_vect, Vinf,i,j, alone = False):
def EVAL_FIELD(x, y, wrapper, gamma_vect, Vinf,i,j, alone = False):
U = 0
W = 0
for i in range(prof.M):
r0 = [prof.x_points[i], prof.y_points[i]]
X_i = utis.RGLOB_LOC(np.array([x,y]), r0, prof.betas[i])
u1, u2, w1, w2 = VOR2D_L(gamma_vect[i], gamma_vect[i+1], X_i, prof.dL[i],0,0,0, x1 = 0)
u1,w1 = utis.rot_vect(utis.MAT_L2G(prof.betas[i]), np.array([u1,w1]))
u2,w2 = utis.rot_vect(utis.MAT_L2G(prof.betas[i]), np.array([u2,w2]))
U += u1+u2
W += w1+w2

twin_j = 0
loc_counter_j = 0
loc_flag_j = wrapper.Ms[loc_counter_j]
for j in range(wrapper.M_total+wrapper.profile_count-1):
if not j == loc_flag_j:
r0 = [wrapper.x_points[j], wrapper.y_points[j]]
X_i = utis.RGLOB_LOC(np.array([x,y]), r0, wrapper.betas[twin_j])
u1, u2, w1, w2 = VOR2D_L(gamma_vect[j], gamma_vect[j+1], X_i, wrapper.dL[twin_j],0,0,0, x1 = 0)
u1, w1 = utis.rot_vect(utis.MAT_L2G(wrapper.betas[twin_j]), np.array([u1,w1]))
u2, w2 = utis.rot_vect(utis.MAT_L2G(wrapper.betas[twin_j]), np.array([u2,w2]))
U += u1+u2
W += w1+w2
twin_j += 1
else:
loc_counter_j += 1
loc_flag_j += wrapper.Ms[loc_counter_j] + 1
if alone:
return(U, W)
else:
return(U+Vinf[0],W+Vinf[1])

def FLOW_FIELD(X, Y, prof, gamma_vect, Vinf, size, alone):
def FLOW_FIELD(X, Y, wrapper, gamma_vect, Vinf, size, alone):
U,W = np.zeros((2, size, size))
Vinf = utis.LST_ARR(Vinf)
for p in range(size):
for q in range(size):
U[p,q],W[p,q] = EVAL_FIELD(X[p,q], Y[p,q], prof, gamma_vect, Vinf,0,0, alone)
U[p,q],W[p,q] = EVAL_FIELD(X[p,q], Y[p,q], wrapper, gamma_vect, Vinf,0,0, alone)
return(U,W)

def GAMMA_PROD(inf_mat, RHS):
return np.matmul(np.linalg.inv(inf_mat), RHS)

def GEN_INF_MATRX(prof, Vinf, sol='VOR2D_L'):
def GEN_INF_MATRX(wrapper, Vinf, sol='VOR2D_L'):
Vinf = utis.LST_ARR(Vinf)
if sol == 'VOR2D_L':
inf_mat = np.zeros((prof.N, prof.N))
tg_mat = np.zeros((prof.N, prof.N))
RHS = np.zeros((prof.N,1))
for i in range(prof.M):
for j in range(prof.M):
r0 = [prof.x_points[j], prof.y_points[j]]
X_i = utis.RGLOB_LOC(np.array([prof.x_mid[i], prof.y_mid[i]]), r0, prof.betas[j])
u1, u2, w1, w2 = VOR2D_L(1, 1, X_i, prof.dL[j], i, j, mode = 'gen', x1 = 0)
u1,w1 = utis.rot_vect(utis.MAT_L2G(prof.betas[j]), np.array([u1,w1]))
u2,w2 = utis.rot_vect(utis.MAT_L2G(prof.betas[j]), np.array([u2,w2]))
inf_mat[i,j] += np.dot([u1,w1], prof.norms[i])
tg_mat[i,j] += np.dot([u1,w1], prof.tgs[i])
tg_mat[i, j+1] += np.dot([u2,w2], prof.tgs[i])
inf_mat[i,j+1] += np.dot([u2, w2], prof.norms[i])
RHS[i] = -np.dot(Vinf, prof.norms[i])
inf_mat[prof.N-1,0] = 1
inf_mat[prof.N-1, -1] = 1
inf_mat = np.zeros((wrapper.N_total, wrapper.N_total))
tg_mat = np.zeros((wrapper.N_total, wrapper.N_total))
RHS = np.zeros((wrapper.N_total,1))
loc_counter_i = 0
loc_counter_j = 0
loc_flag_i = wrapper.Ms[loc_counter_i]
loc_flag_j = wrapper.Ms[loc_counter_j]
twin_i = 0
for i in range(wrapper.M_total+wrapper.profile_count-1):
if not i == loc_flag_i:
twin_j = 0
loc_counter_j = 0
loc_flag_j = wrapper.Ms[loc_counter_j]
for j in range(wrapper.M_total+wrapper.profile_count-1):
if not j == loc_flag_j:
r0 = [wrapper.x_points[j], wrapper.y_points[j]]
X_i = utis.RGLOB_LOC(np.array([wrapper.x_mid[twin_i], wrapper.y_mid[twin_i]]), r0, wrapper.betas[twin_j])
u1, u2, w1, w2 = VOR2D_L(1, 1, X_i, wrapper.dL[twin_j], i, j, mode = 'gen', x1 = 0)
u1,w1 = utis.rot_vect(utis.MAT_L2G(wrapper.betas[twin_j]), np.array([u1,w1]))
u2,w2 = utis.rot_vect(utis.MAT_L2G(wrapper.betas[twin_j]), np.array([u2,w2]))
inf_mat[i,j] += np.dot([u1,w1], wrapper.norms[twin_i])
tg_mat[i,j] += np.dot([u1,w1], wrapper.tgs[twin_i])
tg_mat[i, j+1] += np.dot([u2,w2], wrapper.tgs[twin_i])
inf_mat[i,j+1] += np.dot([u2, w2], wrapper.norms[twin_i])
twin_j += 1
else:
loc_counter_j += 1
loc_flag_j += wrapper.Ms[loc_counter_j] + 1
RHS[i] = -np.dot(Vinf, wrapper.norms[twin_i])
twin_i += 1
else:
loc_counter_i += 1
loc_flag_i += wrapper.Ms[loc_counter_i] + 1
i_counter = wrapper.Ms[0]
j_counter = 0
for i in range(0, wrapper.profile_count):
inf_mat[i_counter, j_counter] = 1
inf_mat[i_counter, j_counter+wrapper.Ms[i]] = 1
j_counter += wrapper.Ms[i]+1
i_counter += wrapper.Ms[i]+1

return inf_mat, RHS, tg_mat

def coef_CL(prof, gamma_vect, Vinf, tg_mat):
CL = 0
V_induc = np.matmul(tg_mat,gamma_vect)
for i in range(prof.M):
CL += (np.dot(V_induc[i]+Vinf, prof.tgs[i]))*prof.dL[i]
return CL

def CPCL(prof, gamma_vect, Vinf, sol = 'VOR2D_L'):
def CPCL(wrapper, gamma_vect, Vinf, sol = 'VOR2D_L'):
if sol == 'VOR2D_L':
dCP = np.zeros(prof.M)
dLift = np.zeros(prof.M)
dCP = np.zeros(wrapper.M_total)
dLift = np.zeros(wrapper.M_total)
Vmod = np.sqrt(Vinf[0]**2 + Vinf[1]**2)
for i in range(prof.M):
loc_V = np.dot(Vinf, prof.tgs[i])/Vmod
U = 0
W = 0
for j in range(prof.M):
r0 = [prof.x_points[j], prof.y_points[j]]
X_i = utis.RGLOB_LOC(np.array([prof.x_mid[i], prof.y_mid[i]]), r0, prof.betas[j])
u1, u2, w1, w2 = VOR2D_L(gamma_vect[j], gamma_vect[j+1], X_i, prof.dL[j],i,j, mode = 'gen', x1 = 0)
u1, w1 = utis.rot_vect(utis.MAT_L2G(prof.betas[j]), np.array([u1,w1]))
u2, w2 = utis.rot_vect(utis.MAT_L2G(prof.betas[j]), np.array([u2,w2]))
U += u1+u2
W += w1+w2
loc_V += np.dot([U,W], prof.tgs[i])/Vmod
dCP[i] = 1-(loc_V**2)
dLift[i] = 1.225*Vmod*(gamma_vect[i]+gamma_vect[i+1])*0.5*prof.dL[i]
loc_counter_i = 0
loc_flag_i = wrapper.Ms[loc_counter_i]
loc_counter_j = 0
loc_flag_j = wrapper.Ms[loc_counter_j]
twin_i = 0
for i in range(wrapper.M_total):
loc_V = np.dot(Vinf, wrapper.tgs[i])/Vmod
U = 0
W = 0
twin_j = 0
loc_counter_j = 0
loc_flag_j = wrapper.Ms[loc_counter_j]
for j in range(wrapper.M_total+wrapper.profile_count-1):
if not j == loc_flag_j:
r0 = [wrapper.x_points[j], wrapper.y_points[j]]
X_i = utis.RGLOB_LOC(np.array([wrapper.x_mid[i], wrapper.y_mid[i]]), r0, wrapper.betas[twin_j])
u1, u2, w1, w2 = VOR2D_L(gamma_vect[j], gamma_vect[j+1], X_i, wrapper.dL[twin_j],i,twin_j, mode = 'gen', x1 = 0)
u1, w1 = utis.rot_vect(utis.MAT_L2G(wrapper.betas[twin_j]), np.array([u1,w1]))
u2, w2 = utis.rot_vect(utis.MAT_L2G(wrapper.betas[twin_j]), np.array([u2,w2]))
U += u1+u2
W += w1+w2
twin_j += 1
else:
loc_counter_j += 1
loc_flag_j += wrapper.Ms[loc_counter_j] + 1
loc_V += np.dot([U,W], wrapper.tgs[i])/Vmod
dCP[i] = 1-(loc_V**2)
dLift[i] = 1.225*Vmod*(gamma_vect[i+twin_i]+gamma_vect[i+1+twin_i])*0.5*wrapper.dL[i]
if i == loc_flag_i:
loc_counter_i += 1
loc_flag_i += wrapper.Ms[loc_counter_i] + 1
twin_i+=1

return dCP, dLift
66 changes: 59 additions & 7 deletions testing.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -133,13 +133,65 @@ def NACA0012_validation():

###Postpro
fig_CLs, ax_CLs = plotter._testing_CLvs([validation_alphas,validation_CLs], CL_lst)

plotter.plot_tester(NACA0012_prof)
#Standalone run

# if __name__ == '__main__':
# import time
# start = time.time()
# NACA0012_validation()
# end = time.time()
# print('Elapsed [s]: ')
# print(end - start)

if __name__ == '__main__':
import time
start = time.time()
NACA0012_validation()
end = time.time()
print('Elapsed [s]: ')
print(end - start)

#NACA 0012 alpha 8 deg validation case

##Set local values
number_of_points = 50
external_flow_Vinf = 1
air_density = 1.225
#CPvsX data
main_alpha = 8

##NACA0012 profile generation
###Cosine geom. discretization
dist_x = np.linspace(np.radians(180), np.radians(0), number_of_points)
dist_x = 0.5*(1-np.cos(dist_x))
dist_y = yp_NACA0012(dist_x)
###Profile obj init
x_profile = np.append(dist_x, np.flip(dist_x[1:-1]))
y_profile = np.append(-dist_y, np.flip(dist_y[1:-1]))
NACA0012_prof = geom_discr.Prof_gen(from_coords = True, x_points = x_profile, y_points = y_profile)
NACA0012_prof.update()

FLAP_prof = geom_discr.Prof_gen(from_coords = True, x_points = x_profile, y_points = y_profile, traslation = [1.01, -0.1])
FLAP_prof.update()

WRAPPER = geom_discr.Assy_wrapper()
WRAPPER.add_profile(NACA0012_prof)
WRAPPER.add_profile(FLAP_prof)
WRAPPER.assy_mtrxs()

plotter.plot_tester([NACA0012_prof, FLAP_prof])

main_alpha_rads = np.radians(main_alpha)
Vinf = [external_flow_Vinf*np.cos(main_alpha_rads), external_flow_Vinf*np.sin(main_alpha_rads)]

inf_mat, RHS, tg_mat = pot_flow.GEN_INF_MATRX(WRAPPER, Vinf)
gamma_vect = pot_flow.GAMMA_PROD(inf_mat, RHS)
CPs, Lifts = pot_flow.CPCL(WRAPPER, gamma_vect,Vinf)
CL = sum(Lifts)*2/(air_density*1*external_flow_Vinf**2)
print(CL)



#Do not use, so slow...
# size = int(2E3)
# X = np.linspace(-0.1,1.3, size)
# Y = np.linspace(-0.4,0.4, size)
# U,W = pot_flow.FLOW_FIELD(X, Y, WRAPPER, gamma_vect, Vinf,size, False)
# X,Y = np.meshgrid(X,Y)
# fig2, ax2 = plt.subplots()
# plotter.plot_Vmap(ax2, X, Y, U, W)