diff --git a/examples/complex/gear/gear.py b/examples/complex/gear/gear.py new file mode 100644 index 0000000..0c4e448 --- /dev/null +++ b/examples/complex/gear/gear.py @@ -0,0 +1,191 @@ +#!/usr/bin/env python +import os +import time + +import numpy as np +import scipy.optimize +from involute_gear import InvoluteGear +from tooth import ToothSketch + +import classy_blocks as cb +from classy_blocks.util import functions as f + +time_start = time.time() + +mesh = cb.Mesh() + +# parameters +RADIUS_EXPANSION = 1.1 + +# create an interpolated curve that represents a gear tooth +fillet = 0.1 # Must be more than zero! +gear = InvoluteGear(fillet=fillet, arc_step_size=0.1, max_steps=1000, teeth=15) +tng_points = gear.generate_tooth_and_gap() + +z_coords = np.zeros(len(tng_points[0])) +tng_points = np.stack((tng_points[0], tng_points[1], z_coords)).T +tng_points = np.flip(tng_points, axis=0) + +# add start and end points exactly on the 2pi/teeth +start_point = f.to_polar(tng_points[0], axis="z") +start_point[1] = np.pi / gear.teeth +start_point = f.to_cartesian(start_point) + +end_point = f.to_polar(tng_points[-1], axis="z") +end_point[1] = -np.pi / gear.teeth +end_point = f.to_cartesian(end_point) + +tng_points = np.concatenate(([start_point], tng_points, [end_point])) + +tooth_curve = cb.LinearInterpolatedCurve(tng_points) +gear_params = np.linspace(0, 1, num=8) + + +# fix points 1 and 3: +# 1 is on radius gear.root_radius + fillet/2 +def frad(t, radius): + return f.norm(tooth_curve.get_point(t)) - radius + + +gear_params[1] = scipy.optimize.brentq(lambda t: frad(t, gear.root_radius + fillet / 2), 0, 0.25) + +# 3 is on radius gear.outer_radius - fillet/2 +gear_params[3] = scipy.optimize.brentq(lambda t: frad(t, gear.outer_radius - fillet / 2), 0.25, 0.5) + +gear_params[6] = 1 - gear_params[1] +gear_params[4] = 1 - gear_params[3] + +gear_params[2] = (gear_params[1] + gear_params[3]) / 2 +gear_params[5] = (gear_params[4] + gear_params[6]) / 2 + +gear_points = np.array([tooth_curve.get_point(t) for t in gear_params]) + +outer_radius = f.norm(gear_points[3] * RADIUS_EXPANSION) +p11_polar = f.to_polar(gear_points[-1], axis="z") +p14_polar = f.to_polar(gear_points[0], axis="z") +angles = np.linspace(p11_polar[1], p14_polar[1], num=4) +tangential_points = np.array([f.to_cartesian([outer_radius, angle, 0]) for angle in angles]) + +radial_points_1 = np.linspace(gear_points[-1], tangential_points[0], axis=0, num=5)[1:-1] +radial_points_2 = np.linspace(tangential_points[-1], gear_points[0], axis=0, num=5)[1:-1] + +outer_points = np.concatenate((gear_points, radial_points_1, tangential_points, radial_points_2)) +inner_points = np.zeros((6, 3)) + +positions = np.concatenate((outer_points, inner_points)) + +# At this point, a smoother would reliably +# produce almost the best blocking if this was a convex sketch. +# Alas, this is a severely concave case so smoothing will produce +# degenerate quads which even optimizers won't be able to fix. +# It's best to manually position points, then optimize the sketch. + + +def mirror(target, source): + # once a position is obtained, the mirrored counterpart is also determined + positions[target] = [positions[source][0], -positions[source][1], 0] + + +# fix points 18 and 23 because optimization doesn't 'see' curved edges +# and produces high non-orthogonality +dir_0 = f.unit_vector(gear_points[0] - gear_points[1]) +dir_2 = f.unit_vector(gear_points[2] - gear_points[1]) +dir_18 = f.unit_vector(dir_0 + dir_2) + +positions[18] = gear_points[1] + dir_18 * f.norm(gear_points[0] - gear_points[1]) / 2 +mirror(23, 18) + +positions[17] = positions[0] + f.unit_vector(positions[17] - positions[0]) * f.norm(positions[18] - positions[1]) +mirror(8, 17) + + +# other points are somewhere between... +def midpoint(target, left, right): + positions[target] = (positions[left] + positions[right]) / 2 + + +midpoint(19, 2, 16) +midpoint(20, 3, 13) + +# and their mirrored counterparts +mirror(22, 19) +mirror(21, 20) + + +sketch = ToothSketch(positions, tooth_curve) + +# Optimize the sketch: +optimizer = cb.SketchOptimizer(sketch) + +# point 2 is on gear curve +optimizer.add_clamp(cb.CurveClamp(positions[2], tooth_curve)) + +# point 13 is movable radially +optimizer.add_clamp(cb.RadialClamp(positions[13], [0, 0, 0], [0, 0, 1])) + +# 15-17 move along a line +for i in (15, 16, 17): + optimizer.add_clamp(cb.LineClamp(positions[i], gear_points[0], tangential_points[-1])) + +# freely movable points (on sketch plane) +# TODO: easier clamp definition for sketch optimization +for i in (19, 20): + optimizer.add_clamp(cb.PlaneClamp(sketch.positions[i], sketch.positions[i], sketch.normal)) + +# Links! +symmetry_pairs = [ + (2, 5), + (19, 22), + (20, 21), + (17, 8), + (16, 9), + (15, 10), +] + +for pair in symmetry_pairs: + optimizer.add_link(cb.SymmetryLink(positions[pair[0]], positions[pair[1]], f.vector(0, 1, 0), f.vector(0, 0, 0))) + +optimizer.optimize() + +stack = cb.TransformedStack( + sketch, + [cb.Translation([0, 0, 4]), cb.Rotation(sketch.normal, 10 * np.pi / 180, [0, 0, 0])], + 2, + [cb.Translation([0, 0, 2]), cb.Rotation(sketch.normal, 5 * np.pi / 180, [0, 0, 0])], +) + +# TODO: this be mighty clumsy; unclumsify +bulk_size = 0.1 +wall_size = 0.01 + +stack.shapes[0].chop(0, start_size=bulk_size) +stack.shapes[0].chop(1, start_size=wall_size, end_size=bulk_size / 2) +stack.shapes[0].operations[10].chop(1, start_size=bulk_size) +stack.chop(count=8) + +# patches +for shape in stack.shapes: + for i in range(7): + shape.operations[i].set_patch("front", "gear") + + for i in (9, 10, 11): + shape.operations[i].set_patch("front", "outer") + +for operation in stack.shapes[0].operations: + operation.set_patch("bottom", "bottom") + +for operation in stack.shapes[-1].operations: + operation.set_patch("top", "top") + +mesh.add(stack) + +for i, angle in enumerate(np.linspace(0, 2 * np.pi, num=gear.teeth, endpoint=False)[1:]): + print(f"Adding tooth {i+2}") + mesh.add(stack.copy().rotate(angle, [0, 0, 1], [0, 0, 0])) + +print("Writing mesh...") +mesh.write(os.path.join("..", "..", "case", "system", "blockMeshDict"), debug_path="debug.vtk") + +time_end = time.time() + +print(f"Elapsed time: {time_end - time_start}s") diff --git a/examples/complex/gear/involute_gear.py b/examples/complex/gear/involute_gear.py new file mode 100644 index 0000000..1934725 --- /dev/null +++ b/examples/complex/gear/involute_gear.py @@ -0,0 +1,280 @@ +"""A copy of py_gear_gen (https://github.com/heartworm/py_gear_gen), +slightly modified to fit in classy_blocks for the gear pump example""" + +import numpy as np + +from classy_blocks.util import functions as f + + +class DimensionError(Exception): + pass + + +def rotation_matrix(theta): + return np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]) + + +def flip_matrix(h, v): + return [[-1 if h else 1, 0], [0, -1 if v else 1]] + + +def polar_to_cart(*coords): + if len(coords) == 1: + coords = coords[0] + r, ang = coords + return r * np.cos(ang), r * np.sin(ang) + + +def cart_to_polar(*coords): + if len(coords) == 1: + coords = coords[0] + x, y = coords + return np.sqrt(x * x + y * y), np.arctan2(y, x) + + +class InvoluteGear: + def __init__( + self, + module=1, + teeth=30, + pressure_angle_deg=20.0, + fillet=0.0, + backlash=0.0, + max_steps=100, + arc_step_size=0.1, + reduction_tolerance_deg=0.0, + dedendum_factor=1.157, + addendum_factor=1.0, + ring=False, + ): + """ + Construct an involute gear, ready for generation using one of the generation methods. + :param module: The 'module' of the gear. (Diameter / Teeth) + :param teeth: How many teeth in the desired gear. + :param pressure_angle_deg: The pressure angle of the gear in DEGREES. + :param fillet: The radius of the fillet connecting a tooth to the root circle. NOT WORKING in ring gear. + :param backlash: The circumfrential play between teeth, + if meshed with another gear of the same backlash held stationary + :param max_steps: Maximum steps allowed to generate the involute profile. Higher is more accurate. + :param arc_step_size: The step size used for generating arcs. + :param ring: True if this is a ring (internal) gear, otherwise False. + """ + + pressure_angle = f.deg2rad(pressure_angle_deg) + self.reduction_tolerance = f.deg2rad(reduction_tolerance_deg) + self.module = module + self.teeth = teeth + self.pressure_angle = pressure_angle + + # Addendum is the height above the pitch circle that the tooth extends to + self.addendum = addendum_factor * module + # Dedendum is the depth below the pitch circle the root extends to. 1.157 is a std value allowing for clearance. + self.dedendum = dedendum_factor * module + + # If the gear is a ring gear, then the clearance needs to be on the other side + if ring: + temp = self.addendum + self.addendum = self.dedendum + self.dedendum = temp + + # The radius of the pitch circle + self.pitch_radius = (module * teeth) / 2 + # The radius of the base circle, used to generate the involute curve + self.base_radius = np.cos(pressure_angle) * self.pitch_radius + # The radius of the gear's extremities + self.outer_radius = self.pitch_radius + self.addendum + # The radius of the gaps between the teeth + self.root_radius = self.pitch_radius - self.dedendum + + # The radius of the fillet circle connecting the tooth to the root circle + self.fillet_radius = fillet if not ring else 0 + + # The angular width of a tooth and a gap. 360 degrees divided by the number of teeth + self.theta_tooth_and_gap = np.pi * 2 / teeth + # Converting the circumfrential backlash into an angle + angular_backlash = backlash / 2 / self.pitch_radius + # The angular width of the tooth at the pitch circle minus backlash, not taking the involute into account + self.theta_tooth = self.theta_tooth_and_gap / 2 + (-angular_backlash if not ring else angular_backlash) + # Where the involute profile intersects the pitch circle, found on iteration. + self.theta_pitch_intersect = None + # The angular width of the full tooth, at the root circle + self.theta_full_tooth = 0 + + self.max_steps = max_steps + self.arc_step_size = arc_step_size + + """ + Reduces a line of many points to less points depending on the allowed angle tolerance + """ + + def reduce_polyline(self, polyline): + vertices = [[], []] + last_vertex = [polyline[0][0], polyline[1][0]] + + # Look through all vertices except start and end vertex + # Calculate by how much the lines before and after the vertex + # deviate from a straight path. + # If the deviation angle exceeds the specification, store it + for vertex_idx in range(1, len(polyline[0]) - 1): + next_slope = np.arctan2( + polyline[1][vertex_idx + 1] - polyline[1][vertex_idx + 0], + polyline[0][vertex_idx + 1] - polyline[0][vertex_idx + 0], + ) + prev_slope = np.arctan2( + polyline[1][vertex_idx - 0] - last_vertex[1], polyline[0][vertex_idx - 0] - last_vertex[0] + ) + + deviation_angle = abs(prev_slope - next_slope) + + if deviation_angle > self.reduction_tolerance: + vertices[0] += [polyline[0][vertex_idx]] + vertices[1] += [polyline[1][vertex_idx]] + last_vertex = [polyline[0][vertex_idx], polyline[1][vertex_idx]] + + # Return vertices along with first and last point of the original polyline + return np.array( + [ + np.concatenate([[polyline[0][0]], vertices[0], [polyline[0][-1]]]), + np.concatenate([[polyline[1][0]], vertices[1], [polyline[1][-1]]]), + ] + ) + + def generate_half_tooth(self): + """ + Generate half an involute profile, ready to be mirrored in order to create one symmetrical involute tooth + :return: A numpy array, of the format [[x1, x2, ... , xn], [y1, y2, ... , yn]] + """ + # Theta is the angle around the circle, however PHI is simply a parameter for iteratively building the involute + phis = np.linspace(0, np.pi, self.max_steps) + points = [] + reached_limit = False + self.theta_pitch_intersect = None + + for phi in phis: + x = (self.base_radius * np.cos(phi)) + (phi * self.base_radius * np.sin(phi)) + y = (self.base_radius * np.sin(phi)) - (phi * self.base_radius * np.cos(phi)) + point = (x, y) + dist, theta = cart_to_polar(point) + + if self.theta_pitch_intersect is None and dist >= self.pitch_radius: + self.theta_pitch_intersect = theta + self.theta_full_tooth = self.theta_pitch_intersect * 2 + self.theta_tooth + elif self.theta_pitch_intersect is not None and theta >= self.theta_full_tooth / 2: + reached_limit = True + break + + if dist >= self.outer_radius: + points.append(polar_to_cart((self.outer_radius, theta))) + elif dist <= self.root_radius: + points.append(polar_to_cart((self.root_radius, theta))) + else: + points.append((x, y)) + + if not reached_limit: + raise Exception("Couldn't complete tooth profile.") + + return np.transpose(points) + + def generate_half_root(self): + """ + Generate half of the gap between teeth, for the first tooth + :return: A numpy array, of the format [[x1, x2, ... , xn], [y1, y2, ... , yn]] + """ + root_arc_length = (self.theta_tooth_and_gap - self.theta_full_tooth) * self.root_radius + + points_root = [] + for theta in np.arange( + self.theta_full_tooth, + self.theta_tooth_and_gap / 2 + self.theta_full_tooth / 2, + self.arc_step_size / self.root_radius, + ): + # The current circumfrential position we are in the root arc, starting from 0 + arc_position = (theta - self.theta_full_tooth) * self.root_radius + # If we are in the extemities of the root arc (defined by fillet_radius), then we are in a fillet + in_fillet = min((root_arc_length - arc_position), arc_position) < self.fillet_radius + + r = self.root_radius + + if in_fillet: + # Add a circular profile onto the normal root radius to form the fillet. + # High near the edges, small towards the centre + # The min() function handles the situation where the fillet size is massive and overlaps itself + circle_pos = min(arc_position, (root_arc_length - arc_position)) + r = r + ( + self.fillet_radius - np.sqrt(pow(self.fillet_radius, 2) - pow(self.fillet_radius - circle_pos, 2)) + ) + points_root.append(polar_to_cart((r, theta))) + + return np.transpose(points_root) + + def generate_roots(self): + """ + Generate both roots on either side of the first tooth + :return: A numpy array, of the format + [ [[x01, x02, ... , x0n], [y01, y02, ... , y0n]], [[x11, x12, ... , x1n], [y11, y12, ... , y1n]] ] + """ + self.half_root = self.generate_half_root() + self.half_root = np.dot(rotation_matrix(-self.theta_full_tooth / 2), self.half_root) + points_second_half = np.dot(flip_matrix(False, True), self.half_root) + points_second_half = np.flip(points_second_half, 1) + self.roots = [points_second_half, self.half_root] + + # Generate a second set of point-reduced root + self.half_root_reduced = self.reduce_polyline(self.half_root) + points_second_half = np.dot(flip_matrix(False, True), self.half_root_reduced) + points_second_half = np.flip(points_second_half, 1) + self.roots_reduced = [points_second_half, self.half_root_reduced] + + return self.roots_reduced + + def generate_tooth(self): + """ + Generate only one involute tooth, without an accompanying tooth gap + :return: A numpy array, of the format [[x1, x2, ... , xn], [y1, y2, ... , yn]] + """ + self.half_tooth = self.generate_half_tooth() + self.half_tooth = np.dot(rotation_matrix(-self.theta_full_tooth / 2), self.half_tooth) + points_second_half = np.dot(flip_matrix(False, True), self.half_tooth) + points_second_half = np.flip(points_second_half, 1) + self.tooth = np.concatenate((self.half_tooth, points_second_half), axis=1) + + # Generate a second set of point-reduced teeth + self.half_tooth_reduced = self.reduce_polyline(self.half_tooth) + points_second_half = np.dot(flip_matrix(False, True), self.half_tooth_reduced) + points_second_half = np.flip(points_second_half, 1) + self.tooth_reduced = np.concatenate((self.half_tooth_reduced, points_second_half), axis=1) + + return self.tooth_reduced + + def generate_tooth_and_gap(self): + """ + Generate only one tooth and one root profile, ready to be duplicated by rotating around the gear center + :return: A numpy array, of the format [[x1, x2, ... , xn], [y1, y2, ... , yn]] + """ + + points_tooth = self.generate_tooth() + points_roots = self.generate_roots() + self.tooth_and_gap = np.concatenate((points_roots[0], points_tooth, points_roots[1]), axis=1) + return self.tooth_and_gap + + def generate_gear(self): + """ + Generate the gear profile, and return a sequence of co-ordinates representing the outline of the gear + :return: A numpy array, of the format [[x1, x2, ... , xn], [y1, y2, ... , yn]] + """ + + points_tooth_and_gap = self.generate_tooth_and_gap() + points_teeth = [ + np.dot(rotation_matrix(self.theta_tooth_and_gap * n), points_tooth_and_gap) for n in range(self.teeth) + ] + points_gear = np.concatenate(points_teeth, axis=1) + return points_gear + + def get_point_list(self): + """ + Generate the gear profile, and return a sequence of co-ordinates representing the outline of the gear + :return: A numpy array, of the format [[x1, y2], [x2, y2], ... , [xn, yn]] + """ + + gear = self.generate_gear() + return np.transpose(gear) diff --git a/examples/complex/gear/tooth-blocking.svg b/examples/complex/gear/tooth-blocking.svg new file mode 100644 index 0000000..7023566 --- /dev/null +++ b/examples/complex/gear/tooth-blocking.svg @@ -0,0 +1,712 @@ + + + + + + + + + + + + + + + + + + + + + + + x + y + O + + + gear + tangential + + radial 1 + + radial 2 + + + + + + + + + + + + + + + + + + + 0 + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + + + 0 + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + + diff --git a/examples/complex/gear/tooth.py b/examples/complex/gear/tooth.py new file mode 100644 index 0000000..afc4a88 --- /dev/null +++ b/examples/complex/gear/tooth.py @@ -0,0 +1,50 @@ +from typing import ClassVar + +import classy_blocks as cb + + +class ToothSketch(cb.MappedSketch): + quads: ClassVar = [ + # layers on tooth wall + [0, 1, 18, 17], # 0 + [1, 2, 19, 18], # 1 + [2, 3, 20, 19], # 2 + [3, 4, 21, 20], # 3 + [4, 5, 22, 21], # 4 + [5, 6, 23, 22], # 5 + [6, 7, 8, 23], # 6 + # surrounding blocks + [8, 9, 22, 23], # 7 + [9, 10, 21, 22], # 8 + [11, 12, 21, 10], # 9 + [12, 13, 20, 21], # 10 + [13, 14, 15, 20], # 11 + [15, 16, 19, 20], # 12 + [16, 17, 18, 19], # 13 + ] + + # In x-direction, only one half needs to be chopped, the other half will be + # copied from teeth on the other side. + + # In y-direction, block 10 should also be in the list + # except that it doesn't need boundary layers that will be created for blocks on the tooth. + # It will be chopped manually. + chops: ClassVar = [ + [1, 2, 9, 10, 11], + [0], + ] + + def __init__(self, positions, curve: cb.LinearInterpolatedCurve): + self.curve = curve + super().__init__(positions, self.quads) + + def add_edges(self) -> None: + for i in range(7): + # If all edges refer to the same curve, it will be + # transformed N-times, N being the number of entities + # where that curve is used. Therefore copy it + # so that each edge will have its own object to work with. + self.faces[i].add_edge(0, cb.OnCurve(self.curve.copy())) + + for i in (9, 10, 11): + self.faces[i].add_edge(0, cb.Origin([0, 0, 0]))