Исследование аппроксимации на примере уравнения теплопроводности

Метод Finite Volume — реализация на примере теплопроводности

Метод Finite Volume (FVM)

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

Термин конечный объем в статье будет часто заменятся на Элемент, будем для удобства считать их эквивалентами (элемент в данной статье не имеет ничего общего с методом конечных элементов).

Есть 2 различных способа решения задачи по FVM:
1) грани контрольного объема совпадают с гранями элемента
2) грани контрольного объема проходят через центры граней элементов(на которые разбита область).Искомые переменные хранятся в вершинах этих элементов.Вокруг каждой вершины строится контрольный объем. Для непрямоугольной сетки этот способ имеет еще 2 подвида.

Мы будем использовать способ 1) с контрольными объемами совпадающими с элементами на которые разбита область.

Некоторые плюсы FVM:

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

Метод FVM реализуем на примере уравнения теплопроводности:

Итак основные шаги при реализации FVM:

1. Перевод дифф уравнения в форму пригодную для FVM — интегрирование по контрольному объему
2. Составление дискретного аналога, выбор способа перевода производных и других подынтегральных выражений в дискретную форму
3. Получение уравнения для каждого из контрольных объемов, на которые разбита область.Составление системы линейных уравнений и ее решение.

Дискретизация по времени.

Немного теории или первый шаг в реализации FVM

FVM на стандартной прямоугольной сетке

На рисунке изображен Элемент С и его соседние элементы справа(E), слева(W), сверху(N) и снизу(S).У элемента С есть 4 грани обозначенные буквами e w n s.Именно эти 4 грани и составляют периметр элемента и по ним производится интегрирование.Для каждого элемента в результате получаем дискретный аналог исходного дифф уравнения.

Составим дискретный аналог для элемента С.Для начала нужно разобраться с интегралом (3).Интеграл это ведь по факту сумма.Поэтому мы и заменяем интеграл по всей поверхности элемента, на сумму по 4-м составляющим этой поверхности, тоесть 4 граням элемента.

Уравнение (7) и есть конечное уравнение для элемента С, из него мы на каждом шаге по времени получаем новое значение температуры (Tnew) в элементе С.

Граничные условия на прямоугольной сетке

Мы рассмотрим только 2 вида граничных условий.

1. Задана температура Tb на границе
2. Задан поток FluxB на границе, рассмотрим только случай когда FluxB=0, т.е. грань e будет теплоизолирована(Insulated)

Случай 2) самый простой, поскольку получается что грань e не потребуется при дискретизации(т.к. все коэффициенты Flux=0) и можно ее просто пропустить.

Теперь рассмотрим случай 1).Дискретизация грани e будет в целом похожа на ту что уже была описана.Будут только 2 изменения — вместо Te будет известное граничное значение Tb и вместо расстояния DXe будет DXe/2.В остальном можно рассматривать значение Te так, как будто это был бы обычный соседний узел E.Теперь подробнее распишем терм для граничного элемента С.

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

FVM в задачах со сложной геометрией

Здесь как раз проявляется преимущество FVM, где также, как и в методе конечных элементов, можно представлять область с круглыми границами через разбиение на треугольники или любые другие полигоны.Но FVM имеет еще 1 плюс — при переходе от треугольников к полигонам с большим числом сторон не требуется абсолютно ничего менять, конечно если код был написан для произвольного треугольника а не равностороннего.Более того, можно без изменения кода использовать смесь разных элементов — треугольники, полигоны, квадраты и тд.

Рассмотрим общий случай, когда вектор соединяющий центры 2-х элементов не совпадает с вектором нормали к общей грани этих элементов.Вычисление потока flux через грань теперь будет состоять из 2-х частей.В первой будет расcчитываться ортогональная составляющая а во второй так называемая «кросс-диффузия».

На картинке изображены 2 элемента, С — текущий рассматриваемый элемент и F — соседний элемент.Опишем дискретизацию для грани, разделяющей эти 2 элемента.Вектор соединяющий центры элементов — DCF.Вектор e — это единичный вектор по направлению DCF.Вектор Sf — направлен по нормали к грани, его длинна равна длине грани.

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

Расчет граничных элементов ничем не отличается от расчетов не на границе, вместо центра соседнего элемента берется центр грани, ну и как обычно подставляется температура на границе.
В моей реализации в итоге получается так:

Примеры и проверка результатов

Описание структуры исходников

Гитхаб с исходниками лежит тут
Основная версия в папке heat2PolyV2.То что относится к вычислительной части лежит в heat2PolyV2\Src\FiniteVolume\.

