regen_tool.py

"""
Programmatically generate rocket engine regenerative cooling geometry using CadQuery 2.
"""

import cadquery as cq
import math
import argparse
from tqdm import tqdm
import time
start_time = time.time()

print(""""
                                                   _  
                                    _             | | 
  ____ _____  ____ _____ ____     _| |_ ___   ___ | | 
 / ___) ___ |/ _  | ___ |  _ \   (_   _) _ \ / _ \| | 
| |   | ____( (_| | ____| | | |    | || |_| | |_| | | 
|_|   |_____)\___ |_____)_| |_|     \__)___/ \___/ \_)
            (_____|                                   

By: Cristian Bicheru
""")

# Read an RPA .txt based contour file
def load_contour(filename):
    with open(filename, "r") as f:
        data = [x for x in f.readlines() if not x.startswith('#')]
    data = [x.strip().replace(' ', '').split('\t') for x in data]
    data = [[float(x[0]), float(x[1])] for x in data if len(x) == 2]
    return data

# Convert a vector into an array, optionally setting the z-component to 0
def unpack_vec(vec, planar=False):
    return vec.x, vec.y, (vec.z if not planar else 0)

# Input parsing
parser = argparse.ArgumentParser(description="Generate regeneratively cooled engine geometry.")
parser.add_argument("-f", type=str, help="The RPA generated contour. Use the uniform txt based contour format from RPA.", required=True)
parser.add_argument("-t1", type=float, help="The CC wall thickness.", required=True)
parser.add_argument("-t2", type=float, help="The outer wall thickness.", required=True)
parser.add_argument("-w", type=float, help="The rib thickness.", required=True)
parser.add_argument("-hc", type=float, help="The height of the channels.", required=True)
parser.add_argument("-n", type=int, help="The number of channels.", required=True)
parser.add_argument("--angle", type=float, default=30.5, help="The overhang angle of the channels in degrees.")
parser.add_argument("--nperp", type=int, default=20, help="Number of interpolation points.")
parser.add_argument("--srib", action='store_true', default=True, help="Whether to scale the rib thickness with the sqrt of the radius.")
parser.add_argument("--full", action='store_true', default=False, help="Whether to export the regen channel negative or to generate the full regen geometry.")
parser.add_argument("out_file", help="The output filename.")
args = parser.parse_args()

# Input parsing
contour_filename = args.f
inner_wall_thickness = args.t1
rib_thickness = args.w
outer_wall_thickness = args.t2
channel_height = args.hc
num_channels = args.n
n_perp = args.nperp

tan_gamma = math.tan( math.pi/2 - math.radians(args.angle) )
theta_per_channel = 2*math.pi/num_channels

if contour_filename.endswith(".dxf"):
    raise RuntimeError("Export nozzle contour as a .txt file.")

print("Loading nozzle contour...")
contour = load_contour(contour_filename)

rmin = min([x[1] for x in contour])

inner_surface_pts = [[] for _ in range(n_perp)]
outer_surface_pts = [[] for _ in range(n_perp)]
theta = 0

