CC BY
77
2017 Математика и механика № 49
УДК 519.6
Б01 10.17223/19988621/49/6
А.А. Семёнова, А.В. Старченко
РАЗНОСТНАЯ СХЕМА ДЛЯ НЕСТАЦИОНАРНОГО УРАВНЕНИЯ
ПЕРЕНОСА, ПОСТРОЕННАЯ С ИСПОЛЬЗОВАНИЕМ ЛОКАЛЬНЫХ ВЕСОВЫХ ИНТЕРПОЛЯЦИОННЫХ КУБИЧЕСКИХ СПЛАЙНОВ1
Представлена новая аппроксимация конвективных членов в нестационарном конвективно-диффузионном уравнении переноса примеси, полученная с использованием локальных весовых кубических сплайнов. На основе сравнительного анализа для двух одномерных тестовых случаев нестационарного переноса примеси с известным аналитическим решением показано её преимущество перед широко используемыми при решении подобных задач мо-нотонизированными схемами второго или третьего порядка аппроксимации.
Ключевые слова: нестационарное конвективно-диффузионное уравнение, локальные весовые кубические сплайны, монотонизированная аппроксимация конвективных членов высокого порядка.
В настоящее время при изучении разнообразных физических процессов широко используются методы математического моделирования, основная идея которых заключается в представлении свойств изучаемого объекта или явления с помощью математических формул, как правило, выражающих тот или иной физический закон. При исследовании процессов в рамках механики сплошных сред применяются законы сохранения массы, импульса, энергии, которые, в общем случае, имеют вид нестационарных интегральных или дифференциальных уравнений, включающих члены, описывающие перенос массы, импульса или энергии как за счет непосредственно самого движения сплошной среды (конвекция), так и за счет молекулярного или турбулентного обмена (диффузия) [1].
Как правило, получение аналитических решений таких уравнений на современном уровне развития теоретической математики представляет определенную проблему, поэтому активно разрабатываются приближенные численные методы решения, обеспечивающие выполнение условий аппроксимации, устойчивости и сходимости [2]. Особое внимание при решении нестационарного уравнения «конвекции - диффузии» уделяется выбору аппроксимационной схемы для конвективного члена уравнения. С одной стороны, требуется, чтобы порядок аппроксимации схемы был выше первого для обеспечения качественного воспроизведения резких изменений переносимой субстанции во времени и пространстве, с другой -используемая аппроксимация (или разностная схема) должна быть монотонной, т.е. обеспечивающей получение монотонного распределения искомой величины из имеющегося монотонного распределения для предыдущего момента времени.
Современные численные методы решения уравнения переноса берут свое начало в том числе и с известной работы Годунова [3], в которой он предложил численный метод для решения локальной проблемы Римана на границе раздела двух объемов сплошной среды. Им представлена очень важная теорема о монотонно-
1 Работа выполнена при финансовой поддержке РФФИ, грант № 16-43-700178.
сти: не существует монотонных линейных разностных схем с порядком аппроксимации по пространству выше первого. Монотонность является очень желательным свойством разностной схемы, ее наличие предотвращает образование вычислительных осцилляций в окрестности сильных градиентов решения. Лакс и Венд-рофф [4] предложили для решения уравнения переноса схему второго порядка аппроксимации, обеспечивающую хорошее предсказание резко меняющегося поведения приближенного решения, но сильно осциллирующую. Для уменьшения вычислительных осцилляций необходимо в эту схему вводить дополнительные дис-сипативные члены.
В настоящее время для того, чтобы обойти жесткое ограничение теоремы Годунова, используется построение нелинейных разностных схем для решения уравнения переноса. В качестве одного из способов построения нелинейных схем является применение кусочно-полиномиальной реконструкции приближенного решения по значениям сеточной функции с введением некоторых ограничений (условий), обеспечивающих уменьшение вычислительных осцилляций при переходе на следующий момент времени. Одной из первых работ этого направления была работа Колгана [5], который предложил применять принцип минимальных значений производных, чтобы получить неосциллирующую схему второго порядка типа схемы Годунова. Другой известный подход предложил ван Лир в своей серии статей [6, 7, 9], в которой развивается теория разностных схем высокого порядка аппроксимации на основе полиномиальной реконструкции и контроля за значениями производных этих полиномов с использованием TVD-свойства (свойства невозрастания полной вариации решения при переходе на новый слой по времени) [10, 11].
На основе полиномиальной TVD-реконструкции ван Лиром, в частности, были предложены схемы MLU (Monotonie Linear Upstream) [7, 8], MUSCL (Monotonie Upstream eentered Seheme for Conservation Laws) [9], которые используют кусочно-линейную интерполяцию и обеспечивают точность второго порядка в случае гладкого монотонного решения. Подход TVD ввел новые усовершенствования в построении схем высокого порядка точности для уравнения переноса и помог избежать ограничений для линейных разностных схем, предоставляемых теоремой Годунова. Но порядок этих схем ограничивается вторым на участке гладкого монотонного изменения решения (см. например, схему Хартена c ограничителем minmod [12] и схему с ограничителем superbee [13]).
Для повышения порядка аппроксимации разностных схем для решения уравнения переноса были введены так называемые существенно неосциллирующие схемы ENO (Essentially Non-Oscillatory sehemes). Схемы ENO для получения ос-редненного по ячейке значения решения содержат алгоритм построения нескольких полиномов высокого порядка на различных сеточных шаблонах, выбирающих наиболее гладкий из полиномов. Такой подход был введен в 1986 г. Harten, Osher, Engquist и Chakravarthy [14].
При численной реализации моделей атмосферного пограничного слоя и переноса примеси важную роль играет выбор качественной разностной схемы для аппроксимации конвективного члена уравнения. Целью данного исследования является разработка построенной на локальных весовых сплайнах разностной схемы для аппроксимации конвективных членов уравнения переноса. Полученная аппроксимация используется в явной разностной схеме, которая сравнивается с другими известными аппроксимациями конвективных слагаемых (MLU, MUSCL,
£N0 и др.) при численном решении двух одномерных нестационарных тестовых задач о распространении примеси.
1. Интерполяционный весовой кубический сплайн
" п—1 "
Введем сетку юИ = |хг; хг+1 = хг + Иг,I = 0, п — 1; х0 = а; ^ Иг = Ь — а| на отрезке
[а,Ь]. На каждом элементарном отрезке [хг,хг+1 ] интерполяционный весовой кубический сплайн есть полином третьей степени (х):
(х) = а0 + аг1 (Х — X ) + а2 (Х — X )2 + аг3 (Х — Х1 )3,
х 6 [ хг, хг+l],
коэффициенты которого надо определить из следующих условий [15]:
8г (хг ) = /г, г = 0 п;
%(хг ) = (хг )Д = 1,п — 1;
М/—1^/— 1 (хг ) = (хг )Д = ит—1.
Здесь м>г > 0, г = 0, п — 1, - весовые коэффициенты.
Построение весового кубического сплайна будет проводиться через наклоны (тг = (хг), г = 0, п ), для получения единственного решения поставленной задачи будем использовать дополнительные условия на границах:
8 ' (х0) = / ,8' (хп) = /'.
Из приведенных выше условий находим формулу для построения сплайнов через наклоны:
8 (х) = (х — хг )3 — (х — хг )3 + 3/±Ц/ (х — ^ )2 —
И? " и} " И2
2т, + т,
г-— (х — xi )2 + mi (х — xi) + /; xi < х < хг+1; г = 0, п — 1.
И
А из условий м{—18г"_ 1(хг) = М'З/Хх, ), г = 1,п — 1, и дополнительных граничных условий получим систему уравнений для определения значений наклонов {тг}:
т0 = /0 ,
тг —1 + 21 Т^ + ~Т Iтг +^тг+1 = 3
Иг —1 I Иг—1 Иг ) Иг
тп = Л .
^/г /г—1 м 1 + /г+1 /г М Л
И2 1 И2
V "г—1 "г
г = 1, п — 1,
Эта система имеет строгое диагональное преобладание, что обеспечивает существование и единственность весового кубического сплайна. Решение полученной системы может быть найдено методом обычной трехточечной прогонки.
Весовые коэффициенты {м>, } можно задать разными способами. Если интерполируемые данные монотонно возрастают, веса можно найти, воспользовавшись следующей теоремой Б.И. Квасова [15].
Теорема. Пусть весовой кубический сплайн (х) е С1 [а,Ь] с краевыми условиями £'(х0) = /0' и 5'(хп) = /П интерполирует монотонные данные {/},г = 0,п . Если выполняются неравенства
0 < /'< 3 / [хо, X ], 0 < /< 3/[ хп_х, хп ], Щг-1 Иг > / [ X , X-+1 ] _ 2 Щ > / [ X-X- ] _ 2 Щ кг-1 /■ [хг-1, X] ' Щ-1 Лг У [Х, X+1] ' то S'(x) > 0 для всех X е [а,Ь], т.е. сплайн 5 монотонен на [а,Ь].
В [16] для монотонно убывающих данных был сформулирован аналог этой теоремы.
Для нахождения весовых коэффициентов можно следовать алгоритму [15]. Пусть весовой коэффициент щг-1 известен. Полагаем щ = щ-1. Проверяем неравенства из теоремы. Если какое-то из них нарушается, то выражаем из соответствующего неравенства щ, заменяя знак неравенства на знак равенства. Алгоритм
начинаем с щ0 = 1 и легко находим все коэффициенты {щ}, обеспечивающие монотонность весового кубического сплайна для произвольных монотонных данных. Если на некотором шаге < е, ^> —^, то полагаем щ = е, ^щ = —^ , где е
- достаточно малое положительное число, позволяющее предотвратить зануление (переполнение).
Рис. 1. Интерполяция таблично заданной функции алгебраическими многочленами Fig. 1. Algebraic polynomials interpolation
Продемонстрируем преимущество весовых кубических сплайн-функций при интерполировании таблично заданных функций. Рассмотрим случай, когда заданы таблично четыре значения функции. Требуется восстановить значение интер-полянты для промежуточных значений X [16]. Для решения этой задачи были
применены два способа: кубический интерполяционный сплайн, построенный через наклоны, и такой же сплайн, но с весами, подобранными по алгоритму, описанному выше. При построении сплайнов на концах выбранного отрезка использовались дополнительные условия равенства нулю производной интерполируемой функции. Заметим, что табличные значения функции выбирались таким образом, чтобы при интерполировании могла возникнуть проблема монотонности построенных интерполяционных многочленов. Данный набор табличных значений интерполируемой функции соответствует профилю аналитического решения задачи «бездиффузионного» переноса примеси, рассмотренной ниже (Тест 1).
Из рис. 1 видно, что интерполяционный кубический сплайн осциллирует на рассматриваемом промежутке между табличными значениями, предсказывая отрицательные значения на отрезке [1,0; 2,0] и больше единицы на отрезке [3,0; 4,0]. Лучший результат демонстрирует интерполяционный весовой кубический сплайн. В данном случае осцилляций визуально не наблюдается и сглаживание в области резкого изменения значений интерполируемой функции меньше, чем в предыдущем случае.
2. Численный метод решения уравнения переноса
Используем полученный выше результат для весового кубического сплайна при построении разностной схемы со сплайн-аппроксимацией конвективного члена уравнения переноса примеси.
Рассмотрим одномерное нестационарное уравнение конвекции - диффузии с постоянными коэффициентами для концентрации примеси:
дФ дФ д 2Ф
-+ и-= В—- + 5Ф , и > 0, В >= 0; 0 < X < X, 0 < Г < Т . (1)
дГ дx дуг
Начальное условие: ( = 0, Ф(0, X) = Ф00 (X).
дФ
Граничные условия: X = 0, Ф(/, 0) = Ф0(/); X = X,-= 0.
дx
Существование и единственность решения рассматриваемой задачи доказаны в [17].
Для отрезка [0,Х] построим равномерную сетку с непересекающимися конечными объемами: х^^,xг-+1/2,г = 0,М - положение граней конечного объема (рис. 2); xi, г = 1,М, - положение узлов в конечном объеме. Введем сеточную
функцию для концентрации примеси Ф^ и Ф(^, xi), tk = xk; k = 0, N; i = 0,M .
W | P | E
-О—I—о—I—о-
i-1 i i+1 i-1/2 i+1/2
Рис. 2. Конечный объем для одномерного случая Fig. 2. Finite volume for the one-dimensional case
Применяя метод конечного объема и явную разностную аппроксимацию по времени первого порядка, получим для уравнения (1) разностную схему следующего вида:
Ф -Фк + и ^-Ф-1/2 = д Ф+1 - 2Ф + + (^ )к ; г ^ к = 0,1,2,... . (2) т И н2
В итоге каждый член этого уравнения выражается через дискретные значения Ф в некоторых соседних узловых точках.
Начальное условие аппроксимируется по формуле: Ф0 =Ф0о(X), ' = 0,M, а граничные: Ф0 = Ф0 (tk, 0), ФМ = Ф^М_i, k = 0, N.
Шаг по времени выбирается из условия устойчивости: т < minI ,.
^ 2\и\ 2D у
Неизвестными в (2) являются значения сеточной функции {фк+1} на k+1-м слое по времени. Однако для их вычисления по явной формуле необходимо первоначально рассчитать значения {фк+1/2|, используя выбранный способ их
представления через значения {фк}. В данной работе для этого применяется
сплайн-интерполяция с использованием весового кубического сплайна, т.е. Фк+1/2 = S (X+1/2), где S, (X) - локальный весовой интерполяционный кубический сплайн, построенный по значениям сеточной функции Фк-1, Фк, Фк+1,
Фк+2,:
m ,+ m q Ф^ -Фк Q Ф^ -Фk
SS, (X) = (X - хг )3 - 2 Ф'+1 3 Ф1 (X - хг )3 + 3 Ф'+1 2 Ф' (X - х- )2 -
h h h
2m, + m,+, ч 2 . „л
--^T^(X - X' )2 + m, (X - X') + Фкк; x, < x < x^.
h
Наклоны mt определяются из решения системы: m,-1 =(3Ф k-1 - 4Фk +Ф k+1) /2h; wj-1mj-1 + 2(-1 + wj )mj + wjmj+1 = '=3 [Wj-1 ((- Фk-1)+wj ((+1 - фk)]/h, j=1,1+1
m,+2 =(k+2 - 4Ф1+1 +Фk )/2h.
Весовые коэффициенты локальных сплайнов на 0 < x < X вычисляются по алгоритму [15].
Приведенный выше алгоритм построения сплайна позволяет монотонно приближать исходные монотонные данные, однако если же данные не были монотонными - возникает проблема, для решения которой к интерполированному значению дополнительно применяется ограничитель, аналогичный ограничителю схемы MUSCL [18].
В конечном итоге сплайн-аппроксимация реализовывалась в разностной схеме в следующем виде (е - достаточно малое положительное число):
, Í Фk + 0,5max [0, min(2©, Ф,2)( (+1 - Фk), и > 0, Ф = ^ ^ ^ (3)
'+1/2 ^ - 0,5max [0, min(2©, Ф, 2)] (+1 - Фk), и < 0,
где
Фк -фк 1 S, (X+1/2) -ф
0 __z_z-1_ ^ __i+1/2 ' .
Фк+1 - Фк + в • sign^h - Фк) ' 0.5 • (<, -Фк) + 8- sign^h - Фк) Ф/+2 -Ф?+1 G _ - Si (X+1/2)
0 _ _i + 2 i + 1_ vp _.
<1 -Фк +Е-sign^h -Фк) ' 0.5 • (<, _Фк) + e-sign^k+1 _Фк) ■
3. Численные эксперименты
Для проведения численного эксперимента рассматривались две задачи переноса примеси с известным аналитическим решением [17, 19]. Кроме построенной в предыдущем разделе разностной схемы были реализованы другие подходы аппроксимации конвективного члена уравнения переноса: противопотоковая схема Upwind, схемы MLU [7, 8], MUSCL [9, 18], ENO-схема третьего порядка точности с ограничителем Годунова [20], схема с ограничителем superbee [13] и схема Хартена [12].
Противопотоковая схема аппроксимирует сеточную функцию с самым низким порядком среди всех рассматриваемых в данной работе схем и действует по правилу
ф Ы, u > а
i+1/2 — I ,, k
1Фы« < 0.
MLU-схема записывается следующим образом:
Фк _|фк + 0,5h * P+1/2, u > 0, i+1/2 R, - 0,5h. ppi+1/2, u < 0,
где p+1/2 _ minmod[(Ф^+j _фf_1)/2h;2minmod((ф£ц _Ф*)/h;^f -Ф^)/h)] .
MUSCL-схема применяется в виде
Г Фп + 0,5max [0,min(2©,(2 + ©)/3,2)]]+j _Фf), u > 0, [фк+j _ 0,5 max[0, min(2©, (2 + ©)/ 3,2)] (+j _ Фf), u < 0,
Фк
i+1/2
где 0,0 имеют такой же вид, как в формуле (3).
Построение ENO-схемы подробно описано в [20]. Из [21] были взяты формулы для построения схем с использованием ограничителя Хартена [12]:
[ Фк + 0,5(1 - Kh )min mod (<, - Фк, Фк - ф^), u > 0, [фh - 0,5(1 + Kh )minmod( - Фк, Фк - ф^), u < 0 и ограничителя superbee [13]:
- 0,5(1 - Kh )max mod [min mod (2 (ф|+ - ф*), ф^ - ф^),
min mod ( - ФJ ,2 (фf - Фf- ))] , u > 0, Фк+1 - 0,5(1 + Kh )max mod [min mod (2 (фf+ 2 - Фf+1), Ф*+, - Фкг ),
minmod((+ 2 -<„2(of+j - фк))], u < 0. Здесь Kh _ ut / h .
ф k
^i+1/2
фк +(
^i+1/2
Тест 1. Решение уравнения переноса в случае «отсутствия» диффузии (О = 0) и 8ф = 0 .
Начальное условие для задачи 1 (X = 2):
,х е [0.25, 0.75],
Фоо(Х) = {0, х г [0.25,0.75]. Граничное: Ф0 (/) = 0 . Аналитическое решение задачи 1 имеет вид [11]
Ф ехас^ Х) = Ф 00( Х - и). Расчеты производились при следующих условиях:
и = 1, М = 100, М = 200, ти /И = 0.2, О = 0.
1-,
0.8-
g 06 Н
a a ft н
Ц 0.4
0.2-
0
0.4
0.8
1.2
"Г
1.6
Рис. 3. Рассчитанные (кружки) по сплайновой схеме значения концентрации примеси при t = Т = 1с Fig. 3. Calculated values (circles) of the impurity concentration calculated at t = T = 1с by the spline scheme
Из рис. 3 видно, что численное решение достаточно хорошо соответствует точному, однако в областях резкого изменения решения наблюдается некоторое
«размытие» профиля {Ф^}.
Данная задача решалась и с использованием перечисленных выше методов аппроксимации конвективного члена уравнения переноса. Для сравнения степени согласования расчетов с точным решением рассматривались нормы погрешности:
Я = max bf -
i=0.....M'
Ф(Т. X )|, Я2
z
-
Ф(Т, X )|2
M
значения которых для рассматриваемой задачи представлены в табл. 1.
0
2
X. м
Таблица 1
Сравнение норм погрешности различных методов аппроксимации конвективных членов нестационарного уравнения переноса для первой тестовой задачи
Метод М = 100 М = 200
Н\ Н 2 Н1 Н 2
1 Иртет! 0,48747 0,17209 0,49108 0,14451
2 мьи 0,32540 0,05414 0,32540 0,03828
3 МШСЬ 0,32989 0,05468 0,32898 0,03867
4 ЕЫО3 0,31528 0,06294 0,31794 0,04479
5 Сплайн 0,26939 0,04202 0,26939 0,02971
6 Схема Хартена 0,43604 0,11087 0,44342 0,08867
7 БирегЪее 0,37625 0,06941 0,38222 0,05024
На основании полученных данных можно сделать вывод, что для рассматриваемой тестовой задачи применение весовой сплайн-функции для аппроксимации конвективных членов имеет преимущество перед остальными схемами. Худший результат показывает противопотоковая схема. Эти результаты также согласуются с теоретическими данными о порядке аппроксимации этих схем. Сходимость схемы проверялась удвоением количества узлов. Нормы погрешности методов Н2 при этом уменьшились.
Тест 2. Распространение примеси от мгновенного точечного источника.
Рассмотрим следующую нестационарную задачу, моделирующую распространение примеси от мгновенного источника, расположенного внутри области [17]:
д 2Ф дх2
дФ дФ
-+ и-+ стФ = В
дt дх
Г = 0: Ф(0,х) = 0, Ф(^0) = 0,
дФ
— ^, X) = 0.
дх
Здесь 5(2) - дельта-функция.
Точное решение данной задачи имеет вид [17]
+ QS(x-х0)8^-0 < х < X, 0 < t < т,
Ф
^, х) =
Q
л/4Пв(7-70)
0, t е [0,t0]
Ехр
-ст(t --
4 - to)
2
t е (to, Т],
Здесь t0 - момент времени, в который произошел выброс примеси, х0 - координата источника, ст - интенсивность осаждения примеси; Q - мощность мгновенного точечного источника. Значения констант модели выбирались следующими: и = 1, В = 0.01, ст = 0, Q = 1, X = 100, Т = 100, ^ = 0.5Т, х0 = 15. Расчеты проводились при М = 100,М = 200, ит / к = 0.2. Результаты вычислений представлены на рис. 4.
Также было проведено сравнение данного решения с решениями полученными методами, описанными выше. Из рис. 4 и табл. 2 видно, что на более подробной сетке наилучший результат дают схемы МЬи и Ми8СЬ. Незначительно им уступает сплайновая схема, но сплайновая схема в данном тестовом примере точнее показывает максимальное значение концентрации, что, безусловно, является преимуществом для воспроизведения пиковых значений концентраций загрязнителей атмосферного воздуха.
0.4-
0.3-
я
a
Э
И
0J
a и
¡2
0.2-
0.1-
50
60
70 x, м
• • • Сплайновая схема
Точное решение
X X X Upwind
о о о MLU
▲ ▲ ▲ MUSCL
SuperBee
★ ★ ★ Схема Хартена
80
90
Рис. 4. Рассчитанные значения концентрации примеси при t = T для M = 200 Fig. 4. Calculated values of the impurity concentration at t = T for M = 200
0
Таблица 2
Сравнение норм погрешности различных методов аппроксимации конвективных членов в задаче о мгновенном точечном источнике
Метод M = 100 M = 200
Hi H 2 H H 2
1 Upwind 0.30272 0.04501 0.30349 0.04164
2 MLU 0.16869 0.02616 0.05377 0.00795
3 MUSCL 0.17177 0.02616 0.06590 0.00841
4 ENO3 0,17091 0,02623 0,14772 0,02143
5 Сплайн 0.15203 0.02370 0.07608 0.00956
6 Схема Хартена 0.25983 0.03921 0.23127 0.03113
7 Superbee 0.20095 0.03104 0.12319 0.01534
Таким образом, можно говорить о том, что использование формосохраняющей сплайновой интерполяции является перспективным направлением при конструировании монотонных разностных схем повышенного порядка аппроксимации.
Заключение
Наиболее популярными методами решения задач конвективно-диффузионного переноса остаются разностные схемы, поскольку в общем случае получить аналитическое решение дифференциальной задачи затруднительно. Построение схем высокого порядка точности является актуальной задачей, поскольку такие схемы дают решение с удовлетворительной погрешностью даже на грубых сетках.
В данной работе для решения нестационарного уравнения «конвекции-диффузии» использовались явные разностные схемы, метод конечного объема и различные способы для аппроксимации конвективного слагаемого уравнения. В частности, рассмотрен вопрос о построении разностных схем высокого порядка точности на основе локальных интерполяционных весовых кубических сплайнов для уравнения «конвекции - диффузии».
Полученная разностная сплайновая схема была применена для решения двух тестовых задач: «бездиффузионного» распространения примеси и распространения примеси от мгновенного точечного источника. Для сравнения использовались разностные схемы с различным способом аппроксимации конвективного члена уравнения: противопотоковая схема, схема Хартена, схема с ограничителем superbee, MLU, MUSCL и схема ENO, которая позволяет получать аппроксимацию по пространству до 6-го порядка. Результаты расчетов, выполненные для различной плотности узлов сетки, показали сеточную сходимость приближенного решения. Для первой тестовой задачи сплайновая схема продемонстрировала преимущество перед остальными схемами по близости приближенного решения к точному. Для второй тестовой задачи, характеризующейся более гладкими пространственными профилями решения, на грубой сетке лучшее качество согласования с точным решением дала сплайновая схема. На более подробной сетке наилучший результат дают схемы MLU и MUSCL. Незначительно им уступает сплайновая схема, но сплайновая схема в данном тестовом примере точнее показывает пиковое значение концентрации, что, безусловно, является преимуществом для воспроизведения пиковых значений концентраций загрязнителей атмосферного воздуха.
ЛИТЕРАТУРА
1. Самарский А.А., Михайлов А.П. Математическое моделирование. М.:Физматлит, 2001. 320 с.
2. Самарский А.А. Теория разностных схем. М.: Наука, 1977. 656 с.
3. Годунов С.К. Разностный метод численного расчета разрывных решений уравнений гидродинамики // Математический сборник. 1959. Т. 47(89). № 3. С. 271-306.
4. Lax P.D., WendroffB. Systems of conservation laws // Communications in Pure and Applied Mathematics. 1960. V. 13. P. 217-237.
5. Колган В.П. Применение принципа минимальных значений производной к построению конечноразностных схем для расчета разрывных решений газовой динамики // Уч. зап. ЦАГИ. 1972. Т. 3. № 6. С. 68-77.
6. van Leer B. Towards the ultimate conservative difference scheme I. The quest for monotonicity // Lecture Notes in Physics. 1973. V. 18. P. 163-168.
7. van Leer B. Towards the ultimate conservative difference scheme II. Monotonicity and conservation combined in a second-order scheme // J. Comp. Phys. 1974. V. 14. P. 361-370.
8. Noll B. Evaluation of a bounded high-resolution scheme for combustor flow computations // AIAA Journal. 1992. V. 30. No. 1. P. 64-68.
9. van Leer B. Towards the ultimate conservative difference scheme V. A Second order sequel to Godunov's method // J. Comput. Phys. 1979. V. 32. P. 101-136.
10. Toro E.F. Riemann solvers and numerical methods for fluid dynamics. 2nd edition. Berlin, Heidelberg: Springer, 1999. 645 p.
11. Лебедев А.С., Черный С.Г. Практикум по численному решению уравнений в частных производных. Новосибирск: НГУ, 2000. 136 с.
12. Harten A. High resolution schemes for hyperbolic conservation laws // J. Comp. Phys. 1983. V. 49. P. 357-393.
13. Roe P.L. Characteristic-based schemes for the Euler equations // Ann. Rev. Fluid Mech. 1986. V. 18. P. 337-365.
14. Harten A., Engquist B., Osher S., and Chakravarthy S.R. Some results on uniformly highorder accurate essentially non-oscillatory schemes // J. Appl. Num. Math. 1986. V. 2. P. 347377.
15. Квасов Б.И. Методы изогеометрической аппроксимации сплайнами. М.: Физматлит, 2006. 360 с.
16. Карпова А.А. Приближение таблично заданных функций с помощью аппроксимацион-ных и интерполяционных весовых сплайнов // Научная конференция студентов механико-математического факультета ТГУ, 24-30 апреля 2014 г.: сб. материалов. Томск, 2014. С. 38-39.
17. Марчук Г.И. Математическое моделирование в проблеме окружающей среды. М.: Наука, 1982. 319 с.: ил.
18. Cada M., Torrilhon M. Compact third-order limiter functions for finite volume methods // J. Computational Physics. 2009. V. 228. No. 11. P. 4118-4145.
19. Starchenko A.V., Danilkin E.A., Semenova A.A., and Bart A.A. Parallel algorithms for a 3D photochemical model of pollutant transport in the atmosphere // Communications in Computer and Information Science. 2016. V. 687. P. 158-171.
20. Chi-Wang Shu. Essentially non-oscillatory schemes for hyperbolic conservation laws // Preprint of Division of Applied Mathematics. Brown University. 1996. 92 p.
21. Ривин Г.С., Воронина П.В. Перенос аэрозоля в атмосфере: выбор конечно-разностной схемы // Оптика атмосферы и океана. 1997. № 6. С. 623-633.
Статья поступила 05.07.2017 г.
Semyonova A.A., Starchenko A.V. (2017) THE FINITE-DIFFERENCE SCHEME FOR THE UNSTEADY CONVECTION-DIFFUSION EQUATION BASED ON WEIGHTED LOCAL CUBIC SPLINE INTERPOLATION. Tomsk State University Journal of Mathematics and Mechanics. 49. pp. 61-74
DOI 10.17223/19988621/49/6
Keywords: unsteady convection-diffusion equation, weighted local cubic splines, monotonized high order approximation for convective terms.
In this paper, special attention is paid to the choice of an approximating scheme for the convective terms of the unsteady convection-diffusion equation. The purpose of this study is to develop a difference scheme for the convection-diffusion equation with weighted local cubic spline approximation for the convective terms.
The advantage of weighted cubic spline functions is shown in comparsion with other methods for interpolating functions that are set by the table of their values for the case when four values of the interpolated function are given. The interpolating local cubic spline oscillates but shows a deviation from the original monotonic distribution. The best result is obtained with the weighted local cubic spline.
The resulting finite difference spline scheme was used to solve two unsteady problems with the known analytical solution: the "diffusionless" propagation of an impurity and the propagation of an impurity from an instantaneous point source. The following finite difference schemes with different approximations for the convective terms of the equation were compared: the upwind scheme, the Harten scheme, the superbee limiter scheme, MLU, MUSCL, and the 3rd order approximating ENO scheme.
The results of the calculations performed for various density of grid nodes show the convergence of the approximate solution to the exact solution. For the first test problem, the spline scheme is at the advantage of the proximity of the calculated solution to the exact one over the other schemes. For the second test problem, which is characterized by smoother spatial solution profiles, on a coarse grid spline scheme gives solution which is in the best agreement with the exact solution. On a more detailed grid, the best results are given by the MLU and MUSCL schemes. The spline proposed is slightly inferior to them, but in this test example the spline scheme predicts the current maximum concentration more accurately, which is certainly an advantage for the representation of peak concentrations of air pollutants.
SEMYONOVA Anastasiya Alexandrovna (Post graduate student, Tomsk State University, Tomsk, Russian Federation). E-mail: semyonova.a@math.tsu.ru
STARCHENKO Alexander Vasil'evich (Doctor of Physics and Mathematics, Tomsk State University, Tomsk, Russian Federation). E-mail: starch@math.tsu.ru
REFERENCES
1. Samarskiy A.A., Mihaylov A.P. (2001) Matematicheskoe modelirovanie [Mathematical modeling]. Moscow: Fizmatlit.
2. Samarskiy A.A. (1977) Teoriya raznostnykh skhem [Theory of finite difference schemes]. Moscow: Nauka.
3. Godunov S.K. (1959) A Difference scheme for numerical solution of discontinuous solution of hydrodynamic equations. Math. Sbornik. 47(89). 3. pp. 271-306.
4. Lax P.D., Wendroff B. (1960) Systems of conservation laws. Communications in Pure and Applied Mathematics. 13. pp. 217-237. DOI: 10.1002/cpa.3160130205.
5. Kolgan V.P. (1972) Primenenie printsipa minimal'nykh znacheniy proizvodnoy k postroeniyu konechnoraznostnykh skhem dlya rascheta razryvnykh resheniy gazovoy dinamiki [Application of the principle of minimum derivatives to the construction of difference schemes for computing discontinuous solutions of gas dynamics]. Uch. Zap. TsaGI. 3(6). pp. 68-77.
6. van Leer B. (1973) Towards the ultimate conservative difference scheme. I. The quest for monotonicity. Lecture Notes in Physics. 18. pp. 163-168. DOI: 10.1007/BFb0118673.
7. van Leer B. (1974) Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme. J. Comp. Phys. 14. pp. 361-370. DOI: 10.1016/0021-9991(74)90019-9.
8. Noll B. (1992) Evaluation of a bounded high-resolution scheme for combustor flow computations. AIAA Journal. 30(1). pp. 64-68. DOI: 10.2514/3.10883.
9. van Leer B. (1979) Towards the ultimate conservative difference scheme. V. A second order sequel to Godunov's method. J. Comp. Physics. 32. pp. 101-136. DOI: 10.1016/0021-9991(79)90145-1.
10. Toro E.F. (1999) Riemann solvers and numerical methods for fluid dynamics. Berlin, Heidelberg: Springer. DOI: 10.1007/978-3-662-03915-1.
11. Lebedev A.S., Chernyj S.G. (2000) Praktikum po chislennomu resheniyu uravnenij v chastnyh proizvodnyh [Practice on the numerical solution of partial differential equations]. Novosibirsk: NSU publ.
12. Harten A. (1983) High resolution schemes for hyperbolic conservation laws. J. Comp. Phys. 49. pp. 357-393. DOI: 10.1016/0021-9991(83)90136-5.
13. Roe P.L. (1986) Characteristic-based schemes for the Euler equations. Ann. Rev. FluidMech. 18. pp. 337-365. DOI: 10.1146/annurev.fl.18.010186.002005.
14. Harten A., Engquist B., Osher S., Chakravarthy S. R. (1986) Some results on uniformly highorder accurate essentially Non-oscillatory schemes. J. Appl. Num. Math. 2. pp. 347-377. DOI: 10.1016/0168-9274(86)90039-5.
15. Kvasov B. I. (2000) Methods of shape-preserving spline approximation. Singapore: World Scientific Publishing Company.
16. Karpova A.A. (2014) Priblizhenie tablichno zadannyh funkcij s pomoshh'yu approksima-cionnyh i interpolyacionnyh vesovyh splajnov [Approximation of tabulated functions by means of approximation and interpolation weighted splines]. Scientific Conference of Students of the Mechanics and Mathematics Department, 24—30 april 2014. Tomsk. pp. 38-39.
17. Marchuk G.I. (1982) Mathematical models in environmental problems. Amsterdam: Elsevier.
18. Cada M., Torrilhon M. (2009) Compact third-order limiter functions for finite volume methods. J. Computational Physics. 228. pp. 4118-4145. DOI: 10.1016/j.jcp.2009. 02.020.
19. Starchenko A.V., Danilkin E.A., Semenova A.A., Bart A.A. (2016) Parallel algorithms for a 3D photochemical model of pollutant transport in the atmosphere. Communications in Computer and Information Science. 687. pp.158-171. DOI: 10.1007/978-3-319-55669-7_13.
20. Chi-Wang Shu (1996) Essentially non-oscillatory schemes for hyperbolic conservation laws. Preprint of Division of Applied Mathematics. Brown University.
21. Rivin G.S., Voronina P.V. (1997) Aerosol transfer in the atmosphere: selection of a finite difference scheme. Atmospheric and Oceanic Optics. 10(6). pp. 386-392.
