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

Введение

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

Установившееся распределение температуры 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]