МОДИФИКАЦИЯ МЕТОДА ПОСТОЯННЫХ ДУГ, ОСНОВАННАЯ НА ИСПОЛЬЗОВАНИИ МАТРИЦЫ СЕКУЩЕЙ
ЖЕСТКОСТИ
В.В. Галишникова
ВолгГАСУ
В настоящее время все более широкое применение в вычислительной механике строительных конструкций получают методы продолжения по параметру, в которых в качестве ведущего параметра используется длина дуги кривой равновесных состояний. В зарубежной литературе эти методы известны под общим названием "Arc Length methods", а в русскоязычных источниках приводятся под различными наименованиями, например "метод продолжения по длине дуги кривой равновесных состояний", "метод дуговых засечек", или "метод длин дуг".
В задачах дискретной нелинейной строительной механики метод продолжения по параметру, в котором в качестве параметра используется длина дуги кривой равновесных состояний, впервые был использован Риксом [1] и Вемпнером [2]. Последующие модификации этого метода были выполнены Крисфилдом [3,4] и рядом других авторов. Наиболее эффективной версией метода длин дуг на настоящее время считается модификация Крисфилда, в которой длина приращения дуги кривой равновесных состояний сохраняется постоянной на всех шагах решения. Для жесткой конструкции приращение дуги в основном зависит от приращения нагрузки.
В подходе Крисфилда при расчете по деформациям используется традиционная совместная линеаризация инкрементальных разрешающих уравнений на шаге нагру-жения. В конце каждого шага вычисляются неуравновешенные силы, возникающие вследствие линейной аппроксимации матрицы жесткости, и используются в виде поправочного вектора сил для итерационного уточнения решения. Эти алгоритмы хорошо работают при выполнении деформационного анализа и позволяют так называемый «обход» критических точек, если эти точки являются предельными. Однако, при необходимости точного вычисления критических точек, в особенности точек бифуркации, а также при анализе закритического поведения конструкции, эти алгоритмы не обеспечивают требуемой точности вычислений.
В представленном исследовании для решения нелинейной задачи деформирования пространственных ферм по методу конечных элементов разработана модификация метода постоянных дуг, в которой используется процедура линеаризации секущей матрицы жесткости конструкции.
Главные особенности предлагаемого подхода состоят в следующем:
• в разрешающих уравнениях метода конечных элементов используется матрица секущей жесткости фермы;
• матрица секущей жесткости фермы формулируется на основе разностных операций;
• в выражениях для коэффициентов матрицы секущей жесткости сохраняются все нелинейные члены исходных разрешающих уравнений;
• матрица секущей жесткости вычисляется итерационно;
• линеаризация уравнений достигается использованием значений переменных, полученных на предыдущем шаге итераций;
• вектор неуравновешенных сил не добавляется к внешней нагрузке, действующей на конструкцию, а используется для вычисления корректирующих перемещений и корректирующих реакций, которые затем вводятся в выражения коэффициентов матрицы секущей жесткости.
В качестве параметра продолжения решения использована постоянная длина дуги кривой равновесных состояний. Размер шага определяется в соответствии с методом постоянных дуг. Шаговый алгоритм решения основных уравнений задачи нелинейного деформирования конструкций описан в п. 4 [5].
Численные исследования точности, устойчивости и надежности метода выполнены при помощи разработанного программного приложения на платформе Java. Полученные результаты сравниваются с результатами точного аналитического решения тестовой задачи, выполненного автором [6].
Equation Section 11. Формирование матрицы секущей жесткости фермы При пошаговом решении задач дискретной нелинейной механики конструкций возникает необходимость дифференцирования понятия «жесткость конструкции». В настоящей работе для описания жесткости используются три термина-определения -полная жесткость, касательная жесткость и секущая жесткость. Для успешной реализации теории в алгоритмах и программах необходим последовательный учет этих определений. На рис.1 наглядно представлены принципы формирования матриц полной, касательной и секущей жесткости.
У/ 'в
//
/ /
//
//
->
A uB
полная жесткость
K[Ub] UB = PB
Pd Pc
"b c
касательная жесткость K,[UB] aU = др
ди = UC - UB
aP = PC - PB
ж
с
/ / / в ->
u
Ur
-"B "C
секущая жесткость К8[ив,ди] ди = др
ди = Uc - UB
ЛР = PC - PB
Р
Р
Р
с
Р
в
в
U
U
Рис. 1. Полная, касательная и секущая жесткость конструкции
Разрешающие уравнения геометрически нелинейного поведения конструкции задают криволинейную траекторию нагружения АВ. Коэффициенты зависимости между перемещениями и нагрузками в этих уравнениях формируют матрицу полной жесткости конструкции. В шаговых методах геометрически нелинейных расчетов используется линеаризация основных уравнений на шаге нагружения путем дифференцирования по направлениям перемещений. В этом случае траектория на шаге нагружения представляет собой касательную ВБ к действительной траектории в начале шага на-
гружения. Коэффициенты зависимости между перемещениями и нагрузками в линеаризованных уравнениях формируют матрицу касательной жесткости конструкции. Линеаризация разрешающих уравнений на основе разностного подхода дает траекторию на шаге в виде секущей СЕ действительной траектории нагружения. Зависимость между перемещениями и нагрузками в этом случае описывается матрицей секущей жесткости.
Матрица секущей жесткости пространственной фермы формулировалась следующим образом.
1.Записываются алгебраические разрешающие уравнения задачи геометрически нелинейного деформирования конструкции для двух последовательных состояний конструкции 5 — 1 И 5 .
2.Определяются разности соответствующих коэффициентов матриц полной жесткости К5 и К51 уравнений состояний 5 и 5 — 1.
3.Формируются коэффициенты матрицы секущей жесткости К $ для данного шага нагружения. При этом производится перегруппировка членов в разностях коэффициентов К5 и К5 1 таким образом, чтобы результирующая матрица секущей жесткости была симметрична.
Матрица жесткости элемента пространственной фермы имеет следующий окончательный вид:
ах а4 а5 -а1 -а 4 ~а5
а4 а2 а6 -а 4 ~а 2 ~а6
АЕ а5 а6 а3 ~а5 ~а6 ~а3
1 -а1 -а4 ~а5 ах а4 а5
-а 4 ~а 2 ~а6 а4 а2 а6
~а5 ~а6 ~а3 а5 а6 а3
— 1 + С1 + С7 + + С14 + С20 , а2 — С2 ^ С7 + С9 + С15 + С20 ,
а3 - С3 + С7 + С10 + С16 + С20 , а4 - С4 + С11 + С17 ,
а5 = С5 + С12 + С18 , а6 = С6 + С13 + С19 , (1.2)
С1 = ^1,1 (2 + ^1,1) , С2 = ^2,1 , С3 = , С4 = ^2,1(1 + ^1,1) ,
С5 = ^3,1(1 + ^1,1) , С6 = ^2,1 ^3,1 , С7 = ^1,1 + 0,50\1 + ^2,1 + ) ,
С8 = д^1,1 (1 + ^1,1) , С9 = д^2,1^2Д , С10 = д^3Д^3,1, Сп = 0,5 (А^2д(1 + УЦ) + ДУцУу) ,
С12 = 0,5(д^3,1(1 + ^1,1) + д^3,1) ,
С13 = 0,5 (д^2,1 ^3,1 +д^3,1 ^2,1) , С14 = 0,5
-2 1,
С15 = 0,5дУ221 , С16 = 0,5 ДУ321, С17 = 0,5 ДУПДУ21,
3
c18 = 0,5 ¿VyAVy, c19 = 0,5 AV2,1AV3,1, c20 = 0,5(д^1,1 +ZVi,1 AVi,1) (L3)
i=1
В выражениях коэффициентов матрицы использованы следующие обозначения: Vi 1 - производные координат перемещений по направлению оси, совпадающей с осью стержня в мгновенной конфигурации, av¡ 1 - инкременты производных. Верхней чертой отмечены члены, аппроксимируемые значениями, взятыми из предыдущей итерации.
Нетрудно убедиться, что матрица жесткости элемента является симметричной. Матрица жесткости конструкции компонуется из матриц составляющих ее элементов по известным правилам. В результате получается следующая система алгебраических разрешающих уравнений относительно инкрементов перемещений и реакций пространственной фермы на шаге нагружения:
Ks s = es + s + *rs . (L4)
Разрешающие уравнения (1.4) являются линеаризованными.
Equation Section 22. Вектор невязки разрешающих уравнений
Инкрементальное разрешающее уравнение (1.4) включает вектор погрешности es , содержащий неуравновешенные силы состояния s фермы. В мгновенной конфигурации конструкции, для которой определен вектор es , этот вектор уравновешен. Составляющие его силы и моменты вызывают перемещения, реакции и внутренние усилия в ферме. Покажем, что неуравновешенные силы не дают вклада в нагрузку, действующую на ферму.
Рассмотрим точные разрешающие уравнения относительно полных перемещений и реакций фермы в состоянии s :
Ks) des) = )+к, где (2.1)
K^es) - точная матрица полной жесткости для состояния s ; As)
de - точный вектор полных перемещении в состоянии s ;
( s)
qe - точный вектор полной нагрузки в состоянии s ;
re - точный вектор полных реакций в состоянии s .
Для того чтобы получить точную матрицу секущей жесткости для шага s , необходимо записать точное разрешающее уравнение для предыдущего шага нагружения фермы s -1:
K^ des-1) = qes-1) + r(s-1). (2.2)
Вычтя уравнение (2.2) из уравнения (2.1), получим
(кes) -Kes-1))d(s-1) + K(s)Ad(s) = Aqes) + Ar(s), где (2.3)
./(s) j(s) j(s-1)
Ad e — d e — d e - точный вектор инкрементов перемещении на шаге s ;
(s) (s) (s-1)
АЧе ~ qe ~ qe - точный вектор инкрементов нагрузки на шаге s ;
(i) (i) (s -1)
аг^' — re ' — re ' - точный вектор инкрементов реакции на шаге s .
Так как разность матриц полной жесткости в состояниях s и s — 1 представляет собой функцию инкрементов перемещений на шаге s , то первый член в уравнении (2.1) может быть записан в следующей форме:
[k{S - к(: -1) yt1 = д к es) ad(s), (2.4)
tr( s)
где дK e - вспомогательная матрица.
Подстановка выражения (2.4) в уравнение (2.3) дает точное разрешающее уравнение относительно инкрементов перемещений на шаге s , которое содержит точную матрицу секущей жесткости:
KAd(s) = Aqis) + аг(s), где (2.5)
K W = Kes) + aK^1 ) - точная матрица секущей жесткости.
Точное разрешающее уравнение (2.5) аппроксимируется разрешающим уравнением, в котором вектор инкрементов нагрузки на шаге является точным, а матрица секущей жесткости и векторы инкрементов перемещений и реакций - приближенными:
K is) Ad(s) = Aq{S + дг(s). (2.6)
Уравнение (2.6) решается относительно инкрементов перемещений Ad(s) и ин-(s)
крементов реакций д^ ' на шаге нагружения s . Так как уравнение (2.6) является лишь приближенным, то состояние s фермы не является равновесным. Неуравновешенные силы равны разности внешних и внутренних сил, действующих в узлах фермы:
е(s) = qis) + г(s) - K) d(s), где (2.7)
(s)
е - вектор погрешности, содержащий неуравновешенные силы в состоянии s .
Предположим, что неуравновешенные силы возникают вследствие использования
в уравнении (2.5) приближенной матрицы секущей жесткости Kis) вместо точной матрицы:
к is) Ad es) = Aqis) + дге( s) + е(s). (2.8)
Подставим вектор инкрементов нагрузки, полученный из уравнения (2.6), в уравнение (2.8):
кis)( Ad es) - Ad(s)) = е(s) + (дг]1) - дг(s)). (2.9)
Уравнение (2.9) решается относительно корректирующих перемещений Ad es) — Ad(s) и корректирующих реакций a^s) — аг (s) в состоянии s . При этом нагрузка в состоянии s не корректируется.
3. Точность и эффективность метода
Исследования точности, устойчивости и надежности метода выполнены при помощи разработанного программного приложения на объектно-ориентированной платформе Java. В качестве тестового примера рассматривалась задача нелинейного деформирования симметричной пространственной трехстержневой фермы, для которой
в работе [6] получено точное аналитическое решение. Шаговый алгоритм численного решения основных уравнений задачи описан в п. 4 [1].
Расчетная модель симметричной пологой пространственной трехстержневой фермы показана на рис. 3.1. Опоры в узлах A, B, и С представляют собой шаровые шарниры. Стержни фермы воспринимают только осевые усилия. Площадь поперечного сечения всех стержней равна 10 кв.см. Модуль упругости материала принят равным
2-105 МПа. В верхнем шарнирном узле Б ферма загружена вертикальной силой Р = -10,0 кН.
Нагрузка прикладывается в 65 ступеней, во время которых ферма проходит полный цикл загружения-разгрузки. Расчет полностью охватывает нелинейную диаграмму работы конструкции. Примеры результатов пошагового расчета приведены в табл. 3.1 с интервалом в 10 шагов. Погрешности в определении перемещений и реакций очень малы и не зависят от величин перемещений.
Таблица 3.1. Погрешность вычисления траектории нагружения фермы
№ Перемещение Нагрузка Точное значение Погрешность
шага М3 , м Рэ, кН Р3, кН е,%
1 -0,003082264 -5,000000000 -4,999854298 0,002914125
5 -0,016185170 -22,858054820 -22,857765439 0,001266007
15 -0,060788321 -48,295855224 -48,239346556 0,117142275
25 -0,111247574 -30,26025916 -30,25931986 0,003104151
35 -0,22123336 46,1693256 46,15869695 0,023026335
45 -0,261347109 41,82901148 41,82691188 0,005019737
55 -0,294715769 8,380986459 8,380659194 0,003905003
65 -0,32230765 -46,06856653 -46,06872352 -0,000340782
Выполненные численные эксперименты показали, что количество итераций на шаге нагружения зависит от принятой предельно допустимой погрешности и началь-
ного коэффициента нагружения, принятых в расчете. Общее количество итераций уменьшается с возрастанием коэффициента погрешности и начального коэффициента нагружения. Среднее количество итераций на шаге не превышает пяти, что свидетельствует о быстрой сходимости процесса. При этом точность численного расчета по деформациям пологой фермы зависит от коэффициента погрешности и мало зависит от значения начального коэффициента нагружения.
4. Заключение
На основе инкрементальных разрешающих уравнений геометрически нелинейного деформирования пространственной фермы получено выражение матрицы секущей жесткости системы. Разработана итерационная процедура учета неуравновешенных сил непосредственно в коэффициентах матрицы секущей жесткости. Сформулированы уравнения метода конечных элементов для задачи, использующие разработанный подход.
На примере задачи деформирования симметричной пространственной трех-стержневой фермы выполнено сравнение результатов численного расчета с результатами аналитического решения, полученного автором. Продемонстрирована высокая точность и хорошая сходимость метода.
Литература
1. Ricks, E. The application of Newton's method to the problem of elastic stability// Trans. ASME. -1972. - E39, N4.-P. 1060-1065. (Рус. пер: Рикс. Применение метода Ньютона к задаче упругой устойчивости // Прикл. механика. - 1972. - №4.-С 204-210).
2. Wempner, G.A., Discrete approximations related to nonlinear theories of solids, Int. J. Solids & Structs., 7, 1971. pp.1581-1599.
3. Crisfield, M.A., Incremental/iterative solution procedures for non-linear structural analysis// Num. Meth. for Non-linear Problems, Vol. 1, ed. C. Taylor et al, Pineridge Press, Swansea, pp.261-290, (1980)
4. Crisfield, M.A., A fast incremental/iterative solution procedure that handles "snap-through"// Computers & Structures, 13, pp.55-62, (1981).
5. Галишникова В.В. Унифицированный и общий подход к геометрически нелинейному расчету строительных конструкций // Вестник Волгогр. гос. архит.-строит. ун-та. Сер.: Техн. науки. 2006. Вып. 6(20). С. 42—66.
6. Галишникова В.В. Аналитическое решение нелинейной задачи устойчивости и исследование закритического поведения трехстержневой фермы // Вестник Волгогр. гос. архит.-строит. ун-та. Сер.: Естеств. науки. 2007. Вып. 6(23). С. 53—64.
Рецензент: д.т.н., проф. И.Г. Овчинников, Саратовский ГТУ