From e4468b3020cb320140945d4521f55f5fcc787379 Mon Sep 17 00:00:00 2001 From: Davide Oberto Date: Wed, 18 Feb 2026 14:32:08 +0100 Subject: [PATCH 1/8] Added possibility to specify thickness and camber max for each radius in the blade __init__ method. To do so, corresponding methods in baseprofile implemented. Added tests too --- bladex/blade.py | 22 +++++++++++++--- bladex/profile/baseprofile.py | 47 ++++++++++++++++++++++++++++++++--- tests/test_profiles.py | 19 ++++++++++++++ 3 files changed, 81 insertions(+), 7 deletions(-) diff --git a/bladex/blade.py b/bladex/blade.py index 74f52f5..fac2567 100644 --- a/bladex/blade.py +++ b/bladex/blade.py @@ -84,6 +84,10 @@ class Blade(object): radial section of the blade. :param array_like skew_angles: 1D array, contains the skew angles (in degrees) for each radial section of the blade. + :param array_like thickness: 1D array, contains the value of the + thickness of each section, specified if sections have not the desired thickness. + :param array_like camber: 1D array, contains the value of the + camber of each section, specified if sections have note the desired camber. Note that, each of the previous array_like parameters must be consistent with the other parameters in terms of the radial ordering of the blade @@ -149,7 +153,7 @@ class Blade(object): """ def __init__(self, sections, radii, chord_lengths, pitch, rake, - skew_angles): + skew_angles, thickness=None, camber=None): # Data are given in absolute values self.sections = sections self.n_sections = len(sections) @@ -158,6 +162,8 @@ def __init__(self, sections, radii, chord_lengths, pitch, rake, self.pitch = pitch self.rake = rake self.skew_angles = skew_angles + self.thickness = thickness + self.camber = camber self._check_params() self.conversion_factor = 1000 # to convert units if necessary @@ -328,12 +334,20 @@ def apply_transformations(self, reflect=True): # Translate reference point into origin self.sections[i].translate(-self.sections[i].reference_point) - if reflect: - self.sections[i].reflect() - # Scale the unit chord to actual length. self.sections[i].scale(self.chord_lengths[i]) + # Setting thickness max is required + if self.thickness is not None: + self.sections[i].set_thickness_max(self.thickness[i]) + + # Setting camber max is required + if self.camber is not None and i < self.n_sections-1: + self.sections[i].set_camber_line_max(self.camber[i]) + + if reflect: + self.sections[i].reflect() + # Rotate according to the pitch angle. # Since the current orientation system is not standard (It is # left-handed Cartesian orientation system, where Y-axis points diff --git a/bladex/profile/baseprofile.py b/bladex/profile/baseprofile.py index 1d08637..b307d8d 100644 --- a/bladex/profile/baseprofile.py +++ b/bladex/profile/baseprofile.py @@ -379,9 +379,6 @@ def deform_camber_line(self, percent_change, n_interpolated_points=None): :param float percent_change: percentage of change of the maximum camber. Default value is None - :param bool interpolate: if True, the interpolated coordinates are - used to compute the camber line and foil's thickness, otherwise - the original discrete coordinates are used. Default value is False. :param int n_interpolated_points: number of points to be used for the equally-spaced sample computations. If None then there is no interpolation, unless the arrays x_up != x_down elementwise which @@ -415,6 +412,50 @@ def deform_camber_line(self, percent_change, n_interpolated_points=None): self.yup_coordinates = self.camber_line[1] + half_thickness self.ydown_coordinates = self.camber_line[1] - half_thickness + def set_camber_line_max(self, camber_max): + """ + Deform camber line according to a given maximum value. + The percentage of camber wrt to the x/c coordinate does not change, i.e., + we rescale the camber value by a scalar + Also reconstructs the deformed airfoil's coordinates. + + Thus, the percentage of change is defined as follows: + + .. math:: + \\frac{\\text{new magnitude of max camber - old magnitude of + maximum \ + camber}}{\\text{old magnitude of maximum camber}} * 100 + + A positive percentage means the new camber is larger than the max + camber value, while a negative percentage indicates the new value + is smaller. + + We note that the method works only for airfoils in the reference + position, i.e. chord line lies on the X-axis and the foil is not + rotated, since the measurements are based on the Y-values of the + airfoil coordinates, hence any measurements or scalings will be + inaccurate for the foils not in their reference position. + + :param float camber_max: maximum camber to be set. + """ + old_camber_max = self.max_camber() + percent_scaling_factor = 100*(camber_max - old_camber_max) / (old_camber_max) + # print(percent_scaling_factor) + self.deform_camber_line(percent_scaling_factor) + + + def set_thickness_max(self, thickness_max): + """ + Deform y_up and y_down coordinates to have the desired thickness max value. + To do so, we compute the ratio between olt thickness and new one + + :param float thickness_max: maximum thickness to be set. + """ + old_thickness = self.max_thickness() + ratio_thickness = thickness_max / old_thickness if old_thickness != 0. else 0. + self.yup_coordinates *= ratio_thickness + self.ydown_coordinates *= ratio_thickness + @property def yup_curve(self): """ diff --git a/tests/test_profiles.py b/tests/test_profiles.py index 1b590c8..3fb309b 100644 --- a/tests/test_profiles.py +++ b/tests/test_profiles.py @@ -14,6 +14,25 @@ def create_custom_profile(): return CustomProfile(xup=xup, yup=yup, xdown=xdown, ydown=ydown) +class TestBaseProfileMethods(TestCase): + def test_max_thickness(self): + profile = NacaProfile(digits='2412') + self.assertAlmostEqual(profile.max_thickness(), .12, delta=1e-4) + + def test_set_thickness_max(self): + profile = NacaProfile(digits='2412') + profile.set_thickness_max(.24) + self.assertAlmostEqual(profile.max_thickness(), .24, delta=1e-4) + + def test_max_camber(self): + profile = NacaProfile(digits='3412') + self.assertAlmostEqual(profile.max_camber(), .03, delta=1e-4) + + def test_set_camber_max(self): + profile = NacaProfile(digits='3412') + profile.set_camber_line_max(.05) + self.assertAlmostEqual(profile.max_camber(), .05, delta=1e-4) + class TestCustomProfile(TestCase): def test_inheritance_custom(self): self.assertTrue(issubclass(CustomProfile, ProfileInterface)) From 412cb5b9ef9afddd3c09b46e8d349b444a9a1ddf Mon Sep 17 00:00:00 2001 From: Davide Oberto Date: Wed, 18 Feb 2026 14:45:17 +0100 Subject: [PATCH 2/8] Fixed uncomplete RunTimeError message --- bladex/profile/customprofile.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bladex/profile/customprofile.py b/bladex/profile/customprofile.py index 72285c3..0ea9fbe 100644 --- a/bladex/profile/customprofile.py +++ b/bladex/profile/customprofile.py @@ -74,8 +74,8 @@ def __init__(self, **kwargs): raise RuntimeError( """Input arguments should be the section coordinates (xup, yup, xdown, ydown) or the section parameters - (camber_perc, thickness_perc, - camber_max, thickness_max, chord_perc).""") + (camber_perc, thickness_perc, camber_max, + thickness_max, chord_perc, chord_len).""") def generate_parameters(self, convention='british'): return super().generate_parameters(convention) From dc9e600b6dd079f4201d4700d68b3bc41c542eb2 Mon Sep 17 00:00:00 2001 From: Davide Oberto Date: Mon, 23 Feb 2026 14:44:28 +0100 Subject: [PATCH 3/8] Fix: generalized reversepropeller code to work with generic blades with z as main direction. Also cleaned a little bit the code --- bladex/reversepropeller.py | 108 ++++++++++++------------------------- 1 file changed, 34 insertions(+), 74 deletions(-) diff --git a/bladex/reversepropeller.py b/bladex/reversepropeller.py index 423b0e3..7c79d2f 100644 --- a/bladex/reversepropeller.py +++ b/bladex/reversepropeller.py @@ -5,63 +5,40 @@ import os, errno import os.path -import matplotlib.pyplot as plt import numpy as np import csv -from bladex import NacaProfile, CustomProfile, Blade, Propeller, Shaft -from OCC.Core.IGESControl import (IGESControl_Reader, IGESControl_Writer, - IGESControl_Controller_Init) +from bladex import CustomProfile, Blade, Propeller, Shaft +from OCC.Core.IGESControl import (IGESControl_Reader) from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_MakeVertex,\ BRepBuilderAPI_MakeEdge, BRepBuilderAPI_MakeWire, BRepBuilderAPI_MakeSolid, \ - BRepBuilderAPI_Sewing, BRepBuilderAPI_Transform, BRepBuilderAPI_NurbsConvert, \ - BRepBuilderAPI_MakeFace -from OCC.Core.BRep import (BRep_Tool, BRep_Builder, BRep_Tool_Curve, - BRep_Tool_CurveOnSurface) + BRepBuilderAPI_Sewing, BRepBuilderAPI_NurbsConvert +from OCC.Core.BRep import BRep_Tool import OCC.Core.TopoDS -from OCC.Core.BRepAlgoAPI import BRepAlgoAPI_Fuse -# from OCC.Core.BRepAlgoAPI import BRepAlgoAPI_Section -from OCC.Core.TopTools import TopTools_ListOfShape, TopTools_MapOfShape +from OCC.Core.TopTools import TopTools_ListOfShape from OCC.Core.TopExp import TopExp_Explorer from OCC.Core.TopAbs import TopAbs_VERTEX, TopAbs_EDGE, TopAbs_FACE, TopAbs_WIRE -from OCC.Core.StlAPI import StlAPI_Writer -from OCC.Core.BRepMesh import BRepMesh_IncrementalMesh -from OCC.Core.gp import gp_Pnt, gp_Pnt2d, gp_Dir, gp_Ax1, gp_Ax2, gp_Trsf, gp_Pln, gp_Vec, gp_Lin -from OCC.Core.TopoDS import TopoDS_Shape -from OCC.Core.TColgp import TColgp_HArray1OfPnt, TColgp_Array1OfPnt -from OCC.Core.GeomAPI import GeomAPI_Interpolate, GeomAPI_IntCS, GeomAPI_ProjectPointOnSurf +from OCC.Core.gp import gp_Pnt, gp_Dir, gp_Ax2, gp_Vec +from OCC.Core.TColgp import TColgp_Array1OfPnt from OCC.Core.BRepAdaptor import BRepAdaptor_Curve from OCC.Core.GCPnts import GCPnts_AbscissaPoint from OCC.Core.BRep import BRep_Tool -from OCC.Core.IntTools import IntTools_FClass2d -from OCC.Core.TopAbs import TopAbs_IN from OCC.Core.BRepExtrema import BRepExtrema_DistShapeShape -from OCC.Core.TopoDS import topods, TopoDS_Edge, TopoDS_Compound -from subprocess import call -from OCC.Core.IntCurvesFace import IntCurvesFace_ShapeIntersector -from OCC.Core.Adaptor3d import Adaptor3d_Curve -from OCC.Core.Geom import Geom_Line +from OCC.Core.TopoDS import topods from OCC.Display.SimpleGui import init_display -from OCC.Core.BRepGProp import (brepgprop_LinearProperties, - brepgprop_SurfaceProperties, - brepgprop_VolumeProperties) -from OCC.Core.GProp import GProp_GProps from OCC.Core.BRepPrimAPI import BRepPrimAPI_MakeCylinder from OCC.Core.GeomLProp import GeomLProp_SLProps from OCC.Core.GCPnts import GCPnts_AbscissaPoint from OCC.Core.BRepAdaptor import (BRepAdaptor_Curve, BRepAdaptor_CompCurve) -from OCC.Core.GCPnts import GCPnts_UniformDeflection from OCC.Core.GeomAPI import GeomAPI_PointsToBSpline -from OCC.Core.GeomAbs import (GeomAbs_C0, GeomAbs_G1, GeomAbs_C1, GeomAbs_G2, - GeomAbs_C2, GeomAbs_C3, GeomAbs_CN) -from OCC.Core.GeomProjLib import geomprojlib -from OCC.Core.TopExp import topexp -from OCC.Core.IntCurvesFace import IntCurvesFace_ShapeIntersector +from OCC.Core.GeomAbs import GeomAbs_C2 +from OCC.Core.Bnd import Bnd_Box +from OCC.Core.BRepBndLib import brepbndlib +from OCC.Core.BRepAlgoAPI import BRepAlgoAPI_Section +from OCC.Display.SimpleGui import init_display -from scipy.spatial import Voronoi, voronoi_plot_2d +from scipy.spatial import Voronoi from scipy.interpolate import interp1d -import matplotlib.pyplot as plt - def distance(P1, P2): """ @@ -249,7 +226,7 @@ def _extract_solid_from_file(self): sewer.Perform() result_sewed_blade = sewer.SewedShape() blade_solid_maker = BRepBuilderAPI_MakeSolid() - blade_solid_maker.Add(OCC.Core.TopoDS.topods_Shell(result_sewed_blade)) + blade_solid_maker.Add(OCC.Core.TopoDS.topods.Shell(result_sewed_blade)) if (blade_solid_maker.IsDone()): self.blade_solid = blade_solid_maker.Solid() else: @@ -263,23 +240,27 @@ def _build_cylinder(self, radius): section taken into account. This method is applied to all the radii in list 'radii_list' given as input. """ - axis = gp_Ax2(gp_Pnt(-0.5 * self.conversion_unit, 0.0, 0.0), + # Base point such that cylinder intersect all section + bbox = Bnd_Box() + brepbndlib.Add(self.blade_solid, bbox) + xmin, ymin, zmin, xmax, ymax, zmax = bbox.Get() + + axis = gp_Ax2(gp_Pnt(xmin-0.2*abs(xmin), 0.0, 0.0), gp_Dir(1.0, 0.0, 0.0)) - self.cylinder = BRepPrimAPI_MakeCylinder(axis, radius, 1 * - self.conversion_unit).Shape() + self.cylinder = BRepPrimAPI_MakeCylinder(axis, radius, 1.2*(xmax-xmin)).Shape() faceCount = 0 faces_explorer = TopExp_Explorer(self.cylinder, TopAbs_FACE) while faces_explorer.More(): - face = OCC.Core.TopoDS.topods_Face(faces_explorer.Current()) + face = OCC.Core.TopoDS.topods.Face(faces_explorer.Current()) faceCount += 1 # Convert the cylinder faces into Non Uniform Rational Basis-Splines geometry (NURBS) nurbs_converter = BRepBuilderAPI_NurbsConvert(face) nurbs_converter.Perform(face) nurbs_face_converter = nurbs_converter.Shape() - nurbs_face = OCC.Core.TopoDS.topods_Face(nurbs_face_converter) + nurbs_face = OCC.Core.TopoDS.topods.Face(nurbs_face_converter) - surface = BRep_Tool.Surface(OCC.Core.TopoDS.topods_Face(nurbs_face)) + surface = BRep_Tool.Surface(OCC.Core.TopoDS.topods.Face(nurbs_face)) self.bounds = 0.0 self.bounds = surface.Bounds() @@ -301,7 +282,7 @@ def _build_intersection_cylinder_blade(self, radius): """ # Construction of the section lines between two shapes (in this case the # blade and the lateral face of the cylinder) - section_builder = BRepAlgo_Section(self.blade_solid, + section_builder = BRepAlgoAPI_Section(self.blade_solid, self.cylinder_lateral_face, False) # Define and build the parametric 2D curve (pcurve) for the section lines defined above section_builder.ComputePCurveOn2(True) @@ -320,9 +301,8 @@ def _build_intersection_cylinder_blade(self, radius): edgeExplorer.Next() wire_maker.Add(edgeList) self.wire = wire_maker.Wire() - # Create a 3D curve from a wire self.curve_adaptor = BRepAdaptor_CompCurve( - OCC.Core.TopoDS.topods_Wire(self.wire)) + OCC.Core.TopoDS.topods.Wire(self.wire)) # Length of the curve section (ascissa curvilinea) self.total_section_length = GCPnts_AbscissaPoint.Length( self.curve_adaptor) @@ -366,8 +346,6 @@ def _voronoi_points(self, radius): # Create the Voronoi diagram for the 2D points vor = Voronoi(self.param_plane_points) vor_points = [] - min_u = 1e8 - self.start = vor.vertices[0] # Find in an iterative way the Voronoi points for i in range(len(vor.vertices)): p = [vor.vertices[i][0], vor.vertices[i][1]] @@ -375,33 +353,15 @@ def _voronoi_points(self, radius): p[1], self.param_plane_points, include_edges=True)): - if (p[0] < min_u): - min_u = p[0] - self.start = p vor_points.append(p) - optimized_vor_points = optimized_path(vor_points, self.start) - us = [] - vs = [] - uss = [] - vss = [] + optimized_vor_points = sorted(vor_points, key=lambda x: x[0]) + #vor_us includes the X coordinates and vor_vs the Z coordinates self.vor_us = [] self.vor_vs = [] - for p in self.orig_param_plane_points: - us.append(p[0]) - vs.append(p[1]) - for p in self.param_plane_points: - uss.append( - p[0] - ) # coordinates of the points of the original parametric curve - vss.append(p[1]) for p in optimized_vor_points: - self.vor_us.append( - p[0]) # coordinates of the points of the Voronoi map - self.vor_vs.append( - p[1] - ) #vor_us includes the X coordinates and vor_vs the Z coordinates - us.append(self.orig_param_plane_points[0][0]) - vs.append(self.orig_param_plane_points[0][1]) + # coordinates of the points of the Voronoi map + self.vor_us.append(p[0]) + self.vor_vs.append(p[1]) def _camber_curve(self, radius): """ @@ -427,7 +387,7 @@ def _camber_curve(self, radius): # BUild the camber line from 3D edges self.edge_3d = BRepBuilderAPI_MakeEdge(camber_curve).Edge() camber_curve_adaptor = BRepAdaptor_Curve( - OCC.Core.TopoDS.topods_Edge(self.edge_3d)) + OCC.Core.TopoDS.topods.Edge(self.edge_3d)) # Compute the length of the camber curve for each section camber_curve_length = GCPnts_AbscissaPoint.Length(camber_curve_adaptor) relative_tolerance = 2.5e-3 @@ -544,7 +504,7 @@ def _camber_curve(self, radius): full_camber_wire = wire_maker.Wire() self.full_camber_curve_adaptor = BRepAdaptor_CompCurve( - OCC.Core.TopoDS.topods_Wire(full_camber_wire)) + OCC.Core.TopoDS.topods.Wire(full_camber_wire)) self.full_camber_length = GCPnts_AbscissaPoint.Length( self.full_camber_curve_adaptor) From fc572f730d78e5c0d51876be2856f60cf47bb436 Mon Sep 17 00:00:00 2001 From: Davide Oberto Date: Tue, 3 Mar 2026 08:28:49 +0100 Subject: [PATCH 4/8] Added case where blade's geometry is made by faces, like for BladeX generated blades --- bladex/reversepropeller.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/bladex/reversepropeller.py b/bladex/reversepropeller.py index 7c79d2f..3579f3b 100644 --- a/bladex/reversepropeller.py +++ b/bladex/reversepropeller.py @@ -59,7 +59,6 @@ def optimized_path(coords, start=None): pass_by.remove(nearest) return path - def point_inside_polygon(x, y, poly, include_edges=True): """ Test if point (x,y) is inside polygon poly. @@ -220,9 +219,18 @@ def _extract_solid_from_file(self): iges_reader = IGESControl_Reader() iges_reader.ReadFile(self.iges_file) iges_reader.TransferRoots() - self.blade_compound = iges_reader.Shape() sewer = BRepBuilderAPI_Sewing(self.tolerance_solid) - sewer.Add(self.blade_compound) + # Case where we have Faces and not closed surface. + # This is the case for BladeX generated blade + if iges_reader.Shape().ShapeType() == 4: + self.blade_compound = iges_reader.OneShape() + exp = TopExp_Explorer(self.blade_compound, TopAbs_FACE) + while exp.More(): + sewer.Add(exp.Current()) + exp.Next() + else: + self.blade_compound = iges_reader.Shape() + sewer.Add(self.blade_compound) sewer.Perform() result_sewed_blade = sewer.SewedShape() blade_solid_maker = BRepBuilderAPI_MakeSolid() From d37f9261442ce729d96bd7e2d4a9a109371f0a7e Mon Sep 17 00:00:00 2001 From: Davide Oberto Date: Wed, 4 Mar 2026 16:22:11 +0100 Subject: [PATCH 5/8] Refactored reversepropeller class by setting a base class and two inheritant classes --- bladex/reversepropeller/__init__.py | 10 + .../reversepropeller/basereversepropeller.py | 392 ++++++++++++++++ .../reversepropeller.py | 304 ++---------- .../reversepropellerbladex.py | 436 ++++++++++++++++++ .../reversepropellerinterface.py | 52 +++ 5 files changed, 921 insertions(+), 273 deletions(-) create mode 100644 bladex/reversepropeller/__init__.py create mode 100644 bladex/reversepropeller/basereversepropeller.py rename bladex/{ => reversepropeller}/reversepropeller.py (67%) create mode 100644 bladex/reversepropeller/reversepropellerbladex.py create mode 100644 bladex/reversepropeller/reversepropellerinterface.py diff --git a/bladex/reversepropeller/__init__.py b/bladex/reversepropeller/__init__.py new file mode 100644 index 0000000..904f92b --- /dev/null +++ b/bladex/reversepropeller/__init__.py @@ -0,0 +1,10 @@ +""" +Profile init +""" +__all__ = ['ReversePropellerInterface', 'BaseReversePropeller', + 'ReversePropeller', 'ReversePropellerBladeX'] + +from .reversepropellerinterface import ReversePropellerInterface +from .basereversepropeller import BaseReversePropeller +from .reversepropeller import ReversePropeller +from .reversepropellerbladex import ReversePropellerBladeX diff --git a/bladex/reversepropeller/basereversepropeller.py b/bladex/reversepropeller/basereversepropeller.py new file mode 100644 index 0000000..9184cd2 --- /dev/null +++ b/bladex/reversepropeller/basereversepropeller.py @@ -0,0 +1,392 @@ +import os, errno +import os.path +import time +import numpy as np +import csv +from bladex import CustomProfile, Blade, Propeller, Shaft +from .reversepropellerinterface import ReversePropellerInterface +from OCC.Core.IGESControl import (IGESControl_Reader) +from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_MakeVertex,\ + BRepBuilderAPI_MakeEdge, BRepBuilderAPI_MakeWire, BRepBuilderAPI_MakeSolid, \ + BRepBuilderAPI_Sewing, BRepBuilderAPI_NurbsConvert +from OCC.Core.BRep import BRep_Tool +import OCC.Core.TopoDS +from OCC.Core.TopTools import TopTools_ListOfShape +from OCC.Core.TopExp import TopExp_Explorer +from OCC.Core.TopAbs import TopAbs_VERTEX, TopAbs_EDGE, TopAbs_FACE, TopAbs_WIRE +from OCC.Core.gp import gp_Pnt, gp_Dir, gp_Ax2, gp_Vec, gp_Pln +from OCC.Core.TColgp import TColgp_Array1OfPnt +from OCC.Core.BRepAdaptor import BRepAdaptor_Curve +from OCC.Core.GCPnts import GCPnts_AbscissaPoint +from OCC.Core.BRep import BRep_Tool +from OCC.Core.TopoDS import topods +from OCC.Display.SimpleGui import init_display +from OCC.Core.BRepPrimAPI import BRepPrimAPI_MakeCylinder +from OCC.Core.GeomLProp import GeomLProp_SLProps +from OCC.Core.GCPnts import GCPnts_AbscissaPoint +from OCC.Core.BRepAdaptor import (BRepAdaptor_Curve, + BRepAdaptor_CompCurve) +from OCC.Core.GeomAPI import GeomAPI_PointsToBSpline +from OCC.Core.GeomAbs import GeomAbs_C2 +from OCC.Core.Bnd import Bnd_Box +from OCC.Core.BRepBndLib import brepbndlib +from OCC.Core.BRepAlgoAPI import BRepAlgoAPI_Section +from OCC.Core.Geom import Geom_Plane +from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_MakeFace +from scipy.interpolate import interp1d +import matplotlib.pyplot as plt + + +class BaseReversePropeller(ReversePropellerInterface): + """ + Base class for the extraction of the parameters of a blade and reconstruction of the parametrized + blade and of the entire propeller. + + Depending on the source code of the generated blade, different child classes can be chosen. + + + Given the IGES file of the blade, the following parameters are extracted: + + - :math:`(X, Y)` coordinates of the blade cylindrical sections after + being expanded in 2D to create airfoils. + + - Chord length, for each cylindrical section. + + - Pitch :math:`(P)`, for each cylindrical section. + + - Rake :math:`(k)`, in distance units, for each cylindrical section. + + - Skew angle :math:`(\\theta_s)`, for each cylindrical section, expressed + in degrees. + + The parameters can be saved in a csv file with the method + 'save_global_parameters' and used to reconstruct the blade, which is saved + in a IGES file (method 'reconstruct_blade'). + The coordinates of each section can be extracted with the method + 'reconstruct_sections'. + Given the shaft, the whole propeller can be reconstructed and saved in a + IGES file with the method 'reconstruct_propeller'. + + -------------------------- + + :param filename: path to the IGES file of the blade. + :param list radii_list: list which contains the radii values of the + sectional profiles. + :param num_points_top_bottom: number of points used to interpolate each + sectional profile. + """ + def __init__(self, filename, radii_list, num_points_top_bottom): + + self.iges_file = filename #filename is the path to the file iges + self.radii_list = radii_list #radii at which we want to measure properties + self.section_points_number = 1200 #number of points for each section profile + self.airfoil_points_number = 1000 + self.num_points_top_bottom = num_points_top_bottom #number of points used to reconstruct the profile, equal for top and lower parts of the profile + # self.num_sections = len(self.radii_list) + # Initialize things we want to compute from the blade IGES file + self.pitch_angles_list = [] + self.pitch_list = [] + self.skew_angles_list = [] + self.skew_list = [] + self.rake_list = [] + self.chord_length_list = [] + self.blade_solid = None + self.blade_compound = None + self.tolerance_solid = 5e-2 + self._extract_solid_from_file() + num_sections = len(radii_list) + self.xup = np.zeros((num_sections, self.num_points_top_bottom)) + self.xdown = np.zeros((num_sections, self.num_points_top_bottom)) + self.yup = np.zeros((num_sections, self.num_points_top_bottom)) + self.ydown = np.zeros((num_sections, self.num_points_top_bottom)) + # Lists of OCC objects used for plotting in display_* methods + self.cylinder_lateral_faces_list = [] + self.section_wires_list = [] + self.camber_wires_list = [] + + def _extract_solid_from_file(self): + """ + Private method that reads the IGES file of the blade and constructs + the correspondent solid object + """ + # Extract the blade from the IGES file; build it as a solid with OCC + iges_reader = IGESControl_Reader() + iges_reader.ReadFile(self.iges_file) + iges_reader.TransferRoots() + sewer = BRepBuilderAPI_Sewing(self.tolerance_solid) + # Case where we have Faces and not closed surface. + # This is the case for BladeX generated blade + if iges_reader.Shape().ShapeType() == 4: + self.blade_compound = iges_reader.OneShape() + exp = TopExp_Explorer(self.blade_compound, TopAbs_FACE) + while exp.More(): + sewer.Add(exp.Current()) + exp.Next() + else: + self.blade_compound = iges_reader.Shape() + sewer.Add(self.blade_compound) + sewer.Perform() + result_sewed_blade = sewer.SewedShape() + blade_solid_maker = BRepBuilderAPI_MakeSolid() + blade_solid_maker.Add(OCC.Core.TopoDS.topods.Shell(result_sewed_blade)) + if (blade_solid_maker.IsDone()): + self.blade_solid = blade_solid_maker.Solid() + else: + self.blade_solid = result_sewed_blade + + def _build_cylinder(self, radius): + """ + Private method that builds the cylinder which intersects the blade + at a specific cylindrical section. + Argument 'radius' is the radius value corresponding to the cylindrical + section taken into account. This method is applied to all the radii in + list 'radii_list' given as input. + """ + # Base point such that cylinder intersect all section + bbox = Bnd_Box() + brepbndlib.Add(self.blade_solid, bbox) + xmin, _, _, xmax, _, _ = bbox.Get() + + axis = gp_Ax2(gp_Pnt(xmin-0.2*abs(xmin), 0.0, 0.0), + gp_Dir(1.0, 0.0, 0.0)) + self.cylinder = BRepPrimAPI_MakeCylinder(axis, radius, 1.2*(xmax-xmin)).Shape() + faceCount = 0 + faces_explorer = TopExp_Explorer(self.cylinder, TopAbs_FACE) + + while faces_explorer.More(): + face = OCC.Core.TopoDS.topods.Face(faces_explorer.Current()) + faceCount += 1 + # Convert the cylinder faces into Non Uniform Rational Basis-Splines geometry (NURBS) + nurbs_converter = BRepBuilderAPI_NurbsConvert(face) + nurbs_converter.Perform(face) + nurbs_face_converter = nurbs_converter.Shape() + nurbs_face = OCC.Core.TopoDS.topods.Face(nurbs_face_converter) + + surface = BRep_Tool.Surface(OCC.Core.TopoDS.topods.Face(nurbs_face)) + self.bounds = 0.0 + self.bounds = surface.Bounds() + + # Compute the normal curve to the surfaces, specifying the bounds considered, + # the maximum order of the derivative we want to compute, the linear tolerance + normal = GeomLProp_SLProps(surface, + (self.bounds[0] + self.bounds[1]) / 2.0, + (self.bounds[2] + self.bounds[3]) / 2.0, + 1, self.linear_tolerance).Normal() + if (normal * axis.Direction() == 0): + self.cylinder_lateral_face = face + self.cylinder_lateral_surface = surface + self.cylinder_lateral_faces_list.append(surface) + faces_explorer.Next() + + def _extract_parameters_and_transform_profile(self, radius): + """ + Private method which extracts the parameters (pitch, rake, skew ,...) related + to a specific section of the blade, and then transforms the camber points + initialized in _initial_camber_points_plane and the prile points initialized + in _initial_airfoil_points_plane according to the parameters found. + This method trasforms the leading and trailing edges initialed in + _initial_leading_trailing_edges and the mid chord point as well. + + Basically all points and edges of the profile section are transformed + to fulfill the properties of the specific sectionwe are considering. + """ + self.pitch_angle = np.arctan2( + self.leading_edge_point_on_plane[1] - + self.trailing_edge_point_on_plane[1], + self.leading_edge_point_on_plane[0] - + self.trailing_edge_point_on_plane[0]) + self.skew = self.mid_chord_point_on_plane[1] + self.rake_induced_by_skew = self.skew / np.tan(self.pitch_angle) + + self.airfoil_points_on_plane[:, 1:2] -= self.skew + self.camber_points_on_plane[:, 1:2] -= self.skew + self.leading_edge_point_on_plane[1] -= self.skew + self.trailing_edge_point_on_plane[1] -= self.skew + self.mid_chord_point_on_plane[1] -= self.skew + + + self.airfoil_points_on_plane[:, 0:1] -= self.rake_induced_by_skew + self.camber_points_on_plane[:, 0:1] -= self.rake_induced_by_skew + self.leading_edge_point_on_plane[0] -= self.rake_induced_by_skew + self.trailing_edge_point_on_plane[0] -= self.rake_induced_by_skew + self.mid_chord_point_on_plane[0] -= self.rake_induced_by_skew + + + self.rake = self.mid_chord_point_on_plane[0] + + self.airfoil_points_on_plane[:, 0:1] -= self.rake + self.camber_points_on_plane[:, 0:1] -= self.rake + self.leading_edge_point_on_plane[0] -= self.rake + self.trailing_edge_point_on_plane[0] -= self.rake + self.mid_chord_point_on_plane[0] -= self.rake + + self.rotation_matrix = np.zeros((2, 2)) + self.rotation_matrix[0][0] = np.cos(-self.pitch_angle) + self.rotation_matrix[0][1] = -np.sin(-self.pitch_angle) + self.rotation_matrix[1][0] = np.sin(-self.pitch_angle) + self.rotation_matrix[1][1] = np.cos(-self.pitch_angle) + + self.airfoil_points_on_plane = self.airfoil_points_on_plane.dot( + self.rotation_matrix.transpose()) + self.camber_points_on_plane = self.camber_points_on_plane.dot( + self.rotation_matrix.transpose()) + self.leading_edge_point_on_plane = self.leading_edge_point_on_plane.dot( + self.rotation_matrix.transpose()) + self.trailing_edge_point_on_plane = self.trailing_edge_point_on_plane.dot( + self.rotation_matrix.transpose()) + self.mid_chord_point_on_plane = self.mid_chord_point_on_plane.dot( + self.rotation_matrix.transpose()) + + + self.chord_length = ((self.leading_edge_point_on_plane[0] - + self.trailing_edge_point_on_plane[0])**2 + + (self.leading_edge_point_on_plane[1] - + self.trailing_edge_point_on_plane[1])**2)**.5 + + self.airfoil_points_on_plane[:, 0:2] /= self.chord_length + self.camber_points_on_plane[:, 0:2] /= self.chord_length + self.leading_edge_point_on_plane[0:1] /= self.chord_length + self.trailing_edge_point_on_plane[0:1] /= self.chord_length + self.mid_chord_point_on_plane[0:1] /= self.chord_length + + self.airfoil_points_on_plane[:, 0:1] *= -1.0 + self.camber_points_on_plane[:, 0:1] *= -1.0 + self.leading_edge_point_on_plane[0] *= -1.0 + self.trailing_edge_point_on_plane[0] *= -1.0 + self.mid_chord_point_on_plane[0] *= -1.0 + + self.airfoil_points_on_plane[:, 0:1] += 0.5 + self.camber_points_on_plane[:, 0:1] += 0.5 + self.leading_edge_point_on_plane[0] += 0.5 + self.trailing_edge_point_on_plane[0] += 0.5 + self.mid_chord_point_on_plane[0] += 0.5 + + self.camber_points_on_plane = np.matrix(self.camber_points_on_plane) + self.camber_points_on_plane[0, 0] = 0.0 + self.camber_points_on_plane[-1, 0] = 1.0 + + # Save the properties we wanted for each section + self.pitch_angles_list.append(self.pitch_angle) + self.pitch_list.append( + abs(2 * np.pi * radius / np.tan(self.pitch_angle)) / 1000.0) + self.skew_angles_list.append((self.skew / radius - np.pi) * 180 / np.pi) + self.skew_list.append((self.skew - np.pi * radius) / 1000.0) + total_rake = self.rake + self.rake_induced_by_skew + rake = total_rake - (self.skew - np.pi * radius) / np.tan( + self.pitch_angle) + rake = -rake / 1000.0 + self.rake_list.append(rake) + self.chord_length_list.append(self.chord_length / 1000.0) + self.camber_points_on_plane = np.sort( + self.camber_points_on_plane.view('float64,float64'), + order=['f0'], + axis=0).view(np.float64) + self.f_camber = interp1d( + np.squeeze(np.asarray(self.camber_points_on_plane[:, 0])), + np.squeeze(np.asarray(self.camber_points_on_plane[:, 1])), + kind='cubic') + + def save_global_parameters(self, filename_csv): + """ + Method that writes all blade properties and points of the sections + in a csv file named filename_csv + """ + + # open the file in the write mode + with open(filename_csv, 'w') as f: + # create the csv writer + writer = csv.writer(f) + # write a row to the csv file + writer.writerow(("Radii list: ", self.radii_list)) + writer.writerow(("Pitch angles list: ", self.pitch_angles_list)) + writer.writerow(("Pitch list: ", self.pitch_list)) + writer.writerow(("Skew angles list: ", self.skew_angles_list)) + writer.writerow(("Rake list: ", self.rake_list)) + writer.writerow(("Chord length list: ", self.chord_length_list)) + for j in range(self.num_sections): + writer.writerow(("x section" + str(j) + " upper coordinates: ", + self.xup[j, :])) + writer.writerow(("x section" + str(j) + " lower coordinates: ", + self.xdown[j, :])) + writer.writerow(("y section" + str(j) + " upper coordinates: ", + self.yup[j, :])) + writer.writerow(("y section" + str(j) + " lower coordinates: ", + self.ydown[j, :])) + + def reconstruct_sections(self): + """ + Method that reconstructs single sections of the blade starting from the + points computed from the original IGES file. + If sections are not enough to reconstruct the blade, one can just export sections + and then reconstruct the blade in a dedicated script. + """ + + for j in range(self.num_sections): + self.recons_sections[j] = CustomProfile( + xup=self.xup[j, :], + yup=self.yup[j, :], + xdown=self.xdown[j, :], + ydown=self.ydown[j, :]) + return self.recons_sections + + def reconstruct_blade(self): + """ + Method that reconstructs the blade starting from the sections + computed from the original IGES file, and then reconstruct the 4 parts + of the blade (upper, lower face, tip, root). + """ + + for j in range(self.num_sections): + self.recons_sections[j] = CustomProfile( + xup=self.xup[j, :], + yup=self.yup[j, :], + xdown=self.xdown[j, :], + ydown=self.ydown[j, :]) + radii = np.array(self.radii_list) / 1000 + self.recons_blade = Blade(sections=self.recons_sections, + radii=radii, + chord_lengths=self.chord_length_list, + pitch=self.pitch_list, + rake=self.rake_list, + skew_angles=self.skew_angles_list) + self.recons_blade.apply_transformations() + self.recons_blade.generate_iges_blade('iges_reconstructed.iges') + + def reconstruct_propeller(self, propeller_iges, shaft_iges, n_blades): + """ + Method which reconstructs the whole propeller with shaft starting from + the iges file of the shaft and the number of blades. The whole propeller + is saved in an IGES file named filename_iges (must be given as input). + """ + shaft_path = shaft_iges + prop_shaft = Shaft(shaft_path) + prop = Propeller(prop_shaft, self.recons_blade, n_blades) + prop.generate_iges(propeller_iges) + prop.display() + + def display_cylinder(self, radius): + """ + Method that displays the cylinder and the blade intersecting each others. + If the radius is not in self.radii_list, error is raised + """ + if radius not in self.radii_list: + raise ValueError("The radius must be among the ones passed in the constructor of the object. The unit is mm.") + index_radius = np.where(self.radii_list == radius)[0][0] + display, start_display, add_menu, add_function_to_menu = init_display() + display.DisplayShape(self.blade_solid, update=True) + display.DisplayShape(self.cylinder_lateral_faces_list[index_radius], update=True) + display.DisplayShape(self.section_wires_list[index_radius], update=True) + start_display() + + def display_section(self, radius): + """ + Method that displays a section profile over the cylinder with its camber line. + """ + if radius not in self.radii_list: + raise ValueError("The radius must be among the ones passed in the constructor of the object. The unit is mm.") + index_radius = np.where(self.radii_list == radius)[0][0] + display, start_display, add_menu, add_function_to_menu = init_display() + display.DisplayShape(self.cylinder_lateral_faces_list[index_radius], update=True) + display.DisplayShape(self.section_wires_list[index_radius], update=True) + display.DisplayShape(self.camber_wires_list[index_radius], color="GREEN", update=True) + start_display() diff --git a/bladex/reversepropeller.py b/bladex/reversepropeller/reversepropeller.py similarity index 67% rename from bladex/reversepropeller.py rename to bladex/reversepropeller/reversepropeller.py index 3579f3b..1404f3b 100644 --- a/bladex/reversepropeller.py +++ b/bladex/reversepropeller/reversepropeller.py @@ -3,11 +3,9 @@ approximated reconstruction of the blade and the whole propeller. """ -import os, errno -import os.path import numpy as np -import csv from bladex import CustomProfile, Blade, Propeller, Shaft +from .basereversepropeller import BaseReversePropeller from OCC.Core.IGESControl import (IGESControl_Reader) from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_MakeVertex,\ BRepBuilderAPI_MakeEdge, BRepBuilderAPI_MakeWire, BRepBuilderAPI_MakeSolid, \ @@ -102,34 +100,8 @@ def point_inside_polygon(x, y, poly, include_edges=True): return inside -class ReversePropeller(object): +class ReversePropeller(BaseReversePropeller): """ - Extraction of the parameters of a blade and reconstruction of the parametrized - blade and of the entire propeller. - - Given the IGES file of the blade, the following parameters are extracted: - - - :math:`(X, Y)` coordinates of the blade cylindrical sections after - being expanded in 2D to create airfoils. - - - Chord length, for each cylindrical section. - - - Pitch :math:`(P)`, for each cylindrical section. - - - Rake :math:`(k)`, in distance units, for each cylindrical section. - - - Skew angle :math:`(\\theta_s)`, for each cylindrical section, expressed - in degrees. - - The parameters can be saved in a csv file with the method - 'save_global_parameters' and used to reconstruct the blade, which is saved - in a IGES file (method 'reconstruct_blade'). - The coordinates of each section can be extracted with the method - 'reconstruct_sections'. - Given the shaft, the whole propeller can be reconstructed and saved in a - IGES file with the method 'reconstruct_propeller'. - - -------------------------- :param filename: path to the IGES file of the blade. :param list radii_list: list which contains the radii values of the @@ -139,36 +111,7 @@ class ReversePropeller(object): """ def __init__(self, filename, radii_list, num_points_top_bottom): - self.iges_file = filename #filename is the path to the file iges - self.radii_list = radii_list #radii at which we want to measure properties - self.section_points_number = 1200 #number of points for each section profile - self.airfoil_points_number = 1000 - self.num_points_top_bottom = num_points_top_bottom #number of points used to reconstruct the profile, equal for top and lower parts of the profile - self.coords = [] - self.num_sections = len(self.radii_list) - self.recons_sections = [0] * self.num_sections - self.start = None - self.x = [] - self.y = [] - self.poly = [] - # Initialize things we want to compute from the blade IGES file - self.pitch_angles_list = [] - self.pitch_list = [] - self.skew_angles_list = [] - self.skew_list = [] - self.rake_list = [] - self.chord_length_list = [] - self.blade_solid = None - self.blade_compound = None - self.tolerance_solid = 5e-2 - self._extract_solid_from_file() - self.xup = np.zeros((self.num_sections, self.num_points_top_bottom)) - self.xdown = np.zeros((self.num_sections, self.num_points_top_bottom)) - self.yup = np.zeros((self.num_sections, self.num_points_top_bottom)) - self.ydown = np.zeros((self.num_sections, self.num_points_top_bottom)) - # At each radius, i.e., for each section, the cylinder with that radius is built - # and the faces of the cylinder are isolated. - + super().__init__(filename, radii_list, num_points_top_bottom) ind_sec = 0 for radius in self.radii_list: self.cylinder = None @@ -178,7 +121,7 @@ def __init__(self, filename, radii_list, num_points_top_bottom): self.linear_tolerance = 1e-3 self.conversion_unit = 1000 self._build_cylinder(radius) - self._build_intersection_cylinder_blade(radius) + self._build_intersection_cylinder_blade() self.leading_edge_point = [] self.trailing_edge_point = [] self.edge_3d = None @@ -194,9 +137,7 @@ def __init__(self, filename, radii_list, num_points_top_bottom): # at a equal relative distance (taken on the curvilinear ascissa) self.param_plane_points = [] self.orig_param_plane_points = [] - for i in range(self.section_points_number): - self.i = i - self._camber_points(radius) + self._camber_points(radius) self._voronoi_points(radius) self._camber_curve(radius) self._initial_leading_trailing_edges_plane(radius) @@ -210,80 +151,8 @@ def __init__(self, filename, radii_list, num_points_top_bottom): self.ydown[ind_sec, :] = self.int_airfoil_bottom ind_sec = ind_sec + 1 - def _extract_solid_from_file(self): - """ - Private method that reads the IGES file of the blade and constructs - the correspondent solid object - """ - # Extract the blade from the IGES file; build it as a solid with OCC - iges_reader = IGESControl_Reader() - iges_reader.ReadFile(self.iges_file) - iges_reader.TransferRoots() - sewer = BRepBuilderAPI_Sewing(self.tolerance_solid) - # Case where we have Faces and not closed surface. - # This is the case for BladeX generated blade - if iges_reader.Shape().ShapeType() == 4: - self.blade_compound = iges_reader.OneShape() - exp = TopExp_Explorer(self.blade_compound, TopAbs_FACE) - while exp.More(): - sewer.Add(exp.Current()) - exp.Next() - else: - self.blade_compound = iges_reader.Shape() - sewer.Add(self.blade_compound) - sewer.Perform() - result_sewed_blade = sewer.SewedShape() - blade_solid_maker = BRepBuilderAPI_MakeSolid() - blade_solid_maker.Add(OCC.Core.TopoDS.topods.Shell(result_sewed_blade)) - if (blade_solid_maker.IsDone()): - self.blade_solid = blade_solid_maker.Solid() - else: - self.blade_solid = result_sewed_blade - - def _build_cylinder(self, radius): - """ - Private method that builds the cylinder which intersects the blade - at a specific cylindrical section. - Argument 'radius' is the radius value corresponding to the cylindrical - section taken into account. This method is applied to all the radii in - list 'radii_list' given as input. - """ - # Base point such that cylinder intersect all section - bbox = Bnd_Box() - brepbndlib.Add(self.blade_solid, bbox) - xmin, ymin, zmin, xmax, ymax, zmax = bbox.Get() - - axis = gp_Ax2(gp_Pnt(xmin-0.2*abs(xmin), 0.0, 0.0), - gp_Dir(1.0, 0.0, 0.0)) - self.cylinder = BRepPrimAPI_MakeCylinder(axis, radius, 1.2*(xmax-xmin)).Shape() - faceCount = 0 - faces_explorer = TopExp_Explorer(self.cylinder, TopAbs_FACE) - - while faces_explorer.More(): - face = OCC.Core.TopoDS.topods.Face(faces_explorer.Current()) - faceCount += 1 - # Convert the cylinder faces into Non Uniform Rational Basis-Splines geometry (NURBS) - nurbs_converter = BRepBuilderAPI_NurbsConvert(face) - nurbs_converter.Perform(face) - nurbs_face_converter = nurbs_converter.Shape() - nurbs_face = OCC.Core.TopoDS.topods.Face(nurbs_face_converter) - - surface = BRep_Tool.Surface(OCC.Core.TopoDS.topods.Face(nurbs_face)) - self.bounds = 0.0 - self.bounds = surface.Bounds() - - # Compute the normal curve to the surfaces, specifying the bounds considered, - # the maximum order of the derivative we want to compute, the linear tolerance - normal = GeomLProp_SLProps(surface, - (self.bounds[0] + self.bounds[1]) / 2.0, - (self.bounds[2] + self.bounds[3]) / 2.0, - 1, self.linear_tolerance).Normal() - if (normal * axis.Direction() == 0): - self.cylinder_lateral_face = face - self.cylinder_lateral_surface = surface - faces_explorer.Next() - - def _build_intersection_cylinder_blade(self, radius): + + def _build_intersection_cylinder_blade(self): """ Private method that constructs the section lines which are the intersections between the cylinder at a fixed radius and the blade, and the camber points. @@ -303,40 +172,42 @@ def _build_intersection_cylinder_blade(self, radius): self.total_section_length = 0.0 edgeList = TopTools_ListOfShape() while edgeExplorer.More(): - edgeCount = edgeCount + 1 + edgeCount = edgeCount + 1 # Numbering from 1 in OCC edge = topods.Edge(edgeExplorer.Current()) edgeList.Append(edge) edgeExplorer.Next() wire_maker.Add(edgeList) self.wire = wire_maker.Wire() + self.section_wires_list.append(self.wire) self.curve_adaptor = BRepAdaptor_CompCurve( OCC.Core.TopoDS.topods.Wire(self.wire)) # Length of the curve section (ascissa curvilinea) self.total_section_length = GCPnts_AbscissaPoint.Length( self.curve_adaptor) + def _camber_points(self, radius): """ Private method which computes the single points of the camber curve and collects them in param_plane_points and in orig_plane_points. - """ firstParam = self.curve_adaptor.FirstParameter() - rel_distance = float(self.i) / float( - self.section_points_number) * self.total_section_length - point_generator = GCPnts_AbscissaPoint(1e-7, self.curve_adaptor, - rel_distance, firstParam) - param = point_generator.Parameter() - self.point = self.curve_adaptor.Value(param) - # Compute the angle for polar reference system - theta = np.arctan2(self.point.Y(), self.point.Z()) - if theta < 0: - theta += 2.0 * np.pi - # Param_plane_points is the list of all points parametrized in cylindric coordinates, - # i.e., the coordinate along X (the axis of the cylinder) and the - # radius multiplying the polar angle. - self.orig_param_plane_points.append([self.point.X(), radius * theta]) - self.param_plane_points.append([self.point.X(), radius * theta]) + for i in range(self.section_points_number): + rel_distance = float(i) / float( + self.section_points_number) * self.total_section_length + point_generator = GCPnts_AbscissaPoint(1e-7, self.curve_adaptor, + rel_distance, firstParam) + param = point_generator.Parameter() + point = self.curve_adaptor.Value(param) + # Compute the angle for polar reference system + theta = np.arctan2(point.Y(), point.Z()) + if theta < 0: + theta += 2.0 * np.pi + # Param_plane_points is the list of all points parametrized in cylindric coordinates, + # i.e., the coordinate along X (the axis of the cylinder) and the + # radius multiplying the polar angle. + self.orig_param_plane_points.append([point.X(), radius * theta]) + self.param_plane_points.append([point.X(), radius * theta]) def _voronoi_points(self, radius): """ @@ -392,10 +263,14 @@ def _camber_curve(self, radius): # with tolerance 1e-1 spline_builder = GeomAPI_PointsToBSpline(Pnts, 3, 3, GeomAbs_C2, 1e-1) camber_curve = spline_builder.Curve() - # BUild the camber line from 3D edges + # Build the camber line from 3D edges self.edge_3d = BRepBuilderAPI_MakeEdge(camber_curve).Edge() camber_curve_adaptor = BRepAdaptor_Curve( OCC.Core.TopoDS.topods.Edge(self.edge_3d)) + # Converting camber edge into wire for plotting purposes + wire_maker = BRepBuilderAPI_MakeWire() + wire_maker.Add(self.edge_3d) + self.camber_wires_list.append(wire_maker.Wire()) # Compute the length of the camber curve for each section camber_curve_length = GCPnts_AbscissaPoint.Length(camber_curve_adaptor) relative_tolerance = 2.5e-3 @@ -756,120 +631,3 @@ def _airfoil_top_and_bottom_points(self): self.int_airfoil_bottom[0] = 0.0 self.int_airfoil_bottom[-1] = 0.0 - def save_global_parameters(self, filename_csv): - """ - Method that writes all blade properties and points of the sections - in a csv file named filename_csv - """ - - # open the file in the write mode - with open(filename_csv, 'w') as f: - # create the csv writer - writer = csv.writer(f) - # write a row to the csv file - writer.writerow(("Radii list: ", self.radii_list)) - writer.writerow(("Pitch angles list: ", self.pitch_angles_list)) - writer.writerow(("Pitch list: ", self.pitch_list)) - writer.writerow(("Skew angles list: ", self.skew_angles_list)) - writer.writerow(("Rake list: ", self.rake_list)) - writer.writerow(("Chord length list: ", self.chord_length_list)) - for j in range(self.num_sections): - writer.writerow(("x section" + str(j) + " upper coordinates: ", - self.xup[j, :])) - writer.writerow(("x section" + str(j) + " lower coordinates: ", - self.xdown[j, :])) - writer.writerow(("y section" + str(j) + " upper coordinates: ", - self.yup[j, :])) - writer.writerow(("y section" + str(j) + " lower coordinates: ", - self.ydown[j, :])) - - def reconstruct_sections(self): - """ - Method that reconstructs single sections of the blade starting from the - points computed from the original IGES file. - If sections are not enough to reconstruct the blade, one can just export sections - and then reconstruct the blade in a dedicated script. - """ - - for j in range(self.num_sections): - self.recons_sections[j] = CustomProfile( - xup=self.xup[j, :], - yup=self.yup[j, :], - xdown=self.xdown[j, :], - ydown=self.ydown[j, :]) - return self.recons_sections - - def reconstruct_blade(self): - """ - Method that reconstructs the blade starting from the sections - computed from the original IGES file, and then reconstruct the 4 parts - of the blade (upper, lower face, tip, root). - """ - - for j in range(self.num_sections): - self.recons_sections[j] = CustomProfile( - xup=self.xup[j, :], - yup=self.yup[j, :], - xdown=self.xdown[j, :], - ydown=self.ydown[j, :]) - radii = np.array(self.radii_list) / 1000 - self.recons_blade = Blade(sections=self.recons_sections, - radii=radii, - chord_lengths=self.chord_length_list, - pitch=self.pitch_list, - rake=self.rake_list, - skew_angles=self.skew_angles_list) - self.recons_blade.apply_transformations() - self.recons_blade.generate_iges_blade('iges_reconstructed.iges') - - def reconstruct_propeller(self, propeller_iges, shaft_iges, n_blades): - """ - Method which reconstructs the whole propeller with shaft starting from - the iges file of the shaft and the number of blades. The whole propeller - is saved in an IGES file named filename_iges (must be given as input). - """ - shaft_path = shaft_iges - prop_shaft = Shaft(shaft_path) - prop = Propeller(prop_shaft, self.recons_blade, n_blades) - prop.generate_iges(propeller_iges) - prop.display() - - def display_cylinders(self): - """ - Method that displays the cylinders and the blade, for each section. - """ - display, start_display, add_menu, add_function_to_menu = init_display() - for radius in self.radii_list: - # Display the blade and the wire corresponding to the intersection with the - # lateral face of the cylinder computed for each section radius - display.DisplayShape(self.blade_solid, update=True) - display.DisplayShape(self.wire, update=True) - - def display_sections(self): - """ - Method that displays the sections profiles. - """ - display, start_display, add_menu, add_function_to_menu = init_display() - for radius in self.radii_list: - - display.DisplayShape(BRepBuilderAPI_MakeVertex( - self.leading_edge_point).Vertex(), - update=True, - color="ORANGE") - display.DisplayShape(BRepBuilderAPI_MakeVertex( - self.trailing_edge_point).Vertex(), - update=True, - color="GREEN") - # Display the camber points, leading and trailing edges projected on plane - display.DisplayShape(self.camber_curve_edge, - update=True, - color="GREEN") - display.DisplayShape(self.edge_3d, update=True, color="RED") - display.DisplayShape(self.first_edge, update=True, color="BLUE1") - display.DisplayShape(self.last_edge, update=True, color="BLUE1") - display.DisplayShape(self.firstSegment, update=True, color="ORANGE") - display.DisplayShape(self.lastSegment, update=True, color="ORANGE") - # Update the blade, cylinder and section when considering a different radius - display.DisplayShape(self.blade_compound, update=True) - display.DisplayShape(self.cylinder_lateral_face, update=True) - display.DisplayShape(self.section, update=True) diff --git a/bladex/reversepropeller/reversepropellerbladex.py b/bladex/reversepropeller/reversepropellerbladex.py new file mode 100644 index 0000000..a11f17d --- /dev/null +++ b/bladex/reversepropeller/reversepropellerbladex.py @@ -0,0 +1,436 @@ +""" +Module for the extraction of the parameters of a blade and for the +approximated reconstruction of the blade and the whole propeller. +""" + +import numpy as np +from bladex import CustomProfile, Blade, Propeller, Shaft +from .basereversepropeller import BaseReversePropeller +from OCC.Core.IGESControl import (IGESControl_Reader) +from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_MakeVertex,\ + BRepBuilderAPI_MakeEdge, BRepBuilderAPI_MakeWire, BRepBuilderAPI_MakeSolid, \ + BRepBuilderAPI_Sewing, BRepBuilderAPI_NurbsConvert +from OCC.Core.BRep import BRep_Tool +import OCC.Core.TopoDS +from OCC.Core.TopTools import TopTools_ListOfShape +from OCC.Core.TopExp import TopExp_Explorer +from OCC.Core.TopAbs import TopAbs_VERTEX, TopAbs_EDGE, TopAbs_FACE, TopAbs_WIRE +from OCC.Core.gp import gp_Pnt, gp_Dir, gp_Ax2, gp_Vec, gp_Pln +from OCC.Core.TColgp import TColgp_Array1OfPnt +from OCC.Core.BRepAdaptor import BRepAdaptor_Curve +from OCC.Core.GCPnts import GCPnts_AbscissaPoint +from OCC.Core.BRep import BRep_Tool +from OCC.Core.TopoDS import topods +from OCC.Display.SimpleGui import init_display +from OCC.Core.BRepPrimAPI import BRepPrimAPI_MakeCylinder +from OCC.Core.GeomLProp import GeomLProp_SLProps +from OCC.Core.GCPnts import GCPnts_AbscissaPoint +from OCC.Core.BRepAdaptor import (BRepAdaptor_Curve, + BRepAdaptor_CompCurve) +from OCC.Core.GeomAPI import GeomAPI_PointsToBSpline +from OCC.Core.GeomAbs import GeomAbs_C2 +from OCC.Core.Bnd import Bnd_Box +from OCC.Core.BRepBndLib import brepbndlib +from OCC.Core.BRepAlgoAPI import BRepAlgoAPI_Section +from OCC.Core.Geom import Geom_Plane +from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_MakeFace +from scipy.interpolate import interp1d +import matplotlib.pyplot as plt + + +class ReversePropellerBladeX(BaseReversePropeller): + """ + + :param filename: path to the IGES file of the blade. + :param list radii_list: list which contains the radii values of the + sectional profiles. + :param num_points_top_bottom: number of points used to interpolate each + sectional profile. + """ + def __init__(self, filename, radii_list, num_points_top_bottom): + + super().__init__(filename, radii_list, num_points_top_bottom) + ind_sec = 0 + for radius in self.radii_list: + self.cylinder = None + self.bounds = None + self.cylinder_lateral_face = None + self.cylinder_lateral_surface = None + self.linear_tolerance = 1e-3 + self.conversion_unit = 1000 + self._build_cylinder(radius) + self._build_intersection_cylinder_blade() + self.leading_edge_point = [] + self.trailing_edge_point = [] + self.edge_3d = None + self.camber_curve_edge = None + self.first_edge = None + self.last_edge = None + self.firstSegment = None + self.lastSegment = None + self.camber_points_on_plane = None + self.leading_edge_point_on_plane = [] + self.trailing_edge_point_on_plane = [] + self.param_plane_points = [] + self.orig_param_plane_points = [] + # points = [] + self._camber_curve(radius) + # for i in range(len(self.vor_us)): + # x, y, z = self.vor_us[i], radius*np.sin(self.vor_vs[i]/radius), radius*np.cos(self.vor_vs[i]/radius) + # pnt = gp_Pnt(float(x), float(y), float(z)) + # points.append(BRepBuilderAPI_MakeVertex(pnt).Vertex()) + # self._camber_curve(radius) + if False: + display, start_display, add_menu, add_function_to_menu = init_display() + # display.DisplayShape(self.blade_solid, update=True) + display.DisplayShape(self.cylinder, update=True) + display.DisplayShape(self.wire_top, color="RED", update=True) + display.DisplayShape(self.wire_bottom, color="BLUE", update=True) + # display.DisplayShape(self.chord_wire, color="ORANGE", update=True) + display.DisplayShape(self.full_camber_wire, color="GREEN", update=True) + display.DisplayShape(self.leading_edge, color="BLACK", update=True) + display.DisplayShape(self.trailing_edge, color="MAGENTA", update=True) + # display.DisplayShape(self.chord_plane, update=True) + display.DisplayShape(self.chord_point, update=True) + display.DisplayShape(self.trimmed_edge, color="BLACK", update=True) + # display.DisplayShape(points, update=True) + start_display() + self._initial_leading_trailing_edges_plane(radius) + self._initial_camber_points_plane(radius) + self._initial_airfoil_points_plane(radius) + self._airfoil_top_and_bottom_points_before_transformations(radius) + self._extract_parameters_and_transform_profile(radius) + self._transform_top_and_bottom() + self._airfoil_top_and_bottom_points_after_transformations() + self.xup[ind_sec, :] = self.ascissa + self.xdown[ind_sec, :] = self.ascissa + self.yup[ind_sec, :] = self.int_airfoil_top + self.ydown[ind_sec, :] = self.int_airfoil_bottom + ind_sec = ind_sec + 1 + + + def _build_intersection_cylinder_blade(self): + """ + Private method that constructs the section lines which are the intersections + between the cylinder at a fixed radius and the blade, and the camber points. + """ + # Construction of the section lines between two shapes (in this case the + # blade and the lateral face of the cylinder) + section_builder = BRepAlgoAPI_Section(self.blade_solid, + self.cylinder_lateral_face, False) + # Define and build the parametric 2D curve (pcurve) for the section lines defined above + section_builder.ComputePCurveOn2(True) + section_builder.Build() + self.section = section_builder.Shape() + + wire_maker = BRepBuilderAPI_MakeWire() + wire_maker_top = BRepBuilderAPI_MakeWire() + wire_maker_bottom = BRepBuilderAPI_MakeWire() + + edgeList = TopTools_ListOfShape() + edgeList_top = TopTools_ListOfShape() + edgeList_bottom = TopTools_ListOfShape() + edgeExplorer = TopExp_Explorer(self.section, TopAbs_EDGE) + edgeCount = 0 + while edgeExplorer.More(): + edgeCount = edgeCount + 1 # Numbering from 1 in OCC + edge = edgeExplorer.Current() + edgeList.Append(edge) + if edgeCount % 2 == 1: + edgeList_top.Append(edge) + else: + edgeList_bottom.Append(edge) + edgeExplorer.Next() + + # Total sectional curve + wire_maker.Add(edgeList) + self.wire = wire_maker.Wire() + self.section_wires_list.append(self.wire) + self.curve_adaptor = BRepAdaptor_CompCurve( + OCC.Core.TopoDS.topods.Wire(self.wire)) + self.total_section_length = GCPnts_AbscissaPoint.Length( + self.curve_adaptor) + # Top part + wire_maker_top.Add(edgeList_top) + self.wire_top = wire_maker_top.Wire() + self.curve_adaptor_top = BRepAdaptor_CompCurve( + OCC.Core.TopoDS.topods.Wire(self.wire_top)) + self.total_section_top_length = GCPnts_AbscissaPoint.Length( + self.curve_adaptor_top) + # Bottom part + wire_maker_bottom.Add(edgeList_bottom) + self.wire_bottom = wire_maker_bottom.Wire() + self.curve_adaptor_bottom = BRepAdaptor_CompCurve( + OCC.Core.TopoDS.topods.Wire(self.wire_bottom)) + self.total_section_bottom_length = GCPnts_AbscissaPoint.Length( + self.curve_adaptor_bottom) + + + def _camber_curve(self, radius): + """ + Computation of the camber points. We get the chord, move along it and fint the intersection + of the othogonal plane to the chord-curvilinear absissa and the top-bottom wires + """ + # self.trailing_edge = self.curve_adaptor_top.Value(0.0) + # self.leading_edge = self.curve_adaptor_top.Value(self.curve_adaptor_top.LastParameter()) + self.trailing_edge = self.curve_adaptor_bottom.Value(0.0) + self.leading_edge = self.curve_adaptor_bottom.Value(self.curve_adaptor_bottom.LastParameter()) + # Equation of generic geodesic from leading to trailing edge: + # x = lambda*s + mu, y = R*sin(a*s+b), z = R*cos(a*s+b) + # with s parametrization of curve and lambda, mu, a, b to be found + # for s=0 we habe leading edge and s=1 we have trailing edge + mu = self.leading_edge.X() + lam = self.trailing_edge.X() - mu + b = np.arcsin(self.leading_edge.Y()/radius) + a = np.arcsin(self.trailing_edge.Y()/radius) - b + # Thus, the new a chord point is, given any s \in [0,1] + # gp_Pnt(lam*s + mu, radius * np.sin(a*s+b), radius * np.cos(a*s+b))) + + # Now, for N chord points, we compute the tangent vector to it, + # define the plane passing to the point and with normal direction equal to the tangent + # and compute the intersecting curve between cylinder and plane. The midpoint is the camber point + self.n_camber_points = 500 + n_plot_point = int (self.n_camber_points/2) + s_values = np.linspace(0, 1, self.n_camber_points) + self.camber_Pnts = TColgp_Array1OfPnt(1, len(s_values)) + for i, s in enumerate(s_values): + if i == 0: + self.camber_Pnts.SetValue(i+1, self.leading_edge) + elif i == len(s_values)-1: + self.camber_Pnts.SetValue(i+1, self.trailing_edge) + else: + chord_point = gp_Pnt(lam*s + mu, radius * np.sin(a*s+b), radius * np.cos(a*s+b)) + # Equation of tangent vector is: x = lam, y = a*R*cos(a*s+b), z = -a*R*sin(a*s+b) + tangent_vector = gp_Vec(lam, a*radius*np.cos(a*s+b), -a*radius*np.sin(a*s+b)) + chord_plane = gp_Pln(chord_point, gp_Dir(tangent_vector)) + chord_plane_face = BRepBuilderAPI_MakeFace(chord_plane, -1e8, 1e8, -1e8, 1e8).Face() + if i == n_plot_point: + self.chord_plane = chord_plane_face + # self.chord_plane = chord_plane + self.chord_point = chord_point + # Now we intersect the plane with the cylinder + section_pc = BRepAlgoAPI_Section(self.cylinder_lateral_face, chord_plane_face) + section_pc.Build() + curve_shape = section_pc.Shape() + exp = TopExp_Explorer(curve_shape, TopAbs_EDGE) + curve_edge = exp.Current() + # Find the intersection of the curve with top and bottom faces of the blade's section + int_top = BRepAlgoAPI_Section(curve_edge, self.wire_top) + int_top.Build() + int_bottom = BRepAlgoAPI_Section(curve_edge, self.wire_bottom) + int_bottom.Build() + exp = TopExp_Explorer(int_top.Shape(), TopAbs_VERTEX) + from OCC.Core.TopoDS import TopoDS_Vertex + while exp.More(): + point_top = exp.Current() + if isinstance(int_top, TopoDS_Vertex): + break + exp.Next() + exp = TopExp_Explorer(int_bottom.Shape(), TopAbs_VERTEX) + while exp.More(): + point_bottom= exp.Current() + if isinstance(int_bottom, TopoDS_Vertex): + break + exp.Next() + # Now we trim the curve to only have inside the two intersection points + curve_handle, u1, u2 = BRep_Tool.Curve(curve_edge) + from OCC.Core.GeomAPI import GeomAPI_ProjectPointOnCurve + proj = GeomAPI_ProjectPointOnCurve(BRep_Tool.Pnt(point_top), curve_handle) + u_top = proj.LowerDistanceParameter() + proj = GeomAPI_ProjectPointOnCurve(BRep_Tool.Pnt(point_bottom), curve_handle) + u_bot = proj.LowerDistanceParameter() + u_min = min(u_top, u_bot) + u_max = max(u_top, u_bot) + trimmed_edge = BRepBuilderAPI_MakeEdge(curve_handle, u_min, u_max).Edge() + if i == n_plot_point: + wire_maker = BRepBuilderAPI_MakeWire() + wire_maker.Add(trimmed_edge) + self.trimmed_edge= wire_maker.Wire() + trimmed_adaptor = BRepAdaptor_Curve(trimmed_edge) + u0 = trimmed_adaptor.FirstParameter() + u1 = trimmed_adaptor.LastParameter() + total_length = GCPnts_AbscissaPoint.Length(trimmed_adaptor, u0, u1) + absc = GCPnts_AbscissaPoint(trimmed_adaptor, 0.5*total_length, u0) + u_mid = absc.Parameter() + self.camber_Pnts.SetValue(i+1, trimmed_adaptor.Value(u_mid)) + # Building an edge from points + spline_builder = GeomAPI_PointsToBSpline(self.camber_Pnts, 3, 3, GeomAbs_C2, 1e-1) + camber_curve = spline_builder.Curve() + # Build the camber line from 3D edges + self.edge_3d = BRepBuilderAPI_MakeEdge(camber_curve).Edge() + # Building wire from edge + wire_maker = BRepBuilderAPI_MakeWire() + wire_maker.Add(self.edge_3d) + + self.full_camber_wire = wire_maker.Wire() + self.full_camber_curve_adaptor = BRepAdaptor_CompCurve( + OCC.Core.TopoDS.topods.Wire(self.full_camber_wire)) + self.full_camber_length = GCPnts_AbscissaPoint.Length( + self.full_camber_curve_adaptor) + self.camber_wires_list.append(self.full_camber_wire) + + def _initial_leading_trailing_edges_plane(self, radius): + """ + Private method which computes the coordinates of the leading and trailing + edges of each section, as were in the initial blade. Thenm transformations + will be applied in order to plot the points of the profiles on a plane + and to map the X coordinates in the interval [0,1]. + """ + + self.trailing_edge_point = self.trailing_edge + self.leading_edge_point = self.leading_edge + leading_edge_theta = np.arctan2(self.leading_edge_point.Y(), + self.leading_edge_point.Z()) + self.leading_edge_point_on_plane = np.array( + [self.leading_edge_point.X(), radius * leading_edge_theta]) + + trailing_edge_theta = np.arctan2(self.trailing_edge_point.Y(), + self.trailing_edge_point.Z()) + self.trailing_edge_point_on_plane = np.array( + [self.trailing_edge_point.X(), radius * trailing_edge_theta]) + + def _initial_camber_points_plane(self, radius): + """ + Private method that defines the points of the camber line projected on + plane, in 2D coordinates. + """ + self.camber_points_on_plane = np.zeros((self.num_points_top_bottom, 2)) + + for i in range(self.num_points_top_bottom): + point = self.camber_Pnts.Value(i+1) + theta = np.arctan2(point.Y(), point.Z()) + if theta < 0: + theta += 2.0 * np.pi + self.camber_points_on_plane[i][0] = point.X() + self.camber_points_on_plane[i][1] = radius * theta + + def _initial_airfoil_points_plane(self, radius): + """ + Private method that evaluates the airfoil points on a plane (the points of + the profile of a single section, which will be then distinguished into + upper and lower part). Those points are defined in a plane, in 2D coordinates. + """ + + self.airfoil_points_on_plane = np.zeros((self.airfoil_points_number, 2)) + + for i in range(self.airfoil_points_number): + firstParam = self.curve_adaptor.FirstParameter() + rel_distance = float(i) / float(self.airfoil_points_number - + 1) * self.total_section_length + param = GCPnts_AbscissaPoint(1e-7, self.curve_adaptor, rel_distance, + firstParam).Parameter() + point = self.curve_adaptor.Value(param) + theta = np.arctan2(point.Y(), point.Z()) + if theta < 0: + theta += 2.0 * np.pi + + self.airfoil_points_on_plane[i][0] = point.X() + self.airfoil_points_on_plane[i][1] = radius * theta + + # Compute the mid point of the chord line + self.mid_chord_point_on_plane = ( + self.leading_edge_point_on_plane + + self.trailing_edge_point_on_plane) / 2.0 + + def _airfoil_top_and_bottom_points_before_transformations(self, radius): + """ + Private method that finds the points of the airfoil belonging to the upper + and lower profiles of each section. + """ + self.airfoil_top = [] + u0_top = self.curve_adaptor_top.FirstParameter() + u1_top = self.curve_adaptor_top.LastParameter() + total_length = GCPnts_AbscissaPoint.Length(self.curve_adaptor_top, u0_top, u1_top) + for i in range(self.num_points_top_bottom): + point_generator = GCPnts_AbscissaPoint(1e-7, self.curve_adaptor_top, + float(i)/(self.num_points_top_bottom-1)*total_length, u0_top) + param = point_generator.Parameter() + point = self.curve_adaptor_top.Value(param) + theta = np.arctan2(point.Y(), point.Z()) + self.airfoil_top.append([point.X(), radius*theta]) + self.airfoil_bottom = [] + u0_bottom = self.curve_adaptor_bottom.FirstParameter() + u1_bottom = self.curve_adaptor_bottom.LastParameter() + total_length = GCPnts_AbscissaPoint.Length(self.curve_adaptor_bottom, u0_bottom, u1_bottom) + for i in range(self.num_points_top_bottom): + point_generator = GCPnts_AbscissaPoint(1e-7, self.curve_adaptor_bottom, + float(i)/(self.num_points_top_bottom-1)*total_length, u0_bottom) + param = point_generator.Parameter() + point = self.curve_adaptor_bottom.Value(param) + theta = np.arctan2(point.Y(), point.Z()) + self.airfoil_bottom.append([point.X(), radius*theta]) + + self.airfoil_top = np.array(self.airfoil_top) + self.airfoil_bottom = np.array(self.airfoil_bottom) + + def _transform_top_and_bottom(self): + """ + Method that transforms the top and bottom coordinates from physical coordinates + to 2D planar with chord aligned to x-axis + skew, pitch, rake etc. have been already compute in _extract_parameters_and_transform_profile() + """ + self.airfoil_top[:, 1] -= self.skew + self.airfoil_bottom[:, 1] -= self.skew + + self.airfoil_top[:, 0] -= self.rake_induced_by_skew + self.airfoil_bottom[:, 0] -= self.rake_induced_by_skew + + self.airfoil_top[:, 0] -= self.rake + self.airfoil_bottom[:, 0] -= self.rake + + self.airfoil_top = self.airfoil_top.dot( + self.rotation_matrix.transpose()) + self.airfoil_bottom = self.airfoil_bottom.dot( + self.rotation_matrix.transpose()) + + self.airfoil_top /= self.chord_length + self.airfoil_bottom /= self.chord_length + + self.airfoil_top[:, 0] *= -1.0 + self.airfoil_bottom[:, 0] *= -1.0 + + self.airfoil_top[:, 0] += 0.5 + self.airfoil_bottom[:, 0] += 0.5 + + # Final check on y components to see which one is really top + if self.airfoil_top[np.abs(self.airfoil_top[:, 1]).argmax(), 1] < 0: + self.airfoil_top[:, 1] *= -1 + self.airfoil_bottom[:, 1] *= -1 + + def _airfoil_top_and_bottom_points_after_transformations(self): + # Now we create two interpolating functions fitting bottom e up points + # In this way we can retrieve as many bottom/up points as we want + self.airfoil_top = np.matrix(self.airfoil_top) + self.airfoil_top = np.sort(self.airfoil_top.view('float64,float64'), + order=['f0'], + axis=0).view(np.float64) + self.airfoil_bottom = np.matrix(self.airfoil_bottom) + self.airfoil_bottom = np.sort( + self.airfoil_bottom.view('float64,float64'), order=['f0'], + axis=0).view(np.float64) + + self.airfoil_top[0, 0] = 0.0 + self.airfoil_top[-1, 0] = 1.0 + self.airfoil_bottom[0, 0] = 0.0 + self.airfoil_bottom[-1, 0] = 1.0 + + self.ascissa = np.linspace(0, + 1, + num=self.num_points_top_bottom, + endpoint=True) + self.f_top = interp1d(np.squeeze(np.asarray(self.airfoil_top[:, 0])), + np.squeeze(np.asarray(self.airfoil_top[:, 1])), + kind='cubic') + self.f_bottom = interp1d( + np.squeeze(np.asarray(self.airfoil_bottom[:, 0])), + np.squeeze(np.asarray(self.airfoil_bottom[:, 1])), + kind='cubic') + # Plot points of camber line, trailing and leading edges, upper and lower profiles + # of each blade section + self.int_airfoil_top = self.f_top(self.ascissa) + self.int_airfoil_top[0] = 0.0 + self.int_airfoil_top[-1] = 0.0 + self.int_airfoil_bottom = self.f_bottom(self.ascissa) + self.int_airfoil_bottom[0] = 0.0 + self.int_airfoil_bottom[-1] = 0.0 + diff --git a/bladex/reversepropeller/reversepropellerinterface.py b/bladex/reversepropeller/reversepropellerinterface.py new file mode 100644 index 0000000..afd6e19 --- /dev/null +++ b/bladex/reversepropeller/reversepropellerinterface.py @@ -0,0 +1,52 @@ +""" +Interface module that provides essential tools for the reversepropeller. +""" + +from abc import ABC, abstractmethod + + +class ReversePropellerInterface(ABC): + """ + Interface for reverse problem propeller blade. + + Given an iges file that has the blade's geometry, we solve an inverse problem + to retrieve the properties of the blades both in terms of: + - sectional properties (chord length, thickness, camber); + - sections location in the space (skew angle, pitch and rake) + """ + @abstractmethod + def _extract_solid_from_file(self): + """ + Abstract method that extract the solid object from the IGES file + """ + pass + + @abstractmethod + def _build_cylinder(self, radius): + """ + Abstract method that, given a radius, it builds the cylinder with specified radius + and main axis along x direction. This cylinder will be intersected with the blade + """ + pass + + @abstractmethod + def _build_intersection_cylinder_blade(self): + """ + Abstract method that computes the intersection between the blade and the cylinder + created by _build_cylinder() method + """ + pass + + @abstractmethod + def _camber_curve(self, radius): + """ + Abstract method that retrieves the camber curve once the intersection is computed + """ + pass + + @abstractmethod + def save_global_parameters(self, filename_csv): + """ + Abstract method to store the retrieved properties and section into a csv file + """ + From cdc592b4bd4254db82cdfc3bb3c1badd3fdd47d7 Mon Sep 17 00:00:00 2001 From: Davide Oberto Date: Thu, 5 Mar 2026 11:41:20 +0100 Subject: [PATCH 6/8] Fixed skew and rake computations --- .../reversepropeller/basereversepropeller.py | 29 +++---- bladex/reversepropeller/reversepropeller.py | 83 +------------------ .../reversepropellerbladex.py | 22 +++-- 3 files changed, 30 insertions(+), 104 deletions(-) diff --git a/bladex/reversepropeller/basereversepropeller.py b/bladex/reversepropeller/basereversepropeller.py index 9184cd2..49cb34a 100644 --- a/bladex/reversepropeller/basereversepropeller.py +++ b/bladex/reversepropeller/basereversepropeller.py @@ -1,6 +1,4 @@ -import os, errno -import os.path -import time +import sys import numpy as np import csv from bladex import CustomProfile, Blade, Propeller, Shaft @@ -13,7 +11,7 @@ import OCC.Core.TopoDS from OCC.Core.TopTools import TopTools_ListOfShape from OCC.Core.TopExp import TopExp_Explorer -from OCC.Core.TopAbs import TopAbs_VERTEX, TopAbs_EDGE, TopAbs_FACE, TopAbs_WIRE +from OCC.Core.TopAbs import TopAbs_VERTEX, TopAbs_EDGE, TopAbs_FACE, TopAbs_WIRE, TopAbs_SHELL from OCC.Core.gp import gp_Pnt, gp_Dir, gp_Ax2, gp_Vec, gp_Pln from OCC.Core.TColgp import TColgp_Array1OfPnt from OCC.Core.BRepAdaptor import BRepAdaptor_Curve @@ -128,7 +126,14 @@ def _extract_solid_from_file(self): sewer.Perform() result_sewed_blade = sewer.SewedShape() blade_solid_maker = BRepBuilderAPI_MakeSolid() - blade_solid_maker.Add(OCC.Core.TopoDS.topods.Shell(result_sewed_blade)) + if result_sewed_blade.ShapeType() == 0: + exp = TopExp_Explorer(result_sewed_blade, TopAbs_SHELL) + while exp.More(): + shell = topods.Shell(exp.Current()) + blade_solid_maker.Add(shell) + exp.Next() + else: + blade_solid_maker.Add(OCC.Core.TopoDS.topods.Shell(result_sewed_blade)) if (blade_solid_maker.IsDone()): self.blade_solid = blade_solid_maker.Solid() else: @@ -265,18 +270,6 @@ def _extract_parameters_and_transform_profile(self, radius): self.camber_points_on_plane[0, 0] = 0.0 self.camber_points_on_plane[-1, 0] = 1.0 - # Save the properties we wanted for each section - self.pitch_angles_list.append(self.pitch_angle) - self.pitch_list.append( - abs(2 * np.pi * radius / np.tan(self.pitch_angle)) / 1000.0) - self.skew_angles_list.append((self.skew / radius - np.pi) * 180 / np.pi) - self.skew_list.append((self.skew - np.pi * radius) / 1000.0) - total_rake = self.rake + self.rake_induced_by_skew - rake = total_rake - (self.skew - np.pi * radius) / np.tan( - self.pitch_angle) - rake = -rake / 1000.0 - self.rake_list.append(rake) - self.chord_length_list.append(self.chord_length / 1000.0) self.camber_points_on_plane = np.sort( self.camber_points_on_plane.view('float64,float64'), order=['f0'], @@ -303,7 +296,7 @@ def save_global_parameters(self, filename_csv): writer.writerow(("Skew angles list: ", self.skew_angles_list)) writer.writerow(("Rake list: ", self.rake_list)) writer.writerow(("Chord length list: ", self.chord_length_list)) - for j in range(self.num_sections): + for j in range(len(self.radii_list)): writer.writerow(("x section" + str(j) + " upper coordinates: ", self.xup[j, :])) writer.writerow(("x section" + str(j) + " lower coordinates: ", diff --git a/bladex/reversepropeller/reversepropeller.py b/bladex/reversepropeller/reversepropeller.py index 1404f3b..1bc98e2 100644 --- a/bladex/reversepropeller/reversepropeller.py +++ b/bladex/reversepropeller/reversepropeller.py @@ -144,6 +144,7 @@ def __init__(self, filename, radii_list, num_points_top_bottom): self._initial_camber_points_plane(radius) self._initial_airfoil_points_plane(radius) self._extract_parameters_and_transform_profile(radius) + self._store_properties(radius) self._airfoil_top_and_bottom_points() self.xup[ind_sec, :] = self.ascissa self.xdown[ind_sec, :] = self.ascissa @@ -463,7 +464,7 @@ def _initial_airfoil_points_plane(self, radius): self.leading_edge_point_on_plane + self.trailing_edge_point_on_plane) / 2.0 - def _extract_parameters_and_transform_profile(self, radius): + def _store_properties(self, radius): """ Private method which extracts the parameters (pitch, rake, skew ,...) related to a specific section of the blade, and then transforms the camber points @@ -475,78 +476,6 @@ def _extract_parameters_and_transform_profile(self, radius): Basically all points and edges of the profile section are transformed to fulfill the properties of the specific sectionwe are considering. """ - self.pitch_angle = np.arctan2( - self.leading_edge_point_on_plane[1] - - self.trailing_edge_point_on_plane[1], - self.leading_edge_point_on_plane[0] - - self.trailing_edge_point_on_plane[0]) - self.skew = self.mid_chord_point_on_plane[1] - self.rake_induced_by_skew = self.skew / np.tan(self.pitch_angle) - - self.airfoil_points_on_plane[:, 1:2] -= self.skew - self.camber_points_on_plane[:, 1:2] -= self.skew - self.leading_edge_point_on_plane[1] -= self.skew - self.trailing_edge_point_on_plane[1] -= self.skew - self.mid_chord_point_on_plane[1] -= self.skew - - self.airfoil_points_on_plane[:, 0:1] -= self.rake_induced_by_skew - self.camber_points_on_plane[:, 0:1] -= self.rake_induced_by_skew - self.leading_edge_point_on_plane[0] -= self.rake_induced_by_skew - self.trailing_edge_point_on_plane[0] -= self.rake_induced_by_skew - self.mid_chord_point_on_plane[0] -= self.rake_induced_by_skew - - self.rake = self.mid_chord_point_on_plane[0] - - self.airfoil_points_on_plane[:, 0:1] -= self.rake - self.camber_points_on_plane[:, 0:1] -= self.rake - self.leading_edge_point_on_plane[0] -= self.rake - self.trailing_edge_point_on_plane[0] -= self.rake - self.mid_chord_point_on_plane[0] -= self.rake - - rotation_matrix = np.zeros((2, 2)) - rotation_matrix[0][0] = np.cos(-self.pitch_angle) - rotation_matrix[0][1] = -np.sin(-self.pitch_angle) - rotation_matrix[1][0] = np.sin(-self.pitch_angle) - rotation_matrix[1][1] = np.cos(-self.pitch_angle) - - self.airfoil_points_on_plane = self.airfoil_points_on_plane.dot( - rotation_matrix.transpose()) - self.camber_points_on_plane = self.camber_points_on_plane.dot( - rotation_matrix.transpose()) - self.leading_edge_point_on_plane = self.leading_edge_point_on_plane.dot( - rotation_matrix.transpose()) - self.trailing_edge_point_on_plane = self.trailing_edge_point_on_plane.dot( - rotation_matrix.transpose()) - self.mid_chord_point_on_plane = self.mid_chord_point_on_plane.dot( - rotation_matrix.transpose()) - - self.chord_length = ((self.leading_edge_point_on_plane[0] - - self.trailing_edge_point_on_plane[0])**2 + - (self.leading_edge_point_on_plane[1] - - self.trailing_edge_point_on_plane[1])**2)**.5 - - self.airfoil_points_on_plane[:, 0:2] /= self.chord_length - self.camber_points_on_plane[:, 0:2] /= self.chord_length - self.leading_edge_point_on_plane[0:1] /= self.chord_length - self.trailing_edge_point_on_plane[0:1] /= self.chord_length - self.mid_chord_point_on_plane[0:1] /= self.chord_length - - self.airfoil_points_on_plane[:, 0:1] *= -1.0 - self.camber_points_on_plane[:, 0:1] *= -1.0 - self.leading_edge_point_on_plane[0] *= -1.0 - self.trailing_edge_point_on_plane[0] *= -1.0 - self.mid_chord_point_on_plane[0] *= -1.0 - - self.airfoil_points_on_plane[:, 0:1] += 0.5 - self.camber_points_on_plane[:, 0:1] += 0.5 - self.leading_edge_point_on_plane[0] += 0.5 - self.trailing_edge_point_on_plane[0] += 0.5 - self.mid_chord_point_on_plane[0] += 0.5 - - self.camber_points_on_plane = np.matrix(self.camber_points_on_plane) - self.camber_points_on_plane[0, 0] = 0.0 - self.camber_points_on_plane[-1, 0] = 1.0 - # Save the properties we wanted for each section self.pitch_angles_list.append(self.pitch_angle) self.pitch_list.append( @@ -559,14 +488,6 @@ def _extract_parameters_and_transform_profile(self, radius): self.rake = -self.rake / 1000.0 self.rake_list.append(self.rake) self.chord_length_list.append(self.chord_length / 1000.0) - self.camber_points_on_plane = np.sort( - self.camber_points_on_plane.view('float64,float64'), - order=['f0'], - axis=0).view(np.float64) - self.f_camber = interp1d( - np.squeeze(np.asarray(self.camber_points_on_plane[:, 0])), - np.squeeze(np.asarray(self.camber_points_on_plane[:, 1])), - kind='cubic') def _airfoil_top_and_bottom_points(self): """ diff --git a/bladex/reversepropeller/reversepropellerbladex.py b/bladex/reversepropeller/reversepropellerbladex.py index a11f17d..08ee028 100644 --- a/bladex/reversepropeller/reversepropellerbladex.py +++ b/bladex/reversepropeller/reversepropellerbladex.py @@ -100,6 +100,7 @@ def __init__(self, filename, radii_list, num_points_top_bottom): self._initial_airfoil_points_plane(radius) self._airfoil_top_and_bottom_points_before_transformations(radius) self._extract_parameters_and_transform_profile(radius) + self._store_properties(radius) self._transform_top_and_bottom() self._airfoil_top_and_bottom_points_after_transformations() self.xup[ind_sec, :] = self.ascissa @@ -122,7 +123,6 @@ def _build_intersection_cylinder_blade(self): section_builder.ComputePCurveOn2(True) section_builder.Build() self.section = section_builder.Shape() - wire_maker = BRepBuilderAPI_MakeWire() wire_maker_top = BRepBuilderAPI_MakeWire() wire_maker_bottom = BRepBuilderAPI_MakeWire() @@ -171,10 +171,8 @@ def _camber_curve(self, radius): Computation of the camber points. We get the chord, move along it and fint the intersection of the othogonal plane to the chord-curvilinear absissa and the top-bottom wires """ - # self.trailing_edge = self.curve_adaptor_top.Value(0.0) - # self.leading_edge = self.curve_adaptor_top.Value(self.curve_adaptor_top.LastParameter()) - self.trailing_edge = self.curve_adaptor_bottom.Value(0.0) - self.leading_edge = self.curve_adaptor_bottom.Value(self.curve_adaptor_bottom.LastParameter()) + self.trailing_edge = self.curve_adaptor_top.Value(0.0) + self.leading_edge = self.curve_adaptor_top.Value(self.curve_adaptor_top.LastParameter()) # Equation of generic geodesic from leading to trailing edge: # x = lambda*s + mu, y = R*sin(a*s+b), z = R*cos(a*s+b) # with s parametrization of curve and lambda, mu, a, b to be found @@ -363,6 +361,20 @@ def _airfoil_top_and_bottom_points_before_transformations(self, radius): self.airfoil_top = np.array(self.airfoil_top) self.airfoil_bottom = np.array(self.airfoil_bottom) + def _store_properties(self, radius): + # Save the properties we wanted for each section + self.pitch_angles_list.append(self.pitch_angle) + self.pitch_list.append( + abs(2 * np.pi * radius / np.tan(self.pitch_angle)) / 1000.0) + self.skew_angles_list.append(-(self.skew / radius) * 180 / np.pi) + self.skew_list.append(self.skew / 1000.0) + total_rake = self.rake + self.rake_induced_by_skew + rake = total_rake - (self.skew) / np.tan( + self.pitch_angle) + rake = -rake / 1000.0 + self.rake_list.append(rake) + self.chord_length_list.append(self.chord_length / 1000.0) + def _transform_top_and_bottom(self): """ Method that transforms the top and bottom coordinates from physical coordinates From b4946c70292bd54156da95be703a5d57d3f77871 Mon Sep 17 00:00:00 2001 From: Davide Oberto Date: Thu, 5 Mar 2026 14:57:42 +0100 Subject: [PATCH 7/8] Added a tutorial on reversepropellerbladex --- bladex/__init__.py | 4 + .../tutorial-7-blades-reversepropeller.ipynb | 145 ++++++++++++++++++ 2 files changed, 149 insertions(+) create mode 100644 tutorials/tutorial-7-blades-reversepropeller.ipynb diff --git a/bladex/__init__.py b/bladex/__init__.py index 948c91d..62566cc 100644 --- a/bladex/__init__.py +++ b/bladex/__init__.py @@ -17,6 +17,10 @@ from .deform import Deformation from .params import ParamFile from .ndinterpolator import RBF, reconstruct_f, scipy_bspline +from .reversepropeller import ReversePropellerInterface +from .reversepropeller import BaseReversePropeller from .reversepropeller import ReversePropeller +from .reversepropeller import ReversePropellerBladeX +from .reversepropeller.reversepropeller import ReversePropeller from .shaft.cylinder_shaft import CylinderShaft from .intepolatedface import InterpolatedFace diff --git a/tutorials/tutorial-7-blades-reversepropeller.ipynb b/tutorials/tutorial-7-blades-reversepropeller.ipynb new file mode 100644 index 0000000..bf27b03 --- /dev/null +++ b/tutorials/tutorial-7-blades-reversepropeller.ipynb @@ -0,0 +1,145 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "c2ec2fd3", + "metadata": {}, + "source": [ + "# BladeX" + ] + }, + { + "cell_type": "markdown", + "id": "ad979fa2", + "metadata": {}, + "source": [ + "## Tutorial 7: Retrieve blade's parameters from iges file using reversepropeller class" + ] + }, + { + "cell_type": "markdown", + "id": "c0ab2647", + "metadata": {}, + "source": [ + "The goal of this tutorial is to show how to start from an iges file containing a blade to obtain its parameters, both for the sections (camber, thickness) and for the whole blade (chord lengths, pitch, rake, skew angles)." + ] + }, + { + "cell_type": "markdown", + "id": "5f21ab8e", + "metadata": {}, + "source": [ + "We start from some useful imports, as usual." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "4362678f", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from bladex import Blade, NacaProfile, ReversePropellerBladeX, ReversePropeller" + ] + }, + { + "cell_type": "markdown", + "id": "45d4e8f4", + "metadata": {}, + "source": [ + "Then, we create the same blade as for Tutorial 6 and store it as iges file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa8f973f", + "metadata": {}, + "outputs": [], + "source": [ + "# Create sections with NACAProfile\n", + "n_sections = 10\n", + "sections = np.asarray([NacaProfile(digits='4412', n_points=50) for i in range(n_sections)])\n", + "\n", + "# Define blade parameters\n", + "radii = np.arange(1.0, 11.0, 1.0)\n", + "chord_lengths = np.array([0.05, 2.5, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7])\n", + "pitch = np.arange(1., 11.)\n", + "rake = np.arange(0.1, 1.1, 0.1)\n", + "skew_angles = np.arange(1., 21, 2.)\n", + "\n", + "# Create the Blade\n", + "blade = Blade(sections=sections,\n", + " radii=radii,\n", + " chord_lengths=chord_lengths,\n", + " pitch=pitch,\n", + " rake=rake,\n", + " skew_angles=skew_angles)\n", + "\n", + "# We build the blade to be able to store its iges version\n", + "# The method build() internally calls apply_transformations()\n", + "blade.build()\n", + "\n", + "# Storing the blade in iges format\n", + "iges_filename='data/blade_reversepropeller.igs'\n", + "blade.export_iges(iges_filename)" + ] + }, + { + "cell_type": "markdown", + "id": "61aabd42", + "metadata": {}, + "source": [ + "Now that we have our starting iges file, we can retrieve the blades parameters. In BladeX you can either use the `ReversePropellerBladeX`, if the blade has been generated using BladeX, or `ReversePropeller` otherwise. The former is robust and more accurate for BladeX blades because it exploits known structure of the iges file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "78b71d18", + "metadata": {}, + "outputs": [], + "source": [ + "# We need to specify: \n", + "# i) the name of the iges file; \n", + "# ii) values of the radii in mm at which the parameters are to be extracted;\n", + "# iii) how many points to use to reconstruct the sections for each radius value\n", + "# In this case we use the same radii values (extrema excluded) of the blade's generation for a better comparison\n", + "radii_reverse = np.arange(2.0, 10., 1.)\n", + "# This operation may take few minutes\n", + "reverse = ReversePropellerBladeX(iges_filename, 1000*radii_reverse, num_points_top_bottom=100)\n", + "\n", + "# Now we can retrieve the blade's parameters\n", + "print(reverse.chord_length_list)\n", + "print(reverse.pitch_list)\n", + "print(reverse.rake_list)\n", + "print(reverse.skew_angles_list)\n", + "\n", + "# Finally, we can store the parameters to a csv file\n", + "reverse.save_global_parameters('data/blade_parameters.csv')" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "marinai", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 69842e920a8781e056733142f0f4d2256a782024 Mon Sep 17 00:00:00 2001 From: Davide Oberto Date: Thu, 5 Mar 2026 15:53:57 +0100 Subject: [PATCH 8/8] Run tutorial and included output cells --- .../tutorial-7-blades-reversepropeller.ipynb | 34 ++++++++++++++++--- 1 file changed, 29 insertions(+), 5 deletions(-) diff --git a/tutorials/tutorial-7-blades-reversepropeller.ipynb b/tutorials/tutorial-7-blades-reversepropeller.ipynb index bf27b03..d1f7660 100644 --- a/tutorials/tutorial-7-blades-reversepropeller.ipynb +++ b/tutorials/tutorial-7-blades-reversepropeller.ipynb @@ -34,7 +34,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 1, "id": "4362678f", "metadata": {}, "outputs": [], @@ -53,10 +53,21 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "aa8f973f", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(10, 3, 50)\n", + "(10, 3, 50)\n", + "(2, 3, 50)\n", + "(2, 3, 50)\n" + ] + } + ], "source": [ "# Create sections with NACAProfile\n", "n_sections = 10\n", @@ -96,10 +107,23 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "id": "78b71d18", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[32;1mTotal number of loaded entities 52.\n", + "\u001b[0m\n", + "[2.4999999998198215, 3.499999999872226, 3.9999999996733515, 4.500000000152133, 4.999999999959586, 5.499999999841911, 6.000000000221554, 6.500000000430857]\n", + "[1.999999999872152, 2.9999999999180984, 3.9999999997822977, 4.9999999997217, 6.000000002259764, 7.000000003635259, 8.000000002720112, 9.000000000963933]\n", + "[0.1999999999957766, 0.30000000000569976, 0.39999999999797053, 0.49999999999262956, 0.6000000000362385, 0.7000000000603949, 0.8000000000447604, 0.9000000000155906]\n", + "[2.9999999959870425, 4.999999997162759, 6.999999998270603, 8.9999999995899, 11.000000000103695, 12.999999998485547, 15.00000000063232, 17.000000001394554]\n" + ] + } + ], "source": [ "# We need to specify: \n", "# i) the name of the iges file; \n",