# 相关讲解视频1：https://www.bilibili.com/video/BV1yx4y1r78v/ # 相关讲解视频2：https://www.bilibili.com/video/BV1bg4y1D7Ya/ # 若使用本文件的任何代码进行科学研究，烦请进行类似致谢： # 中文：该工作使用/借鉴了上海交通大学王有江教授公开的螺旋桨3D几何生成程序，在此表示感谢。 # 英文：We thank Prof. Youjiang Wang from SJTU to share his propeller blade coordinates generating code. import numpy as np import numpy import math import matplotlib.pyplot as plt from scipy import interpolate from scipy import integrate import warnings from datetime import date class foilGeometry: def __init__(self, sharpTE = True): self._sharpTE = sharpTE self.Xp = [] self.Yp = [] self.Zp = [] self.SectionArea = [] self.SectionCoordinate = [] @staticmethod def generate2D(chord, thick2c, camber2c, ctype = "NACA a=0.8(Mod)", ttype = "NACA 66(Mod)" , plot=False): ''' Now only "NACA 66(Mod)" is supported. ''' if (ctype == "NACA a=0.8" or ctype == "NACA a=0.8(Mod)" or ctype == "NACA a=1.0") and ttype == "NACA 66(Mod)": xrate = [0, .005, .0075, .0125, .025, .05, .075, .10, .15, .20, .25, .30, .35, .40, .45, .50, .55, .60, .65, .70, .75, .80, .85, .90, .95, 1.0] if ctype == "NACA a=0.8": yrate = [.0000, .00287, .00404, .00616, .01077, .01841, .02483, .03043, .03985, .04748, .05367, .05863, .06248, .06528, .06709, .06790, .06770, .06644, .06405, .06037, .05514, .04771, .03683, .02435, .01163, .0000] tanf = [0, .48535, .44925, .40359, .34104, .27718, .23868, .21050, .16892, .13734, .11101, .08775, .06634, .04601, .02613, .00620, -.01433, -.03611, -.06010, -.08790, -.12311, -.18412, -.23921, -.25583, -.24904, -.20385] if ctype == "NACA a=0.8(Mod)": yrate = [.0000, .00281, .00396, .00603, .01055, .01803, .02432, .02981, .03903, .04651, .05257, .05742, .06120, .06394, .06571, .06651, .06631, .06508, .06274, .05913, .05401, .04673, .03607, .02452, .01226, .0000] tanf = [0, .47539, .44004, .39531, .33404, .27149, .23378, .20618, .16546, .13452, .10873, .08595, .06498, .04507, .02559, .00607, -.01404, -.03537, -.05887, -.08610, -.12058, -.18034, -.23430, -.24521, -.24521, -.24521] if ctype == "NACA a=1.0": yrate = [.0000, .00250, .00350, .00535, .00930, .01580, .02120, .02585, .03365, .03980, .04475, .04860, .05150, .05355, .05475, .05515, .05475, .05355, .05150, .04860, .04475, .03980, .03365, .02585, 0.1580, .0000] tanf = [0, .42120, .38875, .34770, .29155, .23430, .19995, .17485, .13805, .11030, .08745, .06745, .04925, .03225, .01595, .00000, -.01595, -.03225, -.04925, -.06745, -.08745, -.11030, -.13805, -.17485, -.03430, -.00000] trate = [0, .0665, .0812, .1044, .1466, .2066, .2525, .2907, .3521, .4000, .4363, .4637, .4832, .4952, .5000, .4962, .4846, .4653, .4383, .4035, .3612, .3110, .2532, .1877, .1143, .0333] else: raise ValueError(f"Unsupported cmaberType:{ctype} or thicknessType:{ttype}.") yyrate = numpy.zeros_like(yrate) for k in range(0, len(yrate)): yyrate[k] = yrate[k] / max(yrate) xc = numpy.array(xrate) * chord yc = numpy.array(yyrate) * camber2c * chord yt = numpy.array(trate) * thick2c * chord cosf = 1 / (numpy.array(tanf)**2 + 1) ** 0.5 sinf = cosf * tanf xu = xc - yt*sinf yu = yc + yt*cosf xl = xc + yt*sinf yl = yc - yt*cosf if plot: x2D = xl.tolist()[::-1] + xu.tolist()[1:] y2D = yl.tolist()[::-1] + yu.tolist()[1:] plt.plot(xc, yc, '--') plt.plot(x2D, y2D) plt.axis("equal") plt.grid() plt.show() return xu, yu, xl, yl def saveProfileNX(self, filePath): saveNX_3DCurves(filePath, self.Xp, self.Yp, self.Zp) @staticmethod def CalcSectionArea(x, yu, yd): t = numpy.abs(numpy.array(yu) - numpy.array(yd)) return integrate.simps(t, x) def transform(self, motion): self.Xp = numpy.array(self.Xp) + motion[0] self.Yp = numpy.array(self.Yp) + motion[1] self.Zp = numpy.array(self.Zp) + motion[2] class bladeGeometry(foilGeometry): def __init__(self, sharpTE = True): foilGeometry.__init__(self, sharpTE) self.debug = False self._thickLE = False self._smoothLE = False None def bladeVolume(self, hubR = None): res = 0 if (hubR is None): res = integrate.simps(self.SectionArea, self.SectionCoordinate) else: bladeR = self.SectionCoordinate[-1] r = numpy.linspace(hubR, bladeR, 100) area = interpolate.interp1d( self.SectionCoordinate, self.SectionArea, kind="cubic", fill_value=(self.SectionArea[0], self.SectionArea[-1]))(r) res = integrate.simps(y = area, x = r) return res def fromDesignParameters(self, r2R, chord2D, pitch2D, thick2c, camber2c, rake2D = None, skew = None, D=1, skewIndRake = 1, fName="NACA a=0.8(Mod)", tName="NACA 66(Mod)"): ''' generate 3D coordinate from section profile and section parameters Example: ========= dat4382 = np.array([ [0.20, 0.1740, 1.4550, 0.0000, 0, 0.2494, 0.0430], [0.25, 0.2020, 1.4440, 2.3280, 0, 0.1960, 0.0395], [0.30, 0.2290, 1.4330, 4.6550, 0, 0.1563, 0.0370], [0.40, 0.2750, 1.4120, 9.3630, 0, 0.1069, 0.0344], [0.50, 0.3120, 1.3610, 13.948, 0, 0.0769, 0.0305], [0.60, 0.3370, 1.2850, 18.378, 0, 0.0567, 0.0247], [0.70, 0.3470, 1.2000, 22.747, 0, 0.0421, 0.0199], [0.80, 0.3340, 1.1120, 27.145, 0, 0.0314, 0.0161], [0.90, 0.2800, 1.0270, 31.575, 0, 0.0239, 0.0134], [0.95, 0.2100, 0.9850, 33.788, 0, 0.0229, 0.0140], [1.00, 0.0010, 0.9420, 36.000, 0, 0.0160, 0.0134] ]) geom = bladeGeometry() geom.fromDesignParameters( r2R = dat4382[:,0], chord2D = dat4382[:,1], pitch2D = dat4382[:,2], thick2c = dat4382[:,5], camber2c = dat4382[:,6], rake2D = dat4382[:,4], skew = dat4382[:,3], D=1) geom.visualizeBlade() ''' numR = len(r2R) if (rake2D is None or len(rake2D)==0): rake2D = numpy.zeros_like(r2R) if (skew is None or len(skew)==0): skew = numpy.zeros_like(r2R) Xp = [] Yp = [] Zp = [] SecPos = [] SecArea = [] for ir in range(numR): pitch = pitch2D[ir] * D chord = chord2D[ir] * D rake = rake2D[ir] * D radius = r2R[ir] * D/2 tmp = (pitch**2 + (2 * numpy.pi * radius)**2) ** 0.5 sinPA = pitch/tmp cosPA = 2*numpy.pi*radius/tmp tanPA = pitch/(2*numpy.pi*radius) cskew = numpy.radians(skew[ir]) xu, yu, xl, yl = self.generate2D( chord, thick2c[ir], camber2c[ir], ctype = fName, ttype = tName ) if self._sharpTE: yu[-1] = yl[-1] = 0.5*(yu[-1] + yl[-1]) x2d = numpy.array(xl.tolist()[::-1] + xu.tolist()[1:]) y2d = numpy.array(yl.tolist()[::-1] + yu.tolist()[1:]) x3d = -rake - skewIndRake*cskew*radius*tanPA + (0.5*chord - x2d)*sinPA + y2d*cosPA angle = cskew - ((0.5*chord - x2d)*cosPA - y2d*sinPA) / radius y3d = radius * numpy.sin(angle) z3d = radius * numpy.cos(angle) Xp.append(x3d) Yp.append(y3d) Zp.append(z3d) SecPos.append(radius) SecArea.append(self.CalcSectionArea(0.5*(xu+xl), yu, yl)) self.Xp = Xp self.Yp = Yp self.Zp = Zp self.SectionCoordinate = SecPos self.SectionArea = SecArea return Xp, Yp, Zp def saveNX_3DCurves(filePath, Xp, Yp, Zp): Xp = numpy.array(Xp) Yp = numpy.array(Yp) Zp = numpy.array(Zp) nr = Xp.shape[0] if filePath.endswith(".dat"): filePath = filePath[:-4] for ir in range(nr): fName = filePath + f"_{ir}.dat" xp = Xp[ir,:] yp = Yp[ir,:] zp = Zp[ir,:] with open(fName, "w") as file: for i in range(0, len(xp)): file.write(f"{xp[i]}, {yp[i]}, {zp[i]}\n") if __name__ == "__main__": # r/R, CordLength/D, Pitch/D, skew(deg), Rake/D, MaxThick/C, MaxCamber/C dat4382 = np.array([ [0.20, 0.1740, 1.4550, 0.0000, 0, 0.2494, 0.0430], [0.25, 0.2020, 1.4440, 2.3280, 0, 0.1960, 0.0395], [0.30, 0.2290, 1.4330, 4.6550, 0, 0.1563, 0.0370], [0.40, 0.2750, 1.4120, 9.3630, 0, 0.1069, 0.0344], [0.50, 0.3120, 1.3610, 13.948, 0, 0.0769, 0.0305], [0.60, 0.3370, 1.2850, 18.378, 0, 0.0567, 0.0247], [0.70, 0.3470, 1.2000, 22.747, 0, 0.0421, 0.0199], [0.80, 0.3340, 1.1120, 27.145, 0, 0.0314, 0.0161], [0.90, 0.2800, 1.0270, 31.575, 0, 0.0239, 0.0134], [0.95, 0.2100, 0.9850, 33.788, 0, 0.0229, 0.0140], [1.00, 0.0010, 0.9420, 36.000, 0, 0.0160, 0.0134] ]) geom = bladeGeometry() geom.fromDesignParameters( r2R = dat4382[:,0], chord2D = dat4382[:,1], pitch2D = dat4382[:,2], thick2c = dat4382[:,5], camber2c = dat4382[:,6], rake2D = dat4382[:,4], skew = dat4382[:,3], D=1) geom.visualizeBlade(Z=3) # geom.saveProfileNX("../Data/blade.dat")