Вначале файла Scene2.cs — параметры которые можно менять для отображения в разных цветовых схемах, масштаб, отображение mesh и т.д.Сами примеры хранятся в heat2PolyV2\bin\Debug\Demos\
Выгрузку из Матлаба сделать просто — нужно открыть pde toolbox, открыть m файл (либо создать самому с нуля), зайти в меню Mesh-Экспорт mesh, нажать ОК; перейти в основной Матлаб, в панельке появятся переменные — матрицы p e t, открыть файл savemymesh.m, выполнить его, появится файл p.out, перенести его в папку Demos.
В исходниках для выбора примера необходимо задать имя файла в строке param.file = «p»;(FormParam.cs).Далее необходимо применить граничные условия — для готовых примеров можно просто раскомментировать соответствующие блоки в MainSolver.cs:

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

Разностные аппроксимации

Разностные аппроксимации

1.Примеры разностных аппроксимаций.

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

Пусть u(x) – достаточно гладкая функция, заданная на отрезке [x i-1 , x i+1 ]. Обозначим

называются соответственно правой, левой и центральной разностными производными функции u(x) в точке x i , т.е. при фиксированном x i и при h ® 0 (тем самым при i ® ¥ ) пределом этих отношений является u’(x i ) . Проводя разложение по формуле Тейлора, получим

u x,i – u’(x i ) = 0,5hu’’(x i ) + O(h 2 ),

u x,i – u’(x i ) = -0,5hu’’(x i ) + O(h 2 ),

u x,i – u’(x i ) = O(h 2 ),

Отсюда видно, что левая и правая разностные производные аппроксимируют u’(x) с первым порядком по h , а центральная разностная производная – со вторым порядком. Нетрудно показать, что вторая разностная производная

аппроксимирует u’’(x i ) со вторым порядком по h , причем справедливо разложение

Рассмотрим дифференциальное выражение

(1)

с переменным коэффициентом k(x) . Заменим выражение (1) разностным отношением

(2)

где a=a(x) – функция, определенная на сетке w h . Найдем условия, которым должна удовлетворять функция a(x) для того, чтобы отношение (au x ) x,i аппроксимировало (ku’)’ в точке x i со вторым порядком по h . Подставляя в (2) разложения

где u i ’ = u’(x i ) , получим

С другой стороны, Lu = (ku’)’ = ku’’ + k’u’,

т.е.

Отсюда видно, что L h u–Lu = O(h 2 ) , если выполнены условия

(3)

Условия (3) называются достаточными условиями второго порядка аппроксимации . При их выводе предполагалось, что функция u(x) имеет непрерывную четвертую производную и k(x) – дифференцируемая функция. Нетрудно показать, что условиям (3) удовлетворяют, например, следующие функции:

Заметим, что если положить a i = k(x i ), то получим только первый порядок аппроксимации.

В качестве следующего примера рассмотрим разностную аппроксимацию оператора Лапласа

(4)

Введем на плоскости (x 1 , x 2 ) прямоугольную сетку с шагом h 1 по направлению x 1 и с шагом h 2 по направлению x 2 , т.е. множество точек

Из предыдущих рассуждений следует, что разностное выражение

(5)

аппроксимирует дифференциальное выражение (4) со вторым порядком, т.е. L h u ij – Lu(x i 1 , x j 2 ) = O(h 2 1 ) + O(h 2 2 ). Более того, для функций u(x 1 , x 2 ), обладающих непрерывными шестыми производными, справедливо разложение

Разностное выражение (5) называется пятиточечным разностным оператором Лапласа , так как оно содержит значения функции u(x 1 , x 2 ) в пяти точках сетки, а именно в точках (x 1 i , x 2 j ), (x 1 i ± 1 , x 2 j ), (x 1 i , x 2 j ± 1 ). Указанное множество точек называется шаблоном разностного оператора. Возможны разностные аппроксимации оператора Лапласа и на шаблонах, содержащих большее число точек.

2. Исследование аппроксимации и сходимости

2.1. Аппроксимация дифференциального уравнения.

Ранее рассматривалась краевая задача

(k(x) u’(x))’ – q(x) u(x) + f(x) = 0, 0 (1)

– k(0) u’(0) + b u(0) = m 1 , u(l) = m 2 , (2)

k(x) ³ c 1 > 0, b ³ 0,

для которой интегро-интерполяционным методом была построена разностная схема

(3), (4)

Обозначим через Lu(x) левую часть уравнения (1) и через L h y i – левую часть уравнения (3), т.е.

Пусть u (x) – достаточно гладкая функция и u (x i ) – ее значение в точке x i сетки

Говорят, что разностный оператор L h аппроксимирует дифференциальный оператор L в точке x=x i , если разность L h u i – L h u (x i ) стремится к нулю при h ® 0. В этом случае говорят также, что разностное уравнение (3) аппроксимирует дифференциальное уравнение (1).

