Некоторые вопросы автоматизации проектирования тонкостенных конструкций Текст научной статьи по специальности «Физика»

НЕКОТОРЫЕ ВОПРОСЫ АВТОМАТИЗАЦИИ ПРОЕКТИРОВАНИЯ ТОНКОСТЕННЫХ КОНСТРУКЦИЙ

Трушин С.И. (МГСУ)

При исследовании устойчивости тонкостенных конструкций в нелинейной постановке возникает необходимость определения предельных и бифуркационных точек на кривой равновесных состояний и исследования поведения конструкций в закритической стадии. Для построения кривых равновесных состояний и исследования устойчивости форм равновесия оболочечных конструкций весьма эффективным является класс методов, основная идея которых сводится к нахождению последовательности решений ит=и(рт) (т=0,1,2,...,М) системы нелинейных уравнений вида:

Р(и,р)= 0 (1)

на основе имеющегося начального решения и0=и(р0) при изменении ведущего параметра р от р0 до рт. Эта идея используется в математике давно. Например, метод возмущений или метод малого параметра, в основе которого лежит аналогичный подход, применялся еще в работах У. Леверье и А. Пуанкаре. Для численной процедуры решения нелинейного уравнения данная методика впервые нашло свое место в работе М. Лазя [10]. В дальнейшем различные варианты этого подхода рассматривались И.И. Воровичем и В.Ф. Зипаловой [1], Э.И. Григолюком и В.И. Шалашилиным [2], Д.Ф.Давиденко [3], В.В. Петровым [6], В.И. Феодосьевым [7], Г. Вемпнером [13], М. Крисфилдом [8], Э. Риксом [12] и другими.

Анализ решений нелинейных задач теории пластин и оболочек показывает, что наилучшие результаты дают схемы метода продолжения, основанные на процедурах Ньютона-Рафсона или Рунге-Кутта, использующие в качестве ведущего параметра длину дуги кривой равновесных состояний и включающие в себя вспомогательное уравнение для итераций на сфере. Использование длины дуги как параметра продолжения обеспечивает единый процесс прохождения регулярных, предельных и бифуркационных точек. Отпадает необходимость смены ведущего параметра (например, нагрузки на характерное перемещение) при прохождении предельных точек, да и само понятие предельной точки в такой постановке теряет смысл [2].

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

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

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

Т Т 3 3 и 1 и 0 и 2 и 0

П 1 П 1 т V 1 т И _ И И _ и N

= ' 2 [ / ( Ъат 1« ..... Ъат 6™ , ---, ---) +

2 т = 0 т = 0 и 1 и 2 (2)

3 3 и 3 ц 2 ц 3 ц 1

ЛТ^ 2 ^ V 2 т и и и и м

Ъат 1« >."> Ъат 6™ ,-и-Т-)].

т = 0 т = 0 и 1 и 2

где к\ , Т2 - линейные размеры ячейки в направлении координатных осей; / -подинтегральная функция; ит - векторы искомых функций перемещений в узлах ячейки; ит" - параметры схемы. В предлагаемой схеме (2) порядок аппроксимации производных повышается до 0(к2).

В результате выполненной дискретизации исходная вариационная задача сводится к нахождению экстремума некоторой скалярной функции векторного аргумента с параметром внешней нагрузки р. Условие экстремума дает систему нелинейных уравнений типа (1), которая в развернутом виде записывается следующим образом:

V Ж (и )- ро = 0, (3)

где V Ж (и ) - градиент потенциальной энергии деформации; О - нормированный вектор узловых нагрузок.

