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

Наука к Образование

МГТУ им. Н.Э. Баумана

Сетевое научное издание

Наука и Образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2015. № 03. С. 296-312.

Б01: 10.7463/0315.0759406

Представлена в редакцию: Исправлена:

© МГТУ им. Н.Э. Баумана

УДК 539.3

Конечно-элементное моделирование процессов термоползучести на основе методов Рунге-Кутты

Димитриенко Ю. И.1*, Губарева Е. А.1, Юрин Ю. В.1

26.02.2015 10.03.2015

сИггШ ,ЬтБШ@ gmail.com

1МГТУ им. Н.Э. Баумана, Москва, Россия

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

Ключевые слова: метод конечного элемента, термоползучесть, метод Рунге -Кутты, трехмерная задача термоползучести

Введение

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

настоящее время реализованы в основных коммерческих программных пакетах, в том числе в АКБУБ [4]. Однако, в большинстве случаев для решения нелинейных уравнений ползучести применяется или явные методы или неявные, основанные на методе Эйлера аппроксимации производных по времени. Метод Эйлера достаточно эффективен с точки зрения экономии оперативной памяти при проведении вычислений, однако относительно затратен по времени вычислений и не всегда обеспечивает требуемую точность расчетов деформаций ползучести.

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

Постановка задачи нелинейной теории термоползучести

Рассмотрим краевую трехмерную задачу термоползучести на основе теории течения

[10]:

V-о = 0,

0

О = 4 С • <£ - £° - £)

£С=Р(£С,С)

£ = def(u) = и + (V® и)г) (1)

8° I = 0

1г=0

О - П|8Ь (т) Н!„ = иЬ (т)

Здесь о - тензор напряжений, £ - тензор малых деформаций, £° - тензор деформаций ползучести, £в = а(в-в0)- тензор тепловых деформаций, Р(£°, о) - дважды непрерывно дифференцируемая в некоторой области тензорная функция, описывающая

модель скоростей деформаций термоползучести , а - тензор температурного расширения, и - вектор перемещений, п - вектор внешней нормали, ^ - заданный вектор напряжений

на части поверхности тела ^ , V - набла-оператор [11,12], ® - знак тензорного произведения, (•)- скалярное произведение, 4С - тензор модулей упругости, иь - заданное перемещение.

Температурное поле в будем считать известным, распределенным неравномерно в теле. Среду будем считать изотропной, поэтому для нее

4 С = ЛЕ ® Е +

0

£ = ае Е(в —в0)

Р(£с,о) = [< | Гв-11(в)Е

л

\ат у

У (2) где I (о) = Е--о - первый инвариант тензора напряжений, Е - метрический тензор, <Уи -второй инвариант (интенсивность) тензора напряжений [11], £°и - интенсивность тензора скоростей деформаций, <т ,ш, л, г - константы ползучести, Л,л - упругие константы (параметры Ламе).

Вводя тензор упругих напряжений ое = 4С ••£, тензор напряжений ползучести

0

ос = —4 С • -£С и тензор термонапряжений ов = —4 С •• £, тензор напряжений может быть записан в виде следующего разложения:

е . с . в

о = о + о + о (3)

Численный р -стадийный явный метод Рунге-Кутты

Для численного решения исходной задачи (1),(2) применим р -стадийный явный метод Рунге-Кутты [13-16]. Для этого тензорное соотношение ползучести:

¿с=Р(£с,о) ^

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

Введем сетка по временному параметру: Аы =(гп =0,...,ты =Т) с шагом

8ф(Аы) = (Агп = т, -гп,..., Дгж ! = тм -тМ1). Введем следующие обозначения:

пусть (.)(т) - объект, вычисленный на т -м шаге по временному параметру

(т £ а (.)(тУ - объект, вычисленный на 5 -й (х £ 2р) стадии этого шага. При этом

стадия 5 = 0 является условной и служит для инициализации соответствующих объектов на данном шаге по временному параметру, не приводя к дополнительным вычислениям. Пусть также А = (а)рр - матрица Бутчера рассматриваемого метода Рунге-Кутты, и

w = ^,..., ^ ) - вектор весовых коэффициентов стадий метода. Для явного метода

а) = 0, ) > I.

Для начального шага (т = 0) на основе начального условия £и = 0 имеем следующую систему:

