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

УДК 536.24

Математическое моделирование температурного состояния пространственных стержневых конструкций.

Нестационарные и нелинейные задачи

© И.В. Станкевич

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

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

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

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

Ф^Т Слх) = 1 (Т Слх),л^Л + IV (Л), Слх) е (Лъ Л2) х (0, ; (1)

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

дТ д х

Т (л, 0) = Тс (л), л е [Л1, Л2]; (2)

1Т (Л, х),л1л=т = чш, X ^ 0; (3)

1Т (Л, х),л1Л=Л2 = (Г/ - Т (л, х))|Л=Л2, х ^ 0, (4) где с, р, 1 — удельная теплоемкость, плотность и коэффициент теплопроводности материала стержня соответственно; х — время; Т (л, х) — температура стержня; qy (л) — мощность внутренних источников (стоков) теплоты; Тс (л) — начальная температура стержня (предполагается, что начальное и граничные условия согласованы); — численное значение плотности теплового потока на поверхности Б2, в данном случае принято, что > 0, если теплота отводится от поверхности в2 стержня (Т (л, х),Л| > 0 Ух ^ 0); л1, Л2 — координаты торцевых поверхностей Я2 и стержня соответственно (л1 < Л2); и Tf — коэффициент теплоотдачи и температура внешней среды вблизи поверхности 5з соответственно.

Матричные соотношения МКЭ. При решении нестационарных задач теплопроводности конечно-элементную дискретизацию по пространству целесообразно выполнять на основе метода взвешенных невязок в форме Галёркина [2]. Каждому узлу р (в глобальной нумерации [1, 2]) поставим в соответствие финитную функцию Мр (р = 1, ки, ки — глобальное число узлов сетки конечно-элементной модели). Носителями функций Мр являются объединения областей у(е) тех конечных элементов (е), которые содержат данный узел р. Кроме того, предположим что функции Мр линейно независимы. Тогда температуру и ее производные по координате л и времени х на компакте ¿¡^ можно интерполировать с помощью следующих соотношений

Т = [Ж] {Т} ; (5)

= [^,Л ] {т}; (6)

т

[Л] {т} . (7)

Здесь Т, Т, , Т — интерполированные значения температуры, ее производных по криволинейной координате ^ и времени т соответственно; [Щ — матрица-строка, составленная из финитных функций Мр (р = 1, кц); {Т} — вектор-столбец, составленный из узловых значений температуры Тр (р = 1, ки); |т| — вектор-столбец, составленный из узловых значений производных температуры по времени Тр

(р = 1, ки).

Перенесем все члены уравнения (1) в правую часть, подставим в нее интерполированные значения температуры, ее производных и рассмотрим невязку

К = 1 Т

(т^ ,л + ду — срт.

Теперь в соответствии с методом Галёркина умножим невязку Я последовательно на функции Мр (р = 1, кц), проинтегрируем по области Си и результат приравняем нулю, получим

=

¿Он

1 (т,,ц + Яу — срТ Мр(ц)^ = 0, р =1,ки.

щ!

(8)

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

Перепишем выражение (8) в виде

(т,Л ,л + ду] Ирйю — срТИрйю = 0, р = 1,кц. (9)

¿Он 1 4 ' ^ Зсн

С учетом выражения (7) второй интеграл в формуле (9) представим следующим образом:

[ срТ= [ ср [Ж] Ыр(1у{т\ =

¿Он ¿Си 1 >

= [ ср [МгМр М2МР ...Мки Мр] ¿у\т} , р =1^. (10)

Рассмотрим глобальную матрицу теплоемкости [С], которую на основании выражения (10) можно записать в виде

[С]= [ ср [Ж]т [М] ¿V. (11)

¿Он

Выражение (11) для вычислений и формирования системы линейных алгебраических уравнений является крайне неудобным. Интегрирование следует проводить по объемам конечных элементов, используя

нормированную локальную систему координат 0'%, а результат суммировать. Для этого вместо глобальной матрицы-строки [Щ исполь-

что позволит формирование

V

зуем локальные матрицы-строки глобальной матрицы теплоемкости выполнить по формуле

N1

( )

[С ]

ку /г

§ X „

с(е)р(е)

ж,

( )

У

ж,

( )

У

(е)'

V

(12)

Здесь

( )

— матрица геометрических связей объемного конечного элемента с идентификационной меткой (е) [1,2].

Обозначим интегралы, входящие в выражения (12), следующим образом:

['(?] =

где элемент объема

'у (в)

с(е)р(е)

(е) У

м:

( )

У

¿V,

(13)

¿V = А(е) (X) ¿Л = А(е) (X)

\

г=1

Здесь А(е) (X) — площадь поперечного сечения конечного элемента (е), а производные хг%, г = 1, 3, определяются формулами (19), которые приведены в работе [1].

Площадь поперечного сечения А(е) (X) можно аппроксимировать с помощью функций формы:

А(е)(Х)

К

( )

У

{А(е)}

(14)

где |А(е)} — вектор, составленный из площадей поперечных сечений (^ = 1, Р(е), отнесенных к узлам элемента (е); р(е) — число узлов конечного элемента).

Если материал стержня неоднородный, аналогично можно представить удельную теплоемкость с(е) (X) и плотность р(е)(Х) материала стержня в пределах конечного элемента (е):

с(е) (X) P(e)(X)

м:

( )

У

м:

( )

У

{р( е)}

(15)

(16)

Здесь {с(е)} и {р(е)} — векторы, составленные из значений удельной

теплоемкости ске) и плотности р ке), к =1, р(е), соответственно, отнесенные к узлам к =1, р(е) элемента (е) (р(е) — число узлов конечного элемента).

( )

т

т

Квадратурная формула для вычисления интегралов (13) по объему конечного элемента (е) с учетом соотношений (14)-(16) принимает вид

К?] = Е °{е)(х) р(е)(х)А(е) (х)

к

(е)'

V

К

(е)'

V

1=1

\

=1

н,.

х=х I

где п — число гауссовых точек; I — номер гауссовой точки; X — локальная координата гауссовой точки; Щ — весовой коэффициент квадратурной формулы.

Аналогично можно определить первый интеграл в формуле (9), при преобразовании дивергентной части используют вторую формулу Грина. В результате получены соотношения, полностью совпадающие с формулами (14) и (15), приведенными в работе [1]. Это подтверждает известное утверждение о том, что процедуры метода Ритца и метода Галёркина при решении задач с самосопряженными и положительно определенными операторами приводят к одинаковым по конструкции системам алгебраических уравнений [3].

Таким образом, нестационарная задача (1)-(5) после дискретизации по пространству сводится к решению задачи Коши для линейного матричного дифференциального уравнения первого порядка

[С]{т} + [К] {Т} = {Щ

(17)

с начальным условием

{Т}|т=о = {То} , (18)

где {Т0} — проекция функции Т0 на узлы сетки конечно-элементной модели.

Для ее решения существуют различные подходы [2-4], однако наибольшее распространение получили два. Первый состоит в том, что производную по времени в уравнении (17) заменяют каким-либо конечно-разностным аналогом, а второй заключается в использовании конечных элементов во временной области (метод Галёркина).

Разностный аналог задачи Коши. Рассмотрим основные этапы построения разностного аналога задачи Коши (17), (18) в виде семейства двухслойных разностных схем. В пределах временного шага К = тп — тп-1 векторы узловых температур {Т} и правой части {Д} представим в виде следующих линейных комбинаций:

{Т(т)} = (1 — ш) {Т} п-1 + ш {Т}п ; (19)

{Н(т)} = (1 — ш) {П}п-1 + ш {П}п , (20)

где т е [тп-1,тп] С [0,^; ш — весовой множитель, ш е [0,1]. Здесь

и далее векторы {Т}п-1, {Я}п-1 и {т}п, {я}г времени т = тп-1 и т = тп соответственно.

отнесены к моментам

т

2

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

Т \ = „ {Г}п - {Г}^-1 (21)

Т/= 8х ~ К . (21)

Полагая, что коэффициенты в уравнении (17) постоянны на отрезке [хп-1,хп], и подставляя выражения (19), (20) и (21) в формулу (17), после очевидных преобразований получаем общее выражение для двухслойной схемы с весами:

([С] + КхШ [К]) {Т}га = = ([С] + Мш -1) [К]) {Т}га-1 + Кх ((1 - ш) {Щп-1 + ш {Щп) .

Как известно, фиксированное численное значение параметра ш определяет тип конкретной разностной схемы [4], например: ш = 0 — схема с разностью вперед; ш = 1/2 — схема Кранка — Николсона; ш = 2/3 — схема Галёркина; ш =1 — схема с разностью назад. Вообще говоря, параметр ш может принимать любые численные значения из отрезка [0, 1]. Схема Кранка — Николсона имеет определенные преимущества перед остальными схемами, поскольку аппроксимирует уравнение (17) по переменной х с порядком О (К2), а остальные указанные схемы имеют более низкий порядок аппроксимации — 0(кх). Кроме двухслойных схем существуют и используются трехслойные схемы, построение которых можно найти в работах [2, 4].

Решение нелинейных задач теплопроводности. Нелинейность в задачах теплопроводности возникает, когда коэффициенты в уравнении (1) и граничных условиях (3) и (4) зависят от искомой температуры, например 1 (Т), с (Т), р (Т), qy (л, Т), (Т) и а№ (Т). Предположим, что все эти функции являются измеримыми и ограниченными и, кроме того, имеют ограниченные производные по температуре Т. Тогда, если не применяются линеаризующие процедуры, на каждом шаге по времени необходимо решать систему нелинейных алгебраических уравнений с помощью итерационных методов [4, 5]. Во избежание этого применяют схемы типа предиктор — корректор [6], для которых на каждом временном шаге требуется решать две системы линейных алгебраических уравнений. Существенными недостатками использования схем предиктор — корректор являются: общее усложнение алгоритма решения и дополнительные затраты оперативной памяти. Эти трудности можно исключить, если начально-краевую задачу в каждый момент времени решать методом простых итераций с явным заданием скорректированных значений коэффициентов уравнения теплопроводности и граничных условий, применяя метод Галёркина для построения матричных соотношений МКЭ. Таким образом, в каждой точке временного отрезка решение нелинейной начально-краевой задачи

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

Применение итераций при решении нелинейных уравнений эллиптического и параболического типов является известным и достаточно широко используемым методом [6-12]. Наиболее полные результаты получены для эллиптических уравнений. Особую проблему при использовании итерационного решения составляет сходимость. При анализе сходимости, как правило, рассматривают обобщенное решение в подходящем функциональном классе, а в отдельных случаях анализ проводят на конечномерных подпространствах. Для параболических уравнений характерно рассмотрение сходимости итераций на конечномерных подпространствах, построенных либо с помощью метода Галёркина, либо некоторым разностным методом. Одним из подходов к анализу сходимости последовательных приближений (итераций) является следующий: пусть исходный нелинейный дифференциальный оператор удовлетворяет двум неравенствам, одно из которых является некоторым вариантом условия монотонности, а второе — ограниченной нелинейности в некоторой области; тогда дискретный аналог с помощью процедур какого-либо разностного метода строят так, чтобы сохранились эти два свойства. После этого теоретические исследования сходимости итераций в рамках дискретного аналога проводят в соответствии с некоторыми основными положениями теории монотонных операторов [11-14]. Кроме того, при анализе сходимости итераций как на дифференциальном уровне, так и на конечномерных подпространствах, если это возможно, используют дифференциальные свойства нелинейных операторов — дифференцируемость по Фреше или, в крайнем случае, по Гато.

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

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

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

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

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

На рис. 2 показана исследуемая стержневая конструкция. В качестве конечных элементов выбраны квадратичные трехузловые элементы. Кроме того, на рис. 2 приведена глобальная нумерация узлов соответствующей конечно-элементной модели. Конвективный теплообмен осуществлялся на боковых поверхностях двух конечных элементов е\ и е2: е\ — узлы 24, 25, 26; е2 — узлы 11, 12, 13. В зоне е\ 7} = 900К и аш = 135Вт/(м2 ■ К), а в зоне е2 = 500К и а№ = 60Вт/(м2 ■ К). Начальная температура принята равной 300К. Результаты расчетов представлены на рис. 3-5. Как следует из рисунков, стержневая конструкция имеет неравномерное распределение температуры. Заметно интенсивнее нагрев осуществляется в зоне узла 25, менее нагретыми являются удаленные участки (зоны узлов 3, 9). Кроме того, усиливается влияние зоны наиболее интенсивного нагрева на зону менее интенсивного нагрева. На рис. 5 представлены температурные поля в фиксированные моменты времени, показывающие кинетику роста общего температурного состояния стержневой системы.

к

Рис. 2. Глобальная нумерация узлов конечно-элементной модели стержневой конструкции

Т,К 700 600 500 400 300

Рис. 3. Изменение температуры узлов конечно-элементной модели: 25 (1), 31 (2),

12 (3), 15 (4)

Г,К 1----

600 500 400 300

0 200 400 600 800 т,с

Рис. 4. Изменение температуры узлов конечно-элементной модели: 20 (1), 18 (2),

5 (3), 3 (4), 9 (5)

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

Работа выполнена при финансовой поддержке гранта Президента РФ для государственной поддержки ведущих научных школ (проект НШ-255.2012.8).

Рис. 5. Температурные поля стержневой конструкции в фиксированные моменты времени: 500 с (а), 1 000 с (б), 3 600 с (в) (температура стержня — в кельвинах)

ЛИТЕРАТУРА

[1] Станкевич И.В. Математическое моделирование температурного состояния пространственных стержневых конструкций. Стационарные задачи. Инженерный журнал: наука и инновации, 2013, вып. 8.

URL: http://engiournal.ru/catalGg/mathmodel/technic/893.html

[2] Котович А.В., Станкевич И.В. Решение задач теплопроводности методом конечных элементов. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2010, 84 с.

[3] Марчук Г.И., Агошков В.И. Введение в проекционно-сеточные методы. Москва, Наука, 1981, 416 с.

[4] Галанин М.П., Савенков Е.Б. Методы численного анализа математических моделей. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2010, 591 с.

[5] Ортега Дж., Рейнболдт В. Итерационные методы решения нелинейных систем уравнений со многими неизвестными. Москва, Мир, 1975, 558 с.

[6] Сергиенко И.В., Скопецкий В.В., Дейнека В.С. Математическое моделирование и исследование процессов в неоднородных средах. Киев, Наукова думка, 1991, 432 с.

[7] Кошелев А.И. О сходимости метода последовательных приближений для квазилинейных эллиптических уравнений. Доклады Академии наук СССР, 1962, т. 142, № 5, с. 1007-1010.

[8] Шейбак Т. Построение итерационного процесса для решения квазилинейного эллиптического уравнения с разрывными коэффициентами. Дифференциальные уравнения и их применение, 1985, № 38, с. 61-67.

[9] Арделян Н.В. О сходимости итерационных методов решения нелинейных разностных схем для нелинейного уравнения теплопроводности. Дифференциальные уравнения, 1985, т. XXI, № 12, с. 2131-2137.

[10] Jordan A. Iterative Method of the Analysis of Nonlinear Heat Transfer Problems. Scientific Journal Bialystok University of Technology. Technical Sciences. Electricity, 1992, vol. 83, no. 11, pp. 53-60.

[11] Качуровский Р.И. Нелинейные монотонные операторы в банаховых пространствах. Успехи математических наук, 1968, т. 23, № 2 (140), с. 121-168.

[12] Гаевский Х., Грeгер К., Захариас К. Нелинейные операторные уравнения и операторные дифференциальные уравнения. Москва, Мир, 1978, 336 с.

[13] Станкевич И.В. Сходимость метода простых итераций при решении нелинейных стационарных уравнений теплопроводности. Вестник МГТУ им. Н.Э. Баумана. Машиностроение, 1995, № 3, с. 97-102.

[14] Треногин В.А. Функциональный анализ. Москва, Наука, 1993, 440 с.

Статья поступила в редакцию 20.06.2013

Ссылку на эту статью просим оформлять следующим образом:

Станкевич И.В. Математическое моделирование температурного состояния пространственных стержневых конструкций. Нестационарные и нелинейные задачи. Инженерный журнал: наука и инновации, 2013, вып. 8. URL: http://engjournal.ru/catalog/mathmodel/technic/894.html

Станкевич Игорь Васильевич — д-р техн. наук, проф. кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана. e-mail: aplmex@yandex.ru