Численное моделирование осесимметричных процессов теплопроводности

Введение

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

Установившееся распределение температуры u(x) в прямоугольной области описывается дифференциальным уравнением 

div(grad u) – u = -f(x, y, t),

где x,y,t П, а f(x, y, t) – плотность источников тепла. Для того, чтобы однозначно найти функцию u(x, y, z) необходимо задать два дополнительных условия. Зададим простейшие условия первого рода:

 

u(x, y, t) = при x,y,t границе П,

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

Решением данной краевой  задачи будет являться дважды дифференцируемая в области П функция u(x, y, t), если она является решением дифференциального уравнения и удовлетворяет краевым условиям.

 

Для поиска решения обычно пользуются методом конечных разностей (методом сеток), который является одним из универсальных и широко используемых методов решения краевых задач. Его популярность во многом объясняется относительной простотой подхода к дискретизации дифференциальных уравнений. Область непрерывного изменения аргумента заменяют конечным (дискретным) множеством точек (узлов), называемым сеткой. Вместо функций непрерывного аргумента рассматривают функции, определенные только в узлах сетки – сеточные функции. Производные, которые входят в дифференциальное уравнение и краевые условия заменяют их разностными аналогами – линейными комбинациями значений сеточных функций в некоторых узлах сетки. В результате краевую задачу заменяют дискретной краевой задачей (разностной схемой), представляющей собой систему конечного числа линейных и нелинейных алгебраических уравнений. Решение разностной схемы принимают за приближенное решение краевой задачи.

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

 

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

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

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

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

 

Постановка задачи

в прямоугольных декартовых координатах

 

Пусть R – замкнутый прямоугольный параллелепипед, R’ – полуоткрытый параллелепипед, а  S=R/R’ – множество граничных точек прямоугольного параллелепипеда R, f(x, y, t) – заданная на R достаточно гладкая функция. Требуется найти непрерывную на R функцию u(x, y, t), удовлетворяющую на R’ уравнению теплопроводности (1):

 

 

 

и, кроме того, подчиняющуюся  на S дополнительному условию

u=0    (2)

 

Условие (2) включает в себя как начальное условие 

 

u(x, y, 0) = 0 при t = 0,

 

так и однородные краевые  условия первого рода

 

u(x, y, t) = 0 при x = 0, x = 1, y = 0, y = 1,

 

т.е. на боковых гранях параллелепипеда  R.

Смешанная задача (1), (2) имеет  единственное решение u(x, y, t).

 

 

Дискретизация задачи.

Простейшие  разностные схемы

 

Положим для простоты шаги постоянными , ,

 

 

= -   (3)

 = -   (4)

=   (5)

Явная разностная схема  (6)

,

k,m, = 1,2 … N-1,

проста в применении, но устойчива только при    в равномерных нормах.

 

Неявная разностная схема  (7)

,

k,m, = 1,2 … N-1,

абсолютно устойчива, но на каждом слое по времени требуется  решить систему  уравнений. (Подразумевается, что разностное решение v в точках, расположенных на S равно 0).

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

Метод переменных направлений

 

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

 

   (8)

,   (9)

k,m, = 1,2 … N-1,

где . [1]

 

В разностной схеме (8), (9) шаг  по времени  делится на два полушага (см. рис. 1)

 

Рис.1. Метод переменных направлений.

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

 

   (10)

 

где известно и присоединим к разностному уравнению (10) краевые условия

 

    (11)

в соответствии с условием (2).

 

Видим, что разностная задача (10), (11) распадается на (N-1) независимых трехточечных разностных краевых задач, отвечающих каждому фиксированному m, . Разностная краевая задача (10), (11) решается методом прогонки при каждом m отдельно. Прогонка осуществляется по индексу k, т.е. в направлении оси x. На решение разностной задачи (10), (11) при одном значении затрачивается O(N) арифметических действий, а значит, всего на (N-1) задач расходуется арифметических действий.

 

После того как найдены  все неизвестные   на промежуточном слое с номером переносим их в разностном уравнении (9), соответствующем второму полушагу, вправо. Это разностное уравнение переписываем в виде

 

   (12)

 

где

 

известно и присоединяем к уравнению (12)  в соответствии с условием (2) краевые условия

 

  (13)

