Диссертация (1173093), страница 21
Текст из файла (страница 21)
last ( t )Вычисление производных. dξ1 i ξ i := Dt i , ddξ i dξ i 159Приложение ВПрограмма расчета АЧХ колебаний масс транспортного средства,оснащенного: КВП взамен подвески;стандартной подвеской и традиционными колесамиm⋅d2dtm0 ⋅d2ξ ( t)2dt2ψddξ ( t ) − ψ ( t ) + 2 ⋅ Cb ⋅ ( ξ ( t ) − ψ ( t ) )dt dt+ 2η B ⋅ 0ddddψ ( t) −q ( t) + C s ⋅ ( ψ ( t) − q ( t)) − η B ⋅ ξ ( t) −ψ (t) − C b ⋅ ( ξ ( t) −dtdt dt dt( t ) + η sh ⋅ ψ( t) )0Преобразование Лапласа при нулевых начальных условиях в соответствии с правилами:2m ⋅ p ⋅ ξ ( p ) + 2η B ⋅ ( p ⋅ ξ ( p ) − p ⋅ ψ ( p ) ) + 2 ⋅ C b ⋅ ( ξ ( p ) − ψ ( p ) )02m0 ⋅ p ⋅ ψ ( p ) + η sh ⋅ ( p ⋅ ψ ( p ) − p ⋅ q ( p ) ) + C s ⋅ ( ψ ( p ) − q ( p ) ) − η B ⋅ ( p ⋅ ξ ( p ) − p ⋅ ψ ( p ) ) − C b ⋅ ( ξ ( p ) − ψ ( p ) )0Раскрытие скобок.2m ⋅ p ⋅ ξ ( p ) + 2η B ⋅ p ⋅ ξ ( p ) − 2η B ⋅ p ⋅ ψ ( p ) + 2 ⋅ C b ⋅ ξ ( p ) − 2 ⋅ C b ⋅ ψ ( p )02m0 ⋅ p ⋅ ψ ( p ) + η sh ⋅ p ⋅ ψ ( p ) − η sh ⋅ p ⋅ q ( p ) + C s ⋅ ψ ( p ) − C s ⋅ q ( p ) − η B ⋅ p ⋅ ξ ( p ) + η B ⋅ p ⋅ ψ ( p ) − C b ⋅ ξ ( p ) + C b ⋅ ψ ( p )0Группирование подобных при переменных.ξ (p )(m ⋅ p2()(+ 2η B ⋅ p + 2 ⋅ C b − ψ ( p ) ⋅ 2η B ⋅ p + 2 ⋅ C b())0)2(− ξ ( p ) η B ⋅ p + C b + ψ ( p ) ⋅ m0 ⋅ p + η sh ⋅ p + C s + η B ⋅ p + C b − q ( p ) ⋅ η sh ⋅ p + C s)0Перенос свободных членов в правую часть.ξ (p )(m ⋅ p)− ( p ) ( B ⋅ p + Cb) + ( p ) ⋅ ( m0 ⋅ pξη2(+ 2η B ⋅ p + 2 ⋅ C b − ψ ( p ) ⋅ 2η B ⋅ p + 2 ⋅ C bψ2)0)+ η sh ⋅ p + Cs + η B ⋅ p + Cb()q ( p ) ⋅ η sh ⋅ p + Cs160Запись системы уравнений в матричной форме. m ⋅ p 2 + 2η ⋅ p + 2 ⋅ C−( 2η B ⋅ p + 2 ⋅ Cb)Bb ⋅ ξ ( p ) 2−( η B ⋅ p + Cb)m0 ⋅ p + p ⋅ ( η sh + η B) + Cb + Cs ψ ( p ) 0q (p ) ⋅ η ⋅ p + C ( shs) () m ⋅ p 2 + 2η ⋅ p + 2 ⋅ C ξ (p ) ⋅ m ⋅ p 2 + 2 ⋅ η ⋅ p + 2 ⋅ C − ψ (p ) ⋅ ( 2 ⋅ C + 2 ⋅ p ⋅ η ) − ( 2η B ⋅ p + 2 ⋅ C b)BbBbbB ⋅ ξ ( p ) → 22ψ(p)− ( η B ⋅ p + C b)m0 ⋅ p + p ⋅ ( η sh + η B) + Cb + Cs ψ ( p ) ⋅ m0 ⋅ p + ( η B + η sh) ⋅ p + Cb + Cs − ξ ( p ) ⋅ ( Cb + p ⋅ η B) Нахождение решения в символьном виде. m ⋅ p 2 + 2η ⋅ p + 2 ⋅ C− ( 2η B ⋅ p + 2 ⋅ C b)bB2m0 ⋅ p + p ⋅ ( η sh + η B) + C b + C s−( η B ⋅ p + C b)ξ (p )q ( p ) ⋅ ( C s + p ⋅ η sh ) ⋅ ( 2 ⋅ C b + 2 ⋅ p ⋅ η B) 2 ⋅ C ⋅ C + 2 ⋅ C ⋅ p ⋅ η + 2 ⋅ C ⋅ p ⋅ η + C ⋅ m ⋅ p2 + C ⋅ m ⋅ p2 + 2 ⋅ C ⋅ m ⋅ p2 + m ⋅ m ⋅ p4 + m ⋅ p3 ⋅ η + m ⋅ p3 ⋅ η + 2 ⋅ m ⋅ p3 ⋅ η + 2 ⋅ p2 ⋅ η ⋅ η 0bssBbshbsb00Bsh0BB sh ⋅→2 q ( p ) ⋅ ( η sh ⋅ p + C s ) q ( p ) ⋅ ( C s + p ⋅ η sh ) ⋅ m ⋅ p + 2 ⋅ η B ⋅ p + 2 ⋅ C b 2 ⋅ C ⋅ C + 2 ⋅ C ⋅ p ⋅ η + 2 ⋅ C ⋅ p ⋅ η + C ⋅ m ⋅ p2 + C ⋅ m ⋅ p2 + 2 ⋅ C ⋅ m ⋅ p2 + m ⋅ m ⋅ p4 + m ⋅ p3 ⋅ η + m ⋅ p3 ⋅ η + 2 ⋅ m ⋅ p3 ⋅ η + 2 ⋅ p2 ⋅ η ⋅ η bssbsb000BshbBshBB sh ()( C s + p ⋅ η sh) ⋅ ( 2 ⋅ C b + 2 ⋅ p ⋅ η B):=2ψ (p )−1222433323332⋅ C b ⋅ C s + 2 ⋅ C b ⋅ p ⋅ η sh + 2 ⋅ C s ⋅ p ⋅ η B + C b ⋅ m ⋅ p + C s ⋅ m ⋅ p + 2 ⋅ C b ⋅ m0 ⋅ p + m ⋅ m0 ⋅ p + m ⋅ p ⋅ η B + m ⋅ p ⋅ η sh + 2 ⋅ m0 ⋅ p ⋅ η B + 2 ⋅ p ⋅ η B ⋅ η sh( C s + p ⋅ η sh) ⋅ ( m ⋅ p 2 + 2 ⋅ η B ⋅ p + 2 ⋅ C b):=22224⋅ C b ⋅ C s + 2 ⋅ C b ⋅ p ⋅ η sh + 2 ⋅ C s ⋅ p ⋅ η B + C b ⋅ m ⋅ p + C s ⋅ m ⋅ p + 2 ⋅ C b ⋅ m0 ⋅ p + m ⋅ m0 ⋅ p + m ⋅ p ⋅ η B + m ⋅ p ⋅ η sh + 2 ⋅ m0 ⋅ p ⋅ η B + 2 ⋅ p ⋅ η B ⋅ η sh⋅ q (p )⋅ q (p )Таким образом, передаточные функции будут иметь вид.Wξ ( p ) :=ξ (p )q (p)Wψ ( p ) :=ψ (p)q (p )Передаточные функции, связывающие выходные переменные ξ и ψ, и входное воздействие q:W ξ ( p ) :=( C s + p ⋅ η sh) ⋅ ( 2 ⋅ Cb + 2 ⋅ p ⋅ η B)2⋅222433323332C b ⋅ C s + 2 ⋅ C b ⋅ p ⋅ η sh + 2 ⋅ C s ⋅ p ⋅ η B + C b ⋅ m ⋅ p + C s ⋅ m ⋅ p + 2 ⋅ C b ⋅ m0 ⋅ p + m ⋅ m0 ⋅ p + m ⋅ p ⋅ η B + m ⋅ p ⋅ η sh + 2 ⋅ m0 ⋅ p ⋅ η B + 2 ⋅ p ⋅ η B ⋅ η sh( C s + p ⋅ η sh) ⋅ ( m ⋅ p 2 + 2 ⋅ η B ⋅ p + 2 ⋅ C b)W ψ ( p ) :=22224⋅ C b ⋅ C s + 2 ⋅ C b ⋅ p ⋅ η sh + 2 ⋅ C s ⋅ p ⋅ η B + C b ⋅ m ⋅ p + C s ⋅ m ⋅ p + 2 ⋅ C b ⋅ m0 ⋅ p + m ⋅ m0 ⋅ p + m ⋅ p ⋅ η B + m ⋅ p ⋅ η sh + 2 ⋅ m0 ⋅ p ⋅ η B + 2 ⋅ p ⋅ η B ⋅ η shПреобразование передаточных функции к дробно-рациональному виду.2W ξ ( p ) :=2W ψ ( p ) :=()()⋅ C b ⋅ C s + 2 ⋅ C b ⋅ η sh + C s ⋅ η B ⋅ p + 2 ⋅ η B ⋅ η sh ⋅ p()22()3⋅ C b ⋅ C s + 2 ⋅ C b ⋅ η sh + C s ⋅ η B ⋅ p + C b ⋅ m + C s ⋅ m + 2 ⋅ C b ⋅ m0 + 2 ⋅ η B ⋅ η sh ⋅ p + m ⋅ η B + m ⋅ η sh + 2 ⋅ m0 ⋅ η B ⋅ p + m ⋅ m0 ⋅ p3) 22342 ⋅ C b ⋅ C s + 2 ⋅ ( C b ⋅ η sh + C s ⋅ η B) ⋅ p + ( C b ⋅ m + C s ⋅ m + 2 ⋅ C b ⋅ m0 + 2 ⋅ η B ⋅ η sh ) ⋅ p + ( m ⋅ η B + m ⋅ η sh + 2 ⋅ m0 ⋅ η B) ⋅ p + m ⋅ m0 ⋅ p2()(4⋅ C b ⋅ C s + 2 ⋅ C b ⋅ η sh + C s ⋅ η B ⋅ p + C s ⋅ m + 2 ⋅ η B ⋅ η sh ⋅ p + m ⋅ η sh ⋅ pВвод параметров колебательной системы.m := 550Подрессоренная масса, кгm0 := 10Масса обода колеса, кгŋsh := 1260Коэффициент демпфирования шины, Н·с/мŋВ := 1035Коэффициент демпфирования упругого элемента, Н·с/мCb := 330000Коэффициент нормальной жесткости упругого элемента, Н/мCS := 94000Коэффициент нормальной жесткости шины, Н/мωn := 0.01Начальная частота расчета, рад / c161ωk := 0.01Конечная частота расчета, рад / c∆ω := 0.1Шаг расчета по частоте, рад / cq0 := 0,03Половина высоты неровности, мДиапазон расчета по частоте.ω():= ω n , ∆ω + ω n ..
ω kПередаточные функции, связывающие выходные переменные z, ξ, ψ и входное воздейст-вие q (с учетом амплитуды входного воздействия).2W ξ ( p ) :=2W ψ ( p ) :=()()⋅ C b ⋅ C s + 2 ⋅ C b ⋅ η sh + C s ⋅ η B ⋅ p + 2 ⋅ η B ⋅ η sh ⋅ p()22()3⋅ C b ⋅ C s + 2 ⋅ C b ⋅ η sh + C s ⋅ η B ⋅ p + C b ⋅ m + C s ⋅ m + 2 ⋅ C b ⋅ m0 + 2 ⋅ η B ⋅ η sh ⋅ p + m ⋅ η B + m ⋅ η sh + 2 ⋅ m0 ⋅ η B ⋅ p + m ⋅ m0 ⋅ p)(⋅ q03) 2⋅ q02342 ⋅ C b ⋅ C s + 2 ⋅ ( C b ⋅ η sh + C s ⋅ η B) ⋅ p + ( C b ⋅ m + C s ⋅ m + 2 ⋅ C b ⋅ m0 + 2 ⋅ η B ⋅ η sh ) ⋅ p + ( m ⋅ η B + m ⋅ η sh + 2 ⋅ m0 ⋅ η B) ⋅ p + m ⋅ m0 ⋅ p2(4⋅ C b ⋅ C s + 2 ⋅ C b ⋅ η sh + C s ⋅ η B ⋅ p + C s ⋅ m + 2 ⋅ η B ⋅ η sh ⋅ p + m ⋅ η sh ⋅ pWdξ ( p ) := p ⋅ Wξ ( p )Передаточные функции для первой и второй производной переменных.Wddξ ( p ) := p ⋅ Wdξ ( p )Wdψ ( p ) := p ⋅ Wψ ( p )Wddψ ( p ) := p ⋅ Wdψ ( p )Для расчета частотных характеристик вместо оператора Лапласа р подставляется ω·j.Амплитудно-частотные характеристики.Aξ ( ω ) := Wξ ( ω ⋅ j )Aψ ( ω ) := Wψ ( ω ⋅ j )Adξ ( ω ) := Wdξ ( ω ⋅ j )Adψ ( ω ) := Wdψ ( ω ⋅ j )Addξ ( ω ) := Wddξ ( ω ⋅ j )Addψ ( ω ) := Wddψ ( ω ⋅ j )(Lξ ( ω ) := 20 ⋅ log Wξ ( ω ⋅ j ))Логарифмические амплитудно-частотная характеристики.()Ldξ ( ω ) := 20 ⋅ log( Wdξ ( ω ⋅ j ) )Ldψ ( ω ) := 20 ⋅ log( Wdψ ( ω ⋅ j ) )Lddξ ( ω ) := 20 ⋅ log( Wddξ ( ω ⋅ j ) )Lψ ( ω ) := 20 ⋅ log Wψ ( ω ⋅ j )(Lddψ ( ω ) := 20 ⋅ log Wddψ ( ω ⋅ j )():= arg Wξ ( ω ⋅ j ) ⋅ψ ξ (ω ))162180πЛогарифмические амплитудно-фазовые частотные характеристики.ψ ψ (ω )():= arg Wψ ( ω ⋅ j ) ⋅(180π)180ψ dξ ( ω ):= arg Wdξ ( ω ⋅ j ) ⋅πψ dψ ( ω ):= arg Wdψ ( ω ⋅ j ) ⋅π()(180)180ψ ddξ ( ω ):= arg Wddξ ( ω ⋅ j ) ⋅πψ ddψ ( ω ):= arg Wddψ ( ω ⋅ j ) ⋅π()180Проверка правильности нахождения передаточных функций.Расчет реакции объекта на единичную ступенчатую функцию по его частотной характеристике.Входной сигнал – прямоугольный импульс.
Длительность импульса входного сигнала(половина периода) больше, чем время памяти объекта, т.е. того времени по истечении которого переходные процессы можно считать законченными.Период входного сигнала.T := 6Число расчетных точек по времени.nt := 500Шаг расчета по времени.∆t∆t:=Tnt= 0.012Расчет вектора времени.kt := 0 .. ntt 0 := 0t k := kt ⋅ ∆tt163Функция расчета коэффициентов ряда Фурье для прямоугольного импульсного входногосигнала (по аналитическим выражениям).fun_c_step( n ) :=cn ← 0c0 ←12for k ∈ 1 , 3 .. nck ← −j ⋅1⋅π1kcЧисло расчетных точек по частоте.nf := 5000Расчет коэффициентов ряда Фурье для прямоугольного импульсного входного сигнала.( )c := fun_c_step nf1f0 :=Tω0:= 2 ⋅ π ⋅ f0Расчет вектора частот.kf := 0 .. nff k := kf ⋅ f0fФункция вычисления обратного преобразования Фурье.(nf)Inv_Fourier c , t , nf := c 0 + 2 ⋅∑Re c k ⋅ e(j ⋅ k ⋅ ω 0 ⋅ t)k= 1Проверка. Расчет входного сигнала по коэффициентам ряда Фурье.(step k := Inv_Fourier c , t k , nftt)Расчет частотной характеристики исследуемого объекта.→ω := ( 2 ⋅ π ⋅ f )→jω := ( ω ⋅ j )→Hξ := Wξ ( jω )Расчет реакции объекта на входной сигнал (в частотной области).→Yξ := Hξ ⋅ c()164Расчет реакции объекта на прямоугольный импульсный входной сигнал (во временнойобласти).(Y_step_ ξ k := Inv_Fourier Yξ , t k , nftt)Расчет частотной характеристики исследуемого объекта.→ω := ( 2 ⋅ π ⋅ f )→jω := ( ω ⋅ j )→Hψ := Wψ ( jω )Расчет реакции объекта на входной сигнал (в частотной области).→Yψ := Hψ ⋅ c()Расчет реакции объекта на прямоугольный импульсный входной сигнал (во временнойобласти).(Y_step_ ψ k := Inv_Fourier Y ψ , t k , nftt)165Приложение ГПрограмма расчета АЧХ колебаний масс транспортного средства,оснащенного традиционными колесами без подвескиПриведение уравнений к форме Коши.Ввод дополнительной переменной: dξ .dξdt:= dξПреобразование исходного уравнения:M⋅ddtddt( dξ ) + 2 ⋅ η sh ⋅ dξ −( dξ ) := −M 1ddtq ( t ) + 2 ⋅ Csh ⋅ ( ξ − q ) := 0⋅ 2 ⋅ η sh ⋅ dξ −ddtq ( t ) + 2 ⋅ Csh ⋅ ( ξ − q )q ( t ) - входное воздействие;dq ( t ) - производная входного воздействия.Преобразование Лапласа при нулевых начальных условиях в соответствии с правилами:p ⋅ ξ ( p ) := dξ ( p )p ⋅ dξ ( p ) := −1⋅ 2 ⋅ η sh ⋅ ( dξ ( p ) − p ⋅ q ( p ) ) + 2 ⋅ Csh ⋅ ( ξ ( p ) − q ( p ) )M Нахождение изображения неизвестных переменных:ξ (p ):=1pdξ ( p ) := −⋅ dξ ( p )1p⋅1⋅ 2 ⋅ η sh ⋅ ( dξ ( p ) − p ⋅ q ( p ) ) + 2 ⋅ Csh ⋅ ( ξ ( p ) − q ( p ) )M Преобразование уравнения:dξ ( p ) := −1p⋅1⋅ 2 ⋅ η sh ⋅ p ⋅ ξ ( p ) + 2 ⋅ Csh ⋅ ( ξ ( p ) − q ( p ) ) +M 1M⋅ q (p )166p ⋅ ξ ( p ) := −1p1⋅M 0:= −( p ⋅ ξ ( p ) ) +0:= −( p ⋅ ξ ( p ) ) +0:= −p ⋅ ξ ( p ) +1M121+M1M1⋅ q (p) −⋅ q (p) −M2ξ (p )p⋅ 2 ⋅ η sh ⋅ p ⋅ ξ ( p ) + 2 ⋅ Csh ⋅ ( ξ ( p ) − q ( p ) ) +1p⋅p1M⋅ p ⋅ q (p) −⋅ 2 ⋅ η sh ⋅ p +1M⋅1M⋅ 2 ⋅ η sh ⋅ p ⋅ ξ ( p ) −⋅ 2 ⋅ η sh ⋅ p ⋅ ξ ( p ) −1M1p⋅⋅ 2 ⋅ η sh ⋅ p ⋅ ξ ( p ) −⋅ 2 ⋅ Csh := q ( p ) ⋅ 1M1⋅p1M1M1M⋅ q (p)1⋅ 2 ⋅ Csh ⋅ ( ξ ( p ) − q ( p ) )M ⋅ 2 ⋅ Csh ⋅ ξ ( p ) +1p⋅ 2 ⋅ Csh ⋅ ξ ( p ) +⋅p+1M⋅1M1M⋅ 2 ⋅ Csh ⋅ q ( p )⋅ 2 ⋅ Csh ⋅ q ( p )⋅ 2 ⋅ Csh Таким образом, передаточная функция будет иметь вид:Wξ ( p ) :=ξ (p )q (p)1Wξ ( p ) :=M2p +1M⋅p+1M⋅ 2 ⋅ Csh⋅ 2 ⋅ η sh ⋅ p +1M⋅ 2 ⋅ CshВвод параметров колебательной системы.М := 550Подрессоренная масса, кгŋsh := 1260Коэффициент демпфирования шины, Н·с/мCS := 94000Коэффициент нормальной жесткости шины, Н/мωn := 0.01Начальная частота расчета, рад / cωk := 0.01Конечная частота расчета, рад / c∆ω := 0.1Шаг расчета по частоте, рад / cq0 := 0,025Половина высоты неровности, мДиапазон расчета по частотеω():= ω n , ∆ω + ω n ..