Численное моделирование динамической задачи обобщенной системы с распределенными параметрами

Рассмотрим динамическую модель сооружения, которая представляет собой систему с распределенной массой m(x) и сосредоточенными массами m_j, установленными в некоторых точках (рис. 1). Данное сооружение башенного типа, в котором основными характеристиками являются изгибная жесткость EI(x) и погонная масса m(x), представляет систему с бесконечным числом степеней свободы. Если предположить, что перемещения сооружения ограничены только одной формой колебаний, то такая модель с обобщенной координатой w(t) может быть использована для расчета на сейсмические воздействия в виде заданной акселерограммы землетрясения. Таким образом, предполагается, что движение сооружения, как обобщенная система с одной степенью свободы, происходит только по одной форме в виде безразмерной функции φ(x) с амплитудой колебаний w(t). Следовательно, при такой форме движения перемещение и соответствующие производные по времени выражаются так:

(1)

При заданном перемещении основания в виде функции w₀(t), общие перемещение, скорость и квадрат скорости записываются таким образом:

(2)

Для определения потенциальной энергии исследуемого объекта, который испытывает изгиб и продольное сжатие, необходимо иметь вторую производную от функции W(x,t) по x:

Уравнение движения. Для вывода уравнений движения обобщенной системы воспользуемся уравнением Лагранжа второго рода, соответствующего вариационному принципу Гамильтона [2–4]:

(3)

где T – кинетическая энергия;

U – потенциальная энергия деформации изгиба;

Q_p – обобщенная сила (сила трения или внешняя возмущающая сила, зависящая от времени);

Q_d – обобщенная сила сопротивления.

Кинетическая и потенциальная энергии модели, представленной на рис. 1, с учетом (1), (2) и (2а) записываются в виде:

(4)

где N(x) – нормальная сила от собственного веса сооружения, которая может изменяться по линейному закону.

Суммарная потенциальная энергия от изгибающего момента и сжимающей силы будет равняться

(5)

Здесь предполагается, что опорная часть сооружения во время землетрясения приобретает только горизонтальное перемещение.

Производные от кинетической и потенциальной энергии с учетом (4) и (5) приобретают вид

(6)

(7)

(8)

Обобщенная сила сопротивления определяется по формуле

где Ф – диссипативная функция Рэлея [5][6], которая для системы с распределенными параметрами затухания представляется в виде

(9)

Внося (6)–(9) в (3), получаем уравнение

которое представляется в виде

где m*, k*, c* – обобщенные масса, жесткость и затухания;

m – коэффициент сейсмического возмущения, который характеризует уровень возможных колебаний сооружения по форме φ(x) [7].

Поделив обе части (10) на обобщенную массу m*, получим

(11)

Дифференциальное уравнение (11) описывает динамическое состояние обобщенной системы с одной степенью свободы от сейсмического воздействия в виде заданной акселерограммы. Здесь предполагается, что функции φ(x), c(x), m(x), EI(x), N(x) будут заданы в зависимости от геометрических и физико-механических свойств исследуемого объекта. Коэффициент μ в (10) характеризует различие между одномассовыми системами с сосредоточенной и обобщенной массами.

Обобщенные характеристики. Рассмотрим сооружение башенного типа переменного сечения по высоте. Предполагается, что сооружение круглого сечения в плане состоит из двух центральных колец постоянного сечения, идущих от основания до верхней части объекта и наружных колец, число которых в зависимости от высоты уменьшается. При этом безразмерная функция формы φ(x), распределенная масса m(x), момент инерции I(x) и нормальная сила N(x) в исследуемом объекте изменяются по следующим законам:

(12)

(13)

(14)

где m₁, I₁ – масса и момент инерции опорной части без учета центральной части;

m₀, I₀ – масса и момент инерции центральной части сооружения;

a – длина нижней части сооружения с переменным сечением.

Тогда из (10а) с учетом (12), (13) получим интегральное выражение обобщенной массы:

(15)

Из (10б) с учетом (14) и производных от функции формы

получим интегральное выражение обобщенной жесткости изгибающего момента

(16)

Если предположить, что нормальная сжимающая сила

N(x) = N = γ A_срl, A_ср(A₁ + A₀)/2,

где γ, A₁, A₀ – объемный вес, площади сечений элементов на внешних кольцах и в центральной части объекта, получаем

(17)

Обобщенная масса в правой части (10) выражается следующим образом

(18)

Если принять, что a = l/3, то непосредственным интегрированием [8] из (15)– (18) получим

Полученные обобщенные параметры системы (19), при заданных значениях m₀, m₁, I₀, I₁, A₀, A₁, позволяют перейти к решению дифференциального уравнения (11). Из решения уравнения (11) определяется w(t), а затем с использованием (1) вычисляются перемещения, скорость и ускорения в произвольных точках по оси x, соответствующие моменту времени t. Далее вычисляются внутренние усилия в сечениях исследуемого объекта.

