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

УДК 519.6

А.Ю. Вишняков, Ю.Н. Дерюгин, В. А. Глазунов, И. Н. Чистякова

Российский федеральный ядерный центр — всероссийский НИИ экспериментальной физики

Пакет программ ЛОГОС. Модуль расчета сопряженных и связанных задач теплопереноса

В работе приводится описание методики и счетного модуля, созданного в рамках пакета программ ЛОГОС [1] для расчета задач теплопроводности в твердотельных конструкциях. Расчетные методики построены на основе метода конечного объема и неявной с весами по нижнему и верхнему слою аппроксимации потока тепла через грани ячеек. Приводятся два алгоритма для расчета теплового потока через грани ячеек: метод отложенной коррекции [2], градиентный метод [3]. Вычислительные алгоритмы основаны на записи разностных уравнений в дельта-форме [4]. Эти алгоритмы ориентированы на решение нестационарных и стационарных задач теплопереноса в твердотельных конструкциях с учетом как изотропной, так и анизотропной теплопроводности на неструктурированных сетках. При записи разностных уравнений в дельта-форме правые части уравнений, которые аппроксимируются явно с предыдущей итерации или с нижнего слоя, являются невязками балансных уравнений (погрешностью аппроксимации). Для решения системы разностных уравнений, получающихся при неявной аппроксимации, используются решатели библиотеки РМЬР [5].

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

1. Основные уравнения

Дифференциальное уравнение нестационарной теплопроводности в неоднородной неподвижной среде [6] имеет вид

Здесь р = р (х,у,г,Т) — плотность вещества, Е = Е (х,у,г,Т) — внутренняя энергия, %=% (х,у, г,Т) — коэффициент теплопроводности. В изотропных материалах коэффициент теплопроводности - скаляр, при расчёте теплопроводности в анизотропных материалах — тензор второго ранга.

Частным случаем задачи теплопроводности является задача определения установившегося теплового поля в некоторой области при неизменных граничных условиях. Распределение теплового поля описывается стационарным уравнением теплопроводности:

Для однозначного определения температуры ставятся начальные и граничные условия. В начальный момент времени задается распределение температуры в области:

На границах цилиндра ^ = {(х, у, z, t) : (х, y,z,) € dQ, 0 < t < iend} могут быть поставлены следующие граничные условия:

• задана температура:

^ = div (х • gradT) + pQT.

(1)

div (х • gradT) = -pQT.

задан тепловой поток:

= гс (*0 • е Б (3)

задан конвективный теплообмен:

дТ \

}ш) =ЛТ-те(*0), <4>

Предполагается, что граница области ^ состоит из непересекающихся частей ^ = У На каждой части границы а) может быть поставлено одно из трех гра-

а

ничных условий: (2) - (4).

При решении сопряженных задач теплообмена излучением в прозрачных средах граничное условие на излучающих поверхностях записывается в виде

- )„.="(й'.(). .Ое д®-

■з ,ь) е ^ ^3 ,

где поток д ^ определяется как разность потоков собственного и эффективного излучения [7].

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

1 дТ

[Т] = -Хдт, (х,у, г) ее,

а дп

где а — коэффициент контактного теплообмена.

В задачах плавления и затвердевания (задача Стефана) фазовый переход сопровождается выделением либо поглощением определенного количества тепла. Задача Стефана допускает обобщенную формулировку [8]:

дТ

(рСР + рНё (Т -Т*)) — = ё1у (х ■ gradт) + рЯт,

где Т* — температура фазового перехода, Н — удельная энтальпия фазового перехода, ( Т - Т* ) — дельта-функция Дирака. В задачах нагрева и воспламенения взрывчатых веществ функция внутреннего источника теплоты задается в виде

д С

Ят = ^вв^тт. д

где Явв — теплота разложения ВВ, а С — массовая концентрация ВВ. Кинетика разложения ВВ описывается по закону Аррениуса [9]:

д^=- - {-&}с-

где Е* — энергия активации, К — скорость разложения ВВ (частотный множитель), « — газовая постоянная.

2. Интегральная форма уравнения теплопроводности

Проинтегрировав дифференциальные уравнения (1) по контрольному объёму и применив далее квадратурные формулы и теорему о среднем, получим полудискретный вид уравнения теплопроводности:

(рДУ^^ - £ ^гаёТ ■ п)} Д^ = (р^ДУ) р .

Для описания численных алгоритмов обозначим через Л(Т) разностный оператор, аппроксимирующий потоки тепла через грани ячейки. С учетом введенного обозначения дискретизация нестационарного уравнения теплопроводности запишется в виде

(рАУ)р (В|) - Л(Т) = (р^АУ)р.

Соответственно, разностная аппроксимация стационарного уравнения теплопроводности в операторной форме будет иметь вид

-Л(Т) = (рЯДУ) р.

3. Аппроксимация производной по времени

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

-4Г1+ и

п I гтП— 1

^ (рД У)р (Е^1 - ЕП) - А(Тп+1) =