print("Generating cross-sections...")
for i in tqdm(range(len(contour))):
    r = contour[i][1]
    z = contour[i][0]
    if i != len(contour)-1:
        dr = contour[i+1][1] - r
        dz = contour[i+1][0] - z
    
    dtheta = math.sqrt( dz**2 / tan_gamma**2 - dr**2 ) / r
    r_inner = r + inner_wall_thickness
    r_outer = r_inner + channel_height

    inner_twist_angle = math.atan(r_inner*dtheta/dz)
    outer_twist_angle = math.atan(r_outer*dtheta/dz)
    flare_angle = math.atan(dr/dz)

    adjusted_rib_thickness = rib_thickness * ( (r / rmin) ** 0.5 if args.srib else 1 )

    inner_eff_rib_thickness = adjusted_rib_thickness / ( 
        math.cos(inner_twist_angle)**2 + math.sin(inner_twist_angle)**2 * math.sin(flare_angle)**2 )**0.5
    outer_eff_rib_thickness = adjusted_rib_thickness / ( 
        math.cos(outer_twist_angle)**2 + math.sin(outer_twist_angle)**2 * math.sin(flare_angle)**2 )**0.5

    theta_rib_inner = 2*math.atan(inner_eff_rib_thickness/(2*r_inner))
    theta_rib_outer = 2*math.atan(outer_eff_rib_thickness/(2*r_outer))

    for j in range(n_perp):
        cw_inner = theta_per_channel - theta_rib_inner
        cw_outer = theta_per_channel - theta_rib_outer

        angle_inner = theta - cw_inner/2 + cw_inner * j / (n_perp-1)
        angle_outer = theta - cw_outer/2 + cw_outer * j / (n_perp-1)

        inner_surface_pts[j].append(cq.Vector(
            r_inner*math.sin(angle_inner), r_inner*math.cos(angle_inner), z
        ))

        outer_surface_pts[j].append(cq.Vector(
            r_outer*math.sin(angle_outer), r_outer*math.cos(angle_outer), z
        ))

    theta += dtheta

print("Making surfaces...")
e1 = cq.Workplane("XY").spline(inner_surface_pts[0])
e2 = cq.Workplane("XY").spline(inner_surface_pts[-1])
e3 = cq.Workplane("XY").spline(outer_surface_pts[-1])
e4 = cq.Workplane("XY").spline(outer_surface_pts[0])

s14 = cq.Face.makeRuledSurface(e1.val(), e4.val())
s23 = cq.Face.makeRuledSurface(e2.val(), e3.val())

s12 = cq.Face.makeSplineApprox(points=inner_surface_pts, tol=1e-8)
s34 = cq.Face.makeSplineApprox(points=outer_surface_pts, tol=1e-8)

fb = (
    cq.Workplane("XY")
    .spline([x[0] for x in inner_surface_pts])
    .lineTo(*unpack_vec(outer_surface_pts[-1][0]))
    .spline([x[0] for x in outer_surface_pts][::-1])
    .close()
)

ft = (
    cq.Workplane("XY", origin=(0,0,contour[-1][0]))
    .spline([unpack_vec(x[-1], planar=True) for x in inner_surface_pts])
    .lineTo(*unpack_vec(outer_surface_pts[-1][-1], planar=True))
    .spline([unpack_vec(x[-1], planar=True) for x in outer_surface_pts][::-1])
    .close()
)

sb = cq.Face.makeFromWires(fb.val())
st = cq.Face.makeFromWires(ft.val())

sh = cq.Shell.makeShell([sb, st, s12, s34, s14, s23])

try:
    rv = cq.Solid.makeSolid(sh)
except TypeError:
    raise RuntimeError("Surface stiching failed, try increasing nperp.")

if not args.full:
    print("Exporting regen channel geometry...")
    cq.exporters.export(rv, args.out_file)
else:
    # Circular pattern the channels
    print("Performing circular pattern...")
    bodies = [rv]
    for i in tqdm(range(1, num_channels)):
        bodies.append(rv.rotate((0,0,0), (0,0,1), 360/num_channels*i))

    # Revolve the contour
    inner_pts = [[x[1], x[0]] for x in contour]
    outer_pts = [[x[1]+inner_wall_thickness+channel_height+outer_wall_thickness, x[0]] for x in contour]

    shell = (cq.Workplane("XZ")
            .spline(inner_pts)
            .lineTo(*outer_pts[-1])
            .spline(outer_pts[::-1])
            .close()
            .revolve(360, (0,0,0), (0,1,0))
    )

    # Cut out the channels
    print("Cutting channels...")
    for x in tqdm(bodies):
        shell = shell.cut(x)

    # Save the result
    print("Exporting FULL geometry...")
    cq.exporters.export(shell, args.out_file)


elapsed_time = time.time() - start_time
print()
print("Geometry exported successfully.")
print(f"Regen geometry generated in {elapsed_time:.2f}s.")