Задача (12), (13) тоже распадается  на N – 1 трехточечных разностных краевых задач, отвечающих различным фиксированным k, . Каждая такая задача решается методом прогонки. Прогонка производится теперь уже по индексу m, т.е. в направлении оси y.  На решение задач (12), (13) при всех k расходуется арифметических действий.

 

Таким образом, при переходе от -го слоя к -му слою по схеме (8), (9) затрачивается число арифметических действий порядка числа искомых неизвестных на одном слое, т.е. . Такая разностная схема называется экономичной. [4]

 

Разностная схема (8), (9) является абсолютно устойчивой в некоторых  естественных нормах. Если решение  u(x, y, t) задачи (1), (2) достаточно гладкое на замкнутом параллелепипеде R , то решение v разностной схемы (8), (9) сходится к u в следующем смысле

 

,    (14)

 

где

 

 

 

 

Cтационарное уравнение

в осесимметричном случае

 

Стационарное уравнение  диффузии или теплопроводности

 

div(k grad u) – qu = -f(r, φ, z)

 

в цилиндрической системе  координат (r, φ, z) в случае, когда решение u = u(r) не зависит ни от z, ни от  φ (имеет место осевая симметрия), принимает вид

 

     (20)

 

 

 

При r = 0 ставится условие ограниченности

     (21)

 

При r = 1 ставится обычное краевое условие, например,

    (22)

 

Пусть   и - линейно независимые решения уравнения (20), причем   ограничено при Тогда справедливы свойства:

1. Если q(0) и f(0) конечны, то и
2. Если то производные ограничены при
3. Второе, линейно независимое с решение уравнения (20) имеет при r = 0 логарифмическую особенность.

Условия (21) и (22) выделяют единственное решение уравнения (20). В силу свойства 1) условие (21) можно заменить требованием

    (23) [2]

 

 

Разностные  схемы для стационарного уравнения

в осесимметричном случае

 

Введем равномерную сетку

 .

Разностную схему для  уравнения (20) напишем при помощи метода баланса 

 

    (24)

где

 

 

 

Аппроксимируя поток  выражением   и заменяя интегралы в уравнении баланса (24) выражениями   и   соответственно, получаем разностное уравнение

 

    (25)

 

где    коэффициенты и выбираются так, чтобы

 

.    (26)

В простейшем случае

 

.      (27)

Аппроксимируем краевое  условие при r = 0. Его можно записать как условие равенства нулю потока при . Покажем, что разностное краевое условие

 

   (28)

имеет погрешность аппроксимации  на решении уравнения (20), удовлетворяющему краевому условию (21).

В самом деле, невязка  для (28) равна

    (29)

Подставим сюда

 

получаем

 

    (30)

 

Из уравнения (20) имеем . Так как и

    (31)

Подставляя (31) в формулу (30) и учитывая (23), получаем

Разностное краевое условие (28) будем записывать в виде

 

 

 

Таким образом, задаче (20) - (22) поставим в соответствие разностную схему

 

 

   (32) [2]

 

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

 

     (33)

 

 

с краевыми условиями

 

      (34)

 

где 

Эта задача решается методом  прогонки, условия устойчивости которой выполнены, так как

 

 

 

 

 

 

 

 

 

 

Оценка  точности разностной схемы

для осесимметричного случая

 

Перейдем к оценке точности схемы (32). Подставляя в (32)  y = z + u, где u – решение задачи (20) – (22), а y – решение задачи (32), получим для погрешности z = y – u  задачу

 

 

 

      (35)

 

где и   - погрешности аппроксимации уравнения

 

    (36)

и краевого условия

 

      (37)

Пользуясь уравнением баланса (24), преобразуем, как обычно, в виду

 

 

 

Положим    и получим разложение интегралов, входящих в формулу для , по степеням h:

 

 

 

и, аналогично,

 

 

так что

 

т.е.

 

Очевидно, что функция т.е.

 

 

Сравнивая формулы (29) и (37), находим, что   Отсюда следует, что так как выше мы показали, что [2]

 

Перейдем к оценке погрешности Нам понадобится

 

Лемма. Пусть z - решение задачи (35), а v – решение той же задачи при . Тогда имеет место неравенство

 

 

 