Чтобы установить наличие аппроксимации, достаточно разложить по формуле Тейлора в точке x=x i значения u i ± 1 = u (x i ± h) , входящие в разностное выражение L h u i . Большая часть этой работы проделана в предыдущей главе, где показано, что при условиях

(8)

Если кроме того, докажем, что

d i = q(x i ) + O(h 2 ), j i = f(x i ) + O(h 2 ) (9)

то тем самым будет установлено, что оператор L h аппроксимирует L со вторым порядком по h , т.е.

L h u i – L u (x i ) = O(h 2 ), i = 1, 2,…, N–1 (10)

Итак, доказательство второго порядка аппроксимации сводится к проверке сводится к проверке условий (8), (9) для коэффициентов (5), (6). Проверим сначала выполнение условий (8). Обозначая p(x) = k -1 (x) , получим

т.е. условия (8) выполнены. Условия (9) выполнены в силу того, что замена интегралов (6) значениями q i , f i соответствует приближенному вычислению этих интегралов по формуле прямоугольников с узлом в середине отрезка интегрирования.

2.2. Аппроксимация граничного условия.

Исследуем погрешность аппроксимации разностного граничного условия (4). Обозначим l h u (0) = –a 1 u x, 0 + b u 0 . Если u (x) – произвольная достаточно гладкая функция, то очевидно

l h u (0) = –k(0) u ’(0) + b u (0) + O(h) ,

т.е. имеет место аппроксимация первого порядка по h . Однако если u =u(x) – решение задачи (1), (2), то разностное граничное условие (4) имеет второй порядок аппроксимации, т.е.

Докажем последнее утверждение. Используя разложение

u x, 0 = (u 1 – u 0 )/h = u’(x 1/2 ) + O(h 2 ), x 1/2 = 0,5h,

a 1 = k 1/2 + O(h 2 )

Учитывая граничное условие (2), получаем

l h u(0) = 0,5h [– (ku’)’(0) + d 0 u 0 – j 0 ] + O(h 2 ) .

Выражение, стоящее в квадратных скобках, преобразуем, учитывая уравнение (1), к виду

– (ku’)’(0) + d 0 u 0 – j 0 = – (ku’)’(0) + q(0)u(0) – f(0) +

+ (d 0 – q(0))u 0 – (f(0) – j 0 ) = (d 0 – q(0))u 0 – (f(0) – j 0 ) .

что и требовалось доказать.

Таким образом, при достаточной гладкости коэффициентов k(x), q(x), f(x) и решения u(x) разностная схема (10) аппроксимирует исходную задачу (2) со вторым порядком по h .

При практическом использовании разностной схемы для нахождения ее коэффициентов не обязательно вычислять интегралы (4), (6) точно. Можно воспользоваться коэффициентами, полученными путем замены этих интегралов квадратурными формулами, имеющими точность O(h 2 ) и выше. Например, в результате применения формулы прямоугольников получим следующие коэффициенты: a i = k(x i – 0,5h), d i = q(x i ), j i = f(x i ).

Применяя формулу трапеций, получим

Представление коэффициентов разностной схемы в виде интегралов (4), (6) оказывается полезным при исследовании сходимости в случае разрывных функций k(x), q(x), f(x) .

2.3. Уравнение для погрешности.

Решение y i = y(x i ) разностной задачи (3), (4) зависит от шага h сетки, y(x i ) = y h (x i ) . По существу, мы имеем семейство решений , зависящее от параметра h . Говорят, что решение y h (x) разностной задачи сходится к решению u(x) исходной дифференциальной задачи, если при h ® 0 погрешность y h (x i ) – u(x i ), i = 0, 1,…, N , стремится к нулю в некоторой норме. В настоящем параграфе в качестве такой нормы будем брать норму в сеточном пространстве C( w h ) , т.е. положим

Говорят, что разностная схема имеет m-й порядок точности (или сходится с порядком m ), если

где m>0, M>0 – константы, не зависящие от h .

Выше было установлено, что схема (3), (4) имеет второй порядок аппроксимации. Докажем теперь, что эта схема имеет и второй порядок точности. Для этого прежде всего выпишем уравнение, которому удовлетворяет погрешность z i = y i – u(x i ) . Поставим y i = z i + u(x i ) в уравнения (3), (4). Тогда получим уравнения

(11), (12)

