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

Полный текст

Введение За последние десять лет в нашей стране было возведено достаточно много воздухоопорных сооружений различного назначения. В настоящее время самым распространенным типом воздухоопорных сооружений являются оболочки на прямоугольном плане (рис. 1). Об интересе к подобным конструкциям в отечественной науке свидетельствует появление обзорных статей по данной теме, например [1]. Сложность расчета таких сооружений на действие внешних нагрузок обусловлена рядом факторов: - необходимая жесткость конструкции создается за счет внутреннего давления, являющегося следящей нагрузкой; - форма оболочки может существенно изменяться в процессе нагружения; - материал оболочки является ортотропным и работает только на растяжение, - при ветровых воздействиях проявляются аэроупругие эффекты [2, с. 396]. Исследованиям колебаний оболочек посвящено множество как классических, например, [3-5], так и современных работ [6; 7]. Монография [3] является фундаментальным трудом, посвященным общим вопросам анализа собственных колебаний оболочек с помощью операторных методов. Для изотропных прямоугольных мембран влияние натяжения на частоты собственных колебаний описано в книге [4]. Некоторые вопросы колебаний гибких ортотропных оболочек канонической формы описаны в главе 3 книги [8]. Следует отметить, что аналитические решения удается получить только для оболочек и мембран простой (канонической) формы. В статье [9] изучаются собственные частоты и формы колебаний пневмобалки (круговой цилиндрической оболочки, предварительно напряженной за счет избыточного давления) с помощью метода динамической жесткости (dynamic stiffness method), разработанного в [10]. Пневмобалка моделируется в виде стержня типа Тимошенко, т.е. учитываются деформации поперечного сдвига. Такая модель пневмобалки позволяет учесть ортотропные физикомеханические свойства материала и нелинейные эффекты, вызванные следящей нагрузкой от внутреннего давления, и получить аналитическое решение [9]. Точная матрица динамической жесткости пневмобалки была получена путем линеаризации нелинейных уравнений около предварительно напряженного состояния. На основе проведенного обзора литературы впервые получены аналитические выражения для определения собственных частот колебаний пневмобалки из ортотропного материала при различных граничных условиях. Достоверность результатов подтверждается их сопоставлением с полученным ранее решением задачи для шарнирно опертой пневмобалки из изотропного материала [10], а также с результатами, полученными для пространственной оболочечной модели методом конечных элементов в ПК Abaqus. Рис. 1. Воздухоопорное покрытие футбольного стадиона в Краснодаре, спроектированного и возведенного ООО «ПСБ “Вертеко”» [Figure 1. Air-supported roof of football stadium at Krasnodar, Russia. Project and erection was performed by LLC “PSB Verteco”] 1. Математическая постановка задачи 1.1. Определение частот и форм собственных колебаний для системы с конечным числом степеней свободы В настоящее время для расчета строительных конструкций, как правило, используются численные методы, а поскольку метод конечных элементов сочетает в себе универсальность алгоритмов решения различных краевых задач с эффективностью компьютерной реализации вычислений [11], то он получил наибольшее распространение. Главный принцип, на котором он основан - это рассмотрение реальной конструкции в виде системы с конечным числом степеней свободы. Запишем в матричном виде уравнение движения для системы с конечным числом степеней свободы [12; 13]: Mu Cu Ku p&& &  , (1) где u - вектор узловых перемещений для всего тела; u& и u&& - векторы узловых скоростей и ускорений (точками обозначены производные по времени); M, , C K - «глобальные» матрицы масс, демпфирования и жесткости для всего тела; p - вектор эквивалентных узловых нагрузок для всего тела. Модальный анализ проводится для определения частот и форм (мод) собственных колебаний конструкций в предположении, что система является линейной. Внешние силы и демпфирование полагаются равными нулю. В этом случае уравнение колебаний конструкции в матричной форме (1) принимает вид Mu Ku&&  0. (2) Модальный анализ заключается в нахождении условий, при которых система совершает гармонические колебаний по закону [13]: u t sin t 0, (3) где  - вектор, характеризующий форму собственных колебаний (соотношения между смещениями узлов);  - круговая частота собственных колебаний; 0 - начальная фаза. Подстановка (3) в (2) дает уравнение собственных колебаний: K M  2 0. (4) Тривиальным решением уравнения является нулевое смещение узлов (система остается неподвижной). Нетривиальное решение уравнения (4) существует лишь тогда, когда величины  i i 1, ..., n обращают в нуль детерминант матрицы K M 2. Соответствующие им формы собственных колебаний  вычисляются лишь с точностью до произвольного множителя, значение которого определяется избранным способом нормировки собственных форм [13]. Таким образом, амплитуды собственный колебаний не определены, из решения уравнения (4) могут быть получены только соотношения между перемещениями различных точек системы. Отметим, что число собственных форм совпадает с числом степеней свободы динамической системы n. Наряду с круговой собственной частотой ωi также вводится собственная частота fi, представляющая собой число колебаний по i-той собственной форме, совершаемых системой за одну секунду: i . (5) fi  2 В качестве основного инструмента для проведения расчетов выбран верифицированный в системе РААСН универсальный программный комплекс (ПК) ANSYS Mechanical, в котором реализована возможность учета предварительного напряжения конструкции при проведении модального анализа, а также возможность использования мембранных конечных элементов из ортотропного материала [14]. 1.2. Особенности модального анализа предварительно напряженной системы Собственные частоты и формы колебаний предварительно напряженной системы отличаются от собственных частот и форм колебаний ненапряженной системы. Этот эффект можно наблюдать при натяжении струн в музыкальных инструментах - для систем, работающих преимущественно на растяжение, частоты растут с увеличением предварительного напряжения [4; 12; 13]. Как указано выше, модальный анализ системы имеет смысл только в линейной постановке, в нелинейных задачах можно говорить о разложении движения системы по формам свободных колебаний в окрестности изучаемого равновесного ее положения при линеаризации поведения системы в этой окрестности. Определение напряженно деформированного состояния (НДС) системы при предварительном напряжении должно выполняться с учетом геометрически нелинейных эффектов: расчет следует проводить по деформированной схеме, принимая во внимание следящий характер нагрузки от давления воздуха [9]. Далее необходимо построить динамическую модель так называемой линеаризованной системы, в которой все составляющие ее элементы рассматриваются в линейной постановке, но с касательными (мгновенными) матрицами жесткости [13]. Для этого в уравнении (4) вместо обычной матрицы жесткости K0 вводится касательная (мгновенная) матрица жесткости K : K  K0 KG, (6) где KG  K1 K2 - матрица геометрической жесткости: K1 - матрица начальных напряжений; K2 - матрица начальных поворотов. Детальное построение матрицы геометрической жесткости KG , а также матриц K1 и K2 подробно описано в книге [13, с. 138-140]. 2. Решение тестовой задачи по определению собственных частот и форм колебаний пневмобалки 2.1. Описание расчетных моделей Объектом исследования выбрана пневмобалка - цилиндрическая оболочка диаметром 0,28 м, длиной 3 м и толщиной 1 мм, исследуемая в работе [9]. Очевидно, что наличие торцов в оболочечной модели приводит к дополнительному продольному растяжению оболочки от действия избыточного давления, однако, в работе [9] не указано, было ли учтено это растяжение в построенной стержневой модели. В данной работе были выполнены расчеты для оболочечной модели пневмобалки как с торцами, так и без них. Для лучшей сходимости численного алгоритма решения форма торцов была принята в виде полусфер, расстояние между центрами которых составляло длину пневмобалки, т.е. 3 м. Трехмерная конечно-элементная (КЭ) модель пневмобалки создавалась в ПК ANSYS Mechanical с использованием трехузловых оболочечных конечных элементов типа SHELL181, в котором реализована теория оболочек Миндлина - Рейсснера. Для решения данной задачи была включена мембранная опция (параметр KEYOPT(1)=1), которая позволяет не учитывать изгибную жесткость и игнорировать угловые степени свободы в узлах элементов [14]. Для исследования сеточной сходимости рассмотрены три варианта расчетной сетки с максимальной длиной ребра КЭ d, равной 0,06, 0,03 и 0,015 м. Результаты расчета для двух последних случаев практически совпадают и приведены ниже. Было рассмотрено 2 варианта закреплений пневмобалки: 1) защемленный левый край и свободный от закреплений правый край (консольная балка); 2) защемленный левый край и шарнирно опертый правый край. 2.2. Параметры расчетов Рассмотрены пневмобалки, выполненные из двух различных ортотропных материалов - типа 1 и типа 2. Физико-механические характеристики приведены в табл. 1 и приняты в соответствии с данными из источника [9]. Таблица 1 Физико-механические характеристики материалов [Table 1. Physical and mechanical properties of materials] Для проведения модального анализа в ПК ANSYS Mechanical использовался блочный метод Ланцоша. Решение системы линейных алгебраических уравнений (СЛАУ) выполнялось разреженным методом (решатель Distributed Sparse Matrix Solver). 2.3. Результаты решения тестовой задачи Результаты решения тестовой задачи представлены в следующем виде: на рис. 2 и 3 пока- Характеристика [Property] Тип 1 [Type 1] Тип 2 [Type 2] Плотность ρ, кг/м3 [Density ρ, kg/m3] 1420 1420 Модуль упругости в продольном (осевом) направлении E1, МПа [Modulus of elasticity in longitudinal (axial) direction E1, MPa] 393,13 18370,0 Модуль упругости в кольцевом направлении E2, МПа [Modulus of elasticity in circular direction E2, MPa] 451,59 14120,0 Модуль сдвига G12, МПа [Shear modulus G12, MPa] 103,0 6460,0 Коэффициент Пуассона ν12 [Poisson ratio ν12] 0,07 0,28 Коэффициент Пуассона ν21 [Poisson ratio ν21] 0,08 0,22 заны формы колебаний для пневмобалок из материала типа 1 при различных граничных условиях, а в таблице 2 приведены собственные частоты для всех рассмотренных случаев в сопоставлении с результатами эталонного исследования [9]. Относительные погрешности ε вычислены по формуле: f0  f  100%, (7) f0 где f0 - эталонные результаты [9]; f - результаты расчетов, выполненных авторами статьи. Таблица 2 Результаты решения тестовой задачи [Table 2. Results of test case solution] Номер изгибной формы колебаний [Bending eigenmode number] Собственная частота fi, Гц [Natural frequency fi, Hz] Результаты авторов статьи [Authors’ results] Источник [Source] [9] εmin/εmax,% С торцами [With ends] Без торцов [Without ends] Abaqus Стержневая модель [beam model] d = 0,03 м, 7708 КЭ d = 0,015 м, 30560 КЭ d = 0,03 м, 7046 КЭ d = 0,015 м, 27888 КЭ d = 0,03 м, 2862 КЭ Случай 1.1. Консоль, материал типа 1 [Case 1.1. Cantilever, material type 1] 1 2,90 2,90 3,09 3,09 3,00 3,14 1,5/7,7 2 16,41 16,42 16,52 16,51 16,38 16,41 0,0/0,8 Случай 1.2. Консоль, материал типа 2 [Case 1.2. Cantilever, material type 2] 1 19,82 19,85 21,73 21,77 20,28 22,13 1,6/10,4 2 114,84 115,05 124,46 124,73 117,69 118,47 2,2/6,0 Случай 2.1. Жесткое закрепление с одной стороны и шарнирная опора с другой, материал типа 1 [Case 2.1. Clamped-simply supported ends, material type 1] 1 12,73 12,75 11,28 11,25 12,85 11,89 0,8/12,5 2 36,12 36,14 34,52 34,45 35,76 31,61 1,0/14,3 Случай 2.2. Жесткое закрепление с одной стороны и шарнирная опора с другой, материал типа 2 [Case 2.2. Fixed-simply supported ends, material type 2] 1 89,12 89,32 88,96 89,17 87,70 85,59 1,4/4,4 2 259,62 259,94 258,83 260,12 255,70 231,89 1,2/12,2 а. f1 = 2,898 Гц (Hz) б. f1 = 3,00 Гц (Hz) [9] в. f3 = 16,412 Гц (Hz) г. f2 = 16,38 Гц (Hz) [9] Рис. 2. Первая (а, б) и вторая (в, г) изгибные собственные формы и частоты колебаний консольной пневмобалки: а, в - полученные в данной работе (d = 0,03 м); б, г - данные источника [9] [Figure 2. The first (а, б) and the second (в, г) bending mode shapes and natural frequencies of cantilever airbeam: a, в - calculated by the authors (d = 0.03 m); б, г - from source [9]] а. f1 = 12,727 Гц (Hz) б. f1 = 12,85 Гц (Hz) в. f3 = 36,117 Гц (Hz) г. f2 = 35,76 Гц (Hz) Рис. 3. Первая (а, б) и вторая (в, г) изгибные собственные формы и частоты колебаний пневмобалки, жестко закрепленной слева и шарнирно опертой справа: a, в - полученные в данной работе (d = 0,03 м); б, г - данные источника [9] [Figure 3. The first (а, б) and the second (в, г) bending mode shapes and natural frequencies of fixed-simply supported airbeam: a, в - calculated by the authors (d = 0.03 m); б, г - from source [9]] Поскольку в работе [9] имеется по два эталонных результата (для аналитической стержневой и численной оболочечной моделей), а в данной работе выполнены расчеты для пневмобалок с торцами и без, то для относительной погрешности (невязки) результатов в табл. 1 приведены максимальное εmax и минимальное εmin значения. Для результатов, полученных в данной работе, достигнута сеточная сходимость, поскольку собственные частоты при сгущении расчетной сетки с 0,03 до 0,015 м изменяются не более чем на 0,05-0,5%. Расхождение между результатами авторов статьи и эталонными результатами для стержневой модели [9] достигает 14,3%. Такая разница в результатах может объясняться различием в моделировании торцевых зон и граничных условий, а также особенностями построения стержневой аналитической модели: в частности, была принята гипотеза, что сечения пневмобалки остаются круглыми и после деформации [9], однако, расчеты с использованием оболочечных моделей показывают, что это условие не всегда выполняется (рис. 4, a). Относительная погрешность для оболочечных конечно-элементных моделей не превышает 1%. а. f24 = 259,62 Гц (Hz) б. f13 = 223,06 Гц (Hz) Рис. 4. 24-я (a) - вторая изгибная и 13-я (б) собственные формы и частоты колебаний пневмобалки для случая 2.2 с торцами [Figure 4. The 24th mode shape (а), corresponded to the second bending shape, and 13th mode shape (б) of airbeam with ends (case 2.2)] Кроме того, из одномерной стержневой модели могут быть получены только изгибные формы колебаний, а в трехмерных конечно-элементных моделях наблюдаются оболочечные формы колебаний, обусловленные волнообразованием как в продольном, так и в кольцевом направлении. Поэтому второй изгибной форме колебаний, как правило, соответствует гораздо больший номер собственной формы. В частности, на рис. 4, a для случая 2.2 (с торцами) показана 24-я форма колебаний, которая является второй изгибной. Ее частота f24 = 259,62 Гц выше, чем, например, для 13-й формы колебаний (оболочечной), показанной на рис. 4, б, т.к. f13 = 223,06 Гц. Образование одной полуволны в продольном направлении и восьми полуволн в кольцевом соответствует 13-й форме. С учетом сделанных замечаний можно считать тестовую задачу успешно решенной, что позволяет использовать примененный подход для анализа собственных частот и форм колебаний воздухоопорной оболочки. 3. Исследование собственных частот и форм колебаний воздухоопорной оболочки 3.1. Описание расчетной модели Объектом исследования выбрана воздухоопорная оболочка на прямоугольном плане размером 20×50 м, ее высота в высшей точке составляет 8,6 м, радиус скругления углов оболочки 2 м. Как правило, рабочее внутреннее давление под оболочкой составляет порядка p = 150 Па [15]. Оно может быть повышено в несколько (обычно от двух до четырех) раз в случае неблагоприятных погодных условий - сильного ветра или снегопада. Собственные частоты и формы колебаний зависят от обобщенной жесткости и распределения массы конструкции, поэтому по их изменению можно оценить, как меняется обобщенная жесткость воздухоопорной оболочки при повышении давления. Начальная равновесная форма воздухоопорной оболочки не является канонической и была предварительно определена численно с помощью программы ixForten 4000. Физико-механические параметры линейноупругого ортотропного материала оболочки (условной технической ткани с покрытием) приведены в табл. 3. Таблица 3 Физико-механические характеристики материала [Table 3. Physical and mechanical properties of material] Характеристика [Property] Величина [Value] Плотность ρ, кг/м3 [Density ρ, kg/m3] 1200 Толщина h, м [Thickness h, m] 0,001 Модуль упругости в направлении нитей основы (вдоль короткой стороны) E1, МПа [Modulus of elasticity in warp direction (along the short side) E1, MPa] 600 Модуль упругости в направлении нитей утка (вдоль длинной стороны) E2, МПа [Modulus of elasticity in weft direction (along the long side) E2, MPa] 300 Модуль сдвига G12, МПа [Shear modulus G12, MPa] 30 Коэффициент Пуассона ν12 [Poisson ratio ν12] 0,1 Коэффициент Пуассона ν21 [Poisson ratio ν21] 0,2 Расчет собственных частот и форм колебаний воздухоопорной оболочки произведен по методике, использованной в п. 2 данной работы. Для определения оптимального размера конечного элемента и получения сеточной сходимости расчетная сетка была рассмотрена в трех вариантах, параметры которых приведены в табл. 4. Расчетные сетки для вариантов 2 и 3 можно увидеть на рис. 5. Таблица 4 Параметры расчетной сетки для рассмотренных моделей воздухоопорной оболочки [Table 4. Parameters of finite-element mesh used for the present models of air-supported structure] Модель [Model] 1 2 3 Размер КЭ, d, м [Size of FE, d, m] 1,50 1,00 0,75 Количество КЭ [Number of FE] 1273 2712 4795 3.2. Результаты расчета воздухоопорной оболочки На рис. 5 для двух вариантов расчетной сетки показаны цветовые карты максимальных главных растягивающих напряжений, рассчитанных при действии внутреннего избыточного давления под оболочкой p = 150 Па. На большей части оболочки напряжения меняются в незначительном диапазоне от 1,25 до 1,47 МПа, это говорит о том, что такую оболочку можно приближенно считать «рав- нонапряженной» и, следовательно, рациональной. Величины напряжений и картина их распределения практически совпадают. На рис. 6 показаны величины первых пяти собственных частот для моделей 1, 2 и 3. Влияние размера конечного элемента сильнее сказывается на более высоких частотах, поскольку размер полуволны становится меньше и более грубая сетка хуже аппроксимирует деформированную поверхность оболочки. На рис. 7, a-д показаны первые пять собственных форм колебаний воздухоопорной оболочки для модели 3 (сетка конечных элементов не показана). Формы собственных колебаний охарактеризованы парой параметров (m, n), где m - число полуволн вдоль длинной, а n - число полуволн вдоль короткой стороны оболочки. На рис. 7, е для сравнения показана первая собственная форма, полученная при игнорировании дополнительной жесткости, создаваемой внутренним избыточным давлением: поскольку мембрана не обладает жесткостью на изгиб, то первой собственной форме соответствует локальное выпучивание оболочки вблизи одного из узлов расчетной сетки. а. d = 1,00 м (m), 2712 КЭ (FE) б. d = 0,75 м (m), 4795 КЭ (FE) Рис. 5. Максимальные главные напряжения при действии внутреннего избыточного давления p = 150 Па [Figure 5. Maximum principal stress due to the action of internal pressure p = 150 Pa] Рис. 6. Сеточная сходимость решений для моделей 1, 2 и 3 (зеленый, красный и синий столбцы соответственно) [Figure 6. Mesh convergence for the models 1, 2 and 3 (green, red and blue respectively)] а. (1, 2), f1 = 1,581 Гц (Hz) б. (1, 3), f2 = 2,073 Гц (Hz) в. (2, 2), f3 = 2,384 Гц (Hz) г. (2, 3), f4 = 2,394 Гц (Hz) д. (1, 4), f5 = 2,618 Гц (Hz) е. p = 0 Pa, f1 = 0,351 Гц (Hz) Рис. 7. Собственные частоты и формы колебаний при действии внутреннего давления p = 150 Па для модели 3 (а-д); е - первая собственная частота и форма колебаний при отсутствии внутреннего избыточного давления [Figure 7. Natural frequencies and mode shapes for model 3 when internal pressure p = 150 Pa is applied (а-д); е - the first natural frequency and mode shape when internal pressure p is neglected] 0 1503004506007509001050 0 1503004506007509001050 p, Pa p, Pa а б Рис. 8. Зависимость собственных частот (а) и квадратов собственных частот (б) для форм (1, 2) и (1, 3) от величины избыточного давления p (модель 3, d = 0,75 м) [Figure 8. Dependence of natural frequencies f (а) and squares of natural frequencies f 2 (б) on internal pressure p for the mode shapes (1, 2) and (1, 3). Model 3, d = 0.75 m] На рис. 8 показаны графики зависимости собственных частот (а) и квадратов собственных частот (б) для форм (1, 2) и (1, 3) от величины избыточного давления p в диапазоне от 50 до 1050 Па. Следует отметить, что с ростом давления собственные формы трансформируются, меняется их порядок следования друг за другом в связи с изменением предварительного напряженно-деформированного состояния оболочки. Формы колебаний (1, 2) и (1, 3) были выбраны потому, что их конфигурация оставалось неизменной в рассмотренном диапазоне избыточного давления. В книге [4, с. 437-438] при модальном анализе прямоугольных предварительно растянутых мембран показано, что квадраты собственных частот прямо пропорциональны растягивающему усилию. В воздухоопорных сооружениях растягивающие усилия практически прямо пропорциональны величине избыточного давления (при отсутствии других нагрузок), поэтому в указанном интервале давлений наблюдается почти линейная связь между величиной избыточного давления и квадратами собственных частот, показанная на рис. 8, б. 4. Анализ результатов и выводы На основе проведенного обзора научной литературы была найдена и успешно решена тестовая задача по расчету собственных частот и форм колебаний предварительно напряженной цилиндрической оболочки из ортотропного материала. Получено хорошее соответствие между результатами из источника и результатами авторов статьи. Расхождения по частотам колебаний для оболочечных моделей не превышают 1,8%, расхождение с аналитическим решением для стержневой модели не превышает 15%, что объясняется упрощениями, используемыми при построении аналитической стержневой модели, а также несколько различными подходами к учету влияния торцов оболочки. Методика решения тестовой задачи была применена для численного исследования собственных частот и форм колебаний воздухоопорной оболочки на прямоугольном плане размером 20×50 м. В работе выявлена практически линейная зависимость между величиной внутреннего давления и квадратами собственных частот, что соответствует результатам других исследователей.

Об авторах

Николай Андреевич Мокин

Национальный исследовательский Московский государственный строительный университет

Автор, ответственный за переписку.
Email: mokiavelli@mail.ru

аспирант кафедры строительной и теоретической механики, Национальный исследовательский Московский государственный строительный университет. Область научных интересов: расчеты воздухоопорных сооружений

Ярославское шоссе, 26, Москва, Российская Федерация, 129337

Алексей Андреевич Кустов

Национальный исследовательский Московский государственный строительный университет

Email: alexeykustov@outlook.com

аспирант кафедры металлических и деревянных конструкций, Национальный исследовательский Московский государственный строительный университет. Область научных интересов: исследования сооружений (в том числе воздухоопорных), изготовленных из технических тканей с покрытием

Ярославское шоссе, 26, Москва, Российская Федерация, 129337

Михаил Иоакимович Ганджунцев

Национальный исследовательский Московский государственный строительный университет

Email: oppmgsu2014@yandex.ru

доцент, кандидат технических наук, доцент кафедры строительной и теоретической механики, Национальный исследовательский Московский государственный строительный университет. Область научных интересов: динамика сооружений, нелинейная строительная механика

Ярославское шоссе, 26, Москва, Российская Федерация, 129337

Список литературы