Fixed skew and rake computations

Author: DavideOberto

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: ",

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):
         """

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