V- ое (0) —-V- о6 ое(0) = 4 С ••£(0) е (0) = ёе^и(0))

о

е(0)

• п = S(0) - о6- П

и (0)Е = иГ

|Е" (5)

Из формулы для разложения тензора напряжений и указанного начального условия имеем следующие соотношения для тензоров на начальном шаге:

1ес(0) = 0 ,(0) _ _е(0)

о

+ о6

(6)

Далее на т -м шаге (т £ имеем следующую вычислительную процедуру для

р -стадийного явного метода Рунге-Кутты: 1. Инициализация тензорных полей:

I „с(т) _ С|>-1)

I е(0) — е

с(т) __с(т)

(0) — -С' 'е( *)

(7)

2. На ^ -й стадии (5 £ №р+1):

2.1. Вычисление тензоров деформаций и напряжений ползучести ^ -го

этапа:

2.2. Решение краевой задачи ^ -й стадии и нахождение соответствующих

тензоров упругих напряжений:

\7 е(т) _ гу (с(ш) в\

У' ®(*) (О(*) + О Л

О

е(т) _ 4

- 4С-8

(т)

(*) ^ ) > 8 т - ае^и (т))),

е(т) I _ о(т) ( С(т) , в\ I О(*) ' Пу - ®Ь " (О (*) + О ) ' П

„(т Л — ,,(т)

и (*) I - и ь

11II

2.3.

Вычисление тензора напряжений:

_(т) _ _е(т) , _с(т) , в - + +О

3. Переход к следующему временному шагу:

С(т) _ с(т) 8 8(Р) '

с(т) _ с(т)

О - О(р) ,

_(т) __(т) О - О(т)

(10) (11)

Различные варианты метода Рунге-Кутты

Приведем матрицу Бутчера А и вектор весовых коэффициентов w для конкретных методов Рунге-Кутта, которые будут далее использованы при численном расчете:

1. "Классический" 4-х стадийный метода Рунге-Кутты (далее просто метод Рунге-Кутты), для него А - матрица 4-го порядка:

А -

( 0 0 0 0^ 12 0 0 0

0 У 0

2 0 0 0 1 0

w -

(>3

1/ 1

'3 /6

2. 3-х стадийный метод Рунге-Кутта 3-го порядка:

А -

0

Л

0 0 0 0

К 0,

/3 У

^-(34 о

3. 2-х стадийный метод Рунге-Кутта второго порядка:

(0 0^

А -

V1 0У

, w -(12 1/2)Т

4. Одностадийный метод Эйлера первого порядка:

А-(0), w-(1)Т

Т

Т

Метод конечных элементов с применением метода Рунге-Кутты для численного решения краевых задач термоползучести

Пусть -

\\Т2(1)(П) = {и : П и. е }

- соболевское пространство вектор-функций, определенных на ограниченной области О с липшицевой границей ЗО [1]. Пусть ТтЕ :^Ь2(£) - оператор, сопоставляющий

элементу и е ^()(О) след на £м. Выделим в пространстве ^()(О) подпространство вектор-функций Х(О) = |и е ^(1)(О) Тг^ (и) = о|, а также подпространства вектор-функций, обладающих классическими (не обобщёнными) производными всех порядков: V(О) = {и е ^(1)(О)Щ е С(ш>(О)}, У0Ш(О) = {и е Vш(О)|и|Е = о} . В пространстве V"(О) рассмотрим задачу для - -й стадии метода Рунге-Кутта:

V- ое(т) =-У- (о ^ + ое ),

о

'(-) е(т) _ 4

(-) г,(т) —

= 4 С ••£

(т) (-) '

£ т=аег(и о)),

(т )ч1 (-)

(11)

о£т) • п1 = sm - (оС« + о* ) • п

(т) (т)

и(£ = иЬ )

V Щ

где и(т) е Ь2(£и), 8Ьт) е Ь2(£ет). Будем предполагать, что тензор модулей упругости 4С удовлетворяет условию положительной определенности [10], т.е. для всякого симметричного тензора второго ранга £ справедливо неравенство( у > 0):

£--4С••£ >у£••£

Для вывода слабой формулировки потребуется следующее тождество, справедливое для дифференцируемого симметричного тензорного поля Т и векторного поля а [11]:

а• ^Т) ^^(а• Т)-ёе£(а)••Т

(13)

Умножив первое соотношение рассматриваемой системы скалярно на V е (О),

применяя последнее тождество (для каждого оператора дивергенции тензора) и интегрируя по области О получим:

| ёеад • о^йУ = ^ (V • (о*0 + о-т + ов)Ж-{ ёеОД ••(о2)т) + ов) йУ (14)

О О О

Применяя для первого интеграла в правой части формулу Остроградского и учитывая силовое граничное условие в системе, получим:

| ёеГ(V) • •о^йУ =| у • 8(т)й£ -1 ёеГ(у) • • (о^ + ое ) йУ

Далее подставляя определяющие соотношения, будем иметь: | ёеГ (V) • •4 С • •ёеГ(и(т)) )йУ = | V • Э(т)й£ +1 ёеГ(у) • •4 С

) + £ | йУ

(15)

(16)

О

£

О

О

£

О

Последнее соотношение будем рассматривать теперь в пространстве ^( )(О), предполагая все производные входящее в него - обобщёнными. Тогда под слабым решением краевой задачи на 5 -й стадии метода Рунге-Кутта в пространстве ^2( )(О)

будем понимать такую вектор функцию и(™) е W21)(Q) , что если w(т) е W2(1)(Q) такая вектор-функция, что Тг^ (w(т)) = и(т), то и(т) - w(m) е Z(Q) и и(^ удовлетворяет последнему интегральному соотношению Уу е Z (О). Решение такой задачи существует и единственно [13].

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

Где >МП - пространство столбцов размерности п. Также введем матрицу дифференциальных операторов Б и матрицу модулей упругости С Е ¿(М, 6,6) в следующем виде [16]:

Б =

I4 0 0 1 I с1111 £>1122 с1133 0 0 0 1

0 д 2 0 ^>2222 £>2233 0 0 0

0 д 2 0 д, дз 0 , С = £>3333 0 2323 0 0 0 0

д 3 0 д> сим. с1313 0

10 д з д 2 У V /">1212 с У

(18)

Тогда интегральное соотношение можно записать в следующем виде:

БУ) СБи^ШУ = | уТ8(тЧЕ +1(БУ) С

Г 0

-.с (т)

'( 5 )

+ £ ШУ

(19)

Слабое решение краевой задачи будем искать на основе метода конечного элемента. В качестве конечного элемента будет использован 10-узловой тетраэдр с квадратичной аппроксимацией. Обозначим т(.) Лебеговы меры в М3и М2 соответственно. Пусть построена регулярная триангуляция Т ( Ь - максимальный диаметр конечного элемента в триангуляции) области О с узлами ЫШ(Ть) с М3:

ОХЯ а Ол = У К, /л{К) < 5(К)

к^т,,

где 5 > 0 - точность разбиения. Пусть = {К еТА : т(дК одОА) ^ 0} и ГА (К) = дК п>дОл для К е Гй6. Тогда справедливо разложение: ГЛ = дС\ = ГЛ(Х).

КеТЬ

На границе Г h введем поле нормалей:

nh(x)=

n r (x), x g г (K )\ ar (K) ^ т(Г ( K ))n k (x)

, x g^ (K)

KgT, :xGK

S m(rh ( K))

KGTh :xGK

а также функцию |(x) = min iy : x + ynh (x) g 3Q} . Триангуляцию T будем называть

