Я студент машиностроения и пытаюсь оценить прогиб вала методом упругой линии. Я использую sympy, потому что я уже оценил все внутренние действия символически.
Метод упругих линий включает систему уравнений, которая состоит из ОДУ первого и второго порядка. Сначала я попробовал этот код:
import sympy as sym
import numpy as np
# VARIABLES
# elastic line
v1, v2, v3, v4, v5, r1, r2, r3, r4, r5 = sym.symbols('v1, v2, v3, v4, v5, r1, r2, r3, r4, r5', cls = sym.Function) # deflection function
# General data
L1, L2, tau_1, tau_2, motor_power, rpm_motor, service_factor, impact_energy, spur_gear_efficiency = sym.symbols('L1, L2, tau_1, tau_2, motor_power, rpm_motor, service_factor, impact_energy, spur_gear_efficiency', real = True)
# Transmission analysis
omega_1, omega_2, torque_1, torque_2, user_power = sym.symbols('omega_1, omega_2, torque_1, torque_2, user_power', real = True)
# Gear design variables
dp1, dg1, dp2, dg2 = sym.symbols('dp1, dg1, dp2, dg2', real = True)
pressure_angle, dp1_act, dg1_act, dp2_act, dg2_act, gear_module_act = sym.symbols('pressure_angle, dp1_act, dg1_act, dp2_act, dg2_act, gear_module_act', real = True)
dp1_e_act, dp1_D1_act, dp2_e_act, dp2_D1_act, dg1_e_act, dg1_D1_act, dg2_e_act, dg2_D1_act = sym.symbols('dp1_e_act, dp1_D1_act, dp2_e_act, dp2_D1_act, dg1_e_act, dg1_D1_act, dg2_e_act, dg2_D1_act', real = True)
z_p1, z_g1, z_p2, z_g2, mass_p1, mass_g1, mass_p2, mass_g2, g1_width, g2_width = sym.symbols('z_p1, z_g1, z_p2, z_g2, mass_p1, mass_g1, mass_p2, mass_g2, g1_width, g2_width', real = True)
# Bearings
b1_width, b2_width, reliability, service_life = sym.symbols('b1_width, b2_width, reliability, service_life', real = True)
d,D,C0,C,Pu = sym.symbols('d,D,C0,C,Pu', real = True)
# Shaft data & equilibrium & internal loading actions & steel properties
L1, L2, L3, axial_clearance = sym.symbols('L1, L2, L3, axial_clearance', real = True)
Rax1, Raz1, Rdx1, Rdz1, Fr_g1, Ft_g1 = sym.symbols('Rax1, Raz1, Rdx1, Rdz1, Fr_g1, Ft_g1', real = True)
Rax2, Raz2, Rdx2, Rdz2, Fr_g2, Ft_g2 = sym.symbols('Rax2, Raz2, Rdx2, Rdz2, Fr_g2, Ft_g2', real = True)
N1,Tx1,Tz1,Mx1,Mz1,Mt1 = sym.symbols('N1,Tx1,Tz1,Mx1,Mz1,Mt1', real = True)
N2,Tx2,Tz2,Mx2,Mz2,Mt2 = sym.symbols('N2,Tx2,Tz2,Mx2,Mz2,Mt2', real = True)
xi = sym.Symbol('xi', real = True)
yield_strength,tensile_strength, fatigue_strength, young_modulus, shear_modulus, safe_strength, poisson_modulus,thermal_expansion_coeff = sym.symbols('yield_strength,tensile_strength, fatigue_strength, young_modulus, shear_modulus, safe_strength, poisson_modulus,thermal_expansion_coeff', real = True)
d1,d2,d3,d4,d5 = sym.symbols('d1,d2,d3,d4,d5', real = True)
moment_inertia1, moment_inertia2, moment_inertia3, moment_inertia4, moment_inertia5 = sym.symbols('moment_inertia1, moment_inertia2, moment_inertia3, moment_inertia4, moment_inertia5', real = True)
polar_moment1, polar_moment2, polar_moment3, polar_moment4, polar_moment5 = sym.symbols('polar_moment1, polar_moment2, polar_moment3, polar_moment4, polar_moment5', real = True)
# general data
general_data = {
L1 : 60, # [mm]
L2 : 120, # [mm]
tau_1 : 0.3, # [--]
tau_2 : 0.64, # [--]
motor_power : 800, # [W]
service_factor : 1.6, # [--]
impact_energy : 1.6, # [J]
rpm_motor: 1400, # [rev/min]
}
sinusoidal_results = {
axial_clearance : 6, # [mm]
}
gear_general_data = {
pressure_angle : np.deg2rad(20),
gear_module_act : 1.5, # [mm]
dg1_act : 91.5, # [mm]
dg1_e_act : 94.5, # [mm]
dg1_D1_act : 16.0, # [mm]
dp1_act : 27.0, # [mm]
dp1_e_act : 30.0, # [mm]
dp1_D1_act : 8.0, # [mm]
dg2_act : 72.0, # [mm]
dg2_e_act : 75.0, # [mm]
dg2_D1_act : 14.0, # [mm]
dp2_act : 46.5, # [mm]
dp2_e_act : 49.5, # [mm]
dp2_D1_act : 12.0, # [mm]
z_g1 : 61, # [--]
z_p1 : 18, # [--]
z_g2 : 48,
z_p2 : 31,
mass_g1 : 1.22, # [--]
mass_p1 : 0.1, # [--]
mass_g2 : 0.7,
mass_p2 : 0.3,
g1_width : 17.0,# [mm]
g2_width : 17.0, # [mm]
spur_gear_efficiency : 0.9 # [--]
}
bearings_data = {
b1_width: 8, # [mm]
b2_width: 8, # [mm]
reliability: 99, # [%]
service_life: 1000 # [h]
}
shaft_dimensions = {
L1: 0.5*bearings_data[b1_width] + sinusoidal_results[axial_clearance] + 0.5*gear_general_data[g1_width], # [mm]
L2: general_data[L1] - (2*sinusoidal_results[axial_clearance] + 0.5*gear_general_data[g1_width] + 0.5*gear_general_data[g2_width]), # [mm]
L3: 0.5*bearings_data[b2_width] + sinusoidal_results[axial_clearance] + 0.5*gear_general_data[g2_width] # [mm]
}
shaft_diameters = {
d1:12, # [mm]
d2:16, # [mm]
d3:24, # [mm]
d4:14, # [mm]
d5:12 # [mm]
}
m_inertia = {
moment_inertia1: (np.pi*(shaft_diameters[d1]**4))/64, # [mm^4]
moment_inertia2: (np.pi*(shaft_diameters[d2]**4))/64, # [mm^4]
moment_inertia3: (np.pi*(shaft_diameters[d3]**4))/64, # [mm^4]
moment_inertia4: (np.pi*(shaft_diameters[d4]**4))/64, # [mm^4]
moment_inertia5: (np.pi*(shaft_diameters[d5]**4))/64, # [mm^4]
polar_moment1: 0.5*np.pi*(shaft_diameters[d1]*0.5)**4, # [mm^4]
polar_moment2: 0.5*np.pi*(shaft_diameters[d2]*0.5)**4, # [mm^4]
polar_moment3: 0.5*np.pi*(shaft_diameters[d3]*0.5)**4, # [mm^4]
polar_moment4: 0.5*np.pi*(shaft_diameters[d4]*0.5)**4, # [mm^4]
polar_moment5: 0.5*np.pi*(shaft_diameters[d5]*0.5)**4, # [mm^4]
}
steel_prop = {
yield_strength : 785, # [MPa] 30 [mm] bar
tensile_strength : 1080, # [MPa] 30 [mm] bar
safe_strength : 0.75*1080, # [MPa] at 1×10^4 cycles (rotating bending)
fatigue_strength : 0.45*1080, # [MPa] at 2×10^6 cycles (rotating bending)
young_modulus : 210000, # [MPa]
poisson_modulus: 0.3,
shear_modulus : 80000, # [MPa]
thermal_expansion_coeff : 24*10**-6 # [1/°C] https://www.steel-grades.com/Steel-Grades/Carbon-Steel/AISI-E9314.html
}
tau_1_act = gear_general_data[dp1_act]/gear_general_data[dg1_act]
tau_2_act = gear_general_data[dp2_act]/gear_general_data[dg2_act]
center_distance_1 = (gear_general_data[dg1_act] + gear_general_data[dp1_act])/2
center_distance_2 = (gear_general_data[dg2_act] + gear_general_data[dp2_act])/2
omega_gear_1 = tau_1_act*((2*np.pi*general_data[rpm_motor])/60) # [1/s]
omega_gear_2 = tau_2_act*((2*np.pi*general_data[rpm_motor])/60) # [1/s]
user_power_val = gear_general_data[spur_gear_efficiency]*general_data[motor_power] #[W]
user_torque_1 = (user_power_val/omega_gear_1)*1000 # [N*mm]
user_torque_2 = (user_power_val/omega_gear_2)*1000 # [N*mm]
Ft_1 = user_torque_1/(0.5*gear_general_data[dg1_act])
Ft_2 = user_torque_2/(0.5*gear_general_data[dg2_act])
Fr_1 = Ft_1*np.tan(gear_general_data[pressure_angle])
Fr_2 = Ft_2*np.tan(gear_general_data[pressure_angle])
transmission_data = {
omega_1 : omega_gear_1,
omega_2 : omega_gear_2,
user_power: user_power_val,
torque_1: user_torque_1,
torque_2: user_torque_2,
Ft_g1: Ft_1,
Ft_g2: Ft_2,
Fr_g1: Fr_1,
Fr_g2: Fr_2
}
Eq_x_1 = sym.Eq(-Rax1 + Ft_g1 - Rdx1,0)
Eq_z_1 = sym.Eq(-Raz1 + Fr_g1 - Rdz1,0)
Eq_m_x_1 = sym.Eq(Fr_g1*L1 - Rdz1*(L1+L2+L3),0)
Eq_m_z_1 = sym.Eq(-Ft_g1*L1 + Rdx1*(L1+L2+L3),0)
equilibrium_system_slow = sym.solve([Eq_x_1, Eq_z_1, Eq_m_x_1, Eq_m_z_1],[Rax1, Raz1, Rdx1, Rdz1], dict = True)
# Section AB
N1_act_ab = sym.Eq(N1, 0)
Tx1_act_ab = sym.Eq(Tx1, equilibrium_system_slow[0][Rax1])
Tz1_act_ab = sym.Eq(Tz1, -equilibrium_system_slow[0][Raz1])
Mx1_act_ab = sym.Eq(Mx1, -equilibrium_system_slow[0][Raz1]*xi)
Mz1_act_ab = sym.Eq(Mz1, -equilibrium_system_slow[0][Rax1]*xi)
Mt1_act_ab = sym.Eq(Mt1, -torque_1)
# Section BD
N1_act_bd = sym.Eq(N1, 0)
Tx1_act_bd = sym.Eq(Tx1, equilibrium_system_slow[0][Rax1] - Ft_g1 )
Tz1_act_bd = sym.Eq(Tz1, Fr_g1 - equilibrium_system_slow[0][Raz1])
Mx1_act_bd = sym.Eq(Mx1, Fr_g1*(xi - L1) - equilibrium_system_slow[0][Raz1]*xi)
Mz1_act_bd = sym.Eq(Mz1, Ft_g1*(xi - L1) - equilibrium_system_slow[0][Rax1]*xi)
Mt1_act_bd = sym.Eq(Mt1, 0)
M_ab_resultant = sym.sqrt((Mx1_act_ab.args[1].subs(equilibrium_system_slow[0]).subs(shaft_dimensions).subs(transmission_data))**2 + (Mz1_act_ab.args[1].subs(equilibrium_system_slow[0]).subs(shaft_dimensions).subs(transmission_data))**2)
M_bd_resultant = sym.sqrt((Mx1_act_bd.args[1].subs(equilibrium_system_slow[0]).subs(shaft_dimensions).subs(transmission_data))**2 + (Mz1_act_bd.args[1].subs(equilibrium_system_slow[0]).subs(shaft_dimensions).subs(transmission_data))**2)
rotation1_01 = sym.Eq(young_modulus*moment_inertia1*sym.Derivative(r1(xi), xi, 1), M_ab_resultant).subs(steel_prop).subs(m_inertia)
rotation1_12 = sym.Eq(young_modulus*moment_inertia2*sym.Derivative(r2(xi), xi, 1), M_bd_resultant).subs(steel_prop).subs(m_inertia)
rotation1_23 = sym.Eq(young_modulus*moment_inertia3*sym.Derivative(r3(xi), xi, 1), M_bd_resultant).subs(steel_prop).subs(m_inertia)
rotation1_34 = sym.Eq(young_modulus*moment_inertia4*sym.Derivative(r4(xi), xi, 1), M_bd_resultant).subs(steel_prop).subs(m_inertia)
rotation1_45 = sym.Eq(young_modulus*moment_inertia5*sym.Derivative(r5(xi), xi, 1), M_bd_resultant).subs(steel_prop).subs(m_inertia)
deflection1_01 = sym.Eq(young_modulus*moment_inertia1*sym.Derivative(v1(xi), xi, 2), M_ab_resultant).subs(steel_prop).subs(m_inertia)
deflection1_12 = sym.Eq(young_modulus*moment_inertia2*sym.Derivative(v2(xi), xi, 2), M_bd_resultant).subs(steel_prop).subs(m_inertia)
deflection1_23 = sym.Eq(young_modulus*moment_inertia3*sym.Derivative(v3(xi), xi, 2), M_bd_resultant).subs(steel_prop).subs(m_inertia)
deflection1_34 = sym.Eq(young_modulus*moment_inertia4*sym.Derivative(v4(xi), xi, 2), M_bd_resultant).subs(steel_prop).subs(m_inertia)
deflection1_45 = sym.Eq(young_modulus*moment_inertia5*sym.Derivative(v5(xi), xi, 2), M_bd_resultant).subs(steel_prop).subs(m_inertia)
ODE_system = [rotation1_01, rotation1_12, rotation1_23, rotation1_34, rotation1_45, deflection1_01, deflection1_12, deflection1_23, deflection1_34, deflection1_45]
dist_1 = bearings_data[b1_width]*0.5
dist_2 = dist_1 + 6 + 0.5*gear_general_data[g1_width]
dist_3 = dist_2 + (shaft_dimensions[L2] - gear_general_data[g2_width])
dist_4 = dist_3 + 6 + gear_general_data[g2_width]
dist_5 = dist_4 + bearings_data[b2_width]*0.5
ics = {v1(0):0, v1(dist_1) : v2(dist_1), v2(dist_2):v3(dist_2), v3(dist_3):v4(dist_3), v4(dist_4):v5(dist_4), v5(dist_5):0, r1(dist_1) : r2(dist_1), r2(dist_2):r3(dist_2), r3(dist_3):r4(dist_3), r4(dist_4):r5(dist_4)}
elastic_line_conf_1 = sym.dsolve(ODE_system, [r1(xi), r2(xi), r3(xi), r4(xi), r5(xi), v1(xi), v2(xi), v3(xi), v4(xi), v5(xi)], ics=ics)
Это не работает, и я получил это сообщение об ошибке:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-86-26e73208c1d3> in <module>
27 ics = {v1(0):0, v1(dist_1) : v2(dist_1), v2(dist_2):v3(dist_2), v3(dist_3):v4(dist_3), v4(dist_4):v5(dist_4), v5(dist_5):0, r1(dist_1) : r2(dist_1), r2(dist_2):r3(dist_2), r3(dist_3):r4(dist_3), r4(dist_4):r5(dist_4)}
28
---> 29 elastic_line_conf_1 = sym.dsolve(ODE_system, [r1(xi), r2(xi), r3(xi), r4(xi), r5(xi), v1(xi), v2(xi), v3(xi), v4(xi), v5(xi)], ics=ics)
c:\python37\lib\site-packages\sympy\solvers\ode.py in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)
599 match['eq'] = eq
600 if len(set(order.values()))!=1:
--> 601 raise ValueError("It solves only those systems of equations whose orders are equal")
602 match['order'] = list(order.values())[0]
603 def recur_len(l):
ValueError: It solves only those systems of equations whose orders are equal
Что я могу сделать?