Решение задач ползучести методом конечных элементов Текст научной статьи по специальности «Физика»

" \ '

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Т о м XII 19 8 1 №6

УДК 629.735.33.015.4—977

РЕШЕНИЕ ЗАДАЧ ПОЛЗУЧЕСТИ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ

Г. Н. Замула, В. А. Павлов

С использованием метода конечных элементов решается задача расчета напряженно-деформированного состояния тонкостенных конструкций в условиях ползучести.

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

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

I. Обозначим (см. рис. 1): х = (хи х2, х3) —точка тела (конструкции) с координатами хи х2, х3, г— номер элемента, г — 1, 2, . . . , ./V,

V — объем элемента с номером г, г, /, .. . — номера узлов элемента,

{и (х)}'. — вектор-столбец перемещений элемента с номером г,

{УУ—вектор-столбец перемещений узлов элемента с номером г,

[ ] —матрица,

Рис. 1

t — время,

Т— температура,

1(9)

д<р(?)== ср(9) — <р(?-1) = Г — приращение функции ср на /(9-1)

интервале времени <7 = 1, 2,...,

/0, £—начальное и конечное время,

Т0 — начальная температура.

Предполагая, что кинематические гипотезы метода конечных элементов в форме перемещений сохраняют свой вид и для соответствующих приращений Ли, А У в условиях ползучести, получаем следующую сводку соотношений.

Основные соотношения МКЭ [1]. Соотношения МКЭ в приращениях.

Зависимость перемещений точек элемента от перемещений

его узлов

№)}'=[ | {Ьи(х)У = [У(х)У{АУ\г;

соотношения связи между деформациями и перемещениями {в (*)}'=[Л]'{и (*)'}', ] {Дв (*)}' = [Л] {Ди(*)}г; деформации элемента в зависимости от перемещений его узлов

{«(*)}'=[£]'{Г}% I {Ле(*)К=[£ПД5Т, где [В\'=[А][)/(х)]'.

Физические соотношения в случае упругого тела представляются обобщенным законом Гука.

Так, для трехмерного изотропного тела они имеют вид

+ ч -г — ' о;- (О

Те же соотношения при неустановившейся ползучести материала по теории течения или изотропного упрочнения запишем в форме [2]:

гч [2 а

]+-§-«<»>' т-Р){°ч-ЧЩ +

• \,'МТ Щ. (2)

где р — Ь в случае теории течения, р — £,

Є = а£(а, Т, є)

в случае теории упрочнения, точка сверху обозначает дифференцирование по времени t, а — интенсивность напряжений, е —интенсивность скоростей деформаций ползучести. Отсюда получаем соотношения связи между напряжениями и деформациями элемента

{о-(л:)}г = р]({•(*)}»■-- Ы*)Н)> {ДзС*)}'= [£)1({Де(-*)}г- {Д£о(*)К), где

<а(Т-Т0) а (Г —Т0)

*(Г-Г0)

о о о

