УДК 624.042.8:629.73
А. С. РАСПОПОВ, В. Е. АРТЕМОВ, С. П. РУСУ (ДИИТ)
МОДЕЛИРОВАНИЕ КОЛЕБАНИЙ БАЛОЧНЫХ ЖЕЛЕЗНОДОРОЖНЫХ МОСТОВ В СРЕДЕ ОБЪЕКТНО-ОРИЕНТИРОВАННОГО ПРОГРАММИРОВАНИЯ DELPHI
Описано основш етапи створення програмного комплексу для розрахунку статики та динамши стержне-вих та балочних конструкцш залiзничних моспв на основi об'еднання методiв скшченних елементiв i рiв-нянь динамши твердого тiла.
Описаны основные этапы создания программного комплекса для расчета статики и динамики стержневых и балочных конструкций железнодорожных мостов на основе объединения методов конечных элементов и уравнений динамики твердого тела.
The basic stages of development of the program complex for calculation of statics and dynamics of the rod and beam constructions of railway bridges based on combining the finite-element method and the equations of solid dynamics are described.
Введение
Наиболее часто для расчета мостовых конструкций применяются различные компьютерные программы и расчетные комплексы, которые реализуют метод конечных элементов (Lira, Scad, Robot, Ansys и др.). Однако применение компьютерного моделирования для динамического расчета стержневых систем вызывает определенные трудности [1]. В технической и справочной литературе по расчету стержневых систем мало внимания уделено получению решений в завершенной матричной форме, пригодной для прикладных вычислений. В данной работе рассмотрены вопросы,
связанные с учетом произвольной ориентации стержней в глобальной системе координат, диссипативных свойств конструкции, неопределенностей, возникающих при формировании общей матрицы жесткости системы. Все приведенные в работе алгоритмы реализованы в программном комплексе Belinda [2].
1. Особенности применения МКЭ для динамического расчета
Алгоритм расчета стержневой системы с использованием метода конечных элементов может быть представлен следующей блок-схемой (рис. 1):
Processor
U
Static
U
SetStaticFixations
U
CalculateStage 1
U
CalculateStage2
U
CalculateStage3
U
Postprocessor
Dynamic (time = t0,..., ti,..., t1)
U
I
Time step t = t0
U
SetlnitialConditions
U
ResetSolver
U
SetDynamicFixations
U
CalculateStagel
Time step t = ti
U
SetDrift
U
CalculateStage2
U
CalculateStage3
Time step t = t1
U
Postprocessor
I
Рис. 1. Алгоритм статического и динамического расчета стержневой системы
© Распопов А. С., Артемов В. Е., Русу С. П., 2010
217
На данной блок-схеме (см. рис. 1) использованы оригинальные обозначения процедур из исходного кода модулей FEM.pas, UnitDynam-ics.pas [2]. Для проведения статического расчета (Static) достаточно подготовить информацию о кинематических закреплениях узлов системы (процедура SetStaticFixations) и последовательно выполнить три этапа расчета: CalculateStagel, 2, 3. В первом блоке CalculateStagel вычисляются матрицы поворота каждого стержневого элемента, формируются вектор кинематических закреплений, матрицы жесткости С и податливости L [3]. На втором этапе CalculateStage2 определяются вектор вынужденных перемещений узлов системы и вектор эквивалентных узловых усилий. На третьем этапе CalculateStage3 осуществляется сбор всех приложенных к системе нагрузок, проводится расчет системы методом конечных элементов в форме метода перемещений и строятся матрицы результирующих узловых усилий и усилий в стержнях. По завершению расчета становятся доступными средства просмотра результатов (процедура Postprocessor).
Как видно из рис. 1, для динамических расчетов (Dynamic) нет необходимости проводить первый этап вычислений CalculateStagel на каждом шаге интегрирования. После подготовки начальных условий (SetlnitialConditions) и подготовки решателя системы дифференциальных уравнений (ResetSolver) действия на первом шаге динамического расчета аналогичны действиям в статическом расчете. Этапы расчетов CalculateStage2, 3 выполняются на каждом шаге интегрирования.
В статическом расчете условия, по которым на перемещения узлов накладываются кинематические закрепления, известны [3]. В динамике к этим условиям добавляется еще одно: узел должен быть кинематически закреплен, если в нем сосредоточена масса или момент инерции массы. Последнее обстоятельство определяет участие кинематических параметров узла в общей системе уравнений движения, а также различие процедур SetStaticFixations и SetDynamicFixations.
2. Произвольная ориентация стержня
Локальная матрица жесткости С}- пространственного стержневого элемента [3] приводится к виду (1).
В расчете балочного пролетного строения со сквозными фермами матрица жесткости может быть использована в форме (1) только для продольных балок, элементов нижнего и верхнего пояса, для которых локальные системы коор-
динат совпадают по направлению с глобальной системой координат.
С =
EA
~т
о о о о о
о
12EJZ
о о о
6EJz
о о
12EJ у I3
о
- 6EJ i2
о
о о
0
GJx
1
о
о
_ 6J 12
0
4EJ y
1
о
о
6EJz
о о
0
4EJ z
1
(1)
Для остальных элементов (поперечных балок, раскосов, связей) матрицу жесткости, а также вектор суммарных узловых нагрузок ^ необходимо преобразовать с помощью соответствующей матрицы поворота:
Со = г • C • T; J Fn = T • F • T
(2)
где Т - унитарная изометричная матрица поворота [4, 5], которую можно представить в блочной форме
T =
ft] о о о
о ft] о о
о о ft] о
о о о ft]
(3)
Элементы матрицы Тх описывают ориентацию локальной системы координат стержня в глобальной системе координат и представляют собой косинусы углов между их осями. Матрица Тх может быть получена различными способами: последовательным вращением системы координат вокруг осей х, у, г на соответствующие углы Эйлера фх, фу, фг; вращением
вокруг заданного вектора Уа на угол а ; преобразованиями кватернионов или заданных ранее матриц поворота; прямым составлением таблицы направляющих косинусов [6]. Одним из самых простых методов построения матрицы поворота является вращение системы координат вокруг заданного вектора. Однако для такого построения необходима информация о его координатах. В общем случае, положение произвольно ориентированного стержня в глобальной системе координат известно (это разность координат узлов, соединяемых стержнем). Остается найти координаты вектора поворота стержня. Эта особенность мало освещена в на-
учнои и справочной литературе, где приводятся примеры в основном с плоскими расчетными схемами, а стержни располагаются либо по направлению осей глобальной системы координат, либо в ее плоскостях. Общий пространственный случай требует проведения некоторых дополнительных математических преобразований.
Рассмотрим стержневой элемент, соединяющий узлы конструкции А, В и занимающий произвольное положение в глобальной системе координат (рис. 2). Согласно принятой схеме, продольная ось стержня х' направлена вдоль
вектора АВ . Условимся считать, что плоскость изгиба стержня, образованная его локальной осью х ' и осью г глобальной системы координат, является вертикальной. Тогда, используя векторное произведение, можно получить вектор, перпендикулярный данной вертикальной плоскости:
У = X
х z; z ={0 0 1}.
(4)
Z = x ' х y' .
(5)
y' в (5) необходимо подставить единичный вектор у :
У = у; у = {0 1 0}.
(6)
Далее, записывая последовательно 9 косинусов углов между векторами x', y', У (осями стержня) и векторами x, y, z (осями глобальной системы координат), получим искомую матрицу поворота
cos (x'А x) cos (У а x)
cos (z А x ) cos (У А y ) cos (z А z )
Tx =
cos(x'Aу) cos (y A у)
cos(x'A z ) cos (у A Z )
' A:
. (7)
Рис. 2. Ориентация стержня в глобальной системе координат
Здесь представляет собой проекцию вектора х' на плоскость ху глобальной системы
координат, а сам вектор х' перед использованием в формуле (4) должен быть нормализован.
Полученный вектор У определяет положительное направление одной из главных центральных осей сечения стержня. Повторным векторным произведением находим недостающую ось триэдра локальной системы координат стержня (рис. 2):
Аналогичный подход к представлению поворота описан в справочной системе расчетного комплекса Autodesk Robot Structural Analysis Professional, с решениями которого сравнивались результаты работы комплекса Belinda. Результаты тестовых статических расчетов по обоим комплексам показали точные совпадения при различных видах деформации концов стержня.
3. Конструкционное демпфирование
Рассмотренная модель пролетного строения моста относится к классу консервативных механических систем. В реальной конструкции часть энергии колебаний будет рассеиваться во внешнюю среду за счет сил вязкого трения в материале конструкции и сил сухого трения в ее опорных частях [7, 8]. Для учета диссипа-тивных свойств материала конструкции введем в расчет силу вязкого сопротивления, пропорциональную скорости перемещения узла:
Fd =•
-ßxX -Руу -ßzZ
Фх Ф у -ßTz
(8)
Формула (4) справедлива во всех случаях, кроме одного, когда продольная ось стержня х параллельна оси г глобальной системы координат (карданов подвес). В таком случае вместо
где Р - коэффициент вязкого сопротивления в направлении соответствующей степени свободы.
Сила вязкого сопротивления, характеризуемая вектором (8), линейно зависит от скорости перемещения узла и прикладывается на каждом шаге динамического расчета системы. При этом используется значение скорости /-го узла V с предыдущего шага расчета. Коэффициент в (эквивалентный коэффициент сопротивления, [9]) определяется по формуле
ß =
ус 2nv
(9)
где у = 25 - коэффициент поглощения энергии; 5 - логарифмический декремент колебаний; c - жесткость узла конструкции; v - частота колебаний конструкции, Гц.
Запись формулы (9) предполагает, что параметры вязкого сопротивления определяются на основе характеристик динамического процесса, т. е. «опережают» решение системы уравнений движения узлов. Для расчета системы с неизвестными характеристиками колебаний необходимо сначала провести расчет без учета демпфирования, установить частоту колебаний, вычислить параметры демпфирования и повторить расчет с учетом вычисленных характеристик. Демпфирующие свойства некоторых механических систем представлены в справочных материалах, например [8], для металлических балочных пролетных строений железнодорожных мостов разрезных систем, логарифмический декремент колебаний 5 = о,о8. Это значение использовано далее при динамическом расчете балочного пролетного строения со сквозными фермами расчетным пролетом 11о,о м.
В программном комплексе Belinda модель сухого трения в линейной постановке реализована в виде сосредоточенных динамических сил Fs, зависящих от перемещений узлов конструкции и коэффициента трения ц [7]:
F ={_^ ^ _ц'г}, (1С)
[_Цф*Ф* _Цфу Фу _Цфгфг J
где x, у, z - линейные перемещения узла в направлении соответствующих степеней свободы; Ф*, фу, фz - то же, угловые перемещения.
Силы сухого трения (1о) прикладываются к узлам системы и вычисляются аналогично силам вязкого сопротивления (8). На мостах силы сухого трения возникают в опорных частях при смещении пролетного строения. Величина коэффициента трения может быть найдена по справочным данным.
4. Кинематические закрепления узлов конструкции
При моделировании следует учитывать некоторые особенности определения жесткост-ных характеристик дискретной стержневой системы. Так, обращение общей матрицы жесткости C и получение матрицы податливости L возможно только при использовании редуцированных форм этих матриц. Наличие нулевых строк и столбцов в матрице жесткости влечет за собой появление в матрице податливости эле-
ментов с бесконечно большими значениями. В вычислительном программировании такие значения называются NaN-величины (Not-a-Number). К этому может привести не только отсутствие опорных закреплений в рассматриваемой модели, но и некоторые другие факторы. Например, при компьютерном моделировании жесткостные характеристики стержней хранятся в вещественных переменных, которые инициализируются нулевыми значениями. Современные расчетные комплексы в таком случае либо сигнализируют пользователю о необходимости переназначить жесткость, либо автоматически присваивают параметру новое ненулевое значение. Такая процедура реализована, например, в пакете конечноэлементного анализа Autodesk Robot Structural Analysis Professional. Последний вариант, однако, нельзя считать универсальным, так как значение 10-5 м4 для момента инерции стержня на кручение и 1,0 м4 для момента инерции сечения на изгиб - конечные величины. Предлагаемое альтернативное решение состоит в назначении фиктивного кинематического закрепления в направлении степени свободы с нулевой жесткостью. Данный подход позволяет не только автоматически учитывать обнуленные жесткости стержней, но и вводить временные кинематические закрепления узлов при динамическом расчете. После установки такого закрепления перемещение соответствующего узла рассматривается как заданное.
5. Контактное взаимодействие с подвижной нагрузкой
Базовая теоретическая модель контакта несущих элементов пролетного строения моста и подвижной нагрузки изложена в [10]. В текущей версии программного комплекса Belinda реализована усовершенствованная модель этого контакта, позволяющая учитывать воздействие на конструкцию группы постоянных, гармонических сил и моментов, движущихся поступательно или вращательно с заданной скоростью.
Программный объект «контакт» представляет собой элемент с двумя параметрами - номером контактной группы и номером элемента конструкции, как правило, стержня. Контактная группа может содержать произвольное число контактов. Каждая динамическая нагрузка относится к той или иной контактной группе, ссылаясь на нее по номеру k = 1, 2,..., n . Если нагрузка приблизилась к элементу конструкции из контактной группы на расстояние менее чем
заданный порог (по умолчанию 0,01 мм), то она прикладывается к элементу. При этом все характеристики нагрузки (величина, частота, скорость), в общем случае, зависят от времени t. Такие предположения позволяют существенно увеличить количество динамических нагрузок и контактов, используемых в модели, без заметного снижения быстродействия программного комплекса.
6. Пример динамического расчета пролетного строения моста
В качестве примера рассмотрим динамический расчет пространственной модели балочного пролетного строения железнодорожного моста со сквозными фермами с применением программного комплекса Belinda.
Пролетное строение имеет следующие геометрические размеры: пролет - 110,0 м; высота
главных ферм - 15,0 м; панель главных ферм -11,0 м; панель продольных связей - 5,5 м; расстояние между фермами - 5,8 м. Проезжая часть включена в совместную работу с нижними поясами главных ферм. Общий вес пролетного строения, включая мостовое полотно и другие детали, составляет 519,7 т. Схема опи-рания конструкции - шарнирная, опорные части на одном из концов позволяют системе свободно перемещаться в продольном направлении.
Дискретная модель пролетного строения показана на рис. 3. С каждым стержнем связывается локальная (местная) система координат. Узлы конструкции определены в глобальной (общей) системе координат О, начало которой совмещено с центром тяжести первой поперечной балки. Модель включает 62 узла, 136 стержней, 744 кинематических параметра.
Рис. 3. Модель пролетного строения в программном комплексе
Нагрузка моделируется в виде четырех подвижных гармонических сил, соответствующих центрам крепления тележек электровоза ВЛ8. Статическая составляющая каждой силы моделирует долю собственного веса локомотива, гармоническая - динамические свойства локомотива, полученные на основе натурных испытаний при движении локомотива по прямолинейному участку пути с постоянной скоростью. Скорость движения нагрузки - 72 км/ч.
Общая система разрешающих уравнений движения узлов состоит из 726-ти дифференциальных уравнений (динамика) и 62-х алгебраических (статика). Динамический процесс наблюдается в течение 10 с, шаг интегрирования составляет 10-4 с . Формирование общей неоп-
тимизированной матрицы жесткости заняло 15 мин., полное время динамического расчета на компьютере с 2 GB оперативной памяти и двуядерным процессором Intel Core 2 Duo 2,2 GHz - 6 часов. Были проведены расчеты пространственной динамики пролетного строения при движении локомотива ВЛ8 со скоростями 5 м/с, 15 м/с и 30 м/с. Перемещения середины продольной балки для различных скоростей движения в зависимости от положения начала локомотива показаны на рис. 4.
Из рис. 4 очевидно, что амплитуда вертикальных колебаний центрального узла продольной балки возрастает с увеличением скорости движения локомотива.
UZ(38) (м) 0,00-1
О 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150
Рис. 4. Перемещения центрального узла продольной балки
Для более простых систем, например, балочных металлических пролетных строений со сплошной стенкой, ребристых железобетонных конструкций, а также плоских моделей порядок матриц жесткости и податливости становится существенно меньше. Например [11], динамический расчет плоской дискретной модели железнодорожного пролетного строения длиной 33,6 м, включая формирование основных матриц МКЭ, составил около 15 мин.
В следующих версиях программного комплекса Belinda планируется ввести распределенные и импульсные динамические нагрузки, а также рассмотреть вопросы взаимодействия подвижных нагрузок с мостовыми конструкциями при наличии эксцентриситетов и других неровностей железнодорожного пути.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Распопов, А. С. К вопросу компьютерного моделирования движения поезда по мосту [Текст] / А. С. Распопов, С. П. Русу, В. Е. Артемов // Методи розв'язування приклад. задач мех. де-формiвного твер. тша: зб. наук. пр. / Дшпро-петр. нац. ун-т. - 2007. - Вип. 8. - С. 133-139.
2. Распопов, А. С. Разработка программного комплекса «Belinda» для расчета нелинейных колебаний железнодорожных мостов [Текст] / А. С. Распопов, С. П. Русу, В. Е. Артемов // Дороги i мости: зб. наук. пр.: в 2 т. - Т. 2. - К.: ДерждорНД1, 2007. - Вип. 7. - С. 136-143.
3. Розин, Л. А. Стержневые системы как системы конечных элементов [Текст] / Л. А. Розин. - Л.: Изд-во Ленингр. ун-та., 1975. - 237с.
4. Хорн, Р. Матричный анализ [Текст] / Р. Хорн, Ч. Джонсон : [пер. с англ.] - М.: Мир, 1989. -655 с.
5. Справочник по строительной механике корабля [Текст]: в 3 т. - Т. 2. Пластины. Теория упругости, пластичности и ползучести. Численные методы / Г. В. Бойцов и др. - Л.: Судостроение, 1982. - 464 с.
6. Корн, Г. Справочник по математике для научных работников и инженеров [Текст] / Г. Корн, Т. Корн. - М. : Наука, Гл. ред. физ.-мат. лит., 1968. - 720 с.
7. Динамика железнодорожных мостов [Текст] / Н. Г. Бондарь и др. - М.: Транспорт, 1965. -412 с.
8. Fryba, L. Dynamics of Railway Bridges [Text] / L. Fryba. - Praha: Academia Praha, 1996. - 330 p.
9. Вибрации в технике [Текст]: Справочник. В 6 т. - Т. 6. Защита от вибрации и ударов / ред. совет: В. Н. Челомей; под ред. К. В. Фролова. -М.: Машиностроение, 1981. - 456 с.
10. Распопов, А. С. Моделирование подвижных нагрузок при расчетах динамики дискретных систем «мост-поезд» в программном комплексе «Belinda» [Текст] / А. С. Распопов, С. П. Русу, В. Е. Артемов // Актуальш проблеми мехашки суцшьного середовища i мщносп конструкцш: мгжнар. наук.-техн. конф. пам'яп акад. НАН Украши В. I. Моссаковського (Дшпропет-ровськ, 17-19 жовт. 2007 р.): тези доп. - Д., 2007. - С. 282-283.
11. Распопов О. С. Динамша балкових конструкцш моспв пвд рухомим навантаженням [Текст] / О. С. Распопов, В. £. Артьомов, С. П. Русу // Мехашка i фiзика руйнування будiвельних ма-терiалiв та конструкцш: Зб. наук. ст. / за заг. ред. Лучка Й. Й. - Вип. 8. - Львiв: Каменяр, 2009. - С. 712-721.
Поступила в редколлегию 12.04.2010.
Принята к печати 23.04.2010.