CC BY
34
Красноруцкий Д.А. Методика расчета механики систем связанных стержней по дифференциальной модели // Вестник Пермского национального исследовательского политехнического университета. Механика. - 2016. - № 2. - С. 6988. DOI: 10.15593/ perm.mech/2016.2.06
Krasnorutskiy D.A. The method for analyzing mechanics of thin elastic rod systems using the differential approach. PNRPUMechanics Bulletin. 2016. No. 2. Рр. 69-88. DOI: 10.15593/perm.mech/2016.2.06
вестник пнипу. механика № 2, 2016 PNRPU MECHANICS BULLETIN
http://vestnik.pstu.ru/mechanics/about/inf/
Б01 10.15593/регт.тесЬ/2016.2.06 УДК 539.3
МЕТОДИКА РАСЧЕТА МЕХАНИКИ СИСТЕМ СВЯЗАННЫХ ТОНКИХ УПРУГИХ СТЕРЖНЕЙ ПО ДИФФЕРЕНЦИАЛЬНОЙ МОДЕЛИ
Д.А. Красноруцкий
Новосибирский государственный технический университет, Новосибирск, Россия
О СТАТЬЕ
АННОТАЦИЯ
Получена: 16 марта 2016 г. Принята: 2 июня 2016 г. Опубликована: 30 июня 2016 г.
Ключевые слова:
система связанных стержней, стержневая система, рама, ферма, вантовая конструкция, дифференциальная модель стержня, нелинейная динамика, большие перемещения, метод решения СЛАУ, разреженная матрица
Предложена методика расчета геометрически нелинейного динамического деформирования системы связанных между собой стержней по дифференциальной модели тонкого упругого стержня. Деформирование каждого стержня описывается 12 дифференциальными уравнениями с краевыми условиями. Краевые условия для каждого стержня формируются из условий закрепления, нагружения и условий стыковки стержней между собой. Полученная нелинейная краевая задача решается численно, для этого требуется многократно решать большие СЛАУ, матрицей которых является матрица Якоби совокупности дифференциальных уравнений, краевых условий и условий стыковки для всех стержней, входящих в систему. На главной диагонали такой матрицы стоят блоки от дифференциальных уравнений и краевых условий для каждого отдельного стержня, и эти блоки связаны между собой в остальной части матрицы элементами, полученными из уравнений связи между стержнями в системе. В статье предлагается метод решения таких блочных разреженных СЛАУ, на главной диагонали которых стоят невырожденные блоки, связанные с другими блоками относительно небольшим количеством уравнений. Решение исходной СЛАУ распадается на решение нескольких систем меньшего размера (для каждого блока) и решения СЛАУ размером, равным числу ненулевых строк, связывающих матричные блоки. Поэтому предлагаемый метод эффективен при относительно небольшом по сравнению с размерами матрицы количестве ненулевых строк, связывающих выделенные матричные блоки, стоящие на главной диагонали. В статье приведены тестовые примеры расчета больших перемещений нагруженных рам, результаты расчета по предложенной методике сравниваются с результатами расчета в конечно-элементном пакете ANSYS. Результаты расчетов по методам практически совпали, точность совпадения зависит от дискретизации моделей.
©ПНИПУ
© Красноруцкий Дмитрий Александрович - кандидат технических наук, доцент, e-mail: [email protected] Dmitry A. Krasnorutskiy - PhD in Technical Sciences, Associate Professor, e-mail: [email protected]
THE METHOD FOR ANALYZING MECHANICS OF THIN ELASTIC ROD SYSTEMS USING THE DIFFERENTIAL APPROACH
D.A. Krasnorutskiy
Novosibirsk State Technical University, Novosibirsk, Russian Federation
ARTICLE INFO
ABSTRACT
Received: 16 March 2016 Accepted: 2 June 2016 Published: 30 June 2016
Keywords:
system-related rods, rod system, frame, farm cabling structure, differential model of the rod, nonlinear dynamics, large displacement, method for solving linear systems, sparse matrix
This paper proposes a method for analyzing geometrically nonlinear dynamic deformation of rod system (several interconnected rods) using differential model of a thin elastic rod. The deformation of each rod is defined by 12 differential equations with boundary conditions. Boundary conditions for each rod are formed by kinematics, loads and interconnections between rods. Appropriate nonlinear boundary value problem is solved numerically. It needs many times to solve a large system of linear equations with matrix formed by a Jacobi matrix of differential equations, boundary conditions and conditions of interconnections between the rods. On the main diagonal of this matrix the blocks are formed by differential equations and boundary conditions for each rod. These blocks are linked by elements formed by the equations of connection between the rods. The paper proposes the method for solving sparse linear systems with blocks on the main diagonal, while these blocks are related with each other by a small number of equations. The solution of original sparse system is splitted to few solutions of smaller systems (for each block on the main diagonal) and to the solution of SLAE size of equal to number of nonzero lines linking the matrix blocks. Therefore, the proposed method is effective with a comparatively small (compared with size of the original matrix) count of nonzero lines linked with isolated matrix blocks on the main diagonal. The paper presents the solution of the test problems of large displacements loaded frames. The numerical results of the proposed method are compared with the results of the calculation in the finite element program ANSYS. The results of calculations are practically the same; the accuracy of matches depends on models discretization.
©PNRPU
Введение
Расчету механики стержневых систем, под которыми обычно понимают систему шарнирно или жестко соединенных между собой в узлах стержней (ферменные и рамные конструкции), посвящено очень много научных трудов [1-10]. Под системой связанных между собой стержней будем понимать более широкий класс конструкций, чем просто рамы и фермы. Это могут быть вантовые конструкции, буксируемые системы, висячие тросовые конструкции, гибкие трубопроводы, системы валов, провода участка цепи линии электропередачи, сети, и даже лопасти винта вертолета можно представить как стержни с определенным аэродинамическим профилем.
Обычно для расчета стержневых систем используется метод конечных элементов [1, 2, 4-6, 11-12]. Широкое распространение конечно-элементного подхода можно объяснить универсальностью его применения для стержневых систем с произвольной геометрией. В ряде публикаций подвергается сомнению корректность получаемых в коммерческих конечно-элементных пакетах (таких как NASTRAN, ANSYS, Cosmos, Lira) результатов расчета геометрически нелинейного деформирования стержневых систем [1]. Эти сомнения связаны прежде всего со слабой проработкой создателями этих программных продуктов вопросов устойчивости и бифуркации форм равновесия при геометрически нелинейном деформировании. При одной и той же нагрузке могут существовать несколько равновесных конфигураций, достижение которых зависит от способа приложения этой нагрузки и предыстории деформирования. В [1] отмечается необходимость при рассмотрении задач о больших перемещениях стержневых систем учитывать инерцию и предла-
гается динамический подход к расчету стержневых систем методом конечных элементов, при этом рассматриваются стержневые системы с узловым расположением масс и нагрузок. В [2] развивается геометрически нелинейная теория деформирования стержневых систем с анализом бифуркаций форм равновесия в конечно-элементной формулировке. В работах [11-12] метод конечных элементов используется для анализа нелинейного деформирования и устойчивости пространственных стержней с учетом поперечного сдвига при больших упругих перемещениях. В работе [3] рассчитывают критическую силу Г-образной рамы с помощью дифференциальных уравнений стержней и сравнивают с решением в КЛБТКАК.
Другой аспект, в котором традиционный конечно-элементный подход можно усовершенствовать, - это скорость сходимости итерационного процесса при расчете нелинейных конфигураций. Так, в [2] предлагается метод, позволяющий сохранить в уравнениях метода конечных элементов все нелинейные члены, что, как утверждает автор, ускоряет сходимость итерационной процедуры решения.
В представленной статье предлагается методика расчета нелинейного деформирования системы связанных между собой тонких упругих стержней произвольной формы по дифференциальной модели стержня [13, 14]. Эта модель позволяет включать в расчетную схему распределенные и сосредоточенные нагрузки любой природы. Это могут быть силы веса, инерционные силы (в том числе возникающие при несовпадении центра масс и центра жесткости поперечного сечения), аэродинамические силы, силы от протекающего потока жидкости по шлангу, внутреннее и внешнее демпфирование (точечные демпферы), силы от контактного взаимодействия, электромагнитные силы, температурные деформации и т.д. Кроме того, в расчетную схему можно внести произвольную зависимость любых параметров нагружения и характеристик поперечных сечений от любых других параметров деформирования (перемещений, углов поворота, скоростей, деформаций, времени, текущих координат и т.д.), следовательно, можно задавать следящие нагрузки, а также реактивные силы и моменты.
Используемая дифференциальная модель криволинейного стержня [13, 14] имеет преимущество перед другими известными дифференциальными моделями [15-18] в том, что в уравнения не входит начальная кривизна осевой линии стержня, что позволяет напрямую описывать стержни с изломами осевой линии и со скачками кривизны. На основе работ [19, 20] в работе [14] разработан алгоритм численного решения нелинейной краевой задачи динамического деформирования тонкого стержня. Преимуществом подхода к решению краевой задачи является назначаемая точность вычисления неизвестных функций из решения нелинейной краевой задачи, которая достигается благодаря методу отложенной коррекции, и специальная учитывающая структуру матрицы Якоби «быстрая» процедура факторизации и решения соответствующих систем линейных уравнений, возникающих в итерационном процессе Ньютона [19, 20, 14]. При составлении системы из стержней в матрице Якоби нарушается блочная трехдиагональная структура из-за условий сочленения стержней между собой. В данной статье предлагается метод решения СЛАУ, позволяющий использовать хорошо зарекомендовавший себя подход [19, 20, 14] для факторизации матриц каждого стержня отдельно и получать решение совместной СЛАУ. В качестве примера приведен пример расчета системы четырех стержней, результаты сравниваются с результатами, полученными в конечно-элементном пакете АКБУБ.
1. Дифференциальная модель тонкого стержня
Точки недеформированной осевой линии стержня задаются радиусом-вектором г (£) = xk (£)• ¡к, где Е, - параметр, связанный с естественной координатой - длиной
стержня s; 1123 - орты глобальной (неподвижной) системы координат; x123 (£) - координаты точек осевой линии. Здесь и далее по повторяющимся индексам ведется суммирование от 1 до 3.
Для описания ориентации поперечных сечений стержня до деформирования используется «матрица поворота» (матрица начальной геометрии) р . Направляющие векторы главных осей инерции поперечных сечений недеформированного стержня имеют вид ej (£) = Р^ (£)• ¡к; е3 - направлен по касательной к осевой линии стержня, здесь и далее
неуказанные индексы принимают значения от 1 до 3. Направляющие векторы главных осей поперечного сечения после деформации выражены через компоненты матрицы начальной геометрии р и матрицы поворота при деформировании X: е* = р^ ¡к* = Р^¡п .
Жесткостные характеристики стержня задаются следующими величинами: Е112 -изгибные жесткости; GJ3 - жесткость на кручение; EF - жесткость на растяжение-сжатие. Жесткости также могут зависеть от текущей деформированной конфигурации - текущей кривизны осевой линии или удлинения. Инерционные характеристики стержня: pF - погонная плотность; pJ123 - произведение плотности и моментов
инерции поперечных сечений. Если ось жесткости стержня не проходит через центры инерции поперечных сечений, то задаются координаты центра тяжести поперечного сечения относительно центра жесткости Д12 . Вектор положения центров масс относительно центров жесткости Й = Д1ё* + Д2е* = ^ ¡п , ^ = (Д1р№ + Д2р2к ы .
На стержень могут действовать внешние распределенные нагрузки и моменты. В общем случае они могут зависеть не только от текущей точки на осевой линии, но и от времени г, пространственного положения, текущей кривизны и других параметров:
Я = я(,г,и,x,...), т = тг,и,x,...). Кроме распределенных силовых факторов на
стержень могут действовать сосредоточенные силы и моменты. На концах стержня они задаются через краевые условия. Если сосредоточенные силы и моменты действуют внутри стержня, то они включаются специальным образом в нелинейные уравнения (после применения метода конечных разностей), при этом сохраняется блочный трех-диагональный вид матрицы Якоби, а её размерность увеличивается на 12 с каждой новой точкой приложения сосредоточенных сил или моментов. Такое включение сосредоточенного фактора равносильно разбиению стержня на два стержня, сочлененных в точке его приложения.
Мгновенные угловые ускорения и скорости поворота поперечного сечения стержня относительно главных осей е123 имеют следующий вид:
^ 1 = _^шп Р3ш^кпР2к ~ ^тпР3ш^кпР2к , = е3е2 = шпР3ш^кпР2к ,
^2 = ^шпР3ш^кпР1к + ^шпР3ш'^кпР1к , ^2 = езе1 = ^шпР3ш^кпР1к ,
«3 Хтп^1тХкг$2к + Хтп^1тХкп^2к , «3 в1 в2 Хтп^\тХкп^2к ,
где точка обозначает производную по времени t.
Система уравнений движения тонкого упругого стержня, описывающая любые повороты поперечных сечений при динамическом деформировании, имеет вид [14]
= (Х кг Аю р,,"Х кр Аюг,,) ,
А® ЛЕ = ХтАт^Р»ХкпМт
С =
хг,, ;
ЕЗ" 0 0
0 ЕГ-1 0
0 0 ОГ3-
й-,, =р^ (,) (и +д),, ," ч, (,,...) (1 + е);
м,,, = г (,) « пРпк V, + р^ (,)[(Ц, +Д) Д- (и+Д) Д + (( (,,...) + х;,й - ;р0г ) (1 + 8),
(1)
+
8 =
X ,,Х.,кбк [ Е¥ (,)Ч,] + аАТ, = ( Хк,, Хк,,) ,
- = 1,2,3, р = 2,3,1, г = 3,1,2,
где а - коэффициент температурного удлинения; АТ - изменение температуры; нижний индекс , после запятой означает производную по этой переменной.
Система (1) состоит из 18 нелинейных дифференциальных уравнений первого порядка по , и второго порядка по времени t. Разрешающими функциями системы (1) являются 18 функций: и123 (,, t) - проекции вектора перемещений точек осевой линии стержня на оси глобальной системы координат; Х, (,, t) - компоненты матрицы, описывающей поворот главных осей инерции поперечных сечений при деформировании; й12 3 (,, t) -глобальные проекции вектора внутренних усилий; М123 (,, t) - глобальные проекции внутренних моментов.
На девять компонентов матрицы поворота X (,, t) накладываются дополнительные
зависимости, которые вытекают из общих свойств матриц поворота. Поэтому на практике [14] вместо компонент используются только 3 функции в качестве параметров поворота, что уменьшает число неизвестных функций и число уравнений с 18 до 12. В качестве параметров поворота выступает ш = ю. 1. - вектор конечного поворота:
1 - соб ю
Х,.,. = 1 "Г" (
ю
бШ ю
ю
"юк +"
1 - соб ю
ю
-ю,ю..
бШ ю
К =-ю. +
ю
1 - соб ю
ю
-ю,. юк, ю = ю =
2 2 2 + ю2 +ю3 ,
г = 1,2,3, . = 2,3,1, к = 3,1,2.
Таким образом, имеем 12 дифференциальных уравнений и 12 неизвестных функций
^ 1,2,3 , ю1,2,3 , й1,2,3 , М1,2,3 .
2. Нелинейная краевая задача для системы стержней
Нелинейная краевая задача для системы (1) решается численно [14]. Значения разрешающих функций разыскиваются в дискретных точках по длине стержня, производные функций заменяются конечными разностями, уравнения (1) удовлетворяются в точках посередине между узловыми значениями. Погрешность определения разрешающих функций при такой аппроксимации и метод снижения этой погрешности до желаемого уровня приведены в [14]. Затем полученная система нелинейных уравнений решается по методу Ньютона. На каждой итерации матрица СЛАУ (матрица Якоби системы нелинейных уравнений) имеет блочную структуру [14, 19] и факторизуется следующим образом:
А =
ГI 0
Р2 I
0 Рэ
0 0
V 8 82
0 0 I
0
6з
Г А1 С1 0
Б2 А2 С2
0 Бэ Аз
0 0 0
V Б 0 0
0
0
0
I
(РN+1 + 8N )
0 ¥ а 0 0
0
1
0 0 >
0 0
0 0 =
А N CN
Б N+1 А Г%■N+1 J
а1 С1 0 ... 0 0 1
0 а2 С2 ... 0 0
0 0 аэ ... 0 0
0 0 0 ... ^ CN
0 0 0 ... 0 а N+1 У
(2)
где А - матрица Якоби (матрица, составленная из частных производных системы функций, полученной при применении численной методики [19] к решению нелинейной краевой задачи); А1, аг. - матрицы размером 12*12, матрицы Бг., С., р 1, Б , имеют размер 6*12; N - число разбиений стержня при применении численного алгоритма [14] для решения нелинейной краевой задачи; (N +1) - количество узлов; I - единичная матрица 12х12.
Таким образом, имеется эффективный способ [14, 19, 20] решения СЛАУ двухточечной нелинейной краевой задачи для одного стержня, но он не подходит для решения задачи для системы связанных между собой стержней. В этом случае блочная матрица (2) для каждого стержня сама будет являться блоком, стоящим на главной диагонали еще большей матрицы СЛАУ. Для системы п стержней матрица Якоби будет иметь следующую структуру:
А,, А!
А =
11
V Ам
12
п,2
А
А1,п
п,п J
■ А« С11) 0. . 0 0 Н(21) 0. . К21) 1 1 • 1 • 1 Н1' 0 .. К1 ■
В® А21) С(1). . 0 0 0 0. . 0 1 1 1 • 1 1 • 1 0 0 .. 0
0 в<1) А31) . . 0 0 0 0. . 0 1 1 1 • 1 1 1 1 • 1 1 1 0 0 .. 0
0 0 0. А(1) N1-1 С(1) N1-1 0 0. . 0 1 ■ 1 1 . 1 ■ 1 1 . 1 1 0 0 .. 0
В(1) 0 0. .. В(1) • "N1 А(1) N1 2 _ 0. . Z(21) 1 1 . 1 1 . 1 1 0 7(1) •• Zn
Н( 2) 0 0 . . 0 К2) А((2) С12) . . 0 1 1 • 1 • 1 нп2) 0 .. кп2)
0 0 0. . 0 0 В22) Л22) . . 0 1 1 1 • 1 1 1 1 • 1 1 1 0 0 •• 0
Ь2) 0 0 . . 0 ^2) б(2) 0. А(2) • N2 1 • 1 1 . 1 • 1 1 . 1 1 т( 2 ) Ьп 0 7( 2) •• Zn
1 1 1 1
••• ••• ••• • • ••• ••• ••• • • ••• 1 ф 1 ф 1 1 *•• •••
Н(п) 0 0 . . 0 К(п) Н2п) 0 . . К<2п) 1 1 . 1 1 . 1 1 л(1п) С(п) м •• 0
0 0 0 . . 0 0 0 0 . . 0 1 1 1 • 1 1 1 1 • 1 1 1 В(п) 2 Л2п) •• 0
Ьп) 0 0 . . 0 Z(п) Ь(2п ) 0 . . Z(2n) 1 • 1 1 1 • • 1 1 1 • 1 Б(п) 0 А(п) •• ^N01) _
, (3)
где Лй - матрицы Якоби (2) для г -го стержня (г = 1...п ), внедиагональные блоки Л. (г Ф]), они содержат связи между диагональными блоками Лй; верхний индекс в скобках у матриц указывает на номер стержня; матрицы Н.) и К.) обусловлены связью начала г -го стержня с началом и концом . -го стержня соответственно; матрицы Ь. и Z(г) обусловлены связью конца г -го стержня с началом и концом . -го стержня соответственно; общая размерность
п
матрицы (3) - 12N, где N = ^ Ык , Ык - число разбиений к -го стержня.
к=1
Решение СЛАУ с матрицей (3) можно получить, например, вычислив обратную матрицу. В силу того что размерность матриц (2) и (3) существенно больше (в среднем от одного до 3 порядков) размерности своих блоков Л г, Вг, Сг, нахождение обратной матрицы невыгодно с вычислительной точки зрения и с точки зрения экономии памяти. Здесь можно поступить так же, как и в случае одного стержня - провести ЬЦ-разложение матрицы (3), как в (2), и учесть особенности строения этой матрицы. Такое ЬЦ-разложение матрицы (3) требует, помимо хранения матриц Якоби (2), дополнительных массивов памяти для каждого стержня, что увеличивает требуемый объем памяти пропорционально количеству стержней.
На рис. 1 представлена структура матрицы Якоби (3) и её ЬЦ-разложения для системы 4 стержней. Элементарный блок - матрица 6*12, белым цветом указаны нулевые элементы, разными цветами - разные типы матриц-блоков: например, красным цветом слева указаны блоки, полученные из уравнений деформирования стержня, синим и зеленым указаны блоки краевых условий или условий стыковки между стержнями.
В данной статье предлагается модифицированный алгоритм решения СЛАУ с матрицей (3), который требует значительно меньше памяти, чем при прямом обращении или факторизации матрицы СЛАУ и факторизации с учетом структуры блоков (см. рис. 1). В задаче о деформировании системы стержней внедиагональные блоки Лг. формируются
Рис. 1. Структура матрицы Якоби и её ЬИ-факторизации для системы 4 стержней
из условий связи между стержнями. Каждый стержень может иметь только 12 краевых условий (в том числе и связей с другими стержнями) по 6 на каждом конце. Поэтому число ненулевых строк у внедиагональных блоков ограничено 12 на каждый стержень и, следовательно, чем больше разбиений, тем меньше относительный объем памяти, требуемый для хранения А^. Согласно предлагаемому ниже методу отношение объема памяти для решения СЛАУ к объему памяти для хранения исходной разреженной матрицы имеет обратную пропорциональность числу разбиений каждого стержня и прямую пропорциональность количеству стержней в системе.
Предлагаемый ниже метод решения СЛАУ с матрицей (3) может использовать хорошо отработанный аппарат [14, 19, 20] решения СЛАУ с матрицей (2), при этом каждый блок в (3) связан с каким-либо другим несколькими уравнениями, число которых значительно меньше, чем общее число уравнений.
Для такого алгоритма решения СЛАУ отношение объемов, требуемых для хранения внедиагональных блоков и блоков на главной диагонали (матриц Якоби для каждого стержня), зависит от количества стержней в системе и не зависит от дискретизации (количества узловых точек, в которых разыскиваются неизвестные функции) каждого стержня.
3. Метод решения СЛАУ
Пусть имеется система линейных алгебраических уравнений
Ах = Ь, (4)
где А - квадратная матрица размерностью N х N, | А| Ф 0 .
Исходную матрицу СЛАУ (4) представим в виде суммы двух матриц:
А = Б + С, (5)
где Б - блочно-диагональная матрица, |б| ф 0, а матрица с содержит п ненулевых
строк ( п < N ), которые описывают связи между блоками на диагонали матрицы А . Тогда из СЛАУ (4) получим
(Б + С)х = Ь ; х = [(Е + СБ 1 )б]_1 Ь ; Бх = (е + СБ1 )_1 Ь = Ь, (6)
где Е - единичная матрица размером N х N .
Решение СЛАУ (6) эквивалентно решению СЛАУ (4). Из выражения (6) можно заметить, что если все элементы матрицы С равны нулю, тогда Ь = Ь и система (6) превращается в (4). В другом случае необходимо вычислить Ь по формуле (6). При этом необ-
ходимо обратить или факторизовать матрицу (Е + СБ 1), которая по размеру идентична
исходной матрице А (N х N ). Из выражения (6) неочевидно, что для вычисления Ь вместо (Е + СБ1) можно обращать матрицу меньшего размера, а именно размера ( п х п ), где
п - количество ненулевых строк в матрице С . Чтобы это стало более очевидным, представим Ь в следующем виде:
Ь = Ь + АЬ, (7)
тогда из (6) с учетом (7) получим
АЬ = ( Е + СБ-1 )-1 Ь - Ь; АЬ = (Е + СБ-1 )-1 Гь -(Е + СБ-1) Ь ]; (8)
АЬ = -( Е + СБ-1 )-1 СБ-1Ь.
Выражение (8) определяет слагаемое в (7), которое показывает, как должна измениться правая часть СЛАУ, чтобы осуществить переход от исходной системы (4) к системе (6). Размер вектора-столбца АЬ совпадает с размером вектора-столбца Ь, однако можно заметить, что в выражении (8) есть произведение матрицы С на вектор-столбец Б-1Ь . Если в матрице С есть нулевые строки, то соответствующие элементы АЬ будут всегда равны нулю вне зависимости от значений элементов вектора столбца Б-1Ь или матрицы (е + СБ-1) . Этот факт дает основание рассматривать вместо (8) следующее
выражение, которое позволяет рассчитывать ненулевые элементы АЬ без выполнения «лишних» операций:
АЬ = -
(I+С Б )-
С Б-1Ь
(9)
где I - единичная матрица размером п х п; матрица С получена из С путем вычеркивания нулевых строк (С имеет размерность пхN); вектор-столбец АЬ получен вычеркиванием строк из АЬ , которые в С являются нулевыми, а матрица Б получена из Б-1 вычеркиванием столбцов (I) имеет размерность N х п ), номера которых соответствуют номерам нулевых строк в матрице С. Таким образом, вместо матрицы (Е + СБ-1)
с размерами N х N имеем матрицу (I + С Б) размерами п х п, где п - число ненулевых строк в матрице С .
Выражение (9) можно получить по-другому. Сгруппируем ненулевые строки (8) и введем обозначения:
'X -1 У ^ . (
сБ-1 =( X -1 У), (I + СII) = X,
СБ-1
0 0
Е + СБ-1 =
'X У ^
у0 IN-п ,
(е + СБ-1 )
'X-1 - X-1У ^
0
-п У
тогда из (8) следует
г х-1 -Х-1У " ( X -1 У >
ЛЬ = - Ь = -
V 0 -п у V 0 0 у
Г X-1 ( X -1 У ) ь ^
. . 0 )'
откуда получаем искомую формулу (9).
ЛЬ = -Х-1 (X -1 У )ь = -( I + С Б )-1 С Б-1Ь.
На практике при решении краевой задачи для системы стержней имеет смысл не обращать матрицу Б, а вычислять её ЬЦ-разложение. Чтобы рассчитать вектор-столбец Б-1Ь в (9), имея ЬЦ-разложение Б = Ьи, нужно решить следующую систему уравнений:
Бг = Ь , или
Ьу = Ь,
(10)
иг = у,
где г = Б-1Ь - искомый вектор-столбец.
Аналогично, чтобы вычислить матрицу (I + С Б), нет необходимости вычислять требуемые п столбцов матрицы Б-1, чтобы составить Б и перемножать с С, а можно, совершив аналогичное количество операций, решить п раз систему уравнений, подобную (10). Так как имеет место равенство
Бт (СБ-1 )Т = Ст, уравнения для определения ((СБ) имеют следующий вид:
ит У = Ст
, где Z = (С Б-1 )т
(11)
Ь Z = У,
откуда вычеркиванием строк получим (СБ) . После решения (11) можно получить матрицу (I + С Б ).
Рассмотрим применение (9) к задаче решения СЛАУ с матрицей (3) в методе Ньютона при решении нелинейной краевой задачи о деформировании системы связанных между собой стержней. В выражении (9) производятся следующие вычисления: по имеющемуся ЬЦ-разложению матрицы Б (имеется разложение для каждого стержня) сначала рассчитывается вектор-столбец х0 = Б-1Ь , затем после вычеркивания строк, которые в С являются нулевыми, вычисляется вектор-столбец Ь0 = -Сх0. Матрица С как
часть матрицы Якоби для системы стержней содержит «связи» между стержнями. Например, пусть три стержня жестко соединены в одной точке, первый стержень сочленяется в г-м узле, второй в у'-м узле, а третий в к-м узле. Тогда уравнения связи будут иметь следующий вид:
1. Краевые условия первого стержня в точке сочленения г: равенство перемещений г-у
и(г} - и'] = 0, и(г' - и2у'' = 0 , и3(г' - и3'' = 0,
равенство поворотов i-j
ш(г) -ш^) = 0, ш(2г) -ш^') = 0, ш3г) ) = 0 .
2. Краевые условия второго стержня в точке сочленения j:
равенство перемещений j-k
Uj') - uf k) = 0, U2j') - U2k) = 0, u3 j') - U3k) = 0,
равенство поворотов j-k
ш|j') - ш|к) = 0, ш2') - ш2к) = 0, ш3') - ш3к) = 0.
3. Краевые условия третьего стержня в точке сочленения к:
уравнения равновесия сил
в!')+в!J)+в!k) = 0, в2г)+q2 j')+e2k)=0, в3г)+q3j')+q3 k) = 0,
моментов
м!г) + м!j) + м( k) = 0, м 2г) + м 2j) + м 2k) = 0, м 3г) + м 3j) + м 3 k) = 0.
Если несколько стержней соединены в одной точке и нет дополнительных ограничений на перемещения этого узла, то краевые условия одного из них должны содержать 3 уравнения равновесия по силам и 3 по моментам. Также допускается, если эти уравнения равновесия будут распределены по нескольким стержням.
Таким образом, имеем по 6 уравнений связи для каждого стержня в точке сочленения. Эти же уравнения, по сути, являются краевыми условиями для каждого отдельного стержня. Дифференцируя уравнения связи по узловым неизвестным, будем получать компоненты матрицы C и некоторые элементы матрицы D .
В общем случае, когда стыковка стержней производится не в глобальных координатах i123, а, например, в подвижных координатах e* =Pjk Xknin, связанных со стержнем, то матрицу C целесообразно строить с помощью численного дифференцирования.
4. Примеры расчета систем стержней
Представленная методика была реализована в программе на языке FORTRAN, рабочее название разрабатываемой программы для расчета нелинейной механики систем стержней - KLPALGSYS. Исходные данные задаются в виде текстовых файлов и имеют определенную универсальную структуру, позволяющую гибко задавать расчетные схемы и управлять процессом расчета; кроме того, эта структура открыта к модернизации. В рамках данной статьи не представляется возможным привести алгоритм работы KLPALGSYS по ряду причин. Во-первых, алгоритм далек от оптимальности и находится на стадии тестирования, отладки и постоянной доработки и развития. Во-вторых, даже его краткое описание займет отдельную статью. Тем не менее ядром этой программы является алгоритм расчета одного стержня KLPALG [14], а основные идеи и принципы описаны выше.
В качестве примера расчета по предлагаемой методике рассмотрим задачу о статическом нелинейном деформировании системы из четырех стержней, жестко соединенных в одной точке 5 под прямыми углами (рис. 2).
Длина каждого стержня 1 метр, поперечное сечение - квадрат со стороной 5 мм. Характеристики материала: модуль упругости Е = 2,1-10пПа, коэффициент Пуассона V = 0,3. В точке 1 стержень жестко защемлен, точка 2 может перемещаться только вдоль оси Охх и свободно поворачиваться, точки 3 и 4 могут свободно перемещаться вдоль оси Ох2 и также свободно поворачиваться. В точке 5 в направлении оси Ох3 (из плоскости
рис. 2) приложен момент 100 Н-м и сила 200 Н.
На рис. 3 представлена деформированная конфигурация системы стержней, рассчитанная по предлагаемому алгоритму. Здесь необходимо отметить, что при расчете можно воспользоваться двумя стратегиями: 1) задаём требуемую точность решения краевой задачи, а программа либо достигает её, либо требует увеличить количество разбиений; 2) задаем дискретизацию, а программа достигает максимально возможной точности решения краевой задачи для такого количества разбиений. В данном примере была выбрана вторая стратегия для того, чтобы можно было сравнивать с решением, полученным по методу конечных элементов.
Рис. 3. Деформированная конфигурация. Расчет по предлагаемой методике
Для сравнения была проведена серия расчетов этой системы стержней в геометрически нелинейной постановке в конечно-элементном пакете АКБУБ (текст исполняемого файла для одной дискретизации модели представлен ниже), на рис. 4 представлена полученная деформированная конфигурация системы стержней. Количество конечных элементов на один стержень задается переменной псе1. Результаты расчета, в том числе и время расчета, записываются в файлы на жесткий диск путем добавления записей.
С помощью кода, приведенного в тексте программы расчета АКБУБ можно сгенерировать и рассчитать раму, состоящую из большего числа стержней, количество клеток задается переменной псе1, при этом значения псе1 должны быть четными. Так, при псе1=2 имеем раму из четырех стержней (см. рис. 2). При псе1=4 будет построена рама из 24 стержней, при псе1=6 будет построена рама из 60 стержней, при псе/=8 - рама из 112 стержней, при псе/=10 рама из 180 стержней и т.д. (количество стержней можно рассчитать как 2*(псе/-1)*псе/). Краевые условия для генерируемых рам соответствуют схеме на рис. 2: слева стержни защемлены, справа подвижные вдоль горизонтального направления шарниры, сверху и снизу подвижные вдоль вертикали шарниры.
Текст программы расчета ANSYS
Finish ksel,s,loc,x,(-bx/2)
/CLEAR dk,all,ux,,,,,uy,uz,rotx,roty,rotz
/UNITS, SI ksel,s,loc,x,(bx/2)
/PREP7 dk,all,uy,,,,,uz
ET,1,beam4 ksel,s,loc,y,(-by/2)
мртемр,,,,,,,, dk,all,ux,,,,,uz
мртемр,1,о ksel,s,loc,y,(by/2)
мРБАТА,ЕХ,1,,2.1Е11 dk,all,ux,,,,,uz
мрвата,веж,1,,7850 allsel
мРОАТА,ОХУ,1,,2.1Е11/(2*(1+0.3)) FK,KP(0,0,0), FZ,200
a=0.005 FK,KP(0,0,0),MZ,100
b=0.005 /SOLU
R,1, a*b, a*(b**3)/12, b*(a**3)/12, a, b, 0 ANTYPE,STATIC
ШОШ, ,a*(b**3)/6, , , , NLGEOM,ON
hx=1 *GET, STime, ACTIVE, 0, TIME, CPU
hy=1 SOLVE
ncel=2 *GET, ETime, ACTIVE, 0, TIME, CPU
nel=50 FINISH
bx=(hx*ncel) /REPLOT
by=(hy*ncel) /POST1
ncel2=(ncel*ncel) /SHOW, WIN32C
*do,i,1,ncel-1 /CONT,1,50
k,i*2-1,(hx*i-bx/2),(-by/2) PLNSOL, u,sum, 2
k,i*2,(hx*i-bx/2),(by-by/2) *dim, mU, array, 3 , 1
l,i*2-1,i*2 n num=NODE(0,0,0)
k,i*2-1+ncel2,(-bx/2),(hy*i-by/2) *GET, mU(1), NODE, n num, UX, 1
k,i*2+ncel2,(bx-bx/2),(hy*i-by/2) *GET, mU(2), NODE, n num, UY, 1
l,i*2-1+ncel2, i*2+ncel2 *GET, mU(3), NODE, n num, UZ, 1
*enddo MUU=sqrt(mU(1)*mU(1)+m U(2)*mU(2)+mU(3)*mU(3))
lovlap,all *cfopen,d:\1 output,prn,,APPEND
lplot *vwrite,ncel,nel,mU(1), mU(2), mU(3), MUU, (ETime-STime)
LATT, 1, 1, 1 %i %i %G %G %G %G %G
lesize,all,,,nel,, *cfclos
lmesh,all
Ниже будут приведены результаты расчетов рам для nce/=2,4,6,8,10 для разных дискретизаций nel.
Рис.4. Деформированная конфигурация. Расчет в А№У8
На рис. 5 представлены графики зависимости времени расчета на логарифмической шкале в зависимости от дискретизации моделей для конечно-элементной и дифференциальной модели соответственно. Из рисунка видно, что время расчета КЭ-модели в АКБУБ до двух порядков (при увеличении количества КЭ) больше времени расчета по предлагаемой дифференциальной модели. Такое преимущество весьма условно, время расчета существенно зависит от количества стержней в системе. С ростом числа стержней это преимущество уменьшается и КЭ-расчет в АКБУБ становится быстрее расчета по дифференциальной модели, далее это будет показано. На рис. 6 представлен график разницы в процентах перемещения точки 5 (центральная точка), рассчитанного в АКБУБ и по предлагаемой методике. На рис. 7 представлены графики полных перемещений центральной точки, рассчитанных по МКЭ в АКБУБ и по предлагаемой дифференциальной модели.
Рис. 5. Время расчета в зависимости от дискретизации моделей (4 стержня, псе1=2)
Рис. 6. Разница между рассчитанным максимальным перемещением в и по предлагаемому алгоритму (4 стержня, псе1=2)
Рис. 7. Полное перемещение центральной точки
Рассмотрим расчеты рам с другим количеством стержней. Были проведены расчеты рам для ncel=4, 6, 8, 10 в ANS YS и KLPALGSYS для разных дискретизаций. На рис. 8 приведен пример деформированной конфигурации для системы 180 стержней, полученный в ANSYS.
На рис. 9-12 приведены результаты расчетов для ncel=4,6,8,10 соответственно. На каждом рисунке представлены три вида зависимостей от величины дискретизации модели: время расчета по двум методам, разница полных перемещений в центральной точке и полные перемещения центральной точки по двум методам. По рис. 9-12 можно проследить, как меняется время расчета по методам с ростом числа стержней в системе. Предлагаемая дифференциальная модель имеет безусловное преимущество в скорости перед КЭ-расчетом в ANSYS примерно до 100 стержней в системе. При 112 и 180 стержнях в системе наблюдается сначала преимущество в скорости ANSYS, но с ростом количества КЭ предлагаемая дифференциальная модель вновь опережает ANSYS.
Рис. 8. Расчет рамы из 180 стержней в ANSYS (ncel=10)
50 100 150 200
Количество КЭ, количество разбиений на один стержень
Рис. 9. Результаты расчетов для системы 24 стержней (ncel = 4)
50 100 150 200 0
Количество КЭ, количество разбиений на один стержень
Рис. 10. Результаты расчетов для системы 60 стержней (псе/=6)
Рис. 11. Результаты расчетов для системы 112 стержней (псе/=8)
Количество КЭ, количество разбиений на один стержень Рис. 12. Результаты расчетов для системы 180 стержней (псе/=10)
Продолжая анализ полученных результатов, необходимо отметить, что получаемые перемещения по двум методикам сходятся с увеличением дискретизации, причем отчетливо заметно, что для получения приемлемой погрешности определения перемещений требуется до нескольких десятков КЭ на одном стержне в системе. При этом время расчета по предлагаемой дифференциальной модели медленнее растет с увеличением дискретизации, чем при расчете по КЭ-модели ЛКБУБ.
Заключение
В представленной статье описана методика расчета механики систем связанных стержней по дифференциальной модели. Для использования хорошо отработанного аппарата решения нелинейной краевой задачи одного стержня, для случая системы связанных между собой стержней разработан метод решения СЛАУ. Такой подход позволяет получать решение системы с большой матрицей из решения нескольких СЛАУ
меньшего размера. Методика расчета реализована в программе на языке FORTRAN, в статье приведены тестовые расчеты нелинейной статики систем стержней. Результаты расчета сравниваются с расчетами методом конечных элементов в пакете ANSYS. Результирующие перемещения двумя методами стремятся друг к другу с ростом дискретизации моделей. По результатам расчетов можно сделать вывод, что предлагаемый дифференциальный подход эффективен и имеет преимущество в скорости расчета перед ANSYS при относительно небольшом количестве стержней в системе (условно до 100 стержней). Это объясняется тем, что при реализации предложенного алгоритма решения СЛАУ из условий стыковки стержней между собой формируется матрица, размер которой растет с увеличением числа стыкуемых стержней, и затраты процессорного времени на решение СЛАУ с этой матрицей быстро растут. Так, например, при расчете стержневой системы из ста стержней в общем случае будем иметь матрицу связей размером 1200*1200 элементов. Затраты на решение СЛАУ с такой матрицей становятся сравнимы с затратами на все остальные вычисления, даже притом, что размер исходной матрицы на порядок, а то и на два превышает размеры матрицы связей. Решение данной проблемы видится в применении других численных методов для решения СЛАУ-связей, либо в ней можно выделить блоки и применить тот же метод решения по блокам. Тем не менее предложенный алгоритм решения хорошо подходит для решения достаточно широкого круга задач механики стержневых систем, в которых не требуется рассчитывать одновременно большое количество стержней.
Библиографический список
1. Жиделёв А. В. Динамический подход к расчету геометрически нелинейных стержневых систем : автореф. дис. ... канд. техн. наук. - Волгоград, 2003. - 23 с.
2. Галишникова В. В. Обобщенная геометрически нелинейная теория и численный анализ деформирования и устойчивости пространственных стержневых систем: автореф. дис. ... д-ра техн. наук. - Волгоград, 2011. - 40 с.
3. Raftoyiannis I.G. Nonlinear Stability Analysis of Simple Steel Frames with Load Imperfections // Frontiers of Engineering Mechanics Research. - Feb. 2013. - Vol. 2. - Iss. 1. - P. 15-21.
4. Строительная механика. Стержневые системы: учебник для вузов / А.Ф. Смирнов [и др.] / под ред. А.Ф. Смирнова. - М.: Стройиздат, 1981. - 512 с.
5. Дарков А.В., Шапошников Н.Н., Строительная механика: учебник для строительных спец. вузов - М., 1986. - 607 с.
6. Дьяков И.Ф., Чернов С.А., Черный А.Н. Метод конечных элементов в расчетах стержневых систем: учеб. пособие / Ульян. гос. техн. ун-т. - Ульяновск, 2010. - 133 с.
7. A discrete mechanics approach to Cosserat rod theory. Part 1: Static equilibria / P. Jung, S. Leyendecke, J. Linn, M. Ortiz // Int. J. Numer. Meth. Engng. - 2010. - P. 101-130.
8. Радченко В.П. Введение в механику деформируемых систем: учеб. пособие. - Самара: Са-мар. гос. техн. ун-т, 2009. - 196 с.
9. Petrolito J., Legge K.A. General nonlinear analysis of frames. Proceedings of the 5th Australasian Congress on Applied Mechanics (ACAM 2007). - Brisbane, Australia. 10-12 December, 2007. - Р. 82-88.
10. Чернецов Д. А. Процесс нелинейного деформирования стержневых систем без ограничения степени деформации // Изв. Тул. гос. ун-та. Естественные науки. - 2012. - Вып. 2. - С. 165-176.
11. Levyakov S. V. Formulation of a geometrically nonlinear 3D beam finite element based on kinematic-group approach // Applied Mathematical Modelling. - 2015. - Vol. 39. - Iss. 20. - P. 6207-6222.
12. Левяков С. В. Нелинейный пространственный изгиб криволинейных стержней с учетом поперечного сдвига // ПМТФ. - 2012. - Т. 53, № 2. - С. 128-136.
13. Левин В.Е., Пустовой Н.В. Механика деформирования криволинейных стержней: мо-ногр. - Новосибирск: Изд-во Новосиб. гос. техн. ун-та, 2008. - 208 с.
14. Пустовой Н.В., Левин В.Е., Красноруцкий Д.А. Алгоритм численного решения нелинейной краевой задачи динамического деформирования тонкого стержня // Вестник Пермского национального исследовательского политехнического университета. Механика. - 2014. - № 2. -С.168-199.
15. Жилин П. А. Прикладная механика. Теория тонких упругих стержней: учеб. пособие. -СПб.: Изд-во Политехн. ун-та, 2007. - 101 с.
16. Светлицкий В.А. Строительная механика машин. Механика стержней: в 2 т. Т. 2. Динамика. - М.: Физматлит, 2009. - 383 с.
17. Сорокин Ф.Д. Прямое тензорное представление уравнений больших перемещений гибкого стержня с использованием вектора конечного поворота // Изв. РАН. МТТ. - 1994. - № 1. -С.164-168.
18. Nonlinear dynamic deformation simulation for helical rod like objects / H. Du, W. Xiong, H. Wang, Z. Wang, B. Yuan // Engineering Review. - 2013. - Vol. 33. - Iss. 3. - P. 233-238.
19. Pereyra V. Pasva3: An adaptive finite difference fortran program for first order nonlinear, ordinary boundary problems // Lecture Notes in Computer Science. - 1979. - Vol. 76. - P. 67-88.
20. Keller H.B. Accurate Difference Method For Nonlinear Two-Point Boundary Value Problems // SIAM J. Numer. Anal. April 1974. - Vol. 11. - No. 2. - P. 305-320.
References
1. Zhidelev A.V. Dinamicheskii podkhod k raschetu geometricheski nelineinykh sterzhnevykh sistem [A dynamic approach to the calculation of geometrically nonlinear rod systems]. Abstract of Ph.D. thesis. Volgograd, 2003. 23 р.
2. Galishnikova V.V. Obobshchennaia geometricheski nelineinaia teoriia i chislennyi analiz deformirovaniia i ustoichivosti prostranstvennykh sterzhnevykh sistem [Generalized geometrically nonlinear theory and numerical analysis of deformation and stability of spatial rod systems]. Abstract Doctors degree thesis. Volgograd, 2011. 40 p.
3. Raftoyiannis I.G. Nonlinear Stability Analysis of Simple Steel Frames with Load Imperfections. Frontiers of Engineering Mechanics Research, Feb. 2013, vol. 2, iss. 1, pp. 15-21.
4. Smirnov A.F. et al. Stroitel'naia mekhanika. Sterzhnevye sistemy [Structural Mechanics. Core systems]. Moscow: Stroiizdat, 1981. 512 p.
5. Darkov A.V., Shaposhnikov N.N. Stroitel'naia mekhanika [A Textbook for construction]. Moscow, 1986.607 p.
6. D'iakov I.F., Chernov S.A., Chernyi A.N. Metod konechnykh elementov v raschetakh sterzhnevykh sistem [The finite element method in the calculation of rod systems: a tutorial]. Ul'ianovskii gosudarstvennyi tekhnicheskii universitet, 2010. 133 p.
7. Jung P., Leyendecke S., Linn J., Ortiz M. A discrete mechanics approach to Cosserat rod theory Part 1: Static equilibria. Int. J. Numer. Meth. Engng, 2010, pp. 101-130.
8. Radchenko V.P. Vvedenie v mekhaniku deformiruemykh sistem [Introduction to the mechanics of deformable systems]. Samarskii gosudarstvennyi tekhnicheskii universitet, 2009. 196 p.
9. Petrolito J., Legg K.A. General nonlinear analysis of frames. Proceedings of the 5th Australasian Congress on Applied Mechanics (ACAM 2007). Brisbane, Australia, 10-12 December, 2007, pp. 82-88.
10. Chernetsov D.A. Protsess nelineinogo deformirovaniia sterzhnevykh sistem bez ogranicheniia stepeni deformatsii [The process of nonlinear deformation of rod systems are not limited to the degree of deformation]. Izvestiia Tul'skogo gosudarstvennogo universiteta. Estestvennye nauki, 2012, iss. 2, pp.165-176.
11. Levyakov S.V. Formulation of a geometrically nonlinear 3D beam finite element based on kinematic-group approach. Applied Mathematical Modelling, 2015, vol. 39, iss. 20, pp. 6207-6222.
12. Levyakov S.V. Nonlinear spatial bending of shear-deformable curvilinear rods. Journal of Applied Mechanics and Technical Physics. March 2012, vol. 53, iss. 2, pp 258-265.
13. Levin V.E., Pustovoi N.V. Mekhanika deformirovaniia krivolineinykh sterzhnei [The mechanics of deformation of curved rods]. Novosibirskii gosudarstvennyi tekhnicheskii universitet, 2008. 208 p.
14. Pustovoy N.V., Levin V.E., Krasnorutskiy D.A. The numerical algorithm for solving nonlinear boundary problem of thin rod's dynamic deformations. PNRPU Mechanics Bulletin. 2014, no. 2, pp.168-199.
15. Zhilin P.A. Prikladnaia mekhanika. Teoriia tonkikh uprugikh sterzhnei [Applied mechanics. Theory of thin elastic rods]. Sankt Peterburgskii politekhnicheskii universitet, 2007. 101 p.
16. Svetlitskii V.A. Stroitel'naia mekhanika mashin. Mekhanika sterzhnei. Vol. 2. Dinamika [Construction machinery mechanic. Mechanics rods. Vol. 2. Dynamics]. Moskow: Fizmatlit, 2009. 383 p.
17. Sorokin F.D. Priamoe tenzornoe predstavlenie uravnenii bol'shikh peremeshchenii gibkogo sterzhnia s ispol'zovaniem vektora konechnogo povorota. Izvestiia Rosiiiskoi akademii nauk. MTT, 1994, no. 1, pp. 164-168.
18. Du H., Xiong W., Wang H., Wang Z., Yuan B. Nonlinear dynamic deformation simulation for helical rod like objects. Engineering Review, 2013, vol. 33, iss. 3, pp. 233-238.
19. Pereyra V. Pasva3: An adaptive finite difference fortran program for first order nonlinear, ordinary boundary problems. Lecture Notes in Computer Science, 1979, vol. 76, pp. 67-88.
20. Keller H.B. Accurate Difference Method For Nonlinear Two-Point Boundary Value Problems. SIAM J. Numer. Anal, April 1974, vol. 11, no. 2, pp. 305-320.
