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

ЧИСЛЕННЫЕ СХЕМЫ РЕШЕНИЯ ДВУМЕРНОГО УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ В ЦИЛИНДРИЧЕСКИХ КООРДИНАТАХ

А.А. Чермошенцева1, И.С. Плотникова2

12Камчатский государственный технический университет, г. Петропавловск-Камчатский, 683003

e-mail: allachermoshentseva@mail. ru

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

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

Numeric schemes of solution of two-dimensional equation of thermal conductivity in cylindrival coordinates. A.A. Chermoshentseva1, I.S. Plotnikova2 (u 2 Kamchatka State Technical University, Petropavlovsk-Kamchatskiy, Russia, 683003)

Mathematical problem of heat transfer in cylindrical coordinates was presented. Initial and border conditions were described. The necessity of problem solving for using the numerical methods was proved. Various subtractive schemes, patterns, questions of accuracy and convergence were considered.

Key words: two-dimensional equation of thermal conductivity, heat exchange of the well and surroundings, subtractive schemes.

Задача теплопроводности является одной из типичных нестационарных задач математической физики. Разведка и разработка геотермальных месторождений требует решения такой задачи при рассмотрении теплообмена скважины, окруженной практически неограниченным природным массивом горных пород, с потоком движущегося теплоносителя. Учитывая конструкцию скважины, целесообразно ввести цилиндрические координаты. Тогда для скважины с осесимметричным распределением температур, не имеющей внутренних источников теплоты, задача теплообмена с окружающей горной породой при моделировании пароводяного потока в стволе геотермальной скважины сводится к решению двумерной задачи теплопроводности [1, 2]:

ST

— = а

дх

о Т 1 <ЭТ 3 Т

—г ^-----------------------1-г

д г г dr dz

(1)

где Т - температура; т - время; г - радиус-вектор; г - аппликата;

а - коэффициент температуропроводности, определяемый по формуле

а = — , (2)

ср

где X - коэффициент теплопроводности; р -плотность;

с - удельная теплоемкость (су = ср = с).

Чтобы дать полное математическое описание рассматриваемому процессу, к дифференциальному уравнению необходимо еще добавить краевые условия. Начальные условия необходимы для нестационарных задач и представляют собой задание закона распределения температуры массива горных пород (температурное поле) в начальный момент времени. Граничные условия могут быть заданы несколькими способами. Граничные условия первого рода задают распределение температуры на стенке скважины для каждого момента времени, второго рода - значения плотности теплового потока на стенке скважины для любого момента времени, и граничные условия третьего рода задают температуру массива окружающих горных пород и закон теплообмена между скважиной и породой. Отметим также, что в [3] приведены граничные

условия четвертого рода, отображающие нагрев (охлаждение) системы тел, находящихся в соприкосновении (идеальный тепловой контакт).

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

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

нелинейных задач. И тогда основным способом решения являются численные методы.

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

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

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

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

Рассмотрим наиболее распространенные разностные схемы [5-8, 10, 14]:

1) явная схема треугольника (рис. 1, а):

грк+1 грк грк Г)ГГ^ _|_ Т7^

т т _ т-1 т т+1

где Т/ - температура в і-й точке (в узле) в ]-й момент времени;

2) схема Ричардсона (рис. 1, б):

грк+1 грк-1 грк ОТ7^ _і_ Т7^

т____ т т-1 т т +1 .

3) неявная схема треугольника (рис. 1, в):

грк+1 грк грк+1 '~\грк+1 . грк+1

т______ т _ т-1 т т +1 .

4) схема бегущего счета (рис. 1, г и 1, д):

тк —Тк~1 Тк — Тк + Тк~х — Тк~х

т т т—1 т т+1 т

2

к

рк+1 рк

ГТік ГТік ГТ-гк + І ГТ-гк+І

т—1 т т+1 т

схема Кранка - Николсона (рис. 1, е):

рк+1 рк

грк+1 ‘~угтгк+1 . грк+1 -* тп-\ ~ ^ тп -* т+1

грк гт'к . грк

* тп-\ ~ ^ т "І~ т+1

Лт 2/г2 2/г2 '

По шеститочечному шаблону можно записать более общую двуслойную схему:

1-а

грк+1 грк

т______ т ____ рк+1 __ ^рк+І . рк+1

* 7 2 т-1 т т +1

Ах к

/г2

Г , - 2Т + Т7*

т-1 т т+1

Схема содержит весовой множитель при пространственной производной верхнего слоя ст. Меняя ст, можно добиться улучшения тех или иных свойств схемы. Так при ст = 0 двуслойная схема шеститочечного шаблона становится явной схемой треугольника, при ст = 1 получаем

неявную схему треугольника, а при ст = - схему Кранка - Николсона.

Т ‘

*+і и

рк+1

грк •

т-1 -1

Т'

’к к

грк

1 т+1 -1

рк

рк К. грк

г^к

рк

1-І

рк

грк т+1 1-І

1-1 .1

2

к

б

а

в

Т1

* к

т |_

рк

гтік-І 1 т+1 1

Т

к+1 -1

рк

1-І

рк

д

Т

к+1 -1

ш+1 -1

Т

*+1 -1 т-1 ]_