И

согласованной с границей, если:

1. Функция |(.) определенна на всей границе Гh.

2. Справедливо неравенство: ||(x)| < diam(K) если x g^ (K).

Триангуляция согласованная с границей позволяет всякой функции f (.) c D( f) = Г поставить в соответствие функцию f (x) = f (x + |(x)^(x)), для которой: D( fh) = ГА . Построим таким образом вектор-функции u^), S(m), соответствующие вектор-функциям

(m) о( m) ub , Sb .

Введем далее конечно-элементное пространство Fh (ПА) в следующем виде:

Где Ф^: К ^ L(M, 3,30) - матрица базисных аппроксимирующих функций (функций формы), а Л- = (./;', /к2 ./о • • •.//'1 ./а'" ./о) - значения вектор-функции fh в узлах тетраэдра K g T (называемые степенями свободы КЭ). Матрица Ф^ (.) имеет следующий блочный вид:

ФК (x) =(Ф11 ф22 Ф33 Ф44 Ф12 Ф13 Ф14 Ф23 Ф24 Ф34),

Ф j (x) = Nj (x) E,

f^i (x) (2L (x) -1), i = j [4L, (x)Lj (x), i * j '

где Ц (x) - барицентрические координаты точки в тетраэдре K g T, построенные по его вершинам, i, j G М5, E GL(R, 3,3) - единичная матрица. Пусть далее Suh, -аппроксимация частей Ем и границы Г = сО полученная при триангуляции Th .

В пространстве Fh (QÄ) выделим подпространство

ZA (QÄ ) = | f g F (Q/, )| f (x) = 0, x g Suh ^ ГА} . Тогда МКЭ схемой поиска приближенного решения ослабленной задачи на s -й стадии метода Рунге-Кутты на регулярной триангуляции T будем называть задачу поиска такой вектор-функции u^ g F (Qh), что

uS(x) = u(mЧ*) при x g Nd(T) (т.е. интерполирует вектор-функцию u(m) на границе Suh) и для любой вектор-функции vh g Zh (Qh) справедливо соотношение:

Nj (x) = ■

(21)

j (DVh )T CDutfhdV = j vhS(b'm)dS + j (DVh )r C

^I dV

(22)

Q

S

Q

h

h

Раскладывая интегралы по всей области на интегралы по конечным элементам, и подставляя определение элементов пространства ^ (^ ) получим:

X (^ )" —к ] = 0 (23)

К еГй

AK -J BKTCBKdV, fK - J (ф к /S^dZ + J BKTC [е$? + °е ) dV, BK (x) = ^ (x).

к

SKn!^ K V

>m>30

AG -AK -

( AG )

G G -

J N

f f ^N

(ÇG ) ^ J(^K )' , P -V (K ) V K) |o, p £v(K)

(24)

Матрица AK EL (M, 30,3 0) и вектор fK = M3" представляет локальную матрицу жесткости и локальный вектор нагрузок, соответственно. Пусть введена единая нумерация о: Th —степеней свободы триангуляции (инъективная функция, ставящая в соответствие каждому КЭ K е T кортеж номеров соответствующих ему степеней свободы), число которых в Th обозначим N = 3card(Nd(Th)) . Тогда введем матрицу AK E L ( М, N, N) и векторы /£, v£, и(™)G E > MN

' \(Ak )., p = vt (K), q = v} (K)

0, p, q iu(K)

* = {f,v, }.

Построим далее на основе них следующую матрицу и векторы:

KeTh KeTh

Тогда МКЭ схема может быть записана в виде:

(v°)r[^4r-/°]=0

Далее поскольку vA е ZA(QA), то (vG) = 0,р е / , где / = |/; е d(/Q: У/Г е : оЕнЛ) ^ 0 j. Кроме того, компоненты (z/(("j)G ) зафиксированы значениями u(ь™) в узлах Nd(Th) . Но тогда указанные компоненты

(25)

(26)

могут быть исключены из векторов \>'' и и('")Г' (в результате получим векторы Vй, и

,G ,Xm)G

СО

q = card(I)), а известные компоненты в ü'1"'" могут быть перенесены в правую

СО

часть:

r= \{f'G)={fG)-T{¿G) У

V psi ■ J (27)

Тогда исключив из матрицы строки А° и столбцы, а из вектора f"' строки с

номерами p e I (получив матрицу AG EL ( M, N — q ,N — q) и вектор f

G

>ra?N-<?

) и

учитывая произвольность вектора Vе Е > приходим к разрешающей СЛАУ МКЭ

схемы:

A4m)G - fG.

(28)

N

Семейство МКЭ схем при И ^0 сходится к слабому решению краевой задачи [13]. Решение разрешающей СЛАУ в данной работе осуществляется на основе метода сопряженного градиента.

Тестовая задача ползучести бруса при одноосном растяжении

В качестве тестовой задачи для сравнения точности предложенного метода МКЭ с различными вариантами метода Рунге-Кутты рассмотрим задачу растяжения балки под действием продольной силы, с учетом деформаций ползучести. Пусть 12 -

противоположенные торцы балки, п = е - нормаль к торцу балки.

Рисунок 1 К задаче об одноосном растяжении бруса

Тогда исходная задача (1) принимает вид:

V-о = 0, о = 4 С "(8 - 8е ), ёс=¥(ес,а), е = йе{(и),

8е = 0

1г=0

(29)

о - п и

I = РП

1 2

0, V/и

и1 к=

0.

1х=0 ' 1х=0

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

Г л Л

* (8е , О) =

_ с

— Ш£Х1

Г/

с л Г

с

и

V Ст у

V

о — — /1 ( о ) Е

(30)

у

Это допущение позволит найти аналитическое решение задачи. Также для простоты предположим, что коэффициент Пуассона у = 0 (а, следовательно Я = 0, 2¡л = Е). Тогда тензор напряжений будем искать в виде: о = рег ® е. В таком виде он, очевидно, удовлетворяет уравнению равновесия и силовому граничному условию. Тогда определяющие соотношения в компонентах примут вид:

Р8л8а=Е{е1-ес1]) ¿< =

2р 3/;

Р

(1 )

г

£)7

З/7

кат

(1 -Ше^ )

(31)

т=0

= 0 и

где а Ф /,у = 2,3. Из последнего соотношения с учетом начального условия еС

первого определяющего соотношения будем иметь: еа/? = есар = 0.

Второе и третье соотношения представляют собой задачу Коши для ОДУ относительно компонент тензора деформаций ползучести с указанным начальным условием. Решая эту систему для компоненты аСх получим:

есхх(т) =

1

Ш

2

1 — ехр

2ш

V

3л

р Vе7т

шШ Ф 0

3л

р Vе7т

т,ш = 0

(32)

Для остальных компонент имеем:

(т) = —

Л 3^

р

\°т У

(т — щ(т))

^(т) =

1

Ш

1

т +

3л

г

2шр

\

-Л.

V р У

ехр

2ш р

3л

р

\°т У

—1

,Шф 0

(33)

3л

р

т\ш = 0

Из первого определяющего соотношения получаем:

т (т) = р ЯГ1

(т).

(34)

Найдем компоненты вектора перемещений. Из соотношения Коши для компоненты и и граничного условия на торце ^ имеем: и (х,т) = еп (т)х. Для компонент и2, и3 из соотношения Коши с учетом равенства нулю касательных компонент тензора деформаций, имеем:

и2 (у, 2, т) = е22 (т)у + + С2, и3 (у, т) = е33 (т)2 + А1у + А2

Из системы условий: и| = 0, Ухи|т п = 0 имеем: С = С = А = А = 0 . Таким

1х=0

образом, вектор перемещений можно представить как линейную функцию координат:

(х,т) =(е11(т)X е22(т)У, е33(т)2)т .

и^х, с; = т х, е т у, е т 2 ) . ^^

Сравним полученное аналитическое решение, в частности компоненту тензора еи

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

г

<

г

г

г

Параметр Значение

Модуль Юнга Е, ГПа 212

Предел прочности на растяжение <УТ , ГПа 0.4601

Показатель нелинейности деформаций ползучести г 13.19

Коэффициент вязкости ] 0.225

Коэффициент стабилизации ползучести Ш 105

Модуль приложенной силы р, ГПа 0.1

Число узлов сетки 1419

Число КЭ сетки 704

График аналитического решения приведен на рисунке 2. На этом же рисунке приведены результаты конечно-элементных расчетов кривой ползучести. Вследствие высокой точности численных методов. визуально аналитическое и численные КЭ-решения практически неразличимы.

0.000484

0.000482

0.00048

Й 0.000478 с о

& 0.000476 0.000474 0.000472 0.00047

0123456789 10

х104

Ъ ч

Рис.2 Кривая ползучести бруса при одноосном растяжении, рассчитанная с помощью трехмерного метода КЭ и метода Рунге-Кутты и аналитическая кривая ползучести

Для решения задачи была использована равномерная сетка по времени

Т

Аы = (т0 = 0,...,ты = Т) с шагом 8гр() = , где ^ = ^^ для р-стадийного метода

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

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

Кутты

Т , ч £11. аналитическое значение II со" еп, р = 2 е11. Р = 3 е11. Р = 4

0 4.716981132е-04 4.7169811320е-04 4.7169811320е-04 4.71698113207е-04 4.7169811320е-04

8760 4.7543991640е-04 4.7544709587е-04 4.7543986945е-04 4.75439916795е-04 4.7543991640е-04

17520 4.7778161048е-04 4.77799610503е-04 4.7778137367е-04 4.7778161442е-04 4.7778161041е-04

26280 4.7924708872е-04 4.79272444480е-04 4.7924658483е-04 4.7924710134е-04 4.7924708837е-04

35040 4.8016421385е-04 4.8019240234е-04 4.8016346076е-04 4.8016423910е-04 4.8016421290е-04

43800 4.8073816880е-04 4.8076567963е-04 4.8073724130е-04 4.8073820784е-04 4.8073816695е-04

52560 4.8109736110е-04 4.8112207661е-04 4.8109635034е-04 4.810974123е-04 4.8109735818е-04

61320 4.8132215072е-04 4.8134311384е-04 4.8132113834е-04 4.8132221084е-04 4.8132214669е-04

70080 4.8146282848е-04 4.8147987052е-04 4.8146187514е-04 4.8146289343е-04 4.814628234е-04

78840 4.8155086739е-04 4.8156427638е-04 4.8155001090е-04 4.8155093327е-04 4.8155086167е-04

87600 4.8160596388е-04 4.8161624330е-04 4.8160522238е-04 4.8160602745е-04 4.8160595771е-04

Значения относительных отклонений 5 численного решения от аналитического по каждому из методов приведены в следующей таблице 3. Результаты расчетов показывают, что метод Рунге-Кутты с увеличением порядка р позволяет получать все более точное решение - относительная ошибка с увеличением р на 1 уменьшается примерно на порядок. Метод Рунге-Кутты 4-го порядка дает относительную ошибку примерно на 3 порядка меньшую, чем метод Эйлера.

Таблица 3 Сравнительные значения относительной ошибки расчета деформации ползучести в задаче об одноосном растяжении бруса, полученные с помощью аналитического решения и МКЭ решения для

различных вариантов метода Рунге-Кутты

Т , ч II 5, р = 2 3 II кр 5 Р II 4

0 9.3319449390Е-14 9.3319449390Е-14 9.3319449390Е-14 9.3319449390Е-14

8760 1.5100451119Е-05 9.8757083730Е-08 8.1717363246Е-10 7.6329041303Е-12

17520 3.7672729793Е-05 4.9564170148Е-07 8.2387369449Е-09 1.5401800372Е-10

26280 5.2904672503Е-05 1.0514198680Е-06 2.6329473863Е-08 7.4136498245Е-10

35040 5.8702499488Е-05 1.5684045612Е-06 5.2590903113Е-08 1.9835110073Е-09

43800 5.7222949119Е-05 1.9293327567Е-06 8.1203691278Е-08 3.8458315951Е-09

52560 5.1370546951Е-05 2.1009557010Е-06 1.0654157698Е-07 6.0827583606Е-09

61320 4.3551317332Е-05 2.1033354663Е-06 1.2492279531Е-07 8.3594137534Е-09

70080 3.5395118141Е-05 1.9800987829Е-06 1.3490080938Е-07 1.0364486602Е-08

78840 2.7844654156Е-05 1.7786129656Е-06 1.3679428995Е-07 1.1878982039Е-08

87600 2.1343596583Е-05 1.5396425818Е-06 1.3199535375Е-07 1.2795923527Е-08

Заключение

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

Исследование выполнено за счет гранта Российского научного фонда (проект №1419-00847).

Список литературы

1. Nejad M.Z., Kashkoli M.D. Time-dependent thermo-creep analysis of rotating FGM thick-walled cylindrical pressure vessels under heat flux // International Journal of Engineering Science. 2014. Vol. 82. P. 222-237. DOI: l0.l0l6/j.ijengsci.20l4.06.006

2. Li B., Lin J., Yao X. A novel evolutionary algorithm for determining unified creep damage constitutive equations // International Journal of Mechanical Sciences. 2002. Vol.44, no. 5. P. 987-1002. DOI: l0.l0l6/S0020-7403(02)0002l -8

3. Работнов Ю.Н. Ползучесть элементов конструкций. M.: Наука,1966. 752 с.

4. Implicit Creep // ansys.net : a resource for ansys users : website. Available at: http://ansys.net/ansys/papers/nonlinear/conflong_creep.pdf , accessed 01.02.2015.

5. Димитриенко Ю.И., Яковлев Д.О. Сравнительный анализ решений асимптотической теории многослойных тонких пластин и трехмерной теории упругости // Инженерный журнал: наука и инновации. 2013. № 12. Режим доступа: http://engjournal.ru/catalog/mathmodel/technic/899.html (дата обращения 01.02.2015).

6. Димитриенко Ю.И., Губарева E.A., Сборщиков С.В. Aсимптотическая теория конструктивно-ортотропных пластин с двухпериодической структурой // Mатематическое моделирование и численные методы. 2014. № 1. С. 36-57.

7. Димитриенко Ю.И., Губарева E.A., Яковлев Д.О. Aсимптотическая теория вязкоупругости многослойных тонких композитных пластин // Наука и образование. MF^ им. Н.Э. Баумана. Электрон. журн. 2014. № 10. С. 359-382. DOI: 10.7463/1014.0730105

8. Димитриенко Ю.И., Яковлев Н.О., Ерасов В.С., Федонюк Н.Н., Сборщиков С.В., Губарева E.A., Крылов В.Д., Григорьев M.M., Прозоровский A.A. Разработка многослойного полимерного композиционного материала с дискретным

конструктивно-ортотропным заполнителем // Композиты и наноструктуры. 2014. Т. 6, № 1. С. 32-48.

9. Димитриенко Ю.И., Губарева Е.А., Юрин Ю.В. Асимптотическая теория термоползучести многослойных тонких пластин // Математическое моделирование и численные методы. 2014. № 4. С. 18-36.

10. Димитриенко Ю.И. Механика сплошной среды. В 4 т. Т. 4. Основы механики твердых сред. М.: Изд-во МГТУ им. Н.Э. Баумана, 2013. 624 с.

11. Димитриенко Ю.И. Тензорное исчисление. М.: Высшая школа, 2001. 576 с.

12. Димитриенко Ю.И. Механика сплошной среды. В 4 т. Т. 1. Тензорный анализ. М.: Изд-во МГТУ им. Н.Э. Баумана, 2011. 463 с.

13. Ректорис К. Вариационные методы в математической физике и технике. М.: Мир, 1985. 590 с.

14. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Бином, 2001. С. 363-375.

15. Фалейчик Б.В. Одношаговые методы численного решения задачи Коши. Минск: БГУ, 2010. 42 с.

16. Даутов Р.З., Карчевский М.М. Введение в теорию метода конечных элементов: учеб. пособие. Казань: Изд-во КГУ, 2004. 239 с.

Science and Education of the Bauman MSTU, 2015, no. 03, pp. 296-312.

DOI: 10.7463/0315.0759406

Received: Revised:

26.02.2015 10.03.2015

Science^Education

of the Bauman MSTU

I SS N 1994-0408 © Bauman Moscow State Technical Unversity

Finite Element Modeling of Thermo Creep Processes Using Runge-Kutta Method

Yu.I. Dimitrienko1'*, E.A. Gubareva1, ' damit.wtu:ggmaii-com

Yu.V. Yurin1

bauman Moscow State Technical University, Moscow, Russia

Keywords: finite element method, thermo creep, Runge-Kutta method, three dimensional problem of thermo creep

Thermo creep deformations for most heat-resistant alloys, as a rule, nonlinearly depend on stresses and are practically non- reversible. Therefore, to calculate the properties of these materials the theory of plastic flow is most widely used. Finite-element computations of a stress-strain state of structures with account of thermo creep deformations up to now are performed using main commercial software, including ANSYS package. However, in most cases to solve nonlinear creep equations, one should apply explicit or implicit methods based on the Euler method of approximation of time-derivatives. The Euler method is sufficiently efficient in terms of random access memory in computations, however this method is cumbersome in computation time and does not always provide a required accuracy for creep deformation computations.

The paper offers a finite-element algorithm to solve a three-dimensional problem of thermo creep based on the Runge-Kutta finite-difference schemes of different orders with respect to time. It shows a numerical test example to solve the problem on the thermo creep of a beam under tensile loading. The computed results demonstrate that using the Runge-Kutta method with increasing accuracy order allows us to obtain a more accurate solution (with increasing accuracy order by 1 a relative error decreases, approximately, by an order too). The developed algorithm proves to be efficient enough and can be recommended for solving the more complicated problems of thermo creep of structures.

References

1. Nejad M.Z., Kashkoli M.D. Time-dependent thermo-creep analysis of rotating FGM thick-walled cylindrical pressure vessels under heat flux. International Journal of Engineering Science, 2014, vol. 82, pp. 222-237. DOI: 10.1016/j.ijengsci.2014.06.006

2. Li B., Lin J., Yao X. A novel evolutionary algorithm for determining unified creep damage constitutive equations. International Journal of Mechanical Sciences, 2002, vol.44, no. 5, pp. 987-1002. DOI: 10.1016/S0020-7403(02)00021 -8

3. Rabotnov Yu.N. Polzuchest' elementov konstruktsii [Creep of structural elements]. Moscow, Nauka Publ.,1966. 752 p. (in Russian).

4. Implicit Creep. ansys.net : a resource for ansys users : website. Available at: http://ansys.net/ansys/papers/nonlinear/conflong creep.pdf , accessed 01.02.2015.

5. Dimitrienko Iu.I., Iakovlev D.O. Comparison analysis of asymptotic theory of multilayer composite plates and three-dimentional theory of elasticity. Inzhenernyy zhurnal: nauka i innovatsii = Engineering Journal: Science and Innovation, 2013, no. 12. Available at: http://engjournal.ru/catalog/mathmodel/technic/899.html , accessed 01.02.2015. (in Russian).

6. Dimitrienko Iu.I., Gubareva E.A., Sborshchikov S.V. Asymptotic theory of constructive-orthotropic plates with two-periodic structures. Matematicheskoe modelirovanie i chislennye metody, 2014, no. 1, pp. 36-57. (in Russian).

7. Dimitrienko Yu.I., Gubareva E.A., Yakovlev D.O. Asymptotic theory of viscoelasic multilayer thin composite plates. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2014, no. 10, pp. 359-382. DOI: 10.7463/1014.0730105 (in Russian).

8. Dimitrienko Iu.I., Iakovlev N.O., Erasov V.S., Fedoniuk N.N., Sborshchikov S.V., Gubareva E.A., Krylov V.D., Grigor'ev M.M., Prozorovskii A.A. Development of a multilayer polymer composite material with discrete structural-orthotropic fillers. Kompozity i nanostruktury = Composites andNanostructures, 2014, vol. 6, no. 1, pp. 32-48. (in Russian).

9. Dimitrienko Yu.I., Gubareva E.A., Yurin Yu.V. Asymptotic theory of thermocreep for multilayer thin plates. Matematicheskoe modelirovanie i chislennye metody, 2014, no. 4, pp. 18-36. (in Russian).

10. Dimitrienko Iu.I. Mekhanika sploshnoi sredy. V 4 t. T. 4. Osnovy mekhaniki tverdykh sred [Continuum Mechanics. In 4 vols. Vol. 4. Bases of Mechanics of Solid Medium]. Moscow, Bauman MSTU Publ., 2013. 624 p. (in Russian).

11. Dimitrienko Yu.I. Tenzornoe ischislenie [Tensor Culculus]. Moscow, Vysshaya shkola Publ., 2001. 576 p. (in Russian).

12. Dimitrienko Iu.I. Mekhanika sploshnoi sredy. V 4 t. T. 1. Tenzornyi analiz [Continuum Mechanics. In 4 vols. Vol. 1. Tensor Analysis]. Moscow, Bauman MSTU Publ., 2011. 463 p. (in Russian).

13. Rektorys K. Variational Methods in Mathematics, Science, and Engineering. Reidel, 1980. 567 p. (Russ. ed.: Rektorys K. Variatsionnye metody v matematicheskoi fizike i tekhnike. Moscow, Mir Publ., 1985. 590 p.).

14. Bakhvalov N.S., Zhidkov N.P., Kobel'kov G.M. Chislennye metody [Numerical methods]. Moscow, Binom Publ., 2001, pp. 363-375. (in Russian).

15. Faleichik B.V. Odnoshagovye metody chislennogo resheniya zadachi Koshi [One-step methods for the numerical solution of the Cauchy problem]. Minsk, BSU Publ., 2010. 42 p. (in Russian).

16. Dautov R.Z., Karchevskii M.M. Vvedenie v teoriyu metoda konechnykh elementov [Introduction to the theory of finite element method]. Kazan, KSU Publ., 2004. 239 p. (in Russian).