From a25c5ae17e6ff241c64a0c9d762fb6a7843a5462 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rodrigo=20Ram=C3=B3n?= <79613034+velazquezr35@users.noreply.github.com> Date: Thu, 23 Nov 2023 21:15:30 -0600 Subject: [PATCH] Multiple profiles update Now the code supports N individual PROF objs under the WRAPPER class CPs and dCLs are computed for all of them Check TESTING.py --- geom_discr.py | 60 ++++++++++++++++++- plotter.py | 12 ++-- pot_flow.py | 157 +++++++++++++++++++++++++++++++++----------------- testing.py | 66 ++++++++++++++++++--- 4 files changed, 228 insertions(+), 67 deletions(-) diff --git a/geom_discr.py b/geom_discr.py index 86fc995..b322403 100644 --- a/geom_discr.py +++ b/geom_discr.py @@ -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 @@ -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 = [] diff --git a/plotter.py b/plotter.py index 235647a..a9c06db 100644 --- a/plotter.py +++ b/plotter.py @@ -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 diff --git a/pot_flow.py b/pot_flow.py index 972353a..bbdbcd3 100644 --- a/pot_flow.py +++ b/pot_flow.py @@ -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) @@ -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 \ No newline at end of file diff --git a/testing.py b/testing.py index 825f2de..c1a691f 100644 --- a/testing.py +++ b/testing.py @@ -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) \ No newline at end of file + + #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) \ No newline at end of file