CC BY
317
УДК 536.24 (075.8)
в.в. рындин решение разностных уравнений теплопроводности с использованием системы MATHCAD
Вводная часть. В настоящее время при решении математических задач широко используется программирование в средах Fortran, Turbo Pascal, Delphi и др. При этом для выполнения даже небольших математических расчётов требуется знание основ программирования. При написании формул теряется их наглядность. Например, на языке Pascal 4x записывается
как sqrt(x), степень yx как exp(x * ln(y)) и т. п. Каждый раз при выводе на печать результатов расчёта по какоё-либо формуле требуется давать специальное сообщение, а для выдачи графиков требуется написание самим пользователем специальных программ. Этих недостатков лишена новая математическая система Mathcad.
Mathcad - это интегрированная математическая система, позволяющая наглядно вводить исходные данные, проводить традиционное математическое описание решения задачи и получать результаты вычислений как в аналитическом, так и в численном виде с использованием при необходимости их графического представления. Запись математических выражений производится в традиционном виде с применением общепринятых знаков, таких как квадратный корень, знак деления в виде горизонтальной черты, знак интеграла, дифференциала, суммы и т. д.
Эта система имеет хорошо продуманные встроенные текстовый, формульный и графический редакторы. Они снабжены удобным пользовательским интерфейсом, обладают разнообразными математическими возможностями и ориентированы на нужды большинства пользователей - школьников, студентов, инженеров, экономистов, менеджеров, научных работников.
В Mathcad предусмотрен импорт любых графических изображений - от простых и специальных графиков функций до многокрасочных репродукций художественных произведений. Введены средства анимации рисунков и проигрывания видео файлов со звуковым стереофоническим сопровождением.
Mathcad - настолько мощные и гибкие системы, что вполне заслуживают описания во многих книгах. За рубежом им посвящены уже сотни книг. Наименования некоторых книг на русском языке даны в конце статьи [1 - 5]. Однако до настоящего времени эта система не нашла повсеместного
75
применения в учебной практике вузов при выполнении различных расчётов. Это обусловлено отсутствием как опыта работы преподавателей в системе Mathcad, так и соответствующих примеров расчёта в этой системе. В качестве некоторых примеров применения системы Mathcad в учебном процессе ПГУ им. С. Торайгырова можно привести работы [6 - 9].
Несмотря на то, что эта программа в основном ориентирована на пользователей-непрограммистов, Mathcad также используется в сложных проектах, чтобы визуализировать результаты математического моделирования.
В данной статье в качестве примера научного использования программы Mathcad приводится решение дифференциальных уравнений нестационарной теплопроводности конечно-разностным методом с использованием графической интерпретации получаемых температурных зависимостей для различных моментов времени. Цель данной работы - привлечь внимание преподавателей и аспирантов к использованию системы Mathcad при выполнении научных и исследовательских работ.
Суть метода конечных разностей. Аналитическое решение задач теплопроводности получено лишь для тел простой геометрической формы. Для тел сложной формы или при сложных краевых условиях уравнение теплопроводности не всегда возможно решить аналитически.
В связи с развитием вычислительной техники большое распространение получил конечно-разностный метод решения задач нестационарной теплопроводности, или метод сеток. В основе этого метода лежит допущение о возможности замены бесконечно малых изменений температуры во времени и в пространстве малыми, но конечными её изменениями. Тем самым непрерывный процесс изменения температуры в теле при его нагревании или охлаждении заменяется скачкообразным процессом.
В качестве примера рассмотрим одномерное уравнение теплопроводности Фурье [10]
dT ö2T
= а—ö" dt до2
Суть метода сеток заключается в покрытии расчётной области (x, t) сеткой из M х N точек. Тем самым определяются узлы, в которых будет осуществляться поиск решения. Конфигурацию узлов, используемую для разностной записи уравнений в частных производных на сетке, называют шаблоном, а полученную систему разностных уравнений - разностной схемой. На рисунке 1 изображены шаблоны аппроксимации явной схемы (а) и неявной схемы (б) для уравнения теплопроводности.
76
(1)
п+1
t
^п+2 ^п+1
и
а) а)
Рисунок 1 - Шаблоны аппроксимации явной схемы (а) и неявной схемы (б) для уравнения теплопроводности
Затем дифференциальное уравнение теплопроводности для одномерного нестационарного температурного поля (1) заменяется уравнением в конечных разностях для каждого (т'п) -го узла сетки
А2Т
Щ
—= а-
А Ад
(2)
где Аt е Ах - шаги соответственно по времени и по координате; а - температуропроводность вещества, } 2 /п.
В соответствии с шаблоном для явной схемы (см. рисунок 1, а) частная производная по времени определится выражением
АТ{ = Тт,п+1— Тт,п
А ~ Аt ,
а частные производные по координатам соответственно справа (плюс) и слева (минус)
[ Т 1
Тт+1,п Тп Ах
АТх Ах
Тт,п Тт-1,п Ах
Для второй производной в конечных разностях получим следующее выражение:
А 2ТХ Ад2
Ах
АТх Ах
+
АТх
Ах
1
Ах^
(Тт+1,п 2Тт,п + Тт- 1,п )
После подстановки соответствующих выражений в уравнение (2) получим
Тт,п+1—Тт,п = а 2Т + Т )
Аt
Ах 2
(3)
77
Приведём в разностной схеме (3) подобные слагаемые, перенеся в правую часть значения сеточной функции с индексом п (как часто говорят, с предыдущего слоя по времени), а в левую - с индексом п + 1 (т. е. со следующего временного слоя). Кроме того, введём так называемый коэффициент (число) Куранта к для шагов сетки
г Л А(
* — 2йа7 . (4)
С учётом этих замечаний, разностная схема (3) запишется в виде
Тт,п+1 - * Тт-1,п +(1 - *) Тт,п + * Тт+1,п . (5)
Множители для каждого из значений сеточной функции (температуры Т) в узлах шаблона, соответствующие разностному уравнению (5), приведены рядом с каждой точкой шаблона на рисунке 1.
Для получения замкнутой системы разностных алгебраических уравнений систему (5) необходимо дополнить дискретным представлением начального условия
Т (х,0)=То (х)
и граничных условий
Т (0, г )=о (г), Т (А г )=/! (г).
В соответствии с соотношением (5) каждое неизвестное значение сеточной функции со следующего временного слоя, т. е. левая часть этого соотношения явно выражается через три её значения с предыдущего слоя (правая часть), которые уже известны. Из-за того, что разностная схема сводится к такой явной подстановке, её и называют явной.
В отличие от явной разностной схемы в неявной разностной схеме (см. шаблон на рисунке 1, б) для дискретизации пространственной производной берутся значения сеточной функции с верхнего (неизвестного)
слоя по времени. Таким образом, разностное уравнение для (т'п) -го узла будет отличаться от уравнения для явной схемы (3) только индексами по временной координате в правой части:
Тт,п+1 -Тт,п — а т 2Т Т )
ТГ — 2(Тт-1,п +1 - 2Тт,п +1 + Тт+1,п+1) Аг Ах2
78
Теперь надо выразить искомые Tm,n+1 (на следующем шаге, или слое, по времени) через Tm,n . В отличие от явной схемы Эйлера, неизвестные на новом слое связаны между собой в линейное уравнение Поэтому для нахождения массива Т для каждого нового временного слоя требуется решить систему уравнений. Выпишем её для неявной схемы в привычном виде (приводя подобные члены)
0,5KTm-1,n +1 -(K +1)Tm,n+1 + 0,5KTm+1,n +1 = -Tm,n . (6)
Множители при неизвестных значениях сеточной функции в узлах шаблона показаны на рисунке 1, б в виде подписей, подобно тому, как это было сделано для явной схемы (см. рисунок 1, а).
Пример расчёта в системе Mathcad. В качестве примера применения метода конечных разностей рассмотрим случай остывания неравномерно нагретого стержня, на концах (границах) которого поддерживается постоянная температура T^=const Для упрощения математических выкладок (для установления интервала изменения величин по оси ординат и абсцисс от нуля до единицы и использования стандартных программ) принимаем длину стержня Ь=1м, а в качестве искомой величины возьмём относительную избыточную температуру
О
и.-/-, ?
где Т() - начальная температура в центральных областях стержня.
Подставляя выражение для температуры ^ '< ~ ч ^ ^ ч- в (5) и приводя подобные члены, получим явную разностную схему в таком виде:
®m,n +1 = K®m-1,n +(! - K)®m,n + у ®m+1,n . (7)
Явная схема Эйлера. Ниже приведено решение разностного уравнения (7) в системе Mathcad. Заметим, что всё ниже написанное, включая и комментарии, может составлять содержание программы расчёта - программа сама определяет, где текст, а где математические выражения (последние для наглядности будем выделять жирным шрифтом).
Задаём одинаковые граничные условия постоянства температуры на краях стержня (условие Дирихле) Border(t):=0.
Задаём начальное распределение температуры в виде прямоугольника (нагретая центральная часть стержня, см. рисунок 2) с помощью функции
79
Хевисайда Ф(х) (даёт ноль при х < 0 и единицу в противном случае ), например, в таком виде:
©0(д) ;=Ф(х-0.45}-Ф(х-0.55)
(х<0,45^0-0 = 0; 0,45 < х <0,55->1-0 = 1; 0,55 < х -> 1-1 = 0).
Температуропроводность стержня а:=0.01 м 2/с.
Длина стержня м.
Число точек по оси х М:=20 .
Пространственный шаг А: = М = 0.05 м.
Дискретизация начального условия на сетке - переход от непрерывной функции ©0(х) к дискретной функции ©0т, т. е. создание массива (рисунок 3)
т:=0..М, 00т := ©0(т А)
Шаг по времени определяем из условия, что явная разностная схема является устойчивой, если число Куранта (4) меньше единицы. Принимаем, например, К:=0.4 , тогда шаг по времени в:=0.05 с.
Рисунок 2 - Начальное распределение температуры по длине стержня Рисунок 3 - Дискретизация начального условия на сетке
Возможное программное решение разностной схемы приводится ниже в виде программного модуля, причём функция пользователя 0(п) задаёт вектор распределения искомой относительной температуры в каждый момент времени (иными словами, на каждом временном слое), задаваемом целым числом п.
Основными инструментами работы в Mathcad являются математические выражения и функции, записываемые в одну строку, что позволяет создавать лишь линейные программы, то есть осуществляющие последовательные вычисления от начала к концу программы. Для осуществления циклов, а также задания переменных и функций, записываемых в несколько строк (например, возвращение различных значений в за-
80
висимости от условий), используются специальные программные блоки (модули). Программный модуль в тексте документа выделяется жирной вертикальной линией. Программные элементы, входящие в программный блок, выбираются с помощью панели инструментов Программирование (Programming). В качестве основных программных элементов можно отметить следующие (имена программных операторов не следует вводить с клавиатуры):
Add Line (добавить линию), - создаёт и при необходимости расширяет жирную вертикальную линию, справа от которой в шаблонах (местоза-полнителях) задаётся запись программного блока;
^ - символ локального присваивания (значение величины задаётся и определяется только в теле модуля);
if (если) - оператор условного выражения (если условие выполняется, то выполняется выражение, стоящее перед if );
otherwise (иначе) - оператор иного выбора (обычно применяется с if); for (для) - оператор задания цикла с фиксированным числом повторений (с шагом +1 от указанного начального значения до конечного значения).
Программный модуль решения может быть записан в таком виде:
v0 ^ Border(r) vM ^ Border(r) for m ё 1..M -1
F (в): =
-em-i-у+em-a-K+em+x-у
v
v
(n ): =
eo if n = о
F (e( n-l)) otherwise
Начальное распределение температуры ©0m вдоль расчётной области и решение для трёх моментов времени (для трёх временных слоёв n = 1, 3 и 30) показаны на рисунке 4. Физически такое поведение вполне естественно - с течением времени тепло (хаотическое движение) из более нагретой области перетекает в менее нагретую, а зона изначально высокой температуры остывает и размывается.
81
Рисунок 4 - Начальное распределение относительной избыточной температуры и решение для трёх моментов времени
Как уже отмечалось, явная разностная схема Эйлера для задачи теплопроводности является устойчивой, если соотношение Куранта меньше 1. Если же оно превышает 1, то разностная схема неустойчива, и её ответ неверен. Так, если в начале расчёта принять число Куранта К:=1.6, то решение на первых шагах по времени имеет вид, приведённый на рисунке 5.
Из физических соображений очевидно, что это решение неверно. О неустойчивости схемы говорит характерная «разболтка» сеточного решения.
1 ■. ■
0Л1
а *|
ФП^пЙИ
«Ар 1>
1 -Г
Г,,'
■■ .■ - 5 ' -
■
"!
1
Н М И О Ш VI Ы 1
Я А
Рисунок 5 - «Разболтка» сеточного решения при соотношении Куранта больше 1 (К = 1,6)
Неявная схема Эйлера. В отличие от явной схемы Эйлера, неявная схема является безусловно-устойчивой (т. е. не выдающей «разболтки» ни при каких значениях числа Куранта). Однако ценой устойчивости является необходимость решения на каждом шаге по времени системы алгебраических уравнений (6), которая для относительной избыточной температуры 0 имеет вид
82
05K®m_} п+1 - (K + 1)©щ n+1 + 0,5K©m+J n+1 = n . (8)
Неявная схема (8) приводит к системе линейных уравнений A©n+1 = b с матрицей А и вектором правых частей b = -©n . Для реализации неявной схемы, таким образом, можно использовать комбинацию средств программирования Mathcad и встроенной функции решения системы линейных уравнений lsolve(A,b). Ниже приводится один из возможных способов решения:
Border(t):=0. ®0(х):=Ф(х4^)-Ф(х-4.5э).
a:=0.01 м 2/с; L=1 м; M:=20 . А: = M = 0 05 м;
2
K:=1.6 • в: = = 0.2 c .
' 2a '
m:=0..M ; ©0m := ©0(m А) .
m:=1..M-1; Am m := -K -1; А^.! := 0.5 • K .
Am,m+1 := 0.5'K ; A0,0 := 1; AM,M := 1 .
©(n): =
©0 if n = 0
lsolve(A,-©(n -1) otherwise
В первой строке вводятся начальные и граничные условия; во второй строке - температуропроводность, длина стержня, число точек по оси х, шаг по координате; в третьей строке - число Куранта и шаг по времени; в четвёртой строке - дискретизация начального условия на сетке (задание значений температуры в узлах сетки в начальный момент времени); в пятой и шестой строках - формируется матрица системы уравнений А, которая записывается в подходящем для Mathcad виде. В последней строке определяется функция пользователя, вычисляющая сеточную функцию на каждом временном слое (при помощи встроенной функции решения линейных уравнений Ьо1уе(А,Ь)). Как несложно убедиться, столбец правых частей разностных уравнений выражается вычисленными значениями сеточной функции с предыдущего слоя.
Результаты расчётов по неявной схеме для числа Куранта больше единицы (К = 1,6) показаны на рисунке 6.
83
©Üji
1
05
ÖS flLi
Л»»
вДООи 0 3
0 1
о
J
j__ t
f \
i 07 V3. . о
ли fc-fi-i ■L л
ü ] и üi 0 3 iJ-ч
ш А
Ofi Ü 7 0Й 04 ]
Рисунок 6 - Результаты расчёта по неявной схеме для соотношения Куранта больше единицы (К = 1,6)
Как видим, неявная схема устойчива и при K > 1. Сравнение результатов расчёта по явной и неявной расчётным схемам показывает, что при одинаковых числах Куранта имеется некоторое несовпадение для первых временных слоёв, достигающее 5 % при K = 0,4 и до 30 % при K = 0,8 . С уменьшением числа Куранта K несовпадение уменьшается.
Заключение. На примере расчёта изменения температуры по длине стержня при его остывании и представления температурных зависимостей в виде графиков для произвольных моментов времени показана возможность использования системы Mathcad не только для учебных и инженерных расчётов, но и для программирования различных физических процессов при проведении научных исследований.
СПИСОК ЛИТЕРАТУРЫ
1 Кирьянов Д.В. Mathcad 13. - Спб.: БХВ-Петербург, 2006. - 608 с.: ил.
2 В.Ф. Очков. Mathcad 14 для студентов и инженеров: русская версия. - Спб.: BHV-Петербург, 2009. - 467.
3 Кирьянов Д.В. Mathcad 15 / Mathcad Prime 1.0. В подлиннике. - Спб.: БХВ-Петербург, 2012. - 432 с. : ил.
4 Панферов А. И., Лопарев А. В., Пономарев В. К. Применение Mathcad в инженерных расчетах: Учеб. пособие / СПбГУАП. СПб., 2004. - 88 с.: ил.
5 Половко А. М., Ганичев И. В. Mathcad для студента. - СПб.: БХВ-Петербург, 2006. - 336 с.: ил.
6 Ахметов С.А., Рындина Д.В., Рындин В.В. Применение математической системы Mathcad в курсовом проектировании по теории ДВС. / Материалы научной конференции молодых учёных, студентов 84
и школьников «V Сатпаевские чтения». Павлодар: НИЦ ПГУ им. С.Торайгырова, 2005.
7 Рындин В.В. Расчёт турбокомпрессора с использованием математической системы Mathcad. Методические указания к курсовой работе по газовой динамике и агрегатам наддува для специальности 280440 «Двигатели внутреннего сгорания». Павлодар: НИЦ ПГУ им. С.Торайгырова, 2005. - 50 с.: ил.
8 Рындин В.В., Макушев Ю.П. Газовая динамика и агрегаты наддува. Методические указания к лабораторным работам по курсу «Газовая динамика и агрегаты наддува» для специальности 280440 «Двигатели внутреннего сгорания». Павлодар: Кереку, 2007. - 63 с.: ил.
9 Хайбулина Р.А., Рындин В.В. Автоматизированный расчёт гидропривода с использованием системы Mathcad // Наука и техника Казахстана. - 2010. - № 4. - С. 109 - 118..
10 Рындин В.В. Теплотехника. - Павлодар: Издательство «Кереку», 2007. - 460 с.
Павлодарский государственный университет им. С. Торайгырова, г.Павлодар. Материал поступил в редакцию 27.02.2012.
В.В. РЫНДИН
MATHCAD ЖYЙЕСШ К;ОЛДАНЫ1СПЕН ЖЫ1ЛУ 0ТК1ЗГ1ШТ1КТЩ АЙЫРЫШДЫЩ ТЕНДЕУЛЕРДЩ ШЕШ1М1
V.V. RYNDIN
THE DECISION OF DIFFERENTIAL EQUATIONS OF HEAT CONDUCTIVITY WITH USE OF THE SYSTEM OF MATHCAD
Туйшдеме
MathCAD жуйесш цолданыспен ацырлы айырымдъщ жуыктама odiciMeH стационар емес жылу вткiзгiштiктiц дифференциалдыц тецдеулгрдщ шешiмi 6epinsdi.
Resume
The decision of the differential equations of non-stationary heat conductivity by a method of finite-difference approximation with use of the system of Mathcad is given.
85
