dsolve() definitely needs to be improved here. However, since your system seems to be entirely decoupled, you can just pass each ODE to dsolve separately.
Aaron Meurer On Fri, May 24, 2019 at 9:56 AM pull_over93 <[email protected]> wrote: > > I'm a mechanical engineering student and I'm trying to evaluate the > deflection of a shaft by the elastic line method. I'm using sympy because > i've already evaluate all internal action symbolically. > > Elastic line method involve a system of equation which is composed by first > and second order ODEs. First I've tried this code: > > > 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) > > > > This does not work and I got this error message: > > --------------------------------------------------------------------------- > 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 > > > > > What I can do? > > -- > You received this message because you are subscribed to the Google Groups > "sympy" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to [email protected]. > To post to this group, send email to [email protected]. > Visit this group at https://groups.google.com/group/sympy. > To view this discussion on the web visit > https://groups.google.com/d/msgid/sympy/f6279e74-48a7-4841-846b-752c93ffcaf0%40googlegroups.com. > For more options, visit https://groups.google.com/d/optout. -- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To post to this group, send email to [email protected]. Visit this group at https://groups.google.com/group/sympy. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6%2BDq6yW47DVzB1EE%3DBwjXHcupaJa-2YfwVx9wx5FV%2Bp8g%40mail.gmail.com. For more options, visit https://groups.google.com/d/optout.
