CC BY
18
УДК 519.63+536.2
ОСОБЕННОСТИ РЕАЛИЗАЦИИ ЧИСЛЕННЫХ СХЕМ РЕШЕНИЯ ДВУМЕРНОГО УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ В ЦИЛИНДРИЧЕСКИХ КООРДИНАТАХ
А.А. Чермошенцева1, И.С. Плотникова2, М.О. Карноушенко3
1-3Камчатский государственный технический университет, г. Петропавловск-Камчатский, 683003
1-3e-mail:kafedravmkgtu@mail.ru
Рассмотрены разностные схемы численного решения уравнения теплопроводности в цилиндрических координатах. Проведен сравнительный анализ различных схем и шаблонов. Представлены условия сходимости.
Ключевые слова: разностные схемы, уравнение теплопроводности, цилиндрические координаты.
Implementations of numerical schemes for solving two-dimensional heat transfer equation in cylindrical coordinate. A.A. Chermoshentseva1, I.S. Plotnikova2, M.O. Kamoushenko3 (1-3Kamchatka State technical university, Petropavlovsk-Kamchatsky, 683003)
Difference schemes of numerical decision of the heat transfer equation in cylindrical coordinate are considered. Comparative analysis of schemes and patterns is conducted. Conditions to convergence are presented.
Key words: difference schemes, heat transfer equation, cylindrical coordinate.
Развитие экономики требует увеличения энергетических мощностей, все очевиднее становится ограниченность ресурсов традиционных видов топлива, острее проявляются экологические проблемы топливно-энергетического комплекса. Актуальным остается поиск альтернативных источников энергии и совершенствование технологий нетрадиционной
энергетики. В условиях нашего региона (Камчатка) перспективным представляется использование геотермальных месторождений, теплоноситель которых представляет собой высокотемпературную смесь воды и пара. Математическое моделирование пароводяных течений получило широкое распространение, главным образом в связи с расчетом течений в скважинах. Описание гидродинамических процессов в стволе скважины базируется на фундаментальных физических законах: сохранения массы, энергии и второго закона Ньютона применительно к особенностям движущейся среды [1-3] с учетом специфики газожидкостных сред.
При выборе начала системы координат, совпадающей с устьем скважины и оси z с ее вертикальной осью, в одномерном виде уравнение энергии, учитывающее изменение энтальпии с глубиной, записывается следующим образом:
dh dz dsn dO
— + —- +—— H—— = 0, (1)
dz dz dz dz
где h - удельная энтальпия теплоносителя; sK и sn - удельные кинетическая и потенциальная энергии; dQ - потери тепла в окружающие породы, приходящиеся на единицу массы.
Основная сложность в решении этого дифференциального уравнения связана с определением тепловых потерь в окружающие горные породы. В работах [4-9] предложены различные формулы для учета теплообмена с окружающими породами, полученные в результате решения одномерной задачи теплопроводности. В работах [10, 11] учет теплообмена с окружающей горной породой реализован посредством решения двумерной задачи теплопроводности в цилиндрических координатах. Необходимым условием распространения теплоты является неравномерность распределения температуры в рассматриваемой среде, и задача учета теплообмена сводится к необходимости знать распределение температур в массиве горных пород и на стенке скважины в каждый момент времени.
Предполагая осесимметричное распределение температуры на стенках скважины в [10] рассматривается численное решение двумерной задачи теплопроводности с заданными граничными условиями первого рода:
сТ ( 92Т 1 сТ а2Тл
= а\ —-н-------+
дт I Зг 2 г Зг дz2
где Т - температура; г - радиус-вектор; а - коэффициент температуропроводности; г - аппликата; т - время.
Точное аналитическое решение задач математической физики обычно требует интегрирования дифференциальных уравнений с частными производными, включающих искомые функции. Как правило, аналитическое решение можно найти для весьма ограниченного класса одномерных задач в областях простой формы. Использование численных методов позволяет отказаться от упрощенной трактовки математической модели процесса. Универсальность и наличие хорошо разработанной теории позволяет применять разностные методы для приближенного решения уравнения теплопроводности [12-24].
Математические модели физических процессов предполагают непрерывность распределения искомых величин в пространстве и времени. В методе конечных разностей для получения приближенного представления о пространственном распределении величин и их эволюции достаточно иметь информацию об этих значениях на некотором конечном множестве точек пространства (узлах) в фиксированные моменты времени. Очевидно, что уменьшение временных интервалов и расстояний между выбранными точками позволит приблизить такое дискретное представление к непрерывному. При этом дифференциальное уравнение заменяется эквивалентным соотношением в конечных разностях (разностная схема), решение которого сводится к выполнению алгебраических операций.
Разностные схемы отличаются выбором конфигурации узлов, участвующих в расчетах. Построенная разностная схема может быть явной, когда в каждом уравнении содержится только одно значение функции на следующем слое, и неявной, когда в каждом уравнении несколько значений функции на новом слое. Явные схемы имеют важное достоинство: они просто записываются и легко программируются на ЭВМ, при этом требуется выполнение некоторого условия, обеспечивающего сходимость. Применение неявной схемы приводит на каждом слое к системе линейных уравнений. Матрица полученной системы будет трехдиагональной, а саму систему можно решить алгебраической прогонкой. Запишем обобщенную двуслойную разностную схему на шеститочечном шаблоне:
тк:1-тк„ а >1-а
пк+1 ^г/-гк+1 , г/-гк+1 , £-*к г^грк , грк
~ ^ пг "1~ пг+1 3 7г * пг-\ ~ ^ пг "1~ * пг
* 7 2 * пг-\ пг пг+1 ? 2 * т-\ пг пг+1 '
Ат к к
Схема содержит весовой множитель при пространственной производной верхнего СЛОЯ СТ, меняя который можно добиться улучшения тех или иных свойств схемы. Так при ст = 0 двуслойная схема шеститочечного шаблона становится явной схемой треугольника (рис. 1, а):
г^к+1 г^к г^к
пг пг т-1 пг пг+1
Ат к
при ст = 1 получаем неявную схему треугольника (рис. 1,6):
гтгк+Х гт-1 к гт1 к+\ /угтгклА ^ ут£+1
пг пг пг-1 пг пг+1
а при ст = ^ — схему Кранка - Николсона (рис. 1, в):
г^к+1 г^к г^к+1 ^г^к+\ ^ г^к+1 г^к
пг пг пг—1 пг пг+1 _|_ пг—1 пг пг+1
Ат 2Г 2Г
2
h
г^к
-пк+\ І
у» к
г^к
-1 ‘і
----------Ф
грк гуік
т т+1 1
Т.
к-1 [ т
Т
к+1 -і т^1 ні Т7^1 ні
^ І ^ От от+1 ^
Т 1 Т
т-1 _1 і
грк т+11
а б в
Рис. і. Шаблоны разностных схем: а - явная схема треугольника; б - неявная схема треугольника; в - схема Кранка - Николсона
Наибольшее распространение получила схема Кранка - Николсона, позволяющая по сравнению с неявной треугольной схемой выполнять примерно в 10 раз меньше вычислительной работы для расчетов на одном и том же временном интервале [15]. К тому же она имеет более хорошую аппроксимацию и является абсолютно устойчивой.
Для двумерного случая температуру в узловой точке Тктп сопровождают три индекса: m и п -индексы координат, k - индекс времени и обобщенная схема шеститочечного шаблона имеет вид:
грк+1 ггік
т,п т,п
= + А?%Тк+1 + і-ст'ї’* ",
^ 1 2__~ т,п ^
(3)
где Л,- - разностный оператор, определяемый соотношениями:
Л =Т -2 Т +Т
1 -1 т+\,п ^ га,и ^ * га-1,и
л =тк -2 Тк +Тк
-1 т,п+1 ^ т,п 1 Л т,п-1
Устойчивость схемы обеспечивает выполнение неравенства [15]:
1 1
с >-----------
2 4ат
где к1 и к2 - пространственные шаги по каждой из переменной.
Предложенная схема (3) при ст = 0 становится явной, и значение 7 'кт"п вычисляется
непосредственно по значениям с предыдущего слоя, при а - “ получается абсолютно устойчивый вариант схемы.
Рассмотрим конечно-разностное выражение уравнения (2). Заменим все частные производные дифференциального уравнения на разностные соотношения значений функции в соседних узловых точках шаблона (рис. 2) в соответствии с формулами:
£Т
ог
с а2т^
дг2
пк ' т—1./2
получим:
1 ^к+1 грк ^ ® 2^Гк “Н Тк ~ь ^ ^
Q * т,п _______> о 2 т+1,п т,п т-\,п ^ с ;
0Т о.
а
г5„
-Т +
т+\,п т,п -
Т
к+1 ?. Шаблон для посіиу^ьігмя п,г‘тной схемы (узловые точки)
+, - 2Тк +тк , ,
^2 т,«+1 /л,«
г
где 5Г, 6- и й, - интервалы разбиения в радиальном направлении в глубину и по времени соответственно.
Приближенное значение температуры в точке /'*
определяется соотношением
Тк+1 = аЬг
'T’k
т+1,п
V5;
. 'T^k 1 . 'T^k 1 . 'T^k 1
т-1,п ^2 т,п+1 ^2 т,п-\ ^2
т
к
-1
к
тл
г
т.п
+т.
1-а5т
f 2 1 2 ^ з,2 V+sb
Такой подход дает возможность решать задачи при разнообразных краевых условиях, оценить погрешность перехода от дифференциального уравнения к уравнению в конечных разностях, а также провести анализ условий устойчивости и сходимости решения. Представленная разностная схема является явной и позволяет найти решение краевой двумерной нестационарной задачи теплопроводности в цилиндрических координатах. Условие устойчивости, представленное в [10], позволяет осуществлять выбор шагов пространственной сетки и временного интервала.
Литература
1. БэтчелорДж. Введение в динамику жидкости. - М.: Мир, 1973. - 758 с.
2. Кейс В.М. Конвективный тепло- и массообмен. - М.: Энергия, 1972. - 448 с.
3. Ландау Л.Д., Лифшиц Е.М. Гидродинамика.- М.: Наука, 1986. - 736 с.
4. Palacio-Perez A. A computer code for determining the flow characteristics in a geothermal well // Proceedings of the international conference on numerical methods of thermal problems. - Swansen, 1985. - Part 2. - P. 922-933.
5. ДядькинЮ.Д. Разработка геотермальных месторождений. - М.: Недра, 1989. - 229 с.
6. Кирюхин А.В. Математическое моделирование. - Петропавловск-Камчатский: КГАРФ, 1998. - 52 с.
7. Потапов В.В. Тепломассоперенос в фильтрационном, струйном и закрученном потоках (на примере геотермальной среды): Автореф. дис. ... канд. тех. наук. - М., 2000. - 22 с.
8. Шулюпин А.Н. Течение в геотермальной скважине: модель и эксперимент // Вулканология и сейсмология, 1991. - № 4. - С. 25-31.
9. Pedro Sanchez Upton The wellbore simulator SIMU2000 // Proceeding World Geothermal Congress 2000, Kyushu-Tohoku, Japan, 2000. - P. 2851-2856.
10. Чермошенцева А.А. Течение теплоносителя в геотермальной скважине // Математическое моделирование, 2006. - Т. 18. - № 4. - С. 61-76.
11. Шулюпин А.Н, Чермошенцева А.А. Влияние теплообмена с окружающими породами на эксплуатационные параметры пароводяной скважины // Научный журнал КубГАУ [Электронный ресурс]. - Краснодар: КубГАУ, 2007. - № 23(03).
12. Исаченко В.П., Осипова В.А., Сукомел А.С. Теплопередача. - М.: Энергоиздат, 1981. -416 с.
13.Боглаев Ю.П. Вычислительная математика и программирование. - М.: Высш. шк., 1990. -543 с.
14. Власова Е.А., Зарубин В.С., Кувыркин Г.Н. Приближенные методы математической физики. - М.: МГТУ им Н.Э. Баумана, 2001. - 700 с.
15. Калиткин Н.Н. Численные методы. - М.: Наука, 1978. - 512 с.
16. Лыков А.В. Тепломассообмен: Справочник. - М.: Энергия, 1978. - 479 с.
17.МарчукГ.И. Методы вычислительной математики. - М.: Наука, 1980. - 536 с.
18. Самарский А.А. Теория разностных схем. - М.: Наука, 1989. - 614 с.
19. Самарский А.А., Гулин А.В. Устойчивость разностных схем. - М.: Наука, 1973. - 415 с.
20. Самарский А.А., Гулин А.В. Численные методы. - М.: Наука, 1989. - 429 с.
21. Самарский А.А., Попов Ю.П. Разностные методы решения задач газовой динамики. -М.: Наука, 1992. - 422 с.
22. Созинова Т.Е. Разработка метода расчета и исследования теплового и термонапряженного состояния крепи геотермальных скважин: Автореф. дис. ... канд. тех. наук. - Иваново, 1997. - 24 с.
23. Тихонов А.Н., Самарский А.А. Уравнения математической физики. - М.: Наука, 1966. - 724
с.
24. Чермошенцева А.А. Численные методы: Учеб. пособие. - Петропавловск-Камчатский: КамчатГТУ, 2008. - 110 с.