{є0 (х)У =

вектор начальных температурных деформаций элемента;

{Дво (*)}' =

А а(Т-Т0) Г Деп

Да (Г - Т0) Дв22

Да (Г-Г0) ■ + Дез"

0 2Дє?2

0 2Дг?2

0 2Д£?з

:{Део.(*)Г + (Л£о (*)!'

(4)

вектор приращений начальных деформаций элемента;

{Дг"(л:)}г —вектор приращений начальных деформаций за счет ползучести,

1(9) <(?)

Дву= | ТШк, т> Р)(ац-ьцЦг)М> Дг= і а£(3> г)Л> (5)

М-1)

/С?-1)

Е, V приняты не зависящими от t.

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

[*Г = /[ЯГ[Д] \BYdV,

уг

вектор приращений обобщенных внешних объемных сил Е в уз-

лах

вектор приращений обобщенных внешних поверхностных сил /, заданных на части поверхности тела 02)

{ЬрУ= $[УГМ)с1Я1,

В2

вектор приращений обобщенных объемных сил от приращений начальных деформаций:

{Д/>т К = / [Д]'" [£>1 {Де?

V

за счет температурного расширения,

{др«}.= $[ВУ'[В\№{х))гйУ- (6)

V

за счет ползучести.

Штрих над матрицей вверху справа обозначает операцию транспонирования.

В случае упругого тела уравнения равновесия для конечноэлементной модели имеют вид:

[К] {Г} = (7)

{/} вектор нагрузок, получаемый суммированием по всем элементам обобщенных сил [Р\г, {р}г, {Рт}г и сосредоточенных нагрузок в узлах, [/С] — матрица жесткости конструкции.

Уравнения равновесия в приращениях в условиях ползучести для конечно-элементной модели приобретают вид:

[ЛТ] {А К} = {А/}, (8)

где {Д/} — вектор приращений нагрузок, получаемый суммированием по всем элементам приращений обобщенных сил {АР}Г, {Ар\г, {ДРТГ, {АР"}'', и приращений сосредоточенных нагрузок в узлах. Они полностью идентичны предыдущим уравнениям с учетом дополнительных нагрузок от начальных деформаций ползучести.

Кинематические граничные условия на части поверхности тела 2, и в том и в другом случае учитываются обычным образом [1].

Для решения задачи по времени можно использовать различные шаговые схемы [3].

Для явной схемы переходом с шага на шаг по методу Эйлера алгоритм можно представить следующим образом:

1) на нулевом шаге <7 = 0, t — = решается обычная задача

теории упругости при начальном „мгновенном11 нагружении и нагреве, находятся и запоминаются {К(°)}, декомпозированная матрица [Л'Р> матрицы [/(], при решении системы уравнений (7) методом последовательного исключения неизвестных, напряжения в элементах {а(о>}'=[0] [Ву{Уту- [О]{г0т}';

2) по полученным напряжениям находим приближенно приращения деформаций ползучести в элементах

Ае;7 = 4^а1<?“1)’ Т<9_1)’ ^_1)) Л^> =

= «р(а(«-1), Л*"1», /*9-4) Д*<«>, (9)

при # = 1 и соответствующие им приращения сил (6) в узлах элементов; к

' 3) по {АРп)г и приращениям других сил в узлах элементов

формируется вектор приращений нагрузок {А/} для всей модели и решается (в варианте повторного нагружения обратным ходом метода исключения неизвестных) система линейных алгебраических уравнений (8);

4) по полученному {ДК} вычисляются приращения напряжений в элементах {Да}7" = [О] [В]г{АУу— [О] {Дв0}г, новый вектор перемещений (К<°)}-{-(Д К} и новые напряжения в элементах {а(°)}г 4- {Да}'-, г=\, 2, ...,ЛГ;

5) расчет повторяется с пункта 2) при <7 = 2, 3 ... и т. д. до .

Алгоритм неявной схемы с переходом с шага на шаг по методу Эйлера —Коши с итерациями отличается от предыдущего применением на этапе (2) формул

ае?/=4т [? ыг°, Р(‘‘-гу)+? ы?>, т /?<*)] мы (ю).

Г - ' ' : ' ; , - - ■ - ' ' • . ; ' '• !~ 1 ’

с итерационным уточнением напряженного состояния №(п)(п —

= 0, 1, 2,...) на последующем шаге ^ при начальном приближении из предыдущего шага (<7—1) либо полученном экстраполяцией с предыдущих шагов.

„ ' д/9)

Схема типа предиктор-корректор: делается полушаг с —^— по явной схеме, определяется напряженное состояние 0>. 2 ' при

д^(?)

^ = ---2—> после чего выполняется полный шаг с вычисле-

нием приращений деформаций ползучести по формулам

Дв«. = ?(а[9“А Г'9 (11)

Здесь сделаем два замечания.

Лишь в отдельных случаях можно получить явные выражения для (АР"}'’ через {Дг"К без численного интегрирования в (6). Очевидно, это можно сделать для конечных элементов с постоянными по элементу деформациями, напряжениями и температурой, когда [В], [/3>] и Т элементов не зависят от х. Указанными свойствами обладают изотермические стержень, плоский треугольник, тетраэдр с линейной аппроксимацией перемещений (симплекс-элементы) и постоянными вязко-упругими свойствами. Для них справедливы преобразованные соотношения (6):

{АРпу=1Су{Ьвп0}г, $[Ву’[0]с1У. (12)

V

При использовании описанных выше шаговых схем для всех этих элементов вычисляются постоянные в них величины Да", аг;. и запоминаются ®у. В остальных случаях, когда напряжения и деформации ползучести переменны по элементу, нужна процедура численного интегрирования в (6) с вычислением и запоминанием

необходимых для этого величин Де» <з,у. в определенных точках элементов в соответствии с применяемыми квадратурными формулами (см. ниже).

Соотношения (9)—(11) шаговых схем записаны в виде, характерном для теории течения, когда р = 1, или для частного случая установившейся ползучести g = g(?, Т). При использовании теории изотропного упрочнения (/> = е) они должны быть дополнены такими же формулами для Дз в соответствии с (3), (5). Характеристика накопленной деформации ползучести е должна в этом случае подсчитываться и запоминаться для всех конечных элементов соответственно их виду. Для описания и вычисления функции g(a, Т, е) одноосной ползучести материала можно использовать любые известные приемы, включая многочисленные аппрок-симационные формулы, таблично-интерполяционные методы, а также косвенный способ, опирающийся непосредственно на кривые ползучести е = /(а, Т, Ь) в приращениях (см. [2]).

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

Все три элемента удовлетворяют важному условию непрерывности перемещений на границах элементов.

Элемент стержень предполагаем прямолинейным, произвольно расположенным в пространстве, имеющим постоянное поперечное сечение Т7, выполненным из одного материала, равномерно нагретым до температуры Т и работающим только на растяжение — сжатие (см. рис. 2, а). В этом случае при линейной аппроксимации

Г J ос 4 т - ,

г 0 j

а) <П ,

Рис. 2

перемещений их элемента вдоль его местной продольной координаты х для приведенных в п. 1 основных соотношений МКЭ имеем:

{и(х)У=их(х), {ГК=|Р}, [У(х)]г = [1 ~ х/1 х/1],

{е(х)У=вя, [А]=±, [3]' = [-4- {«(*)}' = «„ \D\-E.

Матрица жесткости элемента

где Р—площадь поперечного сечения, / — длина стержня.

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

{Арту = Р I [__1_ _1_] £Да(Г- Т0) йх = ЕРДа(Т— Г0){ “| } ,

за счет ползучести (явная схема)

, {ДРпУ= Е / [-----г -у-]'ВД Лх = ЕЕДе" {

М)

где Де” = Г g{зx, Т, p)^xdt^g(Jr1), Т(9"1), Д,,)а?"1)Д^)-

*(?-!)

Рассмотрим методический пример установившейся ползучести „мгновенно11 нагруженного в точке / силой Р;- и нагретого до температуры Т пояса, далее условия по времени не меняются (рис. 2, б).\ В точке I наложены кинематические граничные условия в виде заделки (К1 = 0).

Уравнения равновесия (8) в приращениях без учета граничных условий в этом случае имеют вид:

ЕР_

I

1—1 ■1 1

Начальные условия при t = tQ = 0 для задачи ползучести дают решение упругой задачи (с учетом граничных условий)

У(Р = ЕЕа(Т — Т0) + Р;,

откуда

ур=^1+чт-т0)1,

^0)д^.

Далее на шаге д(д = 1, 2, . . . ) имеем ^ДГ)9) = ЕЕ Ае"(9),

, д • \

АУ)9) = Де"<9) I, Д<#> == Е —I-----------------Де” <[я) = I

откуда

ДГ<-9) = Ае1{ч) I. До1?)*= Е ( - 1

А^ = ё№)^А^} .

р.

Таким образом, в этом случае не меняется по

времени, что является точным результатом и следствием статической определимости задачи. Решение этой задачи методом конечных элементов для перемещения

Y<■p = ^l+a(T-T0)l-]-g (а'0)) а£>2 А№

9 = 1

точно совпадает с аналитическим решением, линейно зависящим

от времени ( У\М^) = Ц9)