Функция y i , входящая в правую часть уравнения (11), называется погрешностью аппроксимации дифференциального уравнения (1) разностным уравнением (3) на решении задачи (1), (2). В п.1 было доказано, что y i = O(h 2 ) при h ® 0, i=1, 2,…, N–1 . Аналогично, величина n 1 является по определению погрешностью аппроксимации краевого условия (2) разностным краевым условием (4) на решении задачи (1), (2), причем n 1 =O(h 2 ) . Таким образом, структура уравнений для погрешности (11), (12) та же, что и у разностной схемы (3), (4), отличаются только правые части.

Чтобы доказать сходимость разностной схемы, оценим решение задачи (11), (12) через правые части y i , n 1 , т.е. получим неравенство вида

(13)

с константой M 1 , не зависящей от h . Из этого неравенства и будет следовать, что

Отметим, что неравенства вида (13), называемые априорными оценками, нашли широкое применение в теории разностных схем. Поскольку структура для погрешности (11), (12) та же, что и у разностной схемы (3), (4), а отличаются только правые части, то оценка (13) выполняется одновременно с аналогичной оценкой

для разностной схемы (3), (4) при m 2 = 0 . Последняя оценка выражает устойчивость решения разностной задачи по правым частям j и m 1 .

2.4. Разностные тождества и неравенства.

Для того, чтобы доказать неравенство (13), нам потребуются некоторые разностные тождества и неравенства. Будем рассматривать сеточные функции, заданные на сетке (7). Обозначим

Справедливо следующее разностное утверждение:

(y, u x ) = –( u , y x ) + y N u N – y 0 u 1 . (14)

что и требовалось доказать. Тождество (14) называется формулой суммирования по частям .

Подставляя в (14) вместо u выражение az x и вместо y функцию z, получаем первую разностную формулу Грина

(15)

Здесь

В частности, если z N = 0 (как в задаче (11), (12)), то получим

(16)

и докажем, что для любой сеточной функции z i , удовлетворяющей условию z N = 0 , справедливо неравенство

(17)

Для доказательства воспользуемся тождеством

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

Откуда сразу следует неравенство (17).

2.5. Доказательство сходимости.

Возвращаясь к доказательству сходимости схемы (3), (4), получим тождество, которому удовлетворяет погрешность z i = y i – u(x i ) . Для этого умножим уравнение (11) на hz i и просуммируем по i от 1 до N–1 . Тогда получим

Отсюда, применяя разностную формулу Грина (16), получим

Далее, согласно (12) имеем

следовательно, справедливо тождество

(18)

Из этого тождества и будет сейчас выведено требуемое неравенство вида (13).

Заметим прежде всего, что если

k(x) ³ c 1 > 0, b ³ 0, q(x) ³ 0,

то коэффициенты разностной схемы (3), (4) удовлетворяют неравенствам

a i ³ c 1 > 0, b ³ 0, d i ³ 0. (19)

Это утверждение сразу следует из явного представления коэффициентов (5), (6).

Воспользовавшись (19), оценим слагаемые, входящие в левую часть тождества (18), следующим образом:

Тогда придем к неравенству

(20)

Оценим сверху правую часть этого неравенства. Будем иметь

Подставляя эту оценку в (20) и учитывая неравенство (17), получим

(21)

Посколькуиз неравенства следует,

что погрешность z i = y i – u(x i ) также является величиной O(h 2 ) при h ® 0. Итак, справедливо следующее утверждение.

Пусть k(x) – непрерывно дифференцируемая и q(x), f(x) – непрерывные функции при x Î [0, l], решение u(x) задачи (1), (2) обладает непрерывными четвертыми производными. Пусть коэффициенты разностной схемы (3), (4) удовлетворяют условиям (8), (9), (19). Тогда решение разностной задачи (3), (4) сходится при h ® 0 к решению исходной дифференциальной задачи (1), (2) со вторым порядком по h, так что выполняется оценка

где M – постоянная, не зависящая от h.

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

3.1. Исходная задача.

Будем рассматривать следующую первую краевую задачу для уравнения теплопроводности с постоянными коэффициентами. В области <0 £ T>требуется найти решение уравнения

(1)

удовлетворяющее начальному условию

и граничным условиям

u(0, t) = m 1 (t), u(1, t) = m 2 (t). (3)

Здесь u0(x), m 1 (t), m 2 (t) – заданные функции. Известно, что при определенных предположениях гладкости решение задачи (1)–(3) существует и единственно. В дальнейшем при исследовании аппроксимации разностных схем будем предполагать, что решение u(x, t) обладает необходимым по ходу изложения числом производных по x и по t. Решение задачи (1) – (3) удовлетворяет принципу максимума и тем самым непрерывно зависит от начальных и граничных данных.

3.2. Явная схема.

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

и сетку по переменному t с шагом t , которую обозначим

