ВведениеНаша страна обладает значительными запасами лесных ресурсов, общая доля которых составляет около четвертой части всех ресурсов сосредоточенных на планете. Важным технологическим звеном процесса лесозаготовок, на которое приходится наибольшая энергоемкость, является процесс вывозки лесоматериалов с лесосеки до потребителя лесовозными автомобилями (ЛА), выполняемый по существующим лесовозным дорогам. Удельный вес доли расходов, затрачиваемых на вывозку лесоматериалов ЛАот суммарных затрат, расходуемых на обеспечение функционирования лесного хозяйства, составляет по разным оценкам 30-35 %. Издержки ЛА отрицательно отражаются на конечной себестоимости продукции и экономических результатах деятельности лесозаготовительных предприятий. Поэтому снижение расходов на вывозку лесоматериалов ЛА является важным резервом повышения эффективности лесозаготовок в РФ.Эксплуатируемый в настоящее время на вывозке лесоматериалов лесовозный автомобильный транспорт удовлетворяют почти всем предъявляемым к ним потребителями требованиям. Несмотря на это вывозка лесоматериалов лесовозным автомобильным транспортом, оснащенным традиционными конструкциями зависимых подвесок в условиях лесовозных дорог низкого качества, приводит к существенному ухудшению их эксплуатационных свойств, проявляющихся в следующем: – в ограниченных значениях дорожного просвета и углов поворота управляемых колес ЛА;в снижении надежности, проходимости ЛА в сложных дорожных условиях за счет контакта выступающих узлов ЛА с лесовозной дорогой; в невозможности использования подрамного пространства, располагающегося между колесами ЛА;– в значительной высоте центра тяжести ЛА, приводящей к ухудшению поперечной устойчивости загруженного ЛА; в отсутствии механизмов, позволяющих регулировать высоту дорожного просвета ЛА; в невозможности использования подвесок, выполняемых в виде колесных модулей из-за больших габаритных размеров и сложности конструктивного исполнения;– отсутствие функциональной возможности вывешивания колес ЛА при его движении в порожнем состоянии;не способность обеспечить индивидуальный поворот каждого из колес ЛА.Повышение эксплуатационных свойств подвески играет особенно определяющее значение для ЛА, применяемого для вывозки лесоматериалов в РФ. Точность расчета конструктивных параметров подрессоривания, а также подбор оптимально подходящей для ЛА конструкции подвески дает возможность уменьшить количество денежных средств расходуемых на вывозку лесоматериалов[1].DemyanovD. N., KarabtsevV. S. и MalinovA. I. (2019) в своих исследованиях с целью повышения эффективности функционирования автомобиля с колесной формулой 8 × 8 разработали алгоритмы управления его гидропневматической подвеской в среде Simulink. Выполненная ими двухэтапная проверка позволила выявить, что работа полученных алгоритмов удовлетворяет всем предъявляемым к ним требованиям, опускание и подъем рамы автомобиля в процессе преодоления неровностей дороги выполняется равномерно, с низким коэффициентом погрешности [3].YesilyrtH. Y. и AkalinO. (2018) в своей работе акцентируют внимание на моделировании системы гидропневматической подвески при движении 10 тонного автомобиля по дорожным неровностям. Авторами были скорректированы параметры гидропневматической подвески с учетом коэффициента демпфирования кузова автомобиля. Выявлено, что предельная величина амплитуды вертикального перемещения сиденья водителя может быть нарушена в случае совпадения крейсерской скорости и колесной базы автомобиля с собственной частотой кузова автомобиля на низших скоростях. Также определено, что увеличение средней скорости автомобиля по неровной поверхности приводит к повышению амплитуды вертикального ускорения сиденья водителя [4].В статье LinD., YangF. и др. (2020) рассмотрены динамические свойства компактной гидропневматической амортизационной стойки автомобиля, разработана методика комплексного эксперимента для исследования зависимостей силы трения от рабочего давления в стойке, силы смещения и давления, получена математическая модель амортизационной стойки, включающая процесс политропного газа, нелинейную и гистерезисную силу трения, а также турбулентные потоки через отверстия поршня. Выявлено, что разработанная математическая модель может с высокой точностью выполнять прогнозирование характеристик амортизационной стойки в динамике [5].В работе KoniecznyL. и др. (2016) выполнен анализ конструктивных решений и материалов, применяемых в газовых пружинах гидропневматических подвесок автомобилей. Авторами выявлены основные преимущества использования таких подвесок в легковых автомобилях, ими приведены фундаментальные зависимости, определяющие характеристики газовой пружины, а также описаны свойства материалов, используемых для производства деталей гидропневматической подвески с газовой пружиной [6].В научной статьеYinY., RakhejaS. и др. (2020) приведены результаты исследования зависящих от температуры и частоты демпфирующие свойства компактной гидропневматической амортизационной стойки. Разработана математическая модель, результаты расчета которой при различных температурах стойки и частотах возбуждения показали высокую сходимость с экспериментальными данными [7].IbrahimA. Badwayи др. (2017) в своих исследованиях изучили рабочие характеристики и потенциал гидропневматической подвески, используемой в транспортном средстве, перемещающемся по дороге с неровностями. Разработанные для исследования математическая модель, а также испытательный стенд позволили выявить, что гидропневматическая подвеска обеспечивает комфорт и эффективное демпфирование при движении транспортного средства по дороге с неровностями [8].  аб вга, б, в – кинематическая схема гидропневматической подвески, ее виды сбоку, сзади и сверху; г – гидравлическая схема;1 – колесо; 2 – гидравлический двигатель; 3 – пневмогидравлический  цилиндр; 4 – амортизаторы; 5 – гидроцилиндры; 6 – ползуны; 7 – направляющие; 8 – кронштейны; 9 – рама модуля; 10 – рама автомобиля; 11 – шаговый гидродвигатель; 12 – гофрированный защитный кожух; 13 – подвижный пневмоцилиндр; 14 – поршень; 15 – внутренние шлицы; 16 – внешние шлицы; 17 – шарики;18 –радиальные подшипники; 19 –упорный шарикоподшипник; 20 – зубчатый сектор; 21 – шестерня; 22 – верхний конец штока; 23 – гидроцилиндр; 24 – буфер сжатия; 25-29 – трубопроводы; 30 – гидрораспределитель; 31, 32 – предохранительные клапаны; 33 – гидронасос; 34, 35 – обратные гидравлические клапаны; 36 – пневмогидравлический аккумулятор; 37 – гидробакРис. 1 – Схема гидропневматической подвески модульного типа для ЛА(собственные разработки)  Выполненный анализ существующих работ показал, что в настоящий период времени отсутствуют конструкции гидропневматических подвесок модульного типа для ЛА, обеспечивающие рекуперацию гидравлической энергии при преодолении неровностей и препятствий лесовозной дороги, компенсацию утечек рабочей жидкости, а также дистанционное управление из кабины ЛА величиной дорожного просвета. С целью устранения рассмотренных выше недостатков существующих конструкций подвесок, авторами была предложена перспективная схема модульной гидропневматической подвески для ЛА (рис. 1) [2].Материалы и методыДля поиска наилучших параметров предлагаемой конструкции гидропневматической подвески модульного типа необходимо выполнить двухфакторную оптимизацию, позволяющую исследовать одновременное влияние двух параметров на показатели эффективности, а также выявить благоприятные комбинации двух факторов, которые не могут быть явно определены в результате исследования соответствующих однофакторных зависимостей[9-12].Для проведения исследования использовалась ранее разработанная математическая модель движения лесовозного автомобиля, оснащенного гидропневматической подвеской [1]. Из большого количества конструктивных параметров гидропневматической подвески ЛАвыбраны два параметра, наиболее существенно влияющие на эффективность: диаметр пневматического цилиндра D и коэффициент теплоотдачи kто. В качестве показателей эффективности гидропневматической подвески ЛАвыбраны: максимальное вертикальное ускорение azm подрессоренной массы, максимальное давление Pm в пневмоцилиндре, максимальная амплитуда ΔTm колебаний температуры газа в пневмоцилиндре. Под максимальными понимаются максимальные значения показателей при движении исследуемой механической системы по участку случайной неровной поверхности длиной 60 м на протяжении 20 с горизонтальной скоростью 3 м/с.В этом случае задачу оптимизации конструктивных параметров можно записать следующим образом:   (1) Для решения такой задачи оптимизации выполняется поиск такого интервала изменения исследуемых факторов D и kто, в которой одновременно будут минимальными вертикальное ускорение azm, давление Pm и амплитуда колебания температуры ΔTm газа в пневмоцилиндре.Для решения двухфакторной задачи оптимизации провели девять компьютерных экспериментов, в которых изменяли диаметр пневмоцилиндраD на уровнях 0,2, 0,3, 0,4 м при одновременном изменении коэффициента теплоотдачи kто на уровнях 0, 1, 2 с–1. Исходя из физического смысла, коэффициента теплоотдачи, расширенной единицей измерения является К/(с·К), то есть, на сколько градусов Кельвина изменяется температура газа в пневмоцилиндре за одну секунду при разнице температур газа и окружающей среды в один градус Кельвина.На основе набора данных, полученного в результате проведения девяти оптимизационных компьютерных экспериментов, получены аналитические формулы для функций azm(D, kто), Pm(D, kто) и ΔTm(D, kто). Данные формулы имеют высокую практическую ценность на этапе конструирования гидропневматической подвескиЛА, так как позволяют оценить ее показатели эффективности, не используя длительных компьютерных расчетов. С целью нахождения аналитических зависимостей функции azm(D, kто), Pm(D, kто) и ΔTm(D, kто) выполнили аппроксимацию полиномов второго порядка: K(D, kто) = k1 D2 + k2 kто2 ++k3 D∙kто + k4 D + k5 kто + k6,     (2) где K– искомый критерий оптимизации; k1 ... k6 – значения параметров полинома.Для нахождения значений параметров полиномов использовался метод наименьших квадратов. Данный метод позволяет минимизировать сумму квадратов возможных отклонений аналитической зависимости от полученных результатов компьютерного моделирования:  (3) где NКЭ – общее количествопроводимых в программе компьютерных экспериментов;i– номер проводимого в программе компьютерного эксперимента;Kаналит. – полиноминальная зависимость влияния на критерий K факторов D и kто; KiКЭ – рассчитанное в программе в i-м компьютерном эксперименте с параметрами Di и kтоiзначение критерия K.Выполнение аппроксимации методом наименьших квадратов осуществлялось в программе для инженерных расчетов MathCAD 14. Использование этой программы дало возможность выявить такие аналитические зависимости, которые позволяют определять значения показателей эффективности гидропневматической подвески ЛА: azm(D, kто) = 0,052 kто2 – 0,273kто + 3,15;     (4)Pm(D, kто) = 210,7D2 + 0,057·10–3kто2 + 1,10D∙kто– 171,3D – 0,620kто + 37,80;        (5) ΔTm(D, kто) = 1,00 kто2 – 11,60kто + 179,4;   (6) где D– диаметр пневматического цилиндра, м; kто – коэффициент теплоотдачи, с-1;;azm –минимальное вертикальное ускорение, м/с2;Pm – давление, МПа; ΔTm – амплитуда колебаний температуры, К.Результаты и обсужденияВыявленные зависимости позволяют с высокой скоростью выполнить оценку показателей эффективности гидропневматической подвескиЛА. Для дальнейшего анализа функции azm(D, kто), Pm(D, kто) и ΔTm(D, kто) приведены в виде графиков (рис.2) и картограмм (рис.3). Построенные графики дают возможность визуально оценить характер влияния диаметра пневмоцилиндра и коэффициента теплоотдачи на показатели эффективности гидропневматической подвескиЛА. В частности, обнаружено, что диаметр пневмоцилиндра практически не влияет на azm и ΔTm (рис.2, а, в), а коэффициент теплоотдачи практически не влияет на Pm (рис.2, б).  аб вРис.2 – Графики влияния диаметра пневмоцилиндраD и коэффициента теплоотдачи kтона максимальное вертикальное ускорение подвешенной массы azm (а), максимальное давление газа в пневмоцилиндреPm (б) и амплитуду изменения температуры ΔTm (в)(собственные разработки)  Графики (рис.2) позволяют понять качественный уровень полученных закономерностей, с целью понимания количественного уровня выявленных закономерностей, графики были преобразованы в картограммы, которые являются видом сверху на графики, а также имеют обозначенные поверхности отклика линиями уровня (рис. 3).Если диаметр пневмоцилиндра составляет0,30 м (точка А) и коэффициент теплоотдачи – 1,0 с–1 (точка В), то точка C, являющаяся точкой факторного пространства будет располагаться ближе к линии уровня 2,95. Таким образом, максимальное вертикальное ускорение составит около 2,97 м/с2.На картограммах затемнены области благоприятных (в данном случае минимальных) значений показателей. В результате выполнения разделения на неблагоприятные и благоприятные области пороговые значения показателей эффективности составили:azm = 3,00 м/с2 (рис. 3, а), Pm = 0,5 МПа. (рис. 3, б), ΔTm = 170 К (рис. 3, в), исходя из условий, что благоприятная область будет занимать значительную долю факторного пространства, включать искомые наименьшие значения критерия и граница области будет представлять собой линию уровня на картограмме.  аб вгРис.3 – Картограммы оптимизации конструктивных параметров гидропневматической подвески. Затемнены благоприятные области факторного пространства(собственные разработки) Оптимальная область (рис. 3, г), образуемая в результате пересечения благоприятных областей, позволяет одновременно учитывать требования минимума показателей эффективности.ВыводыНа основе анализа существующих работ в области подвесок транспортных средств, предложена перспективная конструкция модульной гидропневматической подвески для ЛА, позволяющая устранить основные недостатки традиционных конструкций подвесок, а также рекуперировать энергию рабочей жидкости при движении ЛА по лесовозным дорогам низкого качества, компенсировать утечки рабочей жидкости, а также дистанционно из кабины водителя ЛА управлять величиной дорожного просвета.Для оптимизации параметров гидропневматической подвески ЛА в работе была применена двухфакторная оптимизация с анализом наглядных картограмм.В результате двухфакторной оптимизации гидропневматической подвески ЛАобнаружено, что оптимальный диаметр пневмоцилиндра составляет 0,31 ... 0,40 м, оптимальный коэффициент теплоотдачи – 0,9 ... 2,0 с–1. При этом максимальное вертикальное ускорение составит не более 3 м/с2, максимальное давление газа в пневмоцилиндре составит не более 0,5 МПа, максимальная амплитуда колебаний температуры составит не более 170 К.