Введение Низшая частота изгибных колебаний (ЧИК) кузова пассажирского вагона в вертикальной плоскости является одним из основных параметров, характеризующих динамические качества пассажирского вагона [1]. Поэтому ее значение проверяется как теоретически на стадии проектирования [1], так и экспериментально после изготовления [2]. Существует ряд экспериментальных методов определения низшей ЧИК как при стендовых, так и при ходовых испытаниях. Стендовые испытания основаны на возбуждении колебаний с помощью вибромашины или с помощью точечного ударного нагружения. Использование вибромашины дает значительный разброс результатов из-за эффекта Зоммерфельда [3] и требует по этой причине специальных методов обработки. Ударное загружение из-за неопределенности начальных условий соответствующего переходного процесса, а также его массовых и упруго-диссипативных параметров приводит к результатам, зависящим от способов обработки и от объемов выборки, так что репрезентативные выборки требуют для своей обработки использования суперкомпьютерных технологий [4]. В связи со сказанным, в настоящее время начинают распространяться прямые методы, основанные на гибридных технологиях, когда экспериментальные методы дополняются теоретическим модальным анализом – в основном с применением МКЭ [5-8]. При таких оценках динамики исходную КЭ-систему редуцируют, уменьшая (конденсируя) количество степеней свободы (СС) с использованием ряда критериев, которые регламентируются международными требованиями (NSTC 14046 Revision E. Space Shuttle Program. Payload Verification Requirements. USA, NASA, 2000,60 p.) [9]. В основном редуцирование модели проводится на основе МЕМ-критерия (Modal Effective Mass), в соответствии с которым целевым модам должны соответствовать целевые массы, составляющие не менее 10% от соответствующих элементов матрицы инерции, а вторичным (но тоже удерживаемым) – элементы «5%». При этом погрешность (верификации) целевых мод не должна превышать 3…5%. Кроме того, при проведении конденсации учитывают МАС-критерий (Modal Assurance Criterion) модальной достоверности, который по своей сути – просто квадрат косинуса угла между соответствующими собственными векторами расчетных и экспериментальных форм. При этом для сохранения в КЭ-модели i-й СС необходимо, чтобы в МАС-матрице для i-го диагонального элемента выполнялось неравенство , а для остальных (удерживаемых) . Однако в [10] указывается, что полностью полагаться на конденсацию с помощью указанных (и ряда других) критериев нельзя – из-за неполного набора экспериментов это может привести к потере существенных СС. Предлагается задавать существенные СС в узлах-разделителях, удаление которых приводит к делению конструкции на ряд несвязных фрагментов. Описанные приемы отбраковки теоретических мод следует производить не только при гибридном подходе, но и при чисто теоретических оценках (в этом случае МАС-отбраковка сводится к определению взаимной ориентации целевых мод относительно остальных). Это необходимо потому, что МКЭ-подход дает весьма густой спектр, содержащий множество парциальных частот, не представляющих никакого практического интереса. Например в [9] сообщается, что при проведении экспериментов для частот выше 35 Гц не удалось обнаружить ни одной предсказанной МКЭ-частоты: все они оказались парциальными и в эксперименте не верифицировались. Рассматриваемый гибридный подход является весьма дорогостоящим и требует специального технологического оборудования, а также специальных программных средств [11] и по этой причине малоприемлем для проведения экспресс-оценок динамических характеристик кузова, в частности низшей ЧИК, которая является частотой общего изгибания всей конструкции.  Обсуждение возможных подходов В последующих рассуждениях будем исходить из того, что кузов является объектом, протяженным в одном направлении и, по этой причине, для него в качестве модели может использоваться одномерная балочная модель. Следует отметить, что аналогичный подход достаточно давно известен в ракетостроении [12], где он успешно используется для практических расчетов ([13] – ракета-носитель «Циклон-4», [14] – ракета «Ангара» и др.).Теоретически описание линейной вязко-упругой модели в данном случае имеет вид:                                                                                                                  (1)((') означает дифференцирование по времени), где  - матрицы инерции демпфирования и жесткости соответственно. Знание их значений позволяет предсказывать поведение  (или  объекта. Любые параметры из (1) (а также их комбинации) могут служить объектами идентификации. При этом, в качестве данных для идентификации могут выступать обобщенные координаты из , так и некоторые наблюдаемые переменные  (например ускорения), связанные с  в общем случае матричным соотношением [15]:                                                                                                                                          (2)Проектируя  из (2) на , подставляя  в (1), получим уравнение связи наблюдаемых параметров с вектором состояния. Эта связь будет взаимно-однозначной только в случае совпадения порядков  и , при невырожденной квадратной  того же порядка и при точном знании всех величин. В остальных случаях, а их большинство) проектирование по (2) не является однозначным и его следует понимать в некотором обобщенном смысле, так что неизвестные коэффициенты из (1) не могут быть восстановлены однозначно. Такие обратные коэффициентные задачи [16] можно классифицировать по разному: как количественную интерпретацию наблюдений (2), или распознавания, или диагностики – в зависимости от специфики области применения. В случае линейности (2) задача идентификации для (1) сводится к линейной системе обыкновенных дифференциальных уравнений (СОДУ) или интегральных уравнений, которые, в свою очередь, могут быть сведены к системе линейных алгебраических уравнений (СЛАУ):                                .                                                                                                          (3)Так как матрицы  в (1), (2) в общем случае прямоугольные, а количество наблюдений за  или за  произвольно, то идентификационные подходы позволяют конструировать множество моделей, и, следовательно, порождать множество решений. В случае сведения задач к СЛАУ целесообразным является использование сингулярных разложений [17] – тогда количество учитываемых сингулярных собственных чисел можно трактовать как порядок задачи идентификации (аналог структурной идентификации - [15]). В связи с возможной множественностью решений возникает вопрос об адекватности. Наиболее общий подход в этом случае – это рассматривать идентификацию, как задачу оптимизации, сконструировав подходящий экстремальный функционал – скаляр или свертку нескольких критериев [18]. Самый популярный здесь подход – использование принципа минимума взвешенной квадратичной невязки. Это удобно, т.к. для всякой задачи вида (3) существует обобщенное решение  в смысле метода наименьших квадратов (МНК) [19]. Этим решением является уравнение Эйлера для функционала вида                                                                                              (4)((*) – знак транспонирования),где  - вектор весовых множителей.Это приводит к требованию равенства нулю производной Фреше [18]. Получаемое при этом решение является нормальным в норме Фробениуса , однако, оно может быть неустойчивым [20]. Для исключения неустойчивости применяют различные методы, из которых наиболее популярен метод регуляризации Тихонова [21]. Для задачи вида (3) этот метод предполагает разыскание,где - норма Фробениуса,  - оператор дифференцирования порядка , который при  представляет собой единичную матрицу, что соответствует стандартной регуляризации по Тихонову.Если считать справедливой гипотезу Базиля [22], то (1), в соответствии с правилом Виделера примет вид:                          .                                                                                                    (5)Для однородной СОДУ, соответствующей (5), полагая , получим уравнение для определения круговой частоты  в виде:                        ,                                                                                                 (6)что соответствует классической (при диагональной ) или обобщенной проблеме собственных значений. Получение для кузова достаточно точных значений матрицы инерции  не представляет трудностей, чего нельзя сказать о матрице жесткости . Существует большое количество работ, решающих задачу идентификации для (5) или (6), однако, если речь идет о низших модах, то, вместо (5) можно поставить статическую задачу идентификации:                               ,                                                                                                          (7)из которой и определять матрицу .Уравнения (1), (5) и (7) не накладывают никаких ограничений на вид функциональной зависимости  (и на его ее дискретный аналог ), где  - координата вдоль оси вагона, а  - перемещение вдоль некоторой образующей (например, обвязки). В связи с этим, а также, учитывая сходство низших форм изгибных колебаний кузова с балочными модами, можно при решении задачи идентификации считать, что перемещения  могут быть описаны в рамках балочной модели, в частности в рамках ее простейшего варианта – балки Бернулли-Эйлера. Конечно, при этом следует учитывать, что изгибная жесткость  вагона и  модели кузова в виде балки Бернулли-Эйлера на самом деле различны.  модели должно быть таким, чтобы обеспечить совпадение (в определенном смысле) свойств модели (балка) и прототипа (кузова). Выбрав в качестве родовых свойств идентификации совпадение поперечных прогибов прототипа и модели, будем определять матрицу жесткости  в (7) для модели в виде балки на двух жестких опорах (рис. 1).  Рис 1. Расположение опор на схеме идентификации  Из сравнения очевидного равенства                       с (6) следует, что                    ,     ,            (8)                                                                              (8)где  - число точек идентификации, а  - прогибы  в точке  от единичной силы , приложенной в точке . Поскольку для нахождения  достаточно знать геометрию одномерной модели, а также изгибную жесткость , то можно сделать вывод, что для оценки низшей ЧИК достаточно произвести расчетно-экспериментальную оценку . Практически достаточно, чтобы уравнения (6) – (8) позволяли оценивать третью частоту (для схемы на упругих опорах третья частота будет соответствовать низшей частоте колебаний общего изгибания). По самым пессимистическим прогнозам для ее адекватной оценки достаточно рассмотреть систему с 6 СС, т.е. иметь 6 точек идентификации прогибов. При этом единственным разумным предположением об изменении  является предположение о кусочном постоянстве . В [23] приведено решение задачи в случае нагружения кузова сосредоточенными силами с помощью домкратов в ряде точек и замер возникающих при этом перемещений в ряде других точек. Сравнивая теоретические значения перемещений  вычисленных при определенных значениях  с замеренными перемещениями  и используя МНК для взвешенного функционала вида (4), находим значения изгибной жесткости. При поиске экстремума для (4) выяснилось, что из-за требования положительности  методы косвенной оптимизации становятся неэффективными и целесообразно применять прямые методы. При переходе к ним достаточно часто происходили попадания в локальные ямы. Возникновение множества ям объяснялось в основном тем, что приходилось работать с большими выборками зашумленных данных, объединяющих большое количество замеров. Из-за этого приходилось перезапускать процессы поиска. Возникали проблемы с формулировкой критерия останова. Кроме того, по мере приближения к (гипотетическому) глобальному экстремуму возникали колебания – процесс рассыпался. К сожалению эта трудность для рассматриваемой задачи при многомерном поиске оказалась неустранимой. Это связано с предположением о кусочно-постоянном законе изменения , что в нашем случае не только разумно с точки зрения принципа однородной точности, но и согласуется с действительностью. В [24] доказывается, что в случае идентификации негладких зависимостей (удовлетворяющих, однако, условиям Дирикле) обычные методы, даже при малом уровне шумов не только приводят к осцилляциям результатов, но (в случае непринятия специальных мер) и дают отрицательные нефизические значения идентифицируемых параметров. Кроме того, при превышении усилиями от домкратов определенного значения (зависящего от точек его приложения) возможно нагружение контакта кузова с опорами, так что схема рис. 1 перестает соответствовать действительности [25]. В связи с этим схема нагружения для замеров прогиба была изменена и рассматривалась схема весовой догрузки кузова в соответствии с рис. 2. Рис. 2. Нагрузки, используемые при решении задачи идентификации  В эксперименте это соответствует замерам прогибов от некоторого дополнительного нагружения кузова вертикальной нагрузкой. Для регуляризации полученного решения обратной коэффициентной задачи использовалось два подхода. Первый основывался на корректировке элементов  из (8), найденных с использованием экспериментальных прогибов. Эти корректировки сводились к симметризации недиагональных  (они заменялись на ), а затем к проверке положительной определенности получавшейся при этом матрицы. Для этого, в соответствии с критериями Сильвестра [27], проверялись угловые миноры  матрицы  и, если оказывалось, что , то элемент  наращивался до тех пор, пока положительность не восстанавливалась. Более удачным оказался подход, основанный на возможности представления потенциальной энергии в виде суммы произведений слагаемых , где  однозначно определяется геометрией и схемой нагружения модели [26]. В этом случае оказалось возможным свести задачу идентификации к общей задаче линейного программирования. Получаемое здесь решение зависело от разбиения модели на участки кусочно-постоянной жесткости (при прочих равных условиях). Это связано с тем, что значения  находятся с помощью проектирования параметров наблюдения  в пространство , а оценка качества этого проектирования производится с помощью функционала (4).  Основная часть  38Сделаем попытку реально учесть структуру прототипа в нашей идентификационной модели. Специфика изменения жесткостей по длине кузова состоит в том, что на нем обычно присутствует регулярная часть, в которой чередуются два элемента: оконный проем и простенок, что соответствует чередованию двух значений жесткости. Это означает, что количество перемен жесткости для кузова намного превышает количество несовпадающих значений этой жесткости. В связи с этим оказывается целесообразным создание описателей чередования групп жесткости в виде числовой последовательности:                                                                                             (9)где  - индекс группы жесткости для -го участка схемы по рис. 2, а сами изгибные жесткости  образуют массив значений длиной . Заметим, что введение (9) не налагает никаких дополнительных ограничений на решение – просто в общем случае в последовательности (9) все числа оказываются различными.Введем теперь достаточно разумное предположение о том, что фактическое (эффективное значение  с некоторой точностью прямо пропорционально геометрической жесткости поперечного сечения кузова, т.е. его моменту инерции относительно горизонтальной центральной оси:,                                                                                                                (10)где (баз) означает любой ненулевой элемент, принятый в качестве базового. Примем это предположение и перейдем к решению задачи идентификации. Прогиб в любой точке  для схемы по рис.1 может быть вычислен по формуле:,     .                                                                                                   (11)( - число участков перемен жесткости). из (11) в соответствии с теоремой Кастильяно может быть найден в виде интеграла:                                                                                                                    (12) 39 39где  - изгибающий момент от внешней нагрузки, включая моменты от реакции опор , которые находятся так:                                                                            (13)Выражение для  на участке  имеют вид:,                                                                                     (14)где;;                                                                                                   (15).Момент  от единичной силы , приложенной в точке , строится в виде суммы трех линейных слагаемых:,                                                                                 (16)где  вычисляется по формуле:                                                                                       (17)При этом реакции  из (6) находятся так:,     .                                                                                     (18)На основании (14), (17) интеграл по (12) вычисляется следующим образом:,                                                                                        (19)где .;     ;     ;     ;;     .                                                                  (20)Если у нас имеется (s) различных групп жесткости, то формула (11) с учетом (9) приобретает вид,     .                                                                               (21)Сравнивая перемещения ,  из (21) с перемещениями  найденными экспериментально (или иным образом, не совпадающим с (21)) с помощью функционала вида (4), в который дополнительно введены весовые множители :,     ,                                                                                   (22)приходим к СЛАУ, из которой вектор ,  определяется из матричного соотношения:,     ,                                                                           (23)в котором вектор  состоит из  компонентов , вычисляемых по формуле (21).Решение -мерной задачи в виде (23) некорректно по Адамару и может быть регуляризовано так, как это было описано выше. Однако, если воспользоваться разумным предположением (10) и учесть, что формула (23) при  всегда дает качественно (но не количественно) физически обоснованные результаты, то мы придем к одномерной задаче, имеющей решение, выражаемое формулой:,     ,                                                                                         (24)после чего идентифицированные значения изгибных жесткостей участков находятся по формуле, следующей из (10):.                                                                                                                         (25)Решение (24) – (25) имеет ту же точность, что и гипотеза (10). Это означает, что по (24) – (25) восстанавливаются значения  точно, если принять  (проверено на решении модельных задач). 40Предложенная постановка снимает вопрос о возможной множественности решений, т.к. решение (24) не содержит параметра регуляризации. Однако, можно также сказать, что соотношения вида (10) являются естественными регуляризаторами и тогда решение (24) является однозначным с точностью до вектора регуляризации . Предложенное простое решение в некотором смысле близко к нормативной методике 70-х годов. [28, с.135-137], где определение низшей ЧИК сведено к рассмотрению системы с тремя  для уравнения (6). При этом использовалось решение (24) для , которое в этом случае при любом  дает . Некоторая нелогичность методики [28] состоит в том, что теоретический прогиб  в середине кузова определяется от равномерно распределенной нагрузки  под брутто по схеме плоской безраскосной рамы (схема Никольского [29]), но затем интерпретируется как теоретический прогиб середины однородной балки Бернулли-Эйлера по рис.1:,где  - база кузова, ,  - длина консоли.Кроме того, методика [28] предполагает, что , делая тем самым ненужным экспериментальный этап. Имевшийся тогда опыт определения низшей ЧИК, позволяет сделать вывод, что указанная выше нелогичность компенсировалась выбором числового значения . Наверное в настоящее время подход, основанный на использование фиксированных значений даже вектора  не является вполне оправданным из-за быстро меняющейся номенклатуры моделей.После нахождения значений изгибной жесткости несколько ЧИК может быть найдено стандартными методами. При этом следует иметь в виду, что идентификация проводилась для схемы балки Бернулли-Эйлера с кусочно-постоянной жесткостью, имеющей, однако, прямолинейную ось, на которой расположены центры тяжести всех поперечных сечений.Полученные значения прогибов могут быть использованы для предварительных оценок низшей ЧИК без проведения громоздких расчетов. Для этого следует перейти от дискретно-континуального распределения масс (аналогично рис. 2) к дискретному с сосредоточенными массами  ( - количество узлов с заданными прогибами) и принять, что форма колебаний подобна экспериментальной кривой статического деформирования. После этого для схемы опирания по рис. 1, приравнивая работу сил тяжести  кинетической энергии, получим для низшей ЧИК выражения:,,     ,     ,                                                                 (26)где  - абсциссы и экспериментальные ординаты статического прогиба.Если в точках  расположены линейно-упругие опоры с жесткостями , то, представив жесткие (без деформаций) перемещения в виде:,                                                                                                                            (27)получим для частот подпрыгивания галопирования:,                                                                                                             (28)где ,     ,     ,,     ,     . Жесткие статические перемещения для схемы на упругих опорах также имеют вид (27), причем:,     ,где  - ускорение свободного падения.Если теперь представить форму колебаний схемы на упругих опорах в виде:,                                                                                                            (29)где  - нормированные формы колебаний подпрыгивания и галопирования:,а  - коэффициенты распределения амплитуд для жестких форм, удовлетворяющие требованию,( находится по (27)), а также условию ортогональности скелетных схем:,                                                                                             (30) 41откуда получается (например, с помощью МНК):.Если после этого считать, что  в (29) представляет собой экспериментальную статическую форму , то приходим к частотному уравнению вида (6) для матриц 3-го порядка, причем:,а элементы  матрицы  имеют вид:(при этом вследствие (30) будет ).Описанные вычисления соответствуют процедуре Ритца и приводят к завышенным значениям частот. Значительно более точные результаты можно получить, используя метод Граммеля [30]. В этом случае вместо формулы (26) для схемы на жестких опорах получим:,     ,                                                                                            (31)где ,причем значения  вычисляются по формулам (16)-(18) для .В случае упругих опор метод Граммеля вместо (6) приводит к уравнению с матрицами 3-го порядка:,где:     ;     M4=1DQ*ΣQ ;    (32)С =          =        (неуказанные элементы в таблице для матрицы  равны нулю).При этом  вычисляется по (31),  - по (26),  - по (29), в   матрице  первый столбец состоит из единиц, второй – из координат сосредоточенных масс, а третий представляет собой вектор экспериментальных значений перемещений. Матрица  состоит из элементов , которые представляют собой интегралы: ,     ,,причем  вычисляются по формулам (16)-(18) для , . Выводы 1. На основе рассмотрения существующих методов модального анализа различных машиностроительных конструкций сделан вывод о возможности косвенной оценки низшей частоты изгибных колебаний кузова пассажирского вагона в вертикальной плоскости.2. С учетом специфики конструкции предложено для нахождения низшей частоты производить расчетно-экспериментальную оценку изгибной жесткости  кузова, представляя его в виде одномерной балки Бернулли-Эйлера с кусочно-постоянными значениями . Правомерность такого подхода основана на идентичности матричных уравнений, описывающих динамическое поведение различных конструкций. 423. Для упрощения решения предложено использовать тот факт, что для пассажирского вагона количество перемен значений  по длине кузова значительно больше количества различных числовых значений жесткости.4. Использование метода наименьших квадратов при сравнении поперечных перемещений модели (расчетных значений) и прототипа (экспериментальные значения на натурном кузове) сводит задачу идентификации значений  к многомерной обратной коэффициентной задаче.5. Предложена гипотеза о пропорциональности фактического значения  прототипа легко вычисляемому геометрическому значению  модели, что позволило свести многомерную обратную коэффициентную задачу (неустойчивую по Адамару) к элементарной задаче поиска экстремума функции одной переменной.6. Предложенный алгоритм использовался при решении модельных задач, в которых вычислялись перемещения от поперечных нагрузок для схемы с заданными значениями , а затем эти перемещения использовались в качестве исходных данных для алгоритма идентификации. Вычисления показали высокую точность предложенного подхода.7. Предложены простые соотношения для оценки низшей частоты изгибных колебаний, основанные на методиках Ритца и Граммеля, использующие полученные значения кусочно-постоянной жесткости. При этом формула для частоты схемы на жестких опорах (формула (26) вообще не требует решения задачи идентификации – в нее напрямую входят экспериментальные значения прогибов.