Точки (x i , t n ), i = 0, 1,…, N, n = 0, 1,…, K , образуют узлы пространственно-временной сетки w h, t = w h x w t . Узлы (x i , t n ) , принадлежащие отрезкам I 0 = <0 £ x £ 1, t = 0>, I 1 = , I 2 = , называются граничными узлами сетки w h, t , а остальные узлы – внутренними. На рисунке граничные узлы обозначены крестиками, а внутренние – кружочками.

Слоем называется множество всех узлов сетки w h, t , имеющих одну и ту же временную координату. Так, n-м слоем называется множество узлов

(x 0 , t n ), (x 1 , t n ),…, (x N , t n ) .

Для функции y(x, t) , определенной на сетке w h, t , введем обозначения y n i = y(x i , t n ) ,

(4)

Иногда для упрощения записи индексы i и n будем опускать, обозначая

Чтобы аппроксимировать уравнение (1) в точке (x i , t n ), введем шаблон, изображенный на рисунке и состоящий из четырех узлов (x i ± 1 , t n ), (x i , t n ), (x i , t n+1 ). Производную ¶ u/ ¶ t заменим в точке (x i , t n ) разностным отношением y n t, i , а производную ¶ 2 u/ ¶ 2 x – второй разностной производной y n xx, i . Правую часть f(x, t) заменим приближенно сеточной функцией j n i , в качестве j n i можно взять одно из следующих выражений:

В результате получим разносное уравнение

(5)

которое аппроксимирует исходное дифференциальное уравнение в точке (x i , t n ) с первым порядком по t и вторым порядком по h при условии, что разность j n i – f(x i , t n ) имеет тот же порядок малости.

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

(6)

Эта схема представляет собой систему линейных алгебраических уравнений с числом уравнений, равным числу неизвестных. Находить решение такой системы следует по слоям. Решение на нулевом слое задано начальными условиями y 0 i = u 0 (x i ), i = 0, 1,…, N . Если решение y n i , i = 0, 1,…, N , на слое n уже найдено, то решение y i n+1 на слое n+1 находится по явной формуле

(7)

а значениядоопределяются из граничных

условий. По этой причине схема (6) называется явной разностной схемой. Несколько позже мы познакомимся и с неявными схемами, в которых для нахождения y i n+1 при заданных y i n требуется решать систему уравнений.

Погрешность разностной схемы (6) определяется как разность z i n = y i n – u(x i , t n ) между решением задачи (6) и решением исходной задачи (1) – (3). Подставляя в (6) y i n = z i n + u(x i , t n ) , получим уравнение для погрешности

(8)

где – погрешность аппроксимации разностной

схемы (6) на решении задачи (1) – (3), y i n = O( t + h 2 ) . Можно оценить решение z i n уравнения (8) через правую часть y i n и доказать тем самым сходимость разностной схемы (6) с первым порядком по t и вторым – по h. Однако это исследование мы отложим, а сейчас на примере схемы (6) продемонстрируем один распространенный прием исследования разностных схем с постоянными коэффициентами, называемый методом гармоник . Хотя данный метод не является достаточно обоснованным, в частности не учитывает влияния граничных условий и правых частей, он позволяет легко найти необходимые условия устойчивости и сходимости разностных схем. Покажем, например, что явную схему (6) можно применять лишь при условии t £ 0,5h 2 , означающем, что шаг по времени надо брать достаточно малым.

(9)

т.е. однородное уравнение, соответствующее (5). Будем искать частные решения (9), имеющие вид

y j n ( j ) = q n e ijh j ,(10)

где i – мнимая единица, j – любое действительное число и q – число, подлежащее определению. Подставляя (10) в уравнение (9) и сокращая на e ijh j , получим

(11)

Начальные условиясоответствующие решениям вида (10) (их называют гармониками), ограничены. Если для некоторого j множитель q станет по модулю больше единицы, то решение вида (10) будет неограниченно возрастать при n ® ¥ . В этом случае разностное уравнение (9) называется неустойчивым, поскольку нарушается непрерывная зависимость его решения от начальных условий. Если же |q| £ 1 для всех действительных j , то все решения вида (10) ограничены при любом n и разностное уравнение (9) называется устойчивым. В случае неустойчивости найти решение разностной задачи (6) по формулам (7) практически невозможно, так как погрешности (например погрешности округления), внесенные в начальный момент времени, будут неограниченно возрастать при увеличении n. Такие разностные схемы называются неустойчивыми.