Функция находится в явном виде из условий

 

 

 

где определяется формулой (36) при

 

Суммируя уравнения   получим

 

 

 

Подставим в (38)   :

 

 

Просуммируем (39) по

 

 

 

Подставляя в (40)   находим

 

 

 

 

 

так как    Далее имеем

 

 

В результате приходим к  следующему неравенству:

 

 

 

Подставляя в (41) оценки

 

 

 

убеждаемся в том, что и, следовательно,

 

 

 

т.е. схема (32) имеет второй порядок точности в С.

 

Рассмотрим схему второго  типа – «схему на потоковой сетке». Разобьем отрезок на N частей введя узлы (потоковые точки)

 

 

 

Пусть - значения искомой сеточной функции в этих узлах.

Для получения разностной схемы воспользуемся уравнением баланса для (20). Рассматривая уравнение  баланса (аналог(24)) для интервала    получаем

 

 

 

где   выбираются по аналогии с (26) - (27), так что в простейшем случае

 

 

 

Из уравнения баланса  для интервала

 

 

и условия следует разностное уравнение при

 

 

 

где определяются по формулам (43).

 

При ставится обычное условие

 

 

 

В результате получаем разностное уравнение (42) с краевыми условиями (44) и (45).

 

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

 

 

 

 

 

 

 

где - погрешность аппроксимации, равная

 

 

 

 

 

 

 

Отсюда и из уравнения  баланса для интервала   следует, что

 

 

 

 

 

 

 

где

 

Вычисления дают

 

 

 

Для решения задачи (46) с  правой частью

 

 

 

справедлива оценка того же типа, что и для задачи (35). Из этой оценки следует

 

 

 

т.е. схема (42), (44), (45) имеет  второй порядок точности, если

 

 

 

 

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

в полярных координатах

 

Уравнение теплопроводности в полярных координатах

 

 

 

 

 

 

 

Рис.2 Прямоугольная область  в полярных координатах.

 

При r = 1 условие I рода.

При r = 0 особое условие (21).

На верхней и нижней границах области – условие периодичности

 

.

 

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

Циклическая прогонка

 

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

 

Рассмотрим систему уравнений

 

 

   ,         (14)

 

 

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

   

 

Относительно коэффициентов  системы будем предполагать, что   (15)

Приведем получающиеся формулы  решения задачи (14) – формулы циклической  прогонки:

 

 

  (17)

 

,   

 

 

 

Метод циклической прогонки является устойчивым, так как решения  задач (17) ищутся методом прогонки, который устойчив при выполнении условий (15), а знаменатель в выражении для не обращается в нуль. Действительно, из (15), (16) видно, что . .   Предполагая , получаем

 

.  (19)

 

Учитывая (17) и (19), находим . Из всего сказанного следует, что

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Приложение 1

Тестовая  задача

Количество шагов по осям x и y

по t

_шаг

отрезки времени

_{коэффиценты}

 

 

 

 

 

 

 

Заданная функция

_{нижняя граница}

_{передняя граница}

_{задняя граница}

_{левая граница}

_{правая граница}

Заданное решение:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Приложение 2

Циклическая прогонка

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Библиографический список

 

1. Волков Е.А. Численные  методы. М. Наука. 1987. 
2. Самарский А.А. Теория разностных схем М. Наука. 1989. 
3. Самарский А.А. Николаев Е.С. Методы решения сеточных уравнений. М. Наука. 1978.

4. Калиткин Н.Н. Численные  методы. М. Наука. 1978.

5. Амосов А.А., Дубинский  Ю.А., Копченова Н.В. Вычислительные  методы для инженеров. М. МЭИ. 2003.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Оглавление

Введение 3

Постановка задачи в прямоугольных декартовых координатах 5

Дискретизация задачи. Простейшие разностные схемы 6

Метод переменных направлений 7

Cтационарное уравнение в осесимметричном случае 10

Разностные схемы для стационарного уравнения

в осесимметричном случае 12

Оценка точности разностной схемы для осесимметричного случая 15

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

в полярных координатах 21

Циклическая прогонка 22

Приложение 1. Тестовая задача 24

Приложение 2. Циклическая прогонка 28

Библиографический список 29