\?=1 -)

Для плоского треугольного элемента постоянной толщины с узлами г, /, к (рис. 3) имеем

{и(х)}Г =

У

, {е(х)]г = .

[ ах '

■- . {а(*)}г =

2е . алу .

, {У}'=

[У{х)У

[ВУ

1

2/^

V, О V) 0 Ук о о I/, о 1/;. о КА.

61 О 0^0 0 с1 0 с) 0 ск А Ь, с, Ъ} ск Ьк _

, И1 =

дх

0

0 7Г

ду д д .ду дх,

Уг

2Г

г (а, +..Ь1х + с,у), а1 = х]ук —хку],

]

^ = Уу —У». С1 = хк = х}, г 4. Матрица жесткости элемента: [Я]г=^8[5Г [£>]№

Щ

К

и,

*0

и*

Рис. 3

где /Гг, 8 — площадь и толщина элемента, а для плоского напряженного состояния

1 V 0

[О]

1 — №

V 1 0

1 — V

о о

Начальные деформации от температуры:

от ползучести:

Де2'К = Да(Г- Т0)

{Д$о}г =

Де”

д4

2Де"у

где согласно (5)

Д гпху.

Ля)

/(?)

д,;= |

*(?-!> ' ^ '

<(?)

д£;= | 1?(а,т,р)(а У-Ц)<и,

/«7—1) ' ^ 7

