Hutchinson-Tech

import sympy as sp import numpy as np import math from numpy import linalg as la from matplotlib import pyplot as plt NUMERIC = False if NUMERIC: sin = np.sin cos = np.cos sqrt = np.sqrt pi = np.pi Ri = 320 # Internal radius, mm Re = 370 # External radius, mm delta = 19 # Rotation angle applied P = 30 # Pression, MPa ai = 29 ae = 78 # External angle bi = 49.2 be = 55 # Cutting external angle else: sin = sp.sin cos = sp.cos sqrt = sp.sqrt pi = sp.pi Ri = sp.symbols("Ri") Re = sp.symbols("Re") x, y = sp.symbols("x y") # Position of a point, bellow X, Y = sp.symbols("X Y") # Position of a external point delta = sp.symbols("d") # Delta, the angle variaton P = sp.symbols("P") # Pression ae, be = sp.symbols("ae be")

t = sp.symbols("t") # The linear value phi = sp.symbols("phi") C0 = np.zeros(3, dtype="object") Cf = np.zeros(3, dtype="object") R = np.zeros((3, 3), dtype="object") # Rotation Matrix

if NUMERIC: # Transform degrees to radians ai *= pi/180 ae *= pi/180 bi *= pi/180 be *= pi/180 delta *= pi/180 xint = Ri * sin(ai) yint = Ri * cos(ai) xext = Ri * sin(ae) yext = Ri * cos(ae) Xint = Ri*cos(bi)*sin(ai-bi)+sin(bi)*sqrt(Re**2 - Ri**2 * sin(ai-bi)**2) Yint = sqrt(Re**2 - Xint**2) Xext = Ri*cos(be)*sin(ae-be)+sin(be)*sqrt(Re**2 - Ri**2 * sin(ae-be)**2) Yext = sqrt(Re**2-Xext**2) new_ai = np.arcsin(Xint/Re) new_ae = np.arcsin(Xext/Re) x, y = xext, yext X, Y = Xext, Yext

if NUMERIC: # Here we are going to plot theta_int = np.linspace(ai, ae, 1000) theta_ext = np.linspace(new_ai, new_ae, 1000) plot_xint = Ri*np.sin(theta_int) plot_yint = Ri*np.cos(theta_int) plot_xext = Re*np.sin(theta_ext) plot_yext = Re*np.cos(theta_ext) plt.plot(plot_xint, plot_yint, color="b") plt.plot(plot_xext, plot_yext, color="b") plt.plot((plot_xint[0], plot_xext[0]), (plot_yint[0], plot_yext[0]), color="b") plt.plot((plot_xint[-1], plot_xext[-1]), (plot_yint[-1], plot_yext[-1]), color="b") theta_int = np.linspace(-ai, -ae, 1000) theta_ext = np.linspace(-new_ai, -new_ae, 1000) plot_xint = Ri*np.sin(theta_int) plot_yint = Ri*np.cos(theta_int) plot_xext = Re*np.sin(theta_ext) plot_yext = Re*np.cos(theta_ext) plt.plot(plot_xint, plot_yint, color="b") plt.plot(plot_xext, plot_yext, color="b") plt.plot((plot_xint[0], plot_xext[0]), (plot_yint[0], plot_yext[0]), color="b") plt.plot((plot_xint[-1], plot_xext[-1]), (plot_yint[-1], plot_yext[-1]), color="b") plt.ylim(bottom=0) axis = plt.gca() axis.set_aspect('equal', adjustable='box') plt.show()

if NUMERIC: # Here we are going to plot theta_int = np.linspace(ai+delta, ae+delta, 1000) theta_ext = np.linspace(new_ai, new_ae, 1000) plot_xint = Ri*np.sin(theta_int) plot_yint = Ri*np.cos(theta_int) plot_xext = Re*np.sin(theta_ext) plot_yext = Re*np.cos(theta_ext) plt.plot(plot_xint, plot_yint, color="r") plt.plot(plot_xext, plot_yext, color="r") plt.plot((plot_xint[0], plot_xext[0]), (plot_yint[0], plot_yext[0]), color="r") plt.plot((plot_xint[-1], plot_xext[-1]), (plot_yint[-1], plot_yext[-1]), color="r") theta_int = np.linspace(-ai+delta, -ae+delta, 1000) theta_ext = np.linspace(-new_ai, -new_ae, 1000) plot_xint = Ri*np.sin(theta_int) plot_yint = Ri*np.cos(theta_int) plot_xext = Re*np.sin(theta_ext) plot_yext = Re*np.cos(theta_ext) plt.plot(plot_xint, plot_yint, color="r") plt.plot(plot_xext, plot_yext, color="r") plt.plot((plot_xint[0], plot_xext[0]), (plot_yint[0], plot_yext[0]), color="r") plt.plot((plot_xint[-1], plot_xext[-1]), (plot_yint[-1], plot_yext[-1]), color="r") axis = plt.gca() axis.set_aspect('equal', adjustable='box') plt.grid() plt.show()