Для уравнения (9) неравенство |q| £ 1 выполняется согласно (11) при всех j тогда и только тогда, когда g £ 0,5. Таким образом, использование схемы (6) возможно лишь при выполнении условия t £ 0,5h 2 . Разностные схемы, устойчивые лишь при некотором ограничении на отношение шагов по пространству и по времени, называются условно устойчивыми. Следовательно, схема (6) возможно устойчива, причем условие устойчивости имеет вид t /h 2 £ 0,5. Условно устойчивые схемы для уравнений параболического типа используются редко, так как они накладывают слишком сильное ограничение на шаг по времени. Действительно, пусть, например, h = 10 -2 . Тогда шаг t не должен превосходить 0,5 * 10 -4 , и для того чтобы вычислить решение y j n при t = 1 , надо взять число шагов по времени n = t -1 ³ 2 * 10 4 , т.е. провести не менее 2 * 10 4 вычислений по формулам (7).

3.3. Неявные схемы.

Чисто неявной разностной схемой для уравнения теплопроводности теплопроводности (схемой с опережением) называется разностная схема, использующая шаблон (x i , t n ), (x i ± 1 , t n+1 ), (x i , t n+1 ) и имеющая вид

(12)

Здесь j n i = f(x i , t n+1 ) + O( t + h 2 ) . Схема имеет первый порядок аппроксимации по t и второй – по h. Решение системы (12) находится, как и в случае явной схемы, по слоям, начиная с n = 1. Однако, теперь, в отличие от явной схемы, для нахождения y i n+1 по известным y i n требуется решить систему уравнений

(13)

где g = t /h 2 , F i n = y i n + t j i n . Эту систему можно решать методом прогонки, так как условия устойчивости прогонки выполнены.

Для исследования устойчивости разностной схемы (12) будем искать частные решения уравнения

имеющие вид (10). Тогда получим

следовательно, |q| £ 1 при любых j , t , h . Таким образом, схема (12) абсолютно устойчива, т.е. устойчива при любых шагах t и h . Абсолютная устойчивость является основным условием неявных схем. Теперь уже не надо брать шаг t слишком малым, можно взять, например, t = h = 10 -2 . Величина шагов сетки t , h определяются теперь необходимой точностью расчета, а не соображениями устойчивости.

Шеститочечной симметричной схемой называется разностная схема

(14)

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

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

(15)

При s = 0 получим отсюда явную схему, при s = 1 – чисто неявную схему и при s = 0,5 – симметричную схему (14). Исследуем погрешность аппроксимации схемы (15) на решении исходной задачи (1) – (3). Представим решение задачи (15) в виде y i n = u(x i , t n ) + z i n , где u(x i , t n ) – точное решение дифференциальной задачи (1) – (3). Тогда для погрешности получим систему уравнений

(16)

i = 1, 2,…, N – 1, n = 0, 1,…, K – 1,

z 0 n+1 = z N n+1 = 0, n = 0, 1,…, K – 1, z i 0 = 0, i = 0, 1,…, N.

Сеточная функция y i n , входящая в правую часть уравнения (16) и равная

(17)

называется погрешностью аппроксимации схемы (15) на решении задачи (1) – (3). Получим первые члены разложения функции y i n по степеням h и t . Будем разлагать все функции, входящие в выражение для y i n , по формуле Тейлора в точке (x i , t n + 0,5 t ). Учитывая разложения

Отсюда, проводя разложение в точке (x i , t n+1/2 ) и обозначая u = u (x i , t n+1/2 ) , будем иметь

и, перегруппировывая слагаемые, получим, что

Учитывая уравнение (1) u’’ – u = – f и следствие из него u IV – u’’ = –f’’ , окончательно можно записать, что

(18)

Из формулы (18) можно сделать следующие выводы. Если

то схема (15) имеет второй порядок аппроксимации по t и четвертый – по h . Такая схема называется схемой повышенного порядка аппроксимации. Если

то схема (15) имеет второй порядок аппроксимации по t и по h. При остальных значениях s и при j i n º 0 в виде (10), то получим

и |q| £ 1 при всех j , если

(19)

Отсюда видно, в частности, что все схемы с s ³ 0,5 абсолютно устойчивы. Схема повышенного порядка аппроксимации ( s = s * ) также абсолютно устойчива, что проверяется непосредственно.

При s ¹ 0 разностная схема (15) является неявной схемой. Для нахождения решения y i n+1 по заданным y i n требуется решать систему уравнений

(20)

где

Система (20) решается методом прогонки. Условия устойчивости прогонки при s ¹ 0 сводятся к неравенству

|1 + 2 s g | ³ 2 | s | g

и выполнены при s ³ – 1/(4 g ). Последнее неравенство следует из условия устойчивости (19) разностной схемы.

3.4. Уравнения с переменными коэффициентами и линейные уравнения.

Рассмотрим первую краевую задачу для уравнения теплопроводности с переменными коэффициентами