Г I £(3> Т, Р)*хуЫ, а=/с^ + Оу +Зо^у—0ж0у. Дополнительные нагрузки в узлах от начальных деформаций

{ДРт}'=[ОНДгоТК,

{ДРП}-- = [в]" {ДеоТ,

где [0] = [5ПД]-^8:

Четырехугольный плоский Элемент общего вида в пространстве определен точками г, у, / (рис. 4).

Функции формы V (х) в этом случае примем в виде [1]

к,=4-0-9(1-ч). у>=4-о+ Е)о-ч),

V*

4-(1 + «)(1-.+Ч)у ^=4-о-0(1+4).

причем связь координат I, т; и х, у (см. рис. 4) выражается соотношением

х = V, х1 + V] х] + Укхь + Уг Х[,

У — У1У1 + V )У 1 + УкУ и + ^гУг-

Посредством замены переменных интегрирования и с помощью соотношений

йУ = | det [У][ ос/? с1г\, р =■ I, у, I

приходим к следующим выражениям для матрицы жесткости элемента

1 1

[КУ = в | С \ву [О] [В]' I йег [У]| йц

-1 -1

и дополнительных нагрузок в узлах от начальных деформаций

{ДРТ}Г = 51 I [ВУ [О] (Дго Г |<^ [У] I а\ й-п,

-1 -1

{ДР”}' = 8 | I [ВУ [£>] {Д£оп}г^[Урит1.

—1 -1

Интегралы определяются численно с помощью квадратурных формул Гаусса —Лежандра [1] при вычислении и запоминании аХУ ау, о е в точках интегрирования, число и расположение которых в элементе может быть различным в зависимости от степени нелинейности функции формы и закона ползучести.

Так, при наиболее распространенном законе ^ = Л(7)а>'—1 эти точки при ^=1 выбираются в соответствии с рекомендациями [1]

для вычисления [к]г, а при [*>1 на больше в направлении

каждой из координат элемента.

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

3. Для реализации предложенных методов решения задач ползучести создан комплекс программ на алгоритмическом языке РСЖТИА^МУ на ЭВМ серии ЕС. На первом этапе реализована явная шаговая схема для случая, когда нагрузка и нагрев не меняются по времени, а тонкостенная конструкция, плоская или пространственная, представляется набором элементов трех описанных в п. 2 типов.

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

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

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

сегмент 1—считывает данные с магнитного диска и управляющие данные с перфокарт;

сегмент 2 —строит матрицу жесткости конструкции, вызывая, в свою очередь, сегменты второго уровня, которые вычисляют матрицы жесткости отдельных конечных элементов;

сегмент 3 — формирует и решает систему линейных алгебраических уравнений (7), (8) ленточным вариантом метода Гаусса; сегмент 4 — вычисляет напряжения и усилия в элементах;

Рис. 5

сегмент 5 — определяет приращения деформаций ползучести, приращения вектора нагрузки,

сегмент 6 — производит печать результатов.

Комплекс программ имеет минимум ограничений, снабжен средствами управления, диагностики и прерывания решения; обмен данными между сегментами осуществляется через файлы на дисках. По вводу исходных данных, организации и отдельным модулям комплекс преемствен с программой расчета упругого НДС „Отсек“ [4].

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

На рис. 5 сплошными линиями показаны результаты расчета из [5], а крестиком отмечены результаты расчета по МКЭ.

Затраты машинного времени ЭВМ по МКЭ превышают соответствующие затраты при расчетах обычным методом сеток, однако класс решаемых задач значительно шире.

ЛИТЕРАТУРА

1. Сегерленд Л. Применение метода конечных элементов.

М., „Мир“, 1979,

2. Б и р г е р И. А. Термопрочность деталей машин. М., „Маши-ностроение% 1975.

3. 3 а м у л а Г. Н., И в а н о в С. Н. Ползучесть подкрепленных панелей при нестационарном нагреве. „Ученые записки ЦАГИ“, т. VII, № 5, 1976.

4. Д у б и н я В. А., Д у д а р ь к о в Ю. И., 3 а м у л а Г. Н.,

П а в л о в В. А. Программа расчета напряженно-деформированного состояния конструкций летательных аппаратов методом конечных элементов „ОТСЕК-О*. В сб. „Расчеты напряженно-деформированного состояния авиационных конструкций". Труды ЦАГИ, вып. 2063, 1980.

5. За му л а Г. Н. Плоская задача термоползучести в напряжениях. „Ученые записки ЦАГИ", т. IV, № 6, 1973.

7 — .Ученые записки ЦАГИ" № 6.

Рукопись поступила 10\111 1980 г.