C0[0] = x*sp.cos(phi) C0[1] = y C0[2] = x*sp.sin(phi) Cf[0] = X*sp.cos(phi) Cf[1] = Y Cf[2] = X*sp.sin(phi) R[0, 0] = sp.cos(delta) R[0, 1] = sp.sin(delta) R[1, 0] = -sp.sin(delta) R[1, 1] = sp.cos(delta) R[2, 2] = 1 Cr = R @ C0 print(" [", C0[0]) print("C0 = [", C0[1]) print(" [", C0[2]) print("") print(" [", Cf[0]) print("Cf = [", Cf[1]) print(" [", Cf[2]) print("") print(" [", Cr[0]) print("Cr = [", Cr[1]) print(" [", Cr[2])

     [ x*cos(phi)
C0 = [ y
     [ x*sin(phi)

     [ X*cos(phi)
Cf = [ Y
     [ X*sin(phi)

     [ x*cos(d)*cos(phi) + y*sin(d)
Cr = [ -x*sin(d)*cos(phi) + y*cos(d)
     [ x*sin(phi)

C = (1-t)*Cr + t*Cf print(" [" + str(C[0])) print(" C = [" + str(C[1])) print(" [" + str(C[2])) print("") dCdp = sp.diff(sp.Matrix(C), phi) dCdp = np.array(dCdp).reshape(3) dCdt = sp.diff(sp.Matrix(C), t) dCdt = np.array(dCdt).reshape(3) print(" dC [" + str(dCdp[0])) print(" ---- = [" + str(dCdp[1])) print(" dphi [" + str(dCdp[2])) print("") print(" dC [" + str(dCdt[0])) print(" ---- = [" + str(dCdt[1])) print(" dt [" + str(dCdt[2])) print("") N = np.cross(dCdt, dCdp) for i in range(3): N[i] = sp.simplify(N[i]) print(" [" + str(N[0])) print(" N = [" + str(N[1])) print(" [" + str(N[2])) print("")

         [X*t*cos(phi) + (1 - t)*(x*cos(d)*cos(phi) + y*sin(d))
   C  =  [Y*t + (1 - t)*(-x*sin(d)*cos(phi) + y*cos(d))
         [X*t*sin(phi) + x*(1 - t)*sin(phi)

  dC     [-X*t*sin(phi) - x*(1 - t)*sin(phi)*cos(d)
 ---- =  [x*(1 - t)*sin(d)*sin(phi)
 dphi    [X*t*cos(phi) + x*(1 - t)*cos(phi)

  dC     [X*cos(phi) - x*cos(d)*cos(phi) - y*sin(d)
 ---- =  [Y + x*sin(d)*cos(phi) - y*cos(d)
  dt     [X*sin(phi) - x*sin(phi)

         [x*(X - x)*(t - 1)*sin(d)*sin(phi)**2 + (X*t - x*(t - 1))*(Y + x*sin(d)*cos(phi) - y*cos(d))*cos(phi)
   N  =  [(-X + x)*(X*t - x*(t - 1)*cos(d))*sin(phi)**2 + (X*t - x*(t - 1))*(-X*cos(phi) + x*cos(d)*cos(phi) + y*sin(d))*cos(phi)
         [(X*Y*t - X*t*y*cos(d) + X*x*sin(d)*cos(phi) - Y*t*x*cos(d) + Y*x*cos(d) + t*x*y - x*y)*sin(phi)

N = sp.Matrix(N) N = sp.expand(N) Fr = sp.integrate(-P*N, (t, 0, 1)) Fr = sp.integrate(Fr, (phi, 0, sp.pi)) # Fr = Fr.evalf() print(" [" + str(Fr[0])) print(" Fr = [" + str(sp.factor(Fr[1]))) print(" [" + str(sp.factor(Fr[2])))

         [-pi*P*x**2*sin(d)/2
  Fr  =  [-pi*P*(-X**2 + x**2*cos(d))/2
         [P*(-X*Y + X*y*cos(d) - Y*x*cos(d) + x*y)

if NUMERIC: N = np.array(N).reshape(3) for i in range(3): N[i] = N[i].subs(phi, 0) N[i] = N[i].subs(t, 0.5) N[i] = sp.simplify(N[i]) N[i] = float(N[i]) n = -np.copy(N)/la.norm(np.array(N)) print("N = ") print(N) print("n = ") print(n) angle_radians = np.arctan2(n[1], n[0]) angle_degrees = angle_radians * 180/pi print("angle = ", angle_degrees, " degrees")

if NUMERIC: Fr = np.array(Fr).reshape(3) Fr /= 1000 # Transform N into kN print(" [%.3f kN]" % Fr[0]) print("Fr = [%.3f kN]" % Fr[1]) print(" [%.3f kN]" % Fr[2])