(21)

где r (x, t), k(x, t), f(x, t) – достаточно гладкие функции, удовлетворяющие условиям

0 1 £ k(x, t) £ c 2 , r (x, t) ³ c 3 > 0 .(22)

Дифференциальное выражениепри каждом

фиксированном t аппроксимируем в точке (x i , t) так же, как и в стационарном случае, разностным отношением

где разностный коэффициент теплопроводности a(x i , t) должен удовлетворять условиям второго порядка аппроксимации

Наиболее употребительны следующие выражения для a(x i , t) :

Разностная схема с весами для задачи (21) имеет вид

(24)

Здесь в качестве t можно взять любое значение t Î [t n , t n+1 ] , например t = t n + 0,5 t . Если в уравнении (24) t = t n + 0,5 t , s = 0,5 , то схема (24) имеет второй порядок аппроксимации по t и по h . При остальных значениях s и t выполняется первый порядок аппроксимации по t и второй – по h .

При исследовании устойчивости разностных схем с переменными коэффициентами иногда применяется принцип замороженных коэффициентов, сводящий задачу к уравнению с постоянными коэффициентами. Рассмотрим явную схему, соответствующую уравнению (24) с s = 0 и f(x i , t) º 0 , т.е. схему

(25)

Предположим, что коэффициенты r (x i , t), a(x i , t) – постоянные, r (x i , t) º r = const, a(x i , t) º a = const . Тогда уравнение (25) можно записать в виде

или

Из п.2 известно, что последнее уравнение устойчиво при t ’ £ 0,5h 2 , т.е. при

(26)

Принцип замороженных коэффициентов утверждает, что схема (25) устойчива, если условие (26) выполнено при всех допустимых значениях a(x i , t), r (x i , t) , т.е. если при всех x, t выполнены неравенства

(27)

Если известно, что 0 1 £ a(x i , t) £ c 2 , r (x i , t) ³ c 3 > 0 , то неравенство (27) будет выполнено при

Строгое обоснование устойчивости схемы (25) будет дано в примере 2 из главы 2.

Если параметр s ³ 0,5, то из принципа замороженных коэффициентов следует абсолютная устойчивость схемы (24).

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

(28)

В случае нелинейных уравнений, когда заранее неизвестны пределы изменения функции k(u) , избегают пользоваться явными схемами. Чисто неявная схема, линейная относительно y i n+2 , i = 1, 2,…, N – 1 , имеет вид

(29)

где a i = 0,5 (k(y n i ) + k(y n i-1 )) . Эта схема абсолютно устойчива, имеет первый порядок аппроксимации по t и второй – по h . Решение y i n+1 , i = 1, 2,…, N – 1 , находится методом прогонки. Заметим, что схему (29) можно записать в виде

где k i = k(y i n ) .

Часто используется нелинейная схема

(30)

Для реализации этой схемы необходимо применить тот или иной итерационный метод. Например такой:

(31)

Здесь s – номер итерации. Как видим, нелинейные коэффициенты берутся с предыдущей итерации, а в качестве начального приближения для y i n+1 выбирается y i n . Это начальное приближение тем лучше, чем меньше шаг t . Число итераций M задается из соображений точности. В задачах с гладкими коэффициентами при k(u) ³ c 1 > 0 часто бывает достаточно провести две – три итерации. Значения y i (S+1) на новой итерации находятся из системы (31) методом прогонки. При M = 1 итерационный метод (31) совпадает с разностной схемой (29).

Для приближенного решения нелинейного уравнения (28) применяются также схемы предиктор – корректор второго порядка точности, аналогичные методу Рунге – Кутта для обыкновенных дифференциальных уравнений. Здесь переход со слоя n на слой n+1 осуществляется в два этапа. Приведем пример такой схемы. На первом этапе решается неявная линейная система уравнений

из которой находятся промежуточные значения y i n+1/2 , i = 0, 1,…, N . Затем на втором этапе используется симметричная шеститочечная схема для уравнения (28), в которой нелинейные коэффициенты a(y), f(y) вычисляются при y = y i n+1/2 , т.е. схема

Основы численных методов для газодинамики с теплопроводностью (стр. 5 )

Из за большого объема этот материал размещен на нескольких страницах: 1 2 3 4 5 6 7 8

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

Здесь k — некоторый коэффициент. В книге Зельдовича-Райзера [[13]] показано, что это предположение оправдано для случая, когда пробег излучения мал по сравнению с размерами области, в которой ведется расчет. При этом

С учетом указанных упрощений получаются уравнения спектральной неравновесной диффузии:

