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

Уравнение теплопроводности

Ранее (см. разд. 2.1.2, 2.1.3) уже были построены и исследованы разностные схемы решения смешанной задачи для одномерного уравнения теплопроводности:

(2.75)

Были получены две двухслойные схемы — явная (2.3) и неявная (2.4). В явной схеме значения сеточной функции на верхнем (j + 1)-ом слое вычисляли с помощью решения на нижнем слое [соотношение (2.13)]. Для нахождения решения на (j + 1)-м слое по неявной схеме была получена трехдиагональная система линейных алгебраических уравнений (2.22), которую решают методом прогонки.

Неявная схема безусловно устойчива, явная схема устойчива при выполнении условия

Обе схемы сходятся к решению исходной задачи со скоростью .

Схемы (2.3), (2.4) построены для случая, когда значения искомой функции (температуры) Uна границах х = 0, х = 1определяются заданными функциями . Однако граничные условия в смешанной задаче (2.75) могут быть и иными, в них может входить производная искомой функции. Например, если конец стержня х=0 теплоизолирован, то условие имеет вид

В этом случае, как и при решении волнового уравнения, данное условие нужно записывать в схемах (2.3), (2.4) в разностном виде.

Перейдем теперь к построению разностных схем для уравнения теплопроводности с двумя пространственными переменными. Примем для простоты а = 1. Тогда это уравнение можно записать в виде

(2.76)

Пусть при t=0 начальное условие задано в виде

(2.77)

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

Часто задачи теплопроводности или диффузии, описываемые двумерным уравнением (2.76), решаются в ограниченной области. Тогда, кроме начального условия (2.77), нужно формулировать граничные условия. В частности, если расчетная область представляет прямоугольный параллелепипед (рис. 2.24), то нужно задавать граничные условия на его боковых гранях. Начальное условие (2.77) задано на нижнем основании параллелепипеда.

Рис. 2.24. Расчетная область

Введем простейшую сетку с ячейками в виде прямоугольных параллелепипедов, для чего проведем три семейства плоскостей: хi= ih1(i=0,1. I), (j=0,1. J), . Значение сеточной функции в узлах обозначим символом . С помощью этих значений можно построить разностные схемы для уравнения (2.76).

Рассмотренные выше схемы для одномерного уравнения легко обобщаются на двумерный случай.

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

Рис. 2.25. Шаблон двумерной схемы

Отсюда можно найти явное выражение для значения сеточной функции на (k + 1)-ом слое:

(2.78)

Условие устойчивости имеет вид

(2.79)

При получается особенно простой вид схемы (2.78):

(2.80)

Полученная схема сходится со скоростью

Формулы (2.78) или (2.80) представляют собой рекуррентные соотношения для последовательного вычисления сеточной функции во внутренних узлах слоев k = 1,2. К. На нулевом слое используется начальное условие (2.77), которое записывается в виде

(2.81)

Значения в граничных узлах вычисляют с помощью граничных условий.

Алгоритм решения смешанной задачи для двумерного уравнения теплопроводности изображен на рис. 2.26. Здесь решение хранится на двух слоях: нижнем (массив ) и верхнем (массив ). Блоки граничных условий необходимо сформировать в зависимости от конкретного вида этих условий. Результаты выводят на каждом слое, хотя можно ввести шаг выдачи (см. рис. 2.13).

Рис. 2.26. Алгоритм решения двумерного уравнения теплопроводности

Построим теперь абсолютно устойчивую неявную схему для решения уравнения (2.76), аналогичную схеме (2.4) для одномерного уравнения теплопроводности. Аппроксимируя в (2.76) вторые производные по пространственным переменным на (k + 1)-ом слое, получаем следующее разностное уравнение:

(2.82)

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

(2.83)

К этой системе уравнений нужно добавить граничные условия для определения значений сеточной функции в граничных узлах (т.е. при i= 0, I; j = 0, J). На нулевом слое решение находится из начального условия (2.77), представленного в виде (2.81).

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