Система (3) задает в неявном виде кривую равновесных состояний в пространстве X = (и, р У , зависящую от некоторого параметра продолжения, роль которого может выполнять параметр нагрузки р, любая из компонент вектора перемещений и или длина дуги кривой равновесных состояний Ал. Тогда дополнительно к п уравнениям (3) вводится п+ 1-е уравнение вида Р0 (И, р, л )= 0 . Решение основной системы и вспомогательного уравнения выполняется итерационными методами и методами дифференцирования по параметру.

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

V 2 Ж (и " ( )) Д и "+1 = О Дрк+1 + О рк () - V Ж (и " ()) (4) И " +1 = и " + Т " Д и "+1

||*и"+1||2 + (3р"+1 )2 =Ал2 (5)

или

(а 0, ¿и" +1 )+ап +1, Зр"+1 = А л , (6)

где

ди"+1 = ¿и" + Ди"1, ф"+1 = 8р" + Др"+1, а0 =(а1а2...ап)г. Формулы (4) и (5) описывают итерационный процесс нахождения решения на сфере с центром в точке (и(лш-1), р(лш-1)) и радиусом Ал, а соотношения (4) и (6) определяет итерации на плоскости, для которой вектор а является вектором нормали и

ВЕСТНИК

МГСУ

которая расположена на расстоянии Дл от точки (и(5т-1), р(ат-1)).

В методе дифференцирования по параметру система нелинейных уравнений (3) сводится к задаче Коши для системы обыкновенных дифференциальных уравнений в самокорректирующейся форме в совокупности со вспомогательными уравнениями: d и „ dp

Ж (и ( 5 )) — = Q-f- + Z [0 р ( 5 ) ds ds

V 2

и(5о) = и

V Ж (и (5))] (7)

р ( 5 о) = р'

d и 2 Г л

ds [ ^ ds у

= 1

(8)

а

d и

+ а

dp

-= 1

(9)

ds ) " ds

Для решения задачи (7),(8) и (7),(9) разработаны и апробированы численные процедуры на базе метода Эйлера, метода Рунге-Кутта и самокорректирующейся схемы метода Рунге-Кутта [4].

Достаточно универсальный алгоритм вычисления компонент вектора градиента и матрицы Гессе дискретного аналога потенциальной энергии деформации Ж(и) может быть построен на основе приближенных разностных формул, которые были предложены А.Б. Золотовым и использовались при решении ряда задач в работах [5,11]. Эти формулы позволяют организовать гибкий вычислительный процесс, легко перестраиваемый в зависимости от конкретного вида вариационного функционала. Однако, они приводят к многократному вычислению функции Ж(и) при определении каждого коэффициента, что требует значительных затрат машинного времени. Более экономичным и более точным является подход, основанный на вычислении первых и вторых производных Ж(и) по своим явным выражениям. Так например, для упруго-пластических задач в рамках теории течения соответствующие компоненты определяются по формулам [11]:

дЖ(и) _ ггг „ V г» " , ^ 7 1 ''е

0 ' д и /

0и}

0 2 Ж (и)

ди, д и,

= Ш(((«-« о Г О о 1

dV

(10)

о •" •)' О о" - о -)

де.

ди

д2 6

ди, ди,

dV

где е и а - соответственно векторы деформаций и напряжений; Эо -упругопластическая матрица, связывающая приращения напряжений с приращениями деформаций. Из соотношений (10), как частные случаи, следуют зависимости для решения задач теории пластин и оболочек с различной комбинацией нелинейностей. В формулах (10) первая и вторая производные вектора е по узловым перемещениям определяются путем подстановки единичных векторов е, и е,- в линейную и нелинейную части приращения деформаций:

— = а 8 ' (е ) ;

д иг (11)

д 2 е д и д и

= а е

( е,

) - а 8 1,1 (е, ) - а 8 1,1 (е , )

2

Представленные подходы позволяют выполнять расчеты пластин и оболочек в рамках деформационной теории и теории пластического течения с построением кривых равновесных состояний, определением предельных точек и точек ветвления решений, зон пластических деформаций [9].

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

( 2 а )2

к = к = --— = 24 ,

( " Як

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

'"(2г)4 ^

Рис.1. Кривая равновесных состояний пологой сферической оболочки на прямоугольном плане

Анализ последовательности равновесных состояний выполнен с помощью процедур дискретного и непрерывного продолжения (4), (5) и (7), (8).

В общем случае, когда на конструкцию действует система внешних нагрузок, каждая из которых изменяется по своему закону, система нелинейных уравнений, определяющая условие стационарности энергетического функционала, представляется следующим образом:

N

V Г (и ) = X Р.Ъ., (12)

I = 1

где Щи) - дискретный аналог потенциальной энергии деформации; Q1, ... , QN -нормированные векторы внешних нагрузок, соответствующие различным видам нагружения; р1, ... ,pNм параметры нагружения.

Система нелинейных уравнений (12) задает поверхность в пространстве переменных Х=(и, р)т, зависящую от N независимых параметров. Если внешняя нагрузка зависит только от одного параметра, то есть N=1, то система уравнений (12) определяет кривую равновесных состояний деформируемой системы.

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

Для решения нелинейной задачи, описывающей поведение конструкции при действии внешней нагрузки, которая зависит от двух независимых параметров (при N=2), устанавливаются начальные значения параметров нагрузки р1 и р2 и начальные значения компонентов вектора узловых перемещений и. Уравнения (12) решаются с помощью методов дискретного или непрерывного продолжения решения для приращения параметра Др1 при фиксированном значении параметра р2. Таким образом, определяется кривая равновесных состояний, соответствующая заданному значению р2. После этого осуществляется возврат к зафиксированным начальным значениям и выполняется один шаг продолжения решения для приращения параметра Др2 при фиксированном значении параметра р1. Далее устанавливаются новые начальные значения параметров нагрузки и вектора узловых перемещений, после чего вычислительный процесс повторяется.

Движение по первому параметру р1 осуществляется одним из методов дискретного или непрерывного продолжения [9]. При дискретном продолжении функция Е(и, р1, р2) раскладывается в ряд Тейлора с удержанием линейных слагаемых, при этом предполагается, что значениер2 фиксировано:

Р (и + А и , р ! + А р р 2 ) =

д Р д Р

= Р (и , р ,, р 2) + —-А и + --Ар! ,

о и О р 1

тогда, с учетом (12), итерационный процесс запишется в виде:

У2Ж(ик(б'т1))Аик+1 = QlДpl к+1 + QlрД^1) + Q2р2(ли2) - УЩи^1)) (13) 1 2

где к - номер итерации; зт и - значения параметров продолжения; т и п -номера шагов нагружения соответственно по первому и второму параметру; ик(лт') , р1к(^т1) - значения искомых функций на итерации с номером к; Дик+1, Др1к+1 -приращения функций на последующей итерации. Система уравнений (13) дополняется еще одним уравнением:

II 8иШ || 2 = (ДЛт1)2 (14)

(ао, 5иш) +ап+1 8р1к+1 = Д?т' , (15)

5ик+1 = 5ик +Дик+1 , 8р1к+1 = Зр1к + Ар1к+1 ; (16)

ао = (а,1 а2 ... ап)т

Для решения системы уравнений (13) совместно с (14) используется алгоритм, применяемый при однопараметрическом нагружении:

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

У2Ж(нк(зт1))и = б! ; (17)

Ч2Ж(ак(зт1))У = б1 рЛО + б2Р2 (О - УЖУ^т1)).

Зависимость между искомыми функциями Днк+1 , Ар1к+1 и известными векторами и, Vопределяется путем подстановки соотношений (17) в уравнения (13):

Днк+1 = иАр1к+1 +V (18)

2. Подстановка (16) и (18) в дополнительное уравнение (14) преобразует его к виду:

а(Др1к+1)2 + 2ЬАр1к+1 - с = о , (19)

Квадратное уравнение (19) имеет два корня, из которых выбирается тот, для которого угол между предыдущим значением вектора (5ик ,8рк)Т и текущим значением (8иш , §рк+1)Т является острым [8]. Если оба корня обеспечивают выполнение этого условия, то выбирается тот из корней, который ближе к линейному решению.

3. Вектор приращения Диш вычисляется по формуле (18) при известном

к+1

значении Др1 .

4. Определяются величины:

5нк+1 = 5ик +Днк+1 , (2о)

я к+1 я к , Л к+1

8р1 = 8р1 + Др1 . 1. Процесс итераций продолжается до тех пор, пока норма вектора (Диш , Др1к+1) не станет достаточно малой.

При использовании дополнительного уравнения, заданного в форме (15), выполняются следующие действия, аналогичные описанным выше:

1. Находятся решения вспомогательных систем (17), определяющие связь (18)

^ а к+1 1 к+1

для приращении Ди и Др1 .

2. После подстановки (16) и (18) в уравнение (15) и выполнения некоторых >бразований получается следующее выражение для величины Др1к+1 :

А 5т1 - (а о , 5 и к ) - ап +1 8 р 1 к - (а о , V)

А р 1

(а о , и ) + а

3. Вектор Дик+1 находится после подстановки Др1к+1 в уравнение (18).

4. По формулам (2о) определяются величины 5ик+1 и Зр1к+1 .

5. Проверка условия окончания итерационного процесса.

На рис. 2 представлена поверхность равновесных состояний цилиндрической оболочки с параметрами Я=2,54 м, 1=В=о,5о8 м, А=о,ооб35 м, £=3,Ю5*1о9 н/м2, и=о,3 и шарнирно-неподвижно закрепленными прямолинейными краями [9].

При численном построении этой поверхности для каждого фиксированного значения нагрузки р2 с помощью процедуры (17)^(2о) определялась кривая равновесных состояний относительно р1 и

Проекция результирующей поверхности на плоскость р1р2 дает возможность выделить зоны устойчивости и неустойчивости, при этом линии, разделяющие эти зоны, определяют множество катастроф, то есть множество точек, для которых функция полной потенциальной энергии системы имеет несколько сливающихся критических точек. Картина, получаемая на плоскости р1р2, позволяет определять

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

,н

Рис.2. Поверхность равновесных состояний пологой цилиндрической панели

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

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

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

3. В качестве ведущего параметра наиболее целесообразно выбирать длину кривой состояний равновесия рассматриваемой механической системы, при этом использовать вспомогательное уравнение, обеспечивающее итерации на сфере.

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

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

Библиографический список

1. Ворович И.И., Зипалова В.Ф. К решению нелинейных краевых задач теории упругости методом перехода к задаче Коши // ПММ, 1965, т.29, №5, с.894-901.

2. Григолюк Э.И., Шалашилин В.И. Проблемы нелинейного деформирования: Метод продолжения решения по параметру в нелинейных задачах механики твердого деформируемого тела. - М.: Наука, 1988. - 232 с.

3. Давиденко Д.Ф. Об одном новом методе численного решения систем нелинейных уравнений // ДАН СССР, 1953, т.88, №4, с.601-602.

4. Иванов А.С., Трушин С.И. Расчет несущей способности нелинейно деформируемых пологих оболочек с учетом начальных несовершенств / Пространственные конструкции зданий и сооружений: Исследование, расчет, проектирование. Вып.7. - М.: Изд-во ЦНИИСК им.В.А.Кучеренко, 1992, с.24-29.

5. Милейковский И.Е., Трушин С.И. Расчет тонкостенных конструкций. - М.: Стройиздат, 1989. - 200 с.

6. Петров В.В. Метод последовательных нагружений в нелинейной теории пластинок и оболочек. - Саратов: Издательство Саратовского университета, 1975. -119 с.

7. Феодосьев В.И. Об одном способе решения нелинейных задач устойчивости деформируемых систем // ПММ, 1963, т.27, №2, с.265-274.

8. Crisfield M.A. A Fast Incremental/Iterative Solution Procedure that Handles "Snap-Through" // Computers & Structures, 1981, Vol.13, N1, pp.55-62.

9. Ivanov A.S., Trushin S.I. Development of Solution Procedures for Elastic-Plastic Large Deformation Stability Analysis of Shells // Nonlinear Analysis and Design for Shell and Spatial Structures. Proceedings of the SEIKEN-IASS Symposium, Tokyo, 1993, pp. 89-96.

10. Lahaye M.E. Une metode de resolution d'une categorie d'equations transcendentes // Compter Rendus hebdomataires des seances de L'Academie des sciences, 1934, v.198, N21, pp.1840-1842.

11. Mileikovskii I.E., Trushin S.I. Analysis of Thin-Walled Structures. - New Delhi: Oxford & IBH Publishing, 1994. - 187 p.

12. Ricks E. The Application of Newton's Method to the Problems of Elastic Stability // J. Appl. Mech., 1972, 39, pp.1060-1066.

13. Wempner G.A. Discrete Approximations Related to Nonlinear Theories of Solids // Int. J. Solid Structures, 1971, Vol.7, pp.1581-1599.