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])