Дальнейшее упрощение связано с усреднением функций и коэффициентов по частоте (так называемое «серое вещество») и предположением о равновесности процесса, что позволяет исключить из рассмотрения плотность излучения и считать, что поток энергии определяется градиентом плотности равновесного излучения, или с точностью до множителя градиентом температуры. В результате получается уравнение теплопроводности:

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

В таком виде уравнение теплопроводности было использовано во всех основных программах, созданных в 50-60х годах и применявшихся для расчета современных для того времени изделий. Ситуация, как видим, опять напоминает анекдот о пьяном матросе. Однако, в данном случае проблема гораздо серьезнее. При счете детонации приближенная модель макрокинетики при надлежащем подборе параметров достаточно точно описывает штатные режимы. Уточненная двухстадийная модель позволяет рассчитывать процессы возбуждения и затухания детонации. В целом это обеспечивает все основные потребности разработчиков типовых изделий. Иначе обстоит дело с теплопроводностью. Распространение тепла в узком зазоре, ширина которого сравнима с пробегом излучения, является одной из наиболее важных задач и именно она плохо описывается диффузионным приближением. Очень серьезное внимание на это обращал , в свое время он пытался использовать в качестве некоторой поправки так называемый «геометрический пробег». Известно также, что диффузионное приближение плохо передает форму фронта тепловой волны и поток энергии на границе с вакуумом. Усреднение пробегов по частоте заметно влияет на динамику процессов. Заметим, что все эти недостатки типичны также для уравнения переноса нейтронов. Но для нейтронов отказ от диффузионного приближения произошел достаточно быстро, уже в конце 50х годов существовал целый ряд различных методик счета многогруппового (по спектру) уравнения переноса. Для переноса излучения процесс внедрения заметно затянулся, так как задача оказалась более сложной.

Формально диффузионные уравнения для нейтронов и для фотонов выглядят одинаково и относятся к так называемым параболическим уравнениям. Это означает, в частности, что у них нет характеристик и скорость сигнала является бесконечной. Я помню, как при анализе одного из расчетов объяснял мне это таким образом: «Говорят, что если бросить щепотку соли в Одессе, то в Рио-де-Жанейро вода в тот же момент станет чуточку солонее». Тем не менее, есть внутренняя динамика процессов, связанная с величиной коэффициентов уравнений и приводящая к существенному отличию в расчетах. Забегая вперед, скажу, что разностные схемы, используемые для счета диффузии тепла, предъявляют существенно более высокие требования к памяти вычислительных машин и, отчасти, к их быстродействию. Аналогичные различия сохраняются и, более того, многократно возрастают при переходе к уравнению переноса. Это собственно и определило невозможность расчета уравнений переноса излучения на машинах М-20 и БЭСМ-6 с их памятью в несколько десятков килобайт.

Тем не менее, некоторые шаги в этом направлении делались. В частности, мы с в начале 70х годов рассмотрели интегральный метод, обобщающий метод и , реализованный в И-0. Метод оказался достаточно сложным и программно был реализован тогда только для плоского одномерного случая, результаты были изложены в одном из отчетов ВНИИЭФ и в сокращенном виде в моей книге [[14]]. Некоорых результатов удалось добиться также в методе расчета спектральной неравновесной диффузии, результаты того периода также изложены в одном из отчетов и, также сокращенно, в книге [[15]]. В дальнейшем, насколько мне известно, со своими сотрудниками успешно развивает работы по спектральной неравновесной диффузии.

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

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

Через некоторое время матрос попросился обратно на корабль. Возможно, что этому матросу, осознанию им своего места, мы обязаны тем, что имеем сегодня труды , , … Мне, к сожалению, с «матросом» повезло меньше и занятия переносом излучения пришлось оставить.

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

Решение будем искать в виде волны, двигающейся с постоянной скоростью V, для чего введем автомодельную переменную x=x-Vt. Тогда уравнение примет вид

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

Получили уравнение с разделяющимися переменными, при его интегрировании снова возьмем константу равной нулю. Кроме того, выразим T через x и подставим x=x-Vt. Полученное решение имеет вид:

Решение определено в области x£Vt, в точке x=Vt температура обращается в ноль, эту точку будем считать фронтом волны, идущей по нулевому фону с конечной скоростью. Если задача рассматривается в ограниченной области, то на границе левее фронта нужно задать температуру, растущую по определенному закону. При другом граничном условии скорость волны будет конечной, но не постоянной. При k>1, что реально соответствует физике процесса, производная температуры обращается в точке фронта в бесконечность. Это вызывает определенные проблемы при построении разностных схем, однако, такая особенность связана не со свойствами тепловой волны, а с особенностями диффузионного уравнения. При счете уравнения переноса на фронте волны наблюдается острый носик.