= (врду )р+^Т^) (еп -ЕП"1)

где параметр и может принимать значения: 0 ^ и ^ 1. Для схемы первого порядка аппроксимации: и = 0. Для схемы второго порядка аппроксимации: и = 1.

4. Дельта-форма разностных уравнений

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

Проведя линеаризацию потоков тепла и приращения внутренней энергии, получим следующее:

Е1+1 -Еп = Е1 +

/ ВЕ\1 / ВЕ\1

\Вт) (Т7+1 - Т7) -Еп =(^ ДТ + (Е7 -Еп) ,

где ДТ = Т7+1 - Т7.

Линеаризация потоков тепла приводит к выражению:

Л(Т7+1) = Л(Т7+1) - Л(Т7) + Л(Т7) = ^ВЛ) ДТ + Л(Т7).

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

является невязкой, то есть погрешностью разностной аппроксимации уравнения баланса энергии.

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

где ш — релаксационный параметр.

5. Решение разностных уравнений

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

где А = (арр) — разреженная матрица размерности (Ы х N), ДТ = (ДТр) — искомый вектор размерности ( N), В = (Ьр) — вектор правых частей размерности (Ы).

6. Метод отложенной коррекции

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

Данный метод аппроксимации потока тепла на грани является основным в модуле расчета теплопроводности в пакете программ ЛОГОС.

7. Аппроксимация потока тепла для изотропной теплопроводности

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

Для построения аппроксимации потока тепла через внутреннюю грань рассмотрим произвольную грань /, которая является общей для ячеек Р и ^ (рис. 1). Обозначим индексом / центр грани. Проведем через центры Р иF ячеек плоскости, перпендикулярные нормали Пf к грани /. Определим точки Р' и Р' пересечения этих плоскостей с нормальной линией к грани f.

АД Т = В,

дТ

Рис. 1. Шаблон точек для аппроксимации х^— на внутренней грани

дп

Считая, что в точках Р', Р' ив центре грани известны значения температуры, поток тепла в направлении нормали в случае идеального контакта можно определить из условия непрерывности потока:

Тр/ — Тр'

Ту — Тр/ Тр/ — Ту

= ХР- - = ХР-

Д «рЧ + Д « Д « р// Д «р< у

Коэффициент теплопроводности на грани определяется как среднее гармоническое значение.

Значения Тр/ и Тр/ определяются с использованием gradТ, вычисленного в центрах ячеек:

ТР<) = ТР(Р) + )Р(р) ■ д«рр>)) ■

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

Х^Щ) Д5/ ~ {Тр -Тр + (^Т)^ ■ Д«^) - (^Т)р ■ Д«рр^),

где коэффициент Ху для идеального контакта:

ХРХР

Х./

и для неидеального контакта:

XI

Д «

РЧ

Хр +

Д «

■РЧ

-Д5/

ХР

хрхР

Д «

РЧ

ХР +

Д «

РЧ

ХР +

хрхР а

8. Аппроксимация потока тепла для анизотропной теплопроводности

В случае анизотропной среды поток тепла через единицу поверхности в направлении нормали определяется как

( (хртаЛГ) ■ п)/ = ^ ■ п^ = ^ (^Х^ щ^ = (gradТ ■ N)f .

где вектор N имеет следующие компоненты:

N = х ■ п =

£

и=1

3

3

Хг1 ■пг^^Хг2 ■ Щ.^ХъЗ =1 =1

п

Представим вектор N в следующем виде:

N = % • п к.

Здесь % =

N

, а Пк = Щ/

N

- единичный вектор. С учетом введенных обозначений поток

тепла будет выражаться следующим образом: ( (%£гаёТ) • п) / =

£ ЙН ^), = (%к э

/

Окончательное выражение для теплового потока через внутреннюю грань записывается аналогично изотропному случаю.

9. Градиентный метод

Этот метод ориентирован на расчет как изотропной, так и анизотропной теплопроводности. В этом методе определение теплового потока в явном операторе состоит из определения градиента температуры в узлах разностной сетки и проецирования значения градиента из узлов в центр грани. Для определения градиента температуры в узлах строятся среднеме-диальные контрольные объемы относительно узлов разностной сетки. Среднемедиальный контрольный объем А У относительно внутреннего узла г строится таким образом, что геометрические центры ячеек сетки Р^, у которой одна из вершин совпадает с узлом г, соединяются друг с другом через середины разделяющих их граней к и середины ребер Кг (рис. 2). Грани среднемедиального контрольного объема являются треугольниками.

Рис. 2. Часть контрольного объема относительно внутреннего узла г сетки со стороны ячейки Р

Определим в точке Щ (центр ячейки Р)) вектор:

к

Щ = Щ, ^ )г] = ^(nx, пу, пх )г]к АЗг]к,

к=1

где через А Бг ¿к обозначим площадь треугольной грани к контрольного объема узла г. Значение градиента температуры в узлах определяется по формуле Грина-Гаусса:

1 ТпйБ

^гаёТ)г =

4 _ Е^ , N)грТР

III

У

<1У Ауг