При заданной правой части (11) в виде реальной или синтезированной акселерограммы решение осуществляется либо интегралом Дюамеля с последующим численным интегрированием, либо непосредственным численным дифференцированием. Ниже рассматривается реализация обоих методов на примере системы с распределенными параметрами.

Численное дифференцирование. С целью численного решения дифференциального уравнения (11) применим метод последовательных аппроксимаций [9][10], где скорость и ускорения, соответствующие моменту времени t_n, представляются в виде

(20)

(21)

где α_i, β_i – числовые коэффициенты [11];

τ = t_n – t_n–1 – шаг по времени.

Записав уравнение (11) в момент времени t_n и внося (20) и (21) в (11), после некоторых преобразований получим рекуррентную формулу для определения w_n

(22)

Затем с использованием (1) определяются перемещения, скорость и ускорения в произвольной точке координатной оси x, по длине сооружения, соответствующие моменту времени t_n дискретной оси t

(23)

Затем вычисляются значения сил инерции и сил упругости:

(24)

(25)

Поперечная сила и изгибающий момент в точке, расположенной на высоте h от основания сооружения, выражаются следующим образом:

(26)

(27)

Если функции m(x) и φ(x) заданы, то интегралы (26) и (27) вычисляются непосредственным интегрированием.

Численное моделирование с применением интеграла Дюамеля. Решение уравнения (11) можно также представить интегралом Дюамеля [7]:

(28)

Этот интеграл с учетом

sinω(t – τ) sinωt cosωt – cosωτ sinωτ

представляется в виде

w(t) = – λ[a(ω)sinωt – b(ω)cosωt],

(29)

Интегралы a(ω) и b(ω) в (28) вычисляются методом Симпсона [12]. При этом шаг интегрирования выбирается в зависимости от шага оцифровки заданной акселерограммы ẅ₀(t). Затем вычисляются перемещения, скорость и ускорения в произвольной точке сооружения:

(30)

Далее определяются распределенная сила упругости, сдвигающая сила и изгибающий момент в основании сооружения:

На основе изложенных выше алгоритмов численного моделирования были разработаны компьютерные программы и получены результаты от заданного воздействия в виде акселерограммы землетрясения.

Пример. В качестве примера рассматривается сооружение башенного типа (рис. 2), которое состоит из ряда колонн, расположенных во внешних кольцах, а также из двух центральных колец сплошного сечения, соединенных между собой диафрагмами жесткости (рис. 2б). При этом количество внешних колец постепенно уменьшается и на высоте x ≥ h₁ остаются только центральные кольца.

Суммарные площади сечений, распределенные массы и моменты инерции несущих элементов центральных (A₀, m₀, I₀) и пяти внешних колец (A₁, m₁, I₁) равняются:

Следовательно, частота и период свободных колебаний сооружения

ω = = 9,18 рад/с,

T = 2π/ω = 0,68 с.(г)

На основе разработанной программы на языке FORTRAN получены результаты как методом численного дифференцирования, так и методом численного интегрирования. Вначале, на основе (20)–(23), были исследованы свободные колебания сооружения. На рис. 3 показаны графики изменения перемещения w_A и w_L, полученные от действия мгновенного импульса при τ = = 0,005 с и ξ = 0,05.

Из этих графиков следует, что период колебаний, полученный численным решением, совпадает с результатом (г).

На втором этапе на основе метода последовательных аппроксимаций исследуются колебания сооружения от действия землетрясения, представленного в виде заданной акселерограммы El Centro. На рис. 4 сравнивается ускорение, полученное для верхней точки сооружения, с графиком ускорения самого землетрясения. Видно, что относительное ускорение сооружения более чем в два раза превышает ускорение его основания.

Получены также результаты с применением алгоритма (28)–(30) и численного решения интеграла Дюамеля при τ = = 0,01 с
и ξ = 0,05. На рис. 5 для сопоставления приводятся аналогичные рис. 4 графики, полученные численным интегрированием. Можно заметить, что пиковые ускорения второго метода (рис. 5) несколько меньше, чем в первом методе (рис. 4), хотя характер изменения по времени практически совпадает. По-видимому, это объясняется тем, что шаг по времени в первом методе в два раза меньше, чем во втором.

На рис. 6 приведены графики перемещения верхней точки сооружения, полученные по первому (рис. 6а) и второму методам (рис. 6б). Сравнение показывает их практическое совпадение.

Вывод

Разработаны алгоритмы и компьютерные программы, которые позволяют исследовать динамические процессы, связанные с напряженно-деформированным состоянием сооружения, которое представляется в виде обобщенной системы, при различных воздействиях. Достоверность результатов подтверждается сравнением результатов, полученных двумя методами. Разработанные компьютерные программы могут быть использованы для проведения мониторинга уникальных объектов башенного типа.