CC BY
146
УДК 681.5
Методы дискретизации в задачах математического моделирования систем
Д.А. Сидоров
Проведен сравнительный анализ различных методов дискретизации, используемых в задачах математического моделирования двухмерных задач различной физической природы.
The comparative analysis of various methods of the digitization used in problems mathematical modeling of bidimentional problems of the various physical nature is lead.
Постановка задачи
Математическое моделирование какого-либо объекта можно условно разбить на три этапа: модель-алгоритм-программа (рис.1). На первом этапе выбирается (или строится) модель («эквивалент» объекта), отражающая в математической форме важнейшие его свойства - законы, которым он подчиняется. Следующий этап - разработка численного метода решения системы уравнений, описывающих математическую модель. И последний шаг - создание программ, переводящих модель и алгоритм на доступный компьютерный язык. Важное место в решении задачи моделирования играют вопросы дискретизации. Рассмотрим наиболее распространенные методы дискретизации на примере задач гидродинамики.
Рис. 1. Этапы моделирования
Аппроксимация градиента давления
При составлении дискретного аналога уравнения количества движения в направлении оси х для одномерного случая (рис. 2) единственной особенностью является представление члена ёрШх, проинтегрированного по контрольному объему.
В результате интегрирования в дискретный аналог войдет разность давлений (р„ - ре), которая
Рис. 2. Типичный шаблон узловых точек для одномодной залачи
представляет собой силу давления, приложенную к контрольному объему с единичной площадью поперечного сечения. Чтобы выразить (р„ - ре) через значения давления в узловых точках, можно предположить, что давление между узловыми точками изменяется по линейному закону. Если границы контрольного объема выбраны так, что они лежат посередине между соответствующими узловыми точками, то можно записать следующее выражение:
_ Рш + Рр Рр + Ре _ Рш - Ре
■Ре =
■•(1)
2 2 2 Таким образом, дискретный аналог уравнения количества движения будет содержать разность давлений между двумя несоседними точками.
Это означает, что давление берется с сетки более грубой, чем основная, что должно привести к снижению точности решения.
Однако имеются более серьезные недостатки метода. Они лучше видны на рис. 3, на котором поле давления представлено через его значения в узловых точках. Такое зигзагообразное поле нельзя считать физичным. Можно заметить, что для
Р=100 500 100 500 100
Рис. 3. Зигзагообразное поле давления
каждой узловой точки Р соответствующий перепад давления (рш- ре) = 0, так как значения давления в соседних с Р узловых точках равны между собой.
Таким образом, ошибочным следствием данной аппроксимации будет то, что такое волнистое поле давления будет восприниматься в уравнении количества движения как однородное.
Эта трудность еще более усугубляется в двухмерном случае. Так же, как на количество движения в направлении оси х влияет перепад давления (рш - рЕ), на количество движения в направлении оси у влияет перепад давления Р - РЫ), при этом значение давления в точке Р не играет никакой роли. Имея это в виду, можно сделать вывод о том, что показанное на рис. 4 поле давления, образованное из расположенных в шахматном порядке четырех произвольных значений давления, не даст силу давления в направлениях осей х или у.
поле дает нулевой градиент давления, уравнение количества движения «не почувствует» его прибавления. Естественно, что численный метод, который допускает такие абсурдные решения, нежелателен.
Аппроксимация уравнения неразрывности
Аналогичная трудность возникает при построении дискретного аналога уравнения неразрывности. Для стационарного двухмерного течения жидкости с постоянной плотностью уравнение неразрывности имеет вид йи йу
— + —+_ 0. (2)
йх йу
Проинтегрировав это уравнение по изображенному контрольному объему (рис. 5), получим выражение
Пе - + V - Уп _ °. (3)
У
0
Рис. 4. Шахматное поле давления
Таким образом, при рассмотренном способе дискретизации уравнений количества движения сильно неоднородное поле давления будет восприниматься как однородное. Если бы в процессе итерационного решения возникли бы такие поля давления, они бы не сохранились в процессе, так как уравнения количества движения “забудут” об этих полях. Конкретные значения давления на рис. 3 и 4 не имеют какого-либо особого значения, они просто обозначают картину распределения, которую можно было бы составить из любых других значений. Легко представить, что в трехмерном случае может иметь место еще более сложная картина, которая будет восприниматься уравнениями количества движения как однородное поле давления. Если в процессе решения будет получено некоторое гладкое поле давления, можно составить любое количество дополнительных решений путем прибавления к этому решению так называемого шахматного поля давления. Так как это
Рис. 5. Контрольный объем (заштрихованная область) для двухмерного случая
Используя кусочно-линейные профили для и, V и располагая грани контрольного объема посередине между узловыми точками, получаем
(ПР - ПЕ ) - (ПШ - ПР ) + (УР - Уу ) - - УР ) _ 0
2 2 2 2
(4)
ИЛИ
(иЕ - иш ) + К - ) _ °. (5)
Итак, аппроксимация уравнения неразрывности привела к приравниванию скоростей в чередующихся узловых точках, а не в соседних. В результате дискретному аналогу (5) уравнения неразрывности может удовлетворять нефизичное поле скорости (см. рис. 3).
Эти трудности надо исключить до формулировки численного метода решения задачи в переменных, включающих в себя составляющие скорости и давления. В [1,2] описаны методы, в которых этим трудностям не уделяется никакого специального внимания. Возможные нефизичные решения исключаются с помощью некоторой их
специальной обработки на границах, переопределения граничных условий, нижней релаксации относительно гладкого начального приближения и т.п. Но большинство из этих методов воспримут нефизичные поля давления и скорости как удовлетворительные решения; при отсутствии же специальных приемов всегда имеется опасность подобных решений.
Сложности численного анализа связаны с первыми производными давления (в уравнениях количества движения) и скорости (в уравнениях неразрывности).
Шахматная сетка
Описанные выше трудности можно преодолеть, если понять, что не обязательно рассчитывать все переменные в одних и тех же узловых точках. Можно по желанию использовать для каждой зависимой переменной свою сетку. При расчете составляющих скорости значительную выгоду дает определение их на сетке, отличной от такой, которая используется для всех других переменных. В результате такого подхода описанные ранее трудности полностью исчезнут. При расположении сетки в шахматном порядке составляющие скорости рассчитываются для точек, лежащих на гранях, контрольных объемов. Таким образом, составляющая скорости и вдоль оси х рассчитывается на гранях, перпендикулярных направлению оси х. Точки, в которых определяется и, показаны на рис. 6 горизонтальными стрелками, а узловые точки (основные) изображены кружками. Штрихом обозначены грани контрольного объема.
Рис. 6. Расположение и в шахматном порядке: горизонтальные стрелки - места определения и; кружки - места определения других переменных
Точки, в которых определяется и, смещены по отношению к узловым точкам основной сетки, только в направлении оси х. Другими словами, эти точки лежат на отрезках, соединяющих две сосед-
ние (вдоль оси х) основные узловые точки. Находится ли точка, где определяются и, точно посередине между основными узловыми точками, зависит от того, как выделены контрольные объемы. Узловая точка для и должна лежать на грани контрольного объема независимо от того, находится ли последняя посередине между узловыми точками или нет.
Легко выбрать способ размещения узловых точек для составляющих V и и. На рис. 7 показана двухмерная сетка, где узловые точки для и и V помещены на соответствующих гранях контрольных объемов.
Рис. 7. Расположение и и V в шахматном порядке: горизонтальные стрелки - места определения и; вертикальные стрелки - места определения у; кружки - места определения других переменных
Прямым следствием введения шахматной сетки является то, что массовый расход через грани контрольного объема можно теперь определять без интерполяции соответствующей составляющей скорости, что дает некоторые удобства при записи дискретного аналога обобщенного дифференциального уравнения.
Имеются два важных преимущества шахматной сетки. Во-первых, для типичного контрольного объема (заштрихованный участок на рис. 7) дискретный аналог уравнения неразрывности содержит разности составляющих скорости в соседних точках, а это приводит к тому, что волнистое поле скорости (см. рис. 4) не будет удовлетворять уравнению неразрывности. При использовании шахматной сетки только физичные поля скорости могут удовлетворять уравнению неразрывности. Во-вторых, разность давлений между двумя соседними узловыми точками определяет составляющую скорости в точке, расположенной между этими узловыми точками. Поэтому такие поля давления, какие показаны на рис. 3. и 4, не будут восприниматься как равномерные и не могут использоваться как возможные решения.
Трудности, описанные выше, являются результатом расчета всех переменных в одних и тех же узловых точках; они полностью исчезают при использовании шахматных сеток.
Разложение в ряды Тейлора
Основные конечно-разностные формулы для частных производных могут быть получены при помощи разложения в ряды Тейлора. Используемая прямоугольная сетка показана на рис. 8. Нижние индексы / и у относятся к х и у, а верхний индекс п соответствует временному слою. Шаги сетки в направлениях / и у обозначаются через Ах и Ду соответственно. (Для простоты Ах и Ду считаются постоянными.) Переменная / означает какую-либо функцию.
Рис. 8. Прямоугольная конечно-разностная сетка
Формы односторонних разностных представлений для первой производной ё/1ёх можно вывести следующим образом. Предполагая непрерывность производных, раскладываем /+і1 в ряд Тейлора в окрестности точки (і,1). Тогда
г г д/
/г +1,1 - /г, 1 + V
дх
(хі+1„, - Хг, і )+ Т
1 д2 /
2 дх2
Хг+1,, -
)2 = г д/
- Хг,У ) +..= 1 дХ
Л 1 д2 /
Дх + —
2 дх
Дх2 + ЧВП,
(6)
где сокращение «ЧВП» означает «члены высших порядков».
Решая (6) относительно й//йх, получаем
или
дх
/
дх
/г +1,1 - /г,у 1 д7
и 1 ^ 1 Дх
2 дх2
Дх + ЧВП
г, і
/ - /
- г+1,1 г 1 + О(Дх),
Дх
(7)
где запись 0{Дх) читается как «член порядка Дх» и относится к членам, содержащим множители Дх, Дх2 и т. д.
Обозначим конечно-разностный аналог ё/1ёх через ё/1ёх. Тогда для ё/1ёх при разностной аппроксимации вперед получаем выражение
8х
/г +1,, - Л
_ -/і+1,1
Дх
(8)
с ошибкой аппроксимации порядка Дх, т. е. с первым порядком точности.
Раскладывая /¡+1у в окрестности точки (/', у), получаем для ё/1ёх выражение при разностной аппроксимации назад:
А
8х
— 1/1, 1
Дх
(9)
которое также имеет первый порядок точности.
Центральная (симметричная) разностная аппроксимация ё/1ёх получается как разность разложений
г _ г д/
А+1, у _ А у +—
дх
Д 1 д2 /
Дх + —
Дх2 +
1 д3/
+ — у
6 дх3
Л 3 1 д4/
Дх + у
/ _ / - —
Jг-l,1 Jг,J дх
1 д3/
6 дх
Дх3 +
24 дх4 Дх +
^1 1 д4/
2 дх2
г,
Дх4 + О(Дх5)
(10)
1 д2/
Ах2 -
г, і
24 дх4
2 дх2 Дх4 + О(Дх5). (11)
Вычитая (10) из (11), получаем
/ - / _ 2—
/г+1,1 —г-1,1 дх
Д 1 д3/ Дх + -
3 дх3
Дх3 + ЧВП.
Разрешая относительно ё/1ёх, имеем
дх
2Дх
6 дх3
Дх2 + ЧВП. (12)
Таким образом, центральная разностная аппроксимация дает выражение
§х
/г+1,1 /г-к
2Дх
с ошибкой аппроксимации порядка Дх , т.е. со вторым порядком точности.
Интегральный метод
Для построения конечно-разностных аналогов уравнений в частных производных можно использовать методы, основанные на интегрировании этих уравнений (интегральные методы). Рассмотрим уравнение теплопроводности
Эм
э7
:а
Э 2и Эх2
(13)
Построим разностную схему путем интегрирования уравнения теплопроводности по ( и х в окрестности узла (п, у) разностной сетки. Этот узел будем иногда также обозначать как точку (?°, х°). Шаги разностной сетки обозначим Ах и А?1. Так как выбор области интегрирования произволен, то проинтегрируем уравнение (13) от ?° до (?°+ А?) и от (х° - Ах/2) до (х° + Ах/2). Выбор интервала интегрирования от (х° - Ах/2) до (х° + Ах/2) приведет к абсолютно неустойчивой разностной схеме. К сожалению, на этом этапе построения разностной схемы для решения уравнения в частных производных нельзя сказать, какие интервалы интегрирования целесообразно выбрать для обеспечения устойчивости численного метода. Порядок интегрирования в каждой части уравнения выбирается так, чтобы использовать точные дифференциалы:
^ ^ ° +д? /х° +Дх/^2
Г ' Эи
х° +Дх / 2
Эи, и
—йх: Э?
■ а
х° +Дх/2
Э 2
/
(14)
Взяв точно внутренние интегралы, получим
х° +Дх/2
1 [и(?° + Д?, х) - и( ?°, х)йх _
х° +Дх/2
| / (уУу1
[м(?° + Д ?, х°) - м(?°, х° )]Дх:
_а — (?° + Д?, х° +Дх/2)-Эх
-^Т (?° +Д, х°-Дх/2) Эх
Д?.
(18)
Для выражения результата в алгебраической форме выразим производную йи/йх через значения функции и в узлах разностной сетки. Для этого можно использовать уже известные аппроксимации, например центральные разности. С другой стороны, можно придерживаться чисто интегрального метода, тогда на основе теоремы о среднем значении получим
и (?° + Д ?, х° + Дх) _ и (?° + Д ?, х°) +
х° +Дх -ч
+ 1 Э" ? + Д?, х)»
Эи
» и(?° + Д?, х°) +-(?° + Д?, х° + Дх / 2)Дх. (19)
Эх
Из последнего соотношения следует
—(?° + Д ?, х° + Дх/2)@
Эх
и (?° + Д?, х° + Дх) - и (?° + Д?, х°)
Дх
(2°)
_а I —(?, х° + Дх/2)-----------(?, х° -Дх/2) й.
- Эх Эх
?°
(15)
Для вычисления оставшихся интегралов воспользуемся теоремой о среднем значении. Из этой теоремы следует, что для любой непрерывной функции /у) справедливо выражение
У1 +ДУ
7 )Ду, (16)
где у - некоторое значение у из интервала у! < у < (у1+Ау). В соответствии с этой теоремой любое у из указанного интервала позволяет получить приближенное значение интеграла от непрерывной функции:
у +Ду
■■/ (у)Ду, у1 £ у £ у +Ду . (17)
При вычислении интеграла в правой части (19) по теореме о среднем значении было произвольно выбрано, что х =х° + Ах/2 (средняя точка интервала), поэтому интеграл в правой части вычисляется приближенно. Найдя аналогично конечно-разностные аппроксимации остальных первых производных, получим конечно-разностный аналог уравнения теплопроводности [и(?° +Д?,х°) - и(?°, х°)]Дх _
_ —[м(/„ + Д, х° +Дх)-Дх
- 2м (?° + Д?, х°)+м (?° + Д?, х° -Дх)]Д?. (21)
Обозначая значения функций в узлах индексами п, у, где п - номер шага по времени, а у - номер шага по пространственной координате, перепишем (21)в виде
а
— \и".'‘ - 2и" ' ‘ + ),_у-1
п+1 . п
иу +иу
( п+1
\2 К+1
п +1 п +1
2и + и -1 .
Для дальнейшего упрощения интеграла (16) при помощи теоремы о среднем значении возьмем значение подынтегральной функции при х =х°, а при вычислении интеграла в правой части - при ? _ ?° + Д?. Тогда получим
(22)
Д (Дх)2 ' ,т' ' "''
Последнее выражение совпадает с полностью неявной схемой для уравнения теплопроводности. Неявная разностная схема получена благодаря тому, что интеграл в правой части вычислялся по теореме о среднем при значении подынтегральной
и +Д?
функции в момент времени (?°+А?). Если бы при вычислении этих интегралов были использованы значения подынтегральной функции при ? = ?° то получилась бы явная схема.
Полиномиальная аппроксимация
Другой метод получения конечно-разностных выражений основан на применении аппроксимирующей аналитической функции со свободными параметрами, которая строится по значениям в узлах сетки и затем аналитически дифференцируется. Это обычный метод нахождения производных по экспериментальным данным. В идеале вид аппроксимирующей функции должен определяться приближенным аналитическим решением, однако обычно в качестве аппроксимирующих функций используются полиномы. Продемонстрируем настоящий метод на примере параболической аппроксимации.
Предположим, что значения функции / заданы в точках 7 - 1, (и 7 + 1 и проведем параболическую аппроксимацию функции
/ (х) - а + Ьх + ех2
(23)
е
2Дх2
а разрешая их относительно Ь, находим
/7 +1 - /-1
Ь
2Дх
(25)
(26)
будет
8х
= [Ь + 2ех]х-о _ Ь :
а значение второй производной
82/
8х2
: 2е.
(27)
(28)
зуются для определения а и Ь - значения / и /+1 или / и /_!, для 8 // 8 х получаются формулы с разностями вперед или назад соответственно. Очевидно, что при линейной аппроксимации / нельзя получить выражение для 8 /8 х2. Однако, если использовать полином первой степени для построения разностных аналогов первых производ-
ных
8/
8х
і+112
8х
которые соответственно
причем для удобства за начало координат (х = °) примем точку 7. Тогда уравнение (23), записанное в точках 7 - 1,7 и 7 + 1 соответственно, даст
/-1 _ а - ЬДх + сДх2, / _ а,
/+1 _ а + ЬДх + сДх2. (24)
Складывая первое и последнее из этих равенств, получаем
/+1 + /-1 - 2/
представляются разностями вперед и назад, то для 8 / 8 х получится выражение, совпадающее с выражением с центральными разностями.
Разностные формулы для производных более высокого порядка выводятся с использованием полиномов высших порядков. Выражения, полученные при помощи полиномов порядков выше второго, уже не идентичны выражениям, полученным разложениями в ряды Тейлора, и в каждом случае ошибка аппроксимации должна проверяться при помощи разложения в ряд Тейлора. В вычислительной гидродинамике метод полиномиальной аппроксимации, как правило, применяется только для вычислений значений производных вблизи границ.
Отметим недостатки полиномиальных аппроксимаций высшего порядка. С увеличением порядка аппроксимации они становятся чувствительными к «шумам», т. е. к более или менее случайно распределенным малым ошибкам в данных. Так, полином шестой степени, график которого проходит через семь точек, точно расположенных на одной прямой, приводит к аппроксимации в виде прямой, изображенной на рис. 9,а. Однако при добавлении к аппроксимируемым значениям шумовых возмущений коэффициенты полинома будут уже определяться этими искаженными данны-
В точке 7 значение первой производной (23)
Формулы (27) и (28) с учетом (25) и (26) в точности совпадают с формулами (13) второго порядка с центральными разностями, полученными разложением в ряд Тейлора. Если предположить, что/- полином первой степени, т. е. /(х) _ а + Ьх, то в зависимости от того, какие значения исполь-
Рис. 9. Полиномиальная аппроксимация шестого порядка: а - алгебраические идеальные данные; б - данные с добавлением шумовых возмущений
г-112
ми, и тогда аналитическое вычисление производных в точке 7 может привести к абсурдным результатам, что можно усмотреть на рис. 9,6.
Квадратичная аппроксимация не может отразить наличие точки перегиба в рассматриваемых данных, т. е. точки, где 8 //8 х2=°.
Метод контрольного объема
Этот метод принципиально отличается от уже рассмотренных методов построения разностных схем для уравнений в частных производных. При использовании методов, основанных на разложении функций в ряд Тейлора, и интегральных методов предполагалось, что уравнения в частных производных корректны и в соответствующей форме описывают законы сохранения, поэтому при построении разностной схемы применялись математические средства. При таком подходе физические законы сохранения применяются лишь при выводе уравнения в частных производных и никак не используются при построении разностной схемы. В этом смысле разложение в ряд Тейлора и интегральные методы - формальные методы построения разностных схем для уравнений в частных производных.
При использовании метода контрольного объема разностная схема строится на основе физических законов сохранения, следствием которых является рассматриваемое уравнение в частных производных. Сначала этот закон сохранения формулируется словесно для некоторого контрольного объема, окружающего узел разностной сетки, а потом записывается математически с учетом дискретной сетки. Описанная процедура во многом похожа на ту, с помощью которой уравнения в частных производных выводятся из физических законов сохранения, не проводится лишь переход к пределу при стягивании контрольного объема в точку. Если уравнение в частных производных записано в дивергентной форме, то закон сохранения можно получить, интегрируя это уравнение по контрольному объему и используя формулу Гаусса-Остроградского.
Отличительной особенностью метода контрольного объема является то, что он обеспечивает «баланс» физической величины в окрестности узла разностной сетки. Метод контрольного объе-
ма учитывает дискретный характер решения поставленной задачи, поэтому он обеспечивает выполнение законов сохранения в конечной области, а не только в точке при стремлении шага сетки к нулю. Конечно-разностные схемы, построенные методом контрольного объема, почти всегда консервативны.
Таким образом, проведенный сравнительный анализ различных методов дискретизации, используемых в задачах математического моделирования двухмерных задач различной физической природы, показал следующее: использование типичного
или зигзагообразного шаблона узловых точек для одномерных и двухмерных задач ведет к снижению точности дискретизации; возможные нефизичные решения исключаются с помощью некоторой их специальной обработки на границах, переопределения граничных условий, нижней релаксации относительно гладкого начального приближения; сложности подобного численного анализа связаны с первыми производными давления (в уравнениях количества движения) и скорости (в уравнения неразрывности); более предпочтительна дискретизация в виде «шахматной сетки», которая имеет два важных преимущества: во-первых, дискретный аналог
уравнения неразрывности содержит разности составляющих скорости в соседних точках и, во-вторых, разность значений между двумя соседними узловыми точками определяет составляющую в точке, расположенной между этими узловыми точками; наиболее предпочтительны методы построения конечно-разностных аналогов уравнений в частных производных путем использований интегральных методов, а также метод получения конечно-разностных выражений, основанный на применении аппроксимирующей аналитической функции со свободными параметрами, которая строится по значениям в узлах сетки и затем аналитически дифференцируется.
ЛИТЕРАТУРА
1. Полежаев В.П., Бунэ A.B., Верезуб H.A. и др. Математическое моделирование тепломассообмена на основеуравнений Навье-Стокса. - М.: Наука, 1987.
2. Самарский А. А., Михайлов А. П. Математическое моделирование: Идеи. Методы. Примеры. - М.: Наука. Физматлит, 1997.
Поступила 09. 03. 2008 г.