где суммирование ведется по ячейкам Р), у которых узел сетки с номером г является узлом контрольного объема. Объем Д Уг определяется как сумма объемов тетраэдров, у которых одна из вершин совпадает с центром ячейки Р^, вторая вершина является узлом сетки г, третья точка — центр грани Д и четвертая точка — центр ребра Ке.

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

I РСе11

ЕЕ (^^У )г РТР Рс* = г=1 р=1—-1-= N, N2 )рТр ■

Е Д * Р =1 =1

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

10. Тестовые расчеты. Нестационарное нагревание плоской пластины

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

Рассматривается нестационарное нагревание плоской пластины с заданными граничными условиями (ГУ) I и II рода. На одном торце пластины задается условие конвективного теплообмена с коэффициентом теплоотдачи а и с таблично заданной температурой окружающей среды Т^ в зависимости от времени, на остальных поверхностях пластины задаётся условие адиабатичности (рис. 3-4). Данная задача допускает аналитическое решение.

Рис. 3. Однородная плоская стенка

300 -I-1-1-1-1-1-1

0 0.02 0.04 0.06 0.08 0.1 0.12

Расстояние, м

Рис. 4. Сравнение точного и численного решений на момент времени £ = 3600 с

11. Тестовые расчеты. Моделирование тепломассопереноса и переноса излучения

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

Рис. 5. Геометрия расчетной области

Сравнение полученного численного решения с аналитическим решением показано на рис. 6.

Рис. 6. Сравнение точного и численного решений (справа— температура во внутренней сфере, слева - температура в сферической оболочке)

12. Тестовые расчеты. Инициирование детонации внутри дюймового шара из гексагена

Задача состоит в расчете изменения температуры в точке внутри дюймового шара из гексогена. В численном расчёте в силу симметрии задачи рассматривалась половина сферы. В [9] приведены результаты экспериментов для различных значений температуры на внешней границе. Здесь приведены результаты численного расчёта для граничного значения температуры Т1 = 463 К (190°С). На поверхности шара задавалась температура Т1 = 463 К (190°С). Начальная температура То = 298 К (25 °С).

Рис. 7. Изменение температуры со временем в двух точках внутри шара из гексагена (КОХ) диаметром 1 дюйм при начальной температуре 25 °С для различных значений Т

Представлены профили температур из статьи Мейдера при различных начальных температурах Т (рис. 7). В прямоугольник выделен профиль температуры, интересующий нас в данном расчете. Температура отслеживалась в точке, отношение расстояния которой до центра шара к его радиусу равно 0.9, что соответствует штриховой кривой на рис. 7. Сплошная кривая на данном рисунке соответствует изменению температуры в центре шара.

Приводится результат расчета с сенсором, помещенным в точку, отношение расстояния которой до центра шара к его радиусу равно 0.9 (рис. 8). Из этого графика видно, что инициирование ВВ происходит на момент времени ~ 700 с, что соответствует прогнозируемому значению времени согласно [9].

Рис. 8. Изменение температуры на сенсоре

Литература

1. Козелков А.С., Дерюгин Ю.Н., Зеленский Д.К., Глазунов В.А., Голубев А.А., Денисова О.В., Лашкин С.В., Жучков Р.Н., Тарасова Н.В., Сизова М.А. Многофункциональный пакет программ ЛОГОС для расчета задач гидродинамики и тепломассопереноса на многопроцессорных ЭВМ: базовые технологии и алгоритмы // Супервычисления и математическое моделирование: Труды XII международного семинара. — 2010. — С. 215— 230.

2. Schafer M. Computational Engineering — Introduction to Numerical Methods. — Berlin : Springer, 1996.

3. Панов А.И. Методика решения уравнения теплопроводности на неструктурированной сетке // Вопросы атомной науки и техники. Сер. Методики и программы численного решения задач математической физики. - 2004. - Вып. 4. — С. 27-40.

4. Beam R.M., Warming R.F. An implicit finite-difference algorithm for hyperbolic systems in conservation-low form //J. Comput. Phys. — 1976. — V. 22, N 1. — P. 87-110.

5. Артемьев А.Ю., Бартенев Ю.Г., Басалов В.Г., Бондаренко Ю.А., Варгин А.М., Голубев А.А., Ерзунов В.А., Ломтев А.В., Максимов А.С., Панов А.И., Прокофьев А.И., Романова М.Д., Фролова Н.В., Щаникова Е.Б. Библиотека решателей разреженных систем // Труды ВНИИЭФ. Математическое моделирование физических процессов. — 2004.

6. Самарский А.А., Вабищевич П.Н. Вычислительная теплопередача. — М.: ЛИБРОКОМ, 2009.

7. Оцисик М.Н. Сложный теплообмен. — М.: Мир, 1976.

8. Самарский А.А., Моисеенко Б.Д. Экономичная схема сквозного счета для многомерной задачи Стефана. Т. 5 // ЖВММФ. — 1965. — № 5. — С. 816-827.

9. Мейдер Ч. Численное моделирование детонации. — М.: Мир, 1985.

Поступила в редакцию 03.02.2013.