Ф-------

Т

-1-1

уі*+1 Ц „І+1 «І ... . 1 .1

Т'

к к

Тк,

ш+1 і

Рис. 1 Шаблоны разностных схем: а - явная схема треугольника; б - схема Ричардсона; в - неявная схема треугольника; г, д - схема бегущего счета; е - схема Кранка - Николсона

г

е

Наибольшее распространение получила схема Кранка - Николсона, позволяющая по сравнению с неявной треугольной схемой выполнять примерно в 10 раз меньше вычислительной работы для расчетов на одном и том же временном интервале [8]. К тому же она имеет более хорошую аппроксимацию и является абсолютно устойчивой.

Обобщенная схема шеститочечного шаблона для двумерного случая принимает вид:

рк+1 рк

т'п т>п Л і Л Л.ТТ^+1 і 1 Т7^

-----—-----= Л1+Л2 сТт п + 1 — ст Ттп ,

где Лг - разностный оператор:

Л =Тк -2Тк +Тк

1 т+1, п т, п т-1, п

Л, =Тк , -2Тк +Тк

2 т, п+1 т, п т, и-1

Температуру в узловой точке Ткт п сопровождают три индекса: т и п - индексы координат, к - индекс времени.

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

являются асимптотическими при стремлении шага к нулю, но даже быстродействующие современные ЭВМ не позволяют считать сложные реальные многомерные задачи с достаточно малым шагом. Если взять п интервалов по каждой из переменных, то сетка для р-мерной задачи будет содержать пР узлов. Реально может оказаться, что схема первого порядка точности на грубых сетках даст более точный результат, чем схема второго порядка точности, хотя на подробных сетках соотношение будет обратным. Обычно априорные оценки точности схем далеки от оптимальных и бывают завышены в десятки-сотни. Кроме того, не всякая сходящаяся разностная схема гарантирует получение качественного решения. Так, сходимость в гильбертовой норме (среднеквадратическая сходимость) обеспечивает передачу только некоторых интегральных характеристик решения, в чебышевской (сходимость в среднем) - обеспечивает хорошее качество решения только на подробной сетке, а на грубой сетке часто возникает так называемая разболтка, и результаты расчета получаются непригодными. Выбор метрического пространства фактически определяет класс функций, в которых следует искать решение. Чебышевская норма является более сильной, чем гильбертова, но ее использование не всегда удобно и во многих задачах необязательно. Для каждой конкретной задачи выбор пространств определяется в первую очередь физическим смыслом задачи, и только потом принимаются во внимание чисто математические соображения, такие, например, как возможность доказать сходимость. Так, в задачах о нагреве тела потоком тепла даже норма L\ удовлетворительна [8], поскольку температура тела определяется интегралом от потока по времени.

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

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

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

Литература

1. Чермошенцева А.А. Течение теплоносителя в геотермальной скважине // Математическое моделирование. - 2006. - Т. 18, № 4. - С. 61-76.

2. Шулюпин А.Н, Чермошенцева А.А. Влияние теплообмена с окружающими породами на эксплуатационные параметры пароводяной скважины // Научный журнал КубГАУ. - Краснодар: КубГАУ, 2007. - № 23(03).

3. Лыков А.В. Теория теплопроводности. - М.: Высш. шк., 1967. - 599 с.

4. Бронштейн И.Н., Семендяев К.А. Справочник по математике для инженеров. - М.: Наука, 1964. - 608 с.

5. Исаченко В.П., Осипова В.А., Сукомел А.С. Теплопередача. - М.: Энергоиздат, 1981. - 416 с.

6. Боглаев Ю.П. Вычислительная математика и программирование. - М.: Высш. шк., 1990. -543 с.

7. Власова Е.А., Зарубин В.С., Кувыркин Г.Н. Приближенные методы математической физики. - М.: Изд-во МГТУ им Н.Э. Баумана, 2001. - 700 с.

8. Калиткин Н.Н. Численные методы. - М.: Наука, 1978. - 512 с.

9. Лыков А.В. Тепломассообмен: Справочник. - М.: Энергия, 1978. - 479 с.

10. Марчук Г.И. Методы вычислительной математики. - М.: Наука, 1980. - 536 с.

11. Самарский А.А. Теория разностных схем. - М.: Наука, 1989. - 614 с.

12. Самарский А.А., Гулин А.В. Устойчивость разностных схем. - М.: Наука, 1973. - 415 с.

13. Самарский А.А., Гулин А.В. Численные методы. - М.: Наука, 1989. - 429 с.

14. Самарский А.А., Попов Ю.П. Разностные методы решения задач газовой динамики. - М.: Наука, 1992. - 422 с.

15. Созинова Т.Е. Разработка метода расчета и исследования теплового и термонапряженного состояния крепи геотермальных скважин: Автореф. дис ... канд. тех. наук. - Иваново, 1997. - 24 с.

16. Тихонов А.Н., Самарский А.А. Уравнения математической физики. - М.: Наука, 1966. - 724 с.

17. Чермошенцева А.А. Численные методы: Учеб. пособие. (ДВ РУМЦ). - Петропавловск-Камчатский, Изд-во КамчатГТУ, 2008. - 110 с.