1.5. Графические возможности скм

В СКМ Mathematica имеется несколько функций, позволяющих строить различные графические объекты. Простые графики позволяет строить функция Plot. Например команда

Plot [Sin[X], {X, 0, 2*Pi}]

строит график функции sinx в интервале от 0 до 2. Приведем пример создания и отображения функции пользователя. При расчете распределения потенциала в двумерной области разностный аналог уравнения Пуассона решается итерационным методом последовательной верхней релаксации (ПВР). На каждой итерации ПВР рассчитывается максимальная поправка потенциала и критерием завершения итераций является достаточная малость последней. По окончании расчета следует построить график зависимости максимальной поправки от номера итерации, а также трехмерное распределение потенциала.

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

9

deltaPhi = { };

Внутри итерационного цикла в его конце будем сохранять величину поправки, содержащейся в переменной deltaAbs

deltaPhi = Append[ deltaPhi, deltaAbs ];

Отметим, что здесь сохраняется только поправка, но не номер итерации.

По окончании итерационного процесса массив deltaPhi содержит требуемый набор чисел. Однако функции Plot необходимо передать непрерывную (в данном случае кусочно-ломанную) линию. Поэтому сначала проведем линейную интерполяцию имеющегося набора чисел, получим объект (функцию пользователя), который затем передадим в функцию Plot:

fun[x_] = Interpolation[ deltaPhi, InterpolationOrder -> 1] [ x ];

Plot[ fun[x], {x, 1, kolit }];

Параметр InterpolationOrder -> 1 определяет уровень интерполяции (самая простая, линейная), переменная kolit – номер последней итерации ПВР.

Рассчитанное двумерное распределение потенциала отобразим функцией

ListPlot3d[ phi ];

2. Задания для самостоятельной работы

2.1. Расчет двумерного распределения температуры

В нестационарном случае распределение температуры в двумерной области при отсутствии внутри этой области стоков и источников тепла описывается уравнением теплопроводности [2]

(1)

где T – температура в точке (x, y), t – время, – двумерный оператор дифференцирования, a – температуропроводность материла, заполняющего область. Запись означает операцию дивергенции над вектором, заключенным в скобки.

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

(2)

(3)

10

(4)

Граничное условие (2) означает, что в точках области температура известна. Граничное условие (3) описывает плотность температурного потока в точках области . Граничное условие (4) свидетельствует о том, что в точках границы имеют место более сложные процессы взаимодействия с внешней средой. Здесь , , - участки границы , – вектор единичной нормали к ней, – заданные функции.

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

. (5)

Функция , стоящая в правой части (5), описывает распределение температуры в области в момент времени .

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

Перейдем к интегральной форме уравнения теплопроводности. Для этого введем плотность температурного потока

. (6)

Перепишем уравнение (1) в виде

. (7)

Проинтегрируем левую и правую части (7) по произвольной двумерной области :

Применим к правой части теорему Остроградского-Гаусса:

, (8)

11

где – внешняя нормаль к контуру , ограничивающему область . Уравнение (8) является интегральной формой уравнения (1). Оно содержит, как это видно после подстановки (6), лишь первую производную по пространственным координатам.

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

Интегральное тождество (8) справедливо, что очень важно, не только для всех точек области , но и для всех точек границы с неизвестной температурой T. Граничные условия (3), (4) учитываются в (8) путем задания соответствующей плотности температурного потока на участках :

(9)

. (10)

Таким образом, интегральная формулировка обсуждаемой краевой задачи включает интегральное тождество (8), граничные условия (2), (9), (10), начальное условие (5).

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

- содержит лишь первую производную по пространственным переменным, что увеличивает устойчивость разностной схемы;

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

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

Рассмотрим на примере, как на основе интегральной формулировки строится и численно решается разностная схема для краевой задачи (1) ÷ (5).

Требуется рассчитать стационарное распределение температуры в прямоугольной области ABCD (см. рис. 1). Область заполнена однородным материалом, температура его двух граничных участков фиксирована: в точках

12

участка AE температура равна , в точках участка CF температура равна , причем . Остальная часть границы теплоизолирована, т. е. через эти участки отсутствует передача тепла. Последнее означает, что

или

где – проекции вектора на оси и соответственно.

Рис. 1. Пример области моделирования

Сформулированная задача является стационарной, а значит интегральное тождество (8) приобретает вид:

. (11)

Выполним дискретизацию задачи, для чего в области ABCD построим две сетки (см. рис. 2) – основную

с шагом по координате и шагом по координате , при этом

13

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

Рис. 2. Дискретизация области моделирования

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

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

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

14

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

Рис. 3. Ячейка алгебраизации для внутреннего узла сетки

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

Второе допущение позволяет легко аппроксимировать градиент температуры, а значит и плотность температурного потока конечно-разностным соотношением. В частности, плотность температурного потока, например, на левой грани прямоугольной ячейки

Полный поток на этой грани, взятый с обратным знаком (этого требует направление внешней нормали, которая здесь направлена справа налево), соответственно равен

15

. (12 )

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

В итоге разностное уравнение для внутреннего узла сетки можно записать, как

(13)

где

, ,

, , .

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

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

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

16

(его номер показан римской цифрой в левом нижнем углу рисунка), равна a₁, а материала, заполняющего 4-й квадрант – a₂, то полный поток через левую грань ячейки алгебраизации соответственно равен

Этот поток можно записать также в форме (12), с усредненным коэффициентом температуропроводности

.

Пусть узел сетки принадлежит верхней границе прямоугольной области (рис. 4). Ячейка алгебраизации теперь включает лишь два квадранта: II и III. Аппроксимация интегрального тождества (11) на такой ячейке алгебраизации выполняется, как и прежде, путем аппроксимации и суммирования потоков через грани ячейки. При этом температурный поток через верхнюю грань ячейки полагается равным нулю (это следует из граничного условия).

Рис. 4. Ячейка алгебраизации для узла, принадлежащего верхней границе

Конечно-разностная аппроксимация интегрального тождества (11) для таких узлов сетки

17

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

(14)

где

, , ,

.

Разностное уравнение содержит в этом случае температуру в четырех узлах основной сетки. (Коэффициент при оказывается равным нулю.)

Рис. 5. Ячейка алгебраизации для углового граничного узла

Рассмотрим узел сетки, лежащий, в частности, в правом нижнем углу области (см. рис. 5). Его ячейка алгебраизации состоит из одного квадранта. По этой причине конечно-разностная аппроксимация интегрального тождества (11) на такой ячейке заметно упрощается:

.

18

Разностное уравнение , (15)

где

,

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

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

Разностная СЛАУ обладает важным свойством – сильной разреженностью матрицы (в каждом уравнении ненулевыми являются не более пяти коэффициентов). Кроме того, матрица системы, как это следует из расчетных соотношений для , обладает свойством диагонального преобладания. Для решения таких систем можно использовать достаточно простой итерационный метод последовательной верхней релаксации (ПВР) [3].

Перепишем (13) в виде

. (16)

Алгоритм ПВР представляет собой последовательное применение формулы

, (17)

в которой

, (18)

19

 – параметр релаксации, 1   < 2. Формула (18) есть не что иное, как реализация метода Гаусса-Зейделя для разностного уравнения (16).

Для граничных узлов формула (17) сохраняет свой вид, видоизменяется только соотношение (18). В частности, из разностного уравнения (14)

, (19)

из разностного уравнения (15)

. (20)

Аналогичным образом обрабатываются остальные участки границы. Простота формул (17)(20) позволяет напрямую запрограммировать их, не создавая промежуточных массивов для хранения матрицы СЛАУ и вектора правой части. При этом следует организовать обход всех узлов области, уточняя размещение каждого узла и с учетом этого используя подходящую формулу.

Таблица 1

№ варианта	AD	AB	AE	FC	T1	T2
1	10	8	6	6	100	40
2	10	8	6	4	50	150
3	15	8	4	4	140	300
4	15	12	8	10	180	60
5	20	12	10	7	160	350
6	20	12	10	5	320	140
7	25	16	12	8	120	340
8	25	16	12	4	45	15
9	30	16	14	12	80	190
10	30	16	14	16	170	80

Контроль точности решения можно выполнить, вычисляя коэффициент

,

где F₁ – поток температуры, втекающий в область через контакт 1, F₂ – поток температуры, вытекающий из области через контакт 2. Так как внутри области моделирования отсутствуют стоки и источники тепла, алгебраическая сумма потоков тепла, текущих через внешние контакты, должна быть равна нулю. Коэффициент kons показывает, в какой мере решение задачи

20

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

Задание: рассчитать распределение температуры в области, приведенной на рисунке 1, для вариантов исходных данных таблицы 1, a = 25.