CC BY
35
ЧИСЛЕННЫЙ АНАЛИЗ УСТОЙЧИВОСТИ РАВНОВЕСИЯ
ПРОСТРАНСТВЕННЫХ ФЕРМ
В ГЕОМЕТРИЧЕСКИ НЕЛИНЕЙНОЙ ПОСТАНОВКЕ
В.В. ГАЛИШНИКОВА, канд. техн. наук, доцент
Российский университет дружбы народов, Москва
В работе излагаются алгоритмы вычисления критических состояний пространственной фермы и продолжения решения в сингулярных точках, реализующие методики, разработанные автором. Приводится пример расчета сетчатой оболочки на устойчивость равновесия по разработанной программе.
КЛЮЧЕВЫЕ СЛОВА: устойчивость, пространственная ферма, сетчатая оболочка
Вычисление точек бифуркации гибких, подверженных большим деформациям конструкций, а также продолжений траекторий их нагружения после прохождения этих точек остается одной из сложнейших задач вычислительной механики. Очевидно, что задача устойчивости равновесия таких конструкций может быть достоверно решена лишь при учете предшествующей истории нагру-жения, то есть в геометрически нелинейной постановке. Однако, несмотря на то, что численные методы геометрически нелинейного анализа хорошо развиты, в окрестностях точек бифуркации возникают вычислительные трудности, препятствующие созданию надежно работающих алгоритмов продолжения решения в этих точках. Формулировка и проверка достоверности методов геометрически нелинейного анализа конструкций существенно затрудняется недостатком точных аналитических решений задач для конструкций, подвергающихся значительным деформациям перед потерей устойчивости. Достоверность и точность методик, предлагаемых в данной работе, оценивалась на тестовых примерах, для которых автором получены аналитические решения [1].
Для выполнения геометрически нелинейного расчета пространственных ферм на деформации по методу конечных элементов автором разработана методика, основанная на модификации метода постоянных дуг (Arc Length Method), в которой используется процедура итерационного уточнения секущей матрицы жесткости конструкции [2]. В качестве характеристической кривой решения используется нормализованная траектория нагружения. Алгоритмы, описываемые в настоящей работе, являются частью обобщенного метода геометрически нелинейного расчета пространственных ферм, включающего расчет на деформации и потерю устойчивости.
Вычисление сингулярных точек. Конфигурацию конструкции будем называть сингулярной, если определитель матрицы касательной жесткости Kt конструкции в данной конфигурации равен нулю. Конфигурация неустойчива, если коэффициент нагружения на продолжении траектории убывает. Шаговые процедуры решения, использующие обращение матрицы Kt в окрестностях сингулярных точек прекращают работать. Следовательно, необходимо введение дополнительных алгоритмов, решающих следующие задачи:
- выявление конфигураций конструкции, близких к критическим;
- вычисление векторов перемещений и сил системы в критической точке;
- продолжение траекторий нагружения за критическими точками.
Выявление конфигураций с почти сингулярной матрицей касательной жесткости не представляет сложности с теоретической точки зрения. Достаточно лишь сравнить вычисленное значение определителя с некоторым предельным значением. На практике возникают вычислительные сложности, так как значение определителя часто выходит за пределы области типа данных "double" компьютера. Изменение знака определителя также не может служить надежным индикатором, так как для симметричной конструкции в одной и той же конфигурации могут менять знак два собственных значения; при этом знак определителя не меняется. Затруднение удается преодолеть, используя методику отслеживания знаков коэффициентов разложения матрицы касательной жесткости Kt:
detKt = detS detD detST = detD = d0d1...dn_1. (1)
Здесь S и ST - левая и правая треугольные матрицы, а D - диагональная матрица, которая может содержать положительные и отрицательные коэффициенты dk. Определитель матрицы Kt равен произведению коэффициентов диагональной матрицы D. Если количество отрицательных коэффициентов dk на шаге изменилось, это означает, что ферма прошла сингулярную конфигурацию.
В известном методе бисекции [3] для вычисления сингулярных точек используется следующий алгоритм: вычисляются матрицы касательной жесткости конструкции для состояния предшествующего предполагаемой критической точке, и состояния, следующего за ней. Затем интервал, содержащий критическую точку, разбивается пополам. Процедура продолжается до достижения заданной точности. Однако этот подход непригоден для вычисления точек бифуркации, так как в этих точках происходит ветвление решения. Если начальная точка интервала лежит на основной ветви траектории нагружения, а конечная - на второстепенной ветви, то бисекция этого интервала даст точку, не принадлежащую решению.
Автором разработана методика вычисления точки бифуркации, основанная на решении общей проблемы собственных значений. Пусть N - почти сингулярная точка, предшествующая сингулярной точке S траектории нагружения. Перемещения и реакции в точке N могут быть вычислены при помощи обычной процедуры метода постоянных дуг. Инкременты Ad перемещения из точки N в сингулярную точку S записываются как функции неизвестного инкремента коэффициента нагружения АЯ от точки N к точке S:
Ad = АХ dt + f , (2)
где d t - вектор перемещений от модельной нагрузки при жесткости фермы, вычисленной для состояния N, а f - вектор перемещений от неуравновешенных сил того же состояния. Получающееся в результате выражение для матрицы касательной жесткости фермы в сингулярном состоянии является нелинейной функцией неизвестного инкремента коэффициента нагружения АЯ:
KS = K0 + АХ K1 (АХ) (3)
Определитель матрицы касательной жесткости в сингулярной конфигурации полагается равным нулю:
det(K 0 + АХ K1) = 0 (4)
Значение АЯ, приводящее значение определителя матрицы K^ к нулю, вычисляется итерационно. В первой итерации члены, нелинейные относительно АЯ, опускаются. Значение АЯ определяется из решения общей задачи на собственные значения (4), соответствующей уравнению критического равновесия
K 0 x = АА K1 x (5)
Здесь АЯ, х - собственные значения и собственные векторы матриц (K0, Ki). Так как в данном случае требуется отыскать лишь наименьшее собственное значение, то наиболее эффективен метод обратных векторных итераций. Если состояние N, для которого вычисляются матрицы K0 и K1, близко к сингулярному состоянию, то итерации сходятся достаточно быстро.
В последующих циклах итераций на шаге от состояния N к состоянию S нелинейные члены, входящие в матрицу касательной жесткости состояния S аппроксимируются с использованием значения инкремента коэффициента на-гружения АА, полученного в предыдущем цикле итераций. Процедура может быть кратко описана следующим образом:
- цикл 0:
АА = 0 ^ вычисление K(0) ^ решение проблемы СЗ относительно АА(0);
- цикл s :
AA(s-1) ^ вычисление K(s) ^ решение проблемы СЗ относительно AA(s).
Найденное значение AAmjn = min | АА | подставляется в выражение (2) для вычисления инкремента перемещения в сингулярное состояние:
Ad = АА min dt + f (6)
Так как вектор Ad получен с использованием приближенного значения матрицы касательной жесткости, то точка S, вычисленная с его использованием, не будет являться точной сингулярной точкой кривой нагружения. Описанная выше процедура повторяется, используя S как новую почти сингулярную точку. Алгоритм прерывается, когда изменение инкремента коэффициента нагружения становится меньше заданного предельного значения.
Продолжение траекторий нагружения. Поведение траекторий нагружения за сингулярными точками зависит от того, является ли точка предельной или точкой бифуркации. Вычислительные методы для этих двух типов сингулярности также существенно отличаются. Принятая в работе классификация сингулярных точек основана на общей теории бифуркации. Разрешающие уравнения для полного перемещения u конструкции запишем в следующем виде:
g (u) + А q = 0 , (7)
где g - вектор гладких функций некоторого параметра траектории z, A(z) - коэффициент нагружения. Дифференцируя выражение (7) по параметру z, полу-
d u d А
чим K--+ q— = 0, где K =
dz dz
dg dg
дщ дип
(8)
Обозначим значение матрицы К в сингулярной точке как К0, а собственный вектор, соответствующий собственному значению ц = 0 матрицы К0 - как х0. Для собственного состояния (0,хд) частичная проблема СЗ и СВ матрицы К имеет вид:
Кохо = 0 . (9)
Умножив уравнение (8) слева на х0, с учетом (9) получим:
х0я ^ = 0. (10)
dz
Сингулярные точки классифицируются в соответствии с полученным в выражении (10) значением скалярного произведения. Если скалярное произведе-
ние отлично от нуля, то сингулярная точка является предельной. Если произведение равно нулю, то сингулярная точка представляет собой точку бифуркации:
Ф 0 ^ (ид,А,0) - предельная точка; (11)
= 0 ^ (ид,Хд)- точка бифуркации. (12)
Из выражений (10) и (11) видно, что производная коэффициента нагруже-ния по параметру траектории г в предельной точке равна нулю:
dXJ dz = 0. (13)
Продолжение траектории нагружения в предельной точке. Рассмотрим уравнение (9) для частного случая предельной с учетом свойства (11):
К0 ^ = 0 . (14)
dz
Сравнивая уравнения (14) и (9) можно убедиться, что производная вектора перемещений по параметру г пропорциональна собственному вектору Х0 :
^ = 6х0 . (15)
dz
Учитывая это, в первом цикле итераций на первом шаге продолжения решения будем задавать инкремент перемещения в направлении собственного вектора х0:
Лd(0) = с Х0 с заданной длиной
Лd
(0)
(16)
Затем вычисляется матрица секущей жесткости конструкции К8 (Лd(0)). В отличие от матрицы касательной жесткости, эта матрица не является сингулярной, и её определитель не равен нулю. Процедура метода постоянных дуг может быть возобновлена. Однако система уравнений метода конечных элементов в окрестности предельной точки может быть плохо обусловленной. В этом случае для её решения используется специальные методы, не описываемые в данной статье.
Продолжение решения в точках бифуркации. Рассмотрим условие (12) в точке бифуркации. Пусть производная перемещения на шаге, следующем за точкой бифуркации, представляет собой линейную комбинацию собственного вектора Х0 и некоторого вектора zo , ортогонального вектору Х0 :
^ = 01Х0 +е2 zo. (17)
dz
Подставим полученную производную в уравнение (8). С учетом условия (9) получим:
02 К0 zo = ^q . (18)
dz
Принимая коэффициент 02 равным производной X, стоящей в правой части и сокращая уравнение на эти множители, получим следующую систему уравнений для точки бифуркации:
К0 z 0 = q . (19)
Очевидно, что данная система уравнений не может быть решена относительно z0, так как матрица К0 сингулярна.
В точках бифуркации (ветвления решений) поведение конструкции может фундаментально меняться. Если конструкция следует вторичной ветви траектории нагружения, то касательные к траектории слева и справа от точки бифурка-
ции меняют наклон. Это нарушение неразрывности перемещений сопровождается изменениями свойств численной модели конструкции.
Итерационное решение системы уравнений зависит от величины и знака определителя матрицы касательной жесткости в окрестности траектории на-гружения. Численные исследования определителей матриц касательной жесткости трехстержневых ферм [4, 5] показали, что матрица касательной жесткости конструкции в окрестности траектории продолжения решения остается почти сингулярной. На самой траектории продолжения решения эта матрица сингулярна. Исследование аппроксимаций рядами Тейлора показало, что линейная аппроксимация разрешающих уравнений в точке бифуркации не может привести к надежному решению. Поэтому в данной работе предлагается альтернативный подход, который основывается на использовании матрицы полной жесткости конструкции на шаге решения. Разрешающие уравнения записываются для конца шага нагружения:
К Е dE = qE + гЕ, (20)
где КЕ - матрица полной жесткости конструкции, вычисленная для полного перемещения dE в конце шага; гЕ - вектор полных реакций от полной нагрузки qE в конце шага нагружения. Первый шаг продолжения решения в точке бифуркации аналогичен первому шагу в предельной точке:
dE = ds + с Х0, (21)
^ = Г + с Г). (22)
Здесь dS и rS - векторы полных перемещений и полных реакций в начале шага; Х0 - собственный вектор, соответствующий наименьшему собственному значению матрицы касательной жесткости в начале шага; г0 - вектор реакций от перемещения Х0. Последующие шаги алгоритма описаны ниже.
1. Вычисляется матрица полной жесткости конструкции и система уравнений МКЭ решается относительно перемещений dE и реакций г^
К E ^ = qE + гб .
2. Повторно вычисляется матрица полной жесткости конструкции и вычисляется вектор нагрузки: qE = КE dE - rE .
3. Нагрузка qE раскладывается на составляющую, параллельную модельной нагрузке и составляющую, перпендикулярную ей: qE =0 qt + Лqc, где
т
qt Лqc = 0. Система уравнений МКЭ решается относительно перемещений и реакций от перпендикулярной составляющей Лqc : КE Лdc = Лqc + Лгс .
Векторы полных перемещений, реакций и нагрузки корректируются:
^ , rE -Лгс, qE.
4. Для исправленного состояния вычисляются неуравновешенные силы ДqS на шаге, и система уравнений решается относительно корректирующих векторов от неуравновешенных сил: КE Лds = Лqs - Лг5 .
Корректируются векторы полных перемещений, реакций и нагрузок:
^ , гБ -Лг5, qE-Лqs.
5. Шаги с 1 по 4 повторяются до тех пор, пока норма погрешности в конце шага не станет меньше заданного значения.
Анализ устойчивости равновесия сетчатой оболочки. С использованием разработанной на основе описанных алгоритмов тестовой платформы автором выполнено проектирование и численное исследование устойчивости равновесия сетчатой оболочки, предназначенной для строительства в составе торгового 46
комплекса в г. Саратове. Оболочка перекрывает площадь 14 х 14 м в плане и состоит из сферического сетчатого купола и опорной части, имеющих общую вращательную ось симметрии. Форма оболочки предложена архитекторами проекта. Радиус основания купола равен 6м, стрела подъема купола равна 2,7 м. Опорная часть имеет радиус основания 10,0 м и высоту 6,2 м. Геометрическая схема оболочки показана на рисунке 1.
В процессе работы над проектом выполнены следующие задачи:
1) разработка и программная реализация на базе платформы Java генератора модели, выполняющего геометрический расчет конструкции и создающего её расчетную модель;
2) объединение генератора с расчетным приложением (тестовой платформой геометрически нелинейного анализа);
3) расчет конструкции оболочки на действие расчетных сочетаний нагрузок и подбор сечений стержней оболочки в соответствии с требованиями строительных норм;
4) анализ устойчивости равновесия конструкции при действии всех расчетных комбинаций нагрузок: вычисление нагрузок, при которых происходит потеря устойчивости; определение вида потери устойчивости и вычисление форм потери устойчивости конструкции.
Модель сетчатой оболочки, созданная при помощи разработанного генератора, рассчитывалась на прочность и деформации. Расчет выполнялся на девять различных сочетаний нагрузок. В сочетаниях учитывались: нагрузка от веса конструкций, снеговая и ветровая нагрузки, взятые в соответствии с действующими нормами. Сечения стержней проектировались из круглых тонкостенных труб. Были приняты следующие характеристики стали: расчетное сопротивление 230 МПа, модуль упругости 2,05 -105 МПа. Расчет конструкции оболочки на прочность, подбор сечений стержней и расчет на деформации остаются за рамками данной статьи. Все расчеты выполнялись в соответствии с действую-
щими нормами. В результате была разработана конструкция сетчатой оболочки, подкрепленной четырьмя сетчатыми арками.
Схема сетки оболочки показана на рис. 2. Для разработанной конструкции было выполнено численное исследование устойчивости равновесия. Каждое из сочетаний нагрузок задавалось в виде модельной нагрузки. Расчетное значение нагрузки при этом принималось за условную единичную нагрузку.
а) план
Рис. 2. Схема сетки оболочки
Производился пошаговый нелинейный расчет конструкции до наступления критического состояния. Вычислялась сингулярная точка траектории нагруже-ния, и определялся её тип. Далее вычислялось состояние конструкции после наступления потери устойчивости (форма потери устойчивости). Данные расчета для девяти сочетаний нагрузок приведены в табл. 1. Для описания сочетаний нагрузок использованы следующие обозначения: D - постоянная нагрузка от собственного веса конструкций; £ - симметричная снеговая нагрузка; SPx и SPy - несимметричная снеговая нагрузка, размещенная в положительных направлениях осей х и у; и ЗЫу - несимметричная снеговая нагрузка, размещенная в отрицательных направлениях осей х и у; W - ветровая нагрузка. В табл. 1 приведены значения коэффициента нагружения V® в сингулярной точке для каждого из девяти рассмотренных сочетаний нагрузок. Приведены также значения
скалярного произведения Х0 р собственного вектора матрицы жесткости конструкции в сингулярной точке и вектора нагрузки, а также значения соответствующего угла между этими векторами. В последнем столбце приведены типы сингулярных точек для каждого случая.
Таблица 1. Сингулярные состояния сетчатой оболочки
Сочетание ^ тах Т х0 Р Угол Тип
нагрузок сингулярности
С1: 1,1 * D + 5 2,760 0,00000 90,000 т. бифуркации
С 2 1,1* D + 5Рх 3,996 -0,00762 90,437 предельная т.
С3 1,1* D + 5Ру 4,639 -0,00611 90,350 предельная т.
С 4 0,9 * D + Ж 2,318 0,03694 87,883 предельная т.
С5 1,1* D+Ж + 5 2,268 0,03845 87,796 предельная т.
С 6 1,1* D + Ж+5РХ 2,114 0,02665 88,473 предельная т.
С 7 1,1* D+Ж + 5Ру 2,382 0,03447 88,025 предельная т.
С8 1,1* D + Ж+8Ых 2,114 -0,02665 91,527 предельная т.
С9 1,1 * D+Ж+8Ыу 2,382 -0,03447 91,975 предельная т.
На рис. 3 и 4 показаны формы потери устойчивости конструкции для двух характерных расчетных сочетаний нагрузок, одно из которых является симметричным, другое - асимметричным.
Сочетание 1 включает собственный вес несущих и ограждающих конструкций и симметрично распределенную снеговую нагрузку. Потеря устойчивости в этом случае наступает при значении коэффициента нагружения 2,76 (то есть, при значении нагрузки в 2,76 раз превышающей расчетное). Устойчивость теряет первое кольцо сферического купола. На рис. 3 видно, что чередующиеся узлы смещаются внутрь и наружу этого кольца, образуя "волну". Форма потери устойчивости при этом имеет осевую симметрию. Критическая точка классифицирована как точка бифуркации.
Наиболее опасным оказалось сочетание 6, включающее собственный вес конструкции, снеговую нагрузку на половине купола, в положительном направлении оси XI, а также ветровую нагрузку, действующую в положительном направлении оси у2. Потеря устойчивости происходит при значении нагрузки, превышающем расчетное в 2,11 раз. На рис. 5 видно, что зона потери устойчивости расположена в одной из опорных частей оболочки и включает один узел, перемещающийся к оси оболочки и несколько соседних узлов, перемещающихся наружу. В трех остальных опорных частях оболочки потери устойчивости не происходит.
Результаты расчетов показали, что различные виды сочетаний нагрузок приводят к качественно отличающимся формам потери устойчивости. При отсутствии ветровой нагрузки потеря устойчивости происходит в верхней зоне сферического купола. При действии ветровой нагрузки критической зоной является опорная часть оболочки. Сетчатые арки не теряют устойчивости при любых сочетаниях. Минимальная нагрузка, при которой происходит потеря устойчивости, в 2,14 раз превышает расчетную (сочетания 6 и 8). При этом деформации узлов перед моментом потери устойчивости достигают 150 миллиметров.
Заключение. Рассмотренный пример геометрически нелинейного анализа устойчивости равновесия сетчатой оболочки продемонстрировал, что программное приложение, разработанное на основе описанных в настоящей работе алгоритмов, может быть эффективно использовано для расчета пространственных стержневых конструкций на деформации и устойчивость в геометрически нелинейной постановке.
Л и т е р а т у р а
1. Галишникова В.В. Аналитическое решение нелинейной задачи устойчивости и исследование закритического поведения трехстержневой фермы // Вестник Волгогр. гос. архит.-строит. ун-та. сер.: Естеств. науки. - 2006. - Вып. 6(23). - С. 53—64.
2. Галишникова В.В. Модификация метода постоянных дуг, основанная на использовании матрицы секущей жесткости // Вестник МГСУ. -2009. - № 2. - С. 63—69.
3. Belytscho, T. Nonlinear Finite Elements for Continua and Structures / T. Belytscho, W. Liu, B. Moran. J Wiley & Sons, 2000. ISBN 0-948-749-261.
4. Galishnikova V. V. Stability Analysis of Space Trusses // International Journal for Computational Civil and Structural Engineering. - 2009. - Vol. 5, Issue 1&2. -Pp. 35—44.
5. Galishnikova V.V. Solving the Unsolvable: Unusual Formulations in Computational Mechanics // Proceedings of EG-ICE Conference "Computing in Engineering". Berlin, 2009. Pp. 113—123.
GEOMETRICALLY NONLINEAR NUMERICAL STABILITY ANALYSIS
OF SPACE TRUSSES
V. V. Galishnikova
Algorithms for direct computing of singular states of space trusses at singular states and subsequent continuation of load paths are described in the paper. An example of geometrically nonlinear stability analysis of lattice shell is presented.