Недостатком явной схемы (2.78) является жесткое ограничение на шаг по времени τ, вытекающее из условия (2.79). Существуют абсолютно устойчивые экономичные разностные схемы, позволяющие вести расчет со сравнительно большим значением шага по времени и требующие меньшего объема вычислений. Две из них будут рассмотрены в разд. 2.3.3.

Методы решения УЧП

МЕТОДЫ РЕШЕНИЯ УЧП

1. Уравнение теплопроводности

Рассмотрим численное решение уравнений с частными производными (УЧП) на примере решения уравнения теплопроводности (диффузии)

— плотность тепловых источников, — коэффициент теплопроводности

Если , то мы имеем стационарное уравнение теплопроводности

— уравнение Пуассона.

Если плотность тепловых источников тоже равна 0, то получаемое уравнение называется уравнением Лапласа.

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

Это могут быть или значения искомой функции на границе (задача Дирихле).

Или значения потока тепла на границе (задача Неймана)

Либо смешанные условия.

Если необходимо решать УЧП в области, имеющей круговую симметрию, оператор Лапласа удобна записать в полярных координатах:

2. РЕШЕНИЕ СТАЦИОНАРНОГО УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ

Будем решать задачу Дирихле для уравнения Пуассона.

в прямоугольной области , .

Для численного решения данной задачи применим метод SOR (метод последовательной верхней релаксации). Вначале используем метод конечных разностей. Для этого разобьем отрезок [a, b] на K равных интервалов длиной , а отрезок [c, d] — на L интервалов той же длины. Пусть при этом , . Тогда

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

Запишем вначале в одномерном случае рамках метода конечных разностей производную функции U по х в точке . Получим

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

Аналогичным образом, для второй производной в точке имеем

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

,

Подставив полученные таким образом вторые производные в уравнение Пуассона, имеем

Коэффициенты перед матричными элементами в данном случае равны 1 и 4. Однако в общем случае (полярные координаты, например) мы должны записать

и вычислить A, B, C, D, E.

Перепишем полученное нами уравнение в следующем виде:

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

Подставляя в правую часть (3) уравнение (1) имеем

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

Это и есть основное уравнение метода SOR. Параметр обычно выбирается не тождественно равным 1, изменяющимся параметром от 1 до 1.5. Это делается для лучшей сходимости метода. В качестве начального приближения берется функция, тождественно равная 0 внутри области определения и равная граничным условиям на границе области.

3. РЕШЕНИЕ НЕСТАЦИОНАРНОГО УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ

Будем искать решение этой задачи методом Кранка-Николсона. Рассмотрим вначале одномерное нестационарное уравнение диффузии.

Запишем производную по времени в k-й точке х в n-й момент времени в виде:

Здесь — шаг по времени. В правой части вторую производную по координате можно взять как в момент времени ,

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

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

Для двумерного уравнения диффузии

в рамках метода Кранка-Николсона запишем вторые производные по х и по у

как среднее арифметическое от производных в точках и ,

перенесем все значения функции в момент времени влево, а в — й вправо. Получим.

,

(1)

(2)

Схема расчета по методу Кранка-Николсона такова:

1. В начальный момент времени из начальных и граничных условий методом SOR (или каким-либо другим методом) находим значения функции U.

2. Вычисляем с ее помощью в этот же момент времени по формуле (2) функцию F.

3. По формуле (1) с помощью уже вычисленной функции F и граничных условий находим значение функции U в следующий момент времени.

4. Вычисляем с ее помощью по формуле (2) функцию F.

5. По формуле (1) находим значение функции U в следующий момент времени.

МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ

4. Решение одномерного уравнения теплопроводности

Базисные функции метода конечных элементов

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

Пусть отрезок [a, b], на котором определяется искомая функция, разбит N точками на (N-1) равных отрезков (конечных элементов) длиной . На рисунке показан случай, когда N = 4.

Значения функции в этих точках равны , , . , . на каждом из конечных элементов функцию u(x) можно представить с помощью базисных функций :