The finite-difference scheme for the unsteady convection-diffusion equation based on weighted local cubic spline interpo.pdf В настоящее время при изучении разнообразных физических процессов широко используются методы математического моделирования, основная идея которых заключается в представлении свойств изучаемого объекта или явления с помощью математических формул, как правило, выражающих тот или иной физический закон. При исследовании процессов в рамках механики сплошных сред применяются законы сохранения массы, импульса, энергии, которые, в общем случае, имеют вид нестационарных интегральных или дифференциальных уравнений, включающих члены, описывающие перенос массы, импульса или энергии как за счет непосредственно самого движения сплошной среды (конвекция), так и за счет молекулярного или турбулентного обмена (диффузия) [1]. Как правило, получение аналитических решений таких уравнений на современном уровне развития теоретической математики представляет определенную проблему, поэтому активно разрабатываются приближенные численные методы решения, обеспечивающие выполнение условий аппроксимации, устойчивости и сходимости [2]. Особое внимание при решении нестационарного уравнения «конвекции - диффузии» уделяется выбору аппроксимационной схемы для конвективного члена уравнения. С одной стороны, требуется, чтобы порядок аппроксимации схемы был выше первого для обеспечения качественного воспроизведения резких изменений переносимой субстанции во времени и пространстве, с другой -используемая аппроксимация (или разностная схема) должна быть монотонной, т.е. обеспечивающей получение монотонного распределения искомой величины из имеющегося монотонного распределения для предыдущего момента времени. Современные численные методы решения уравнения переноса берут свое начало в том числе и с известной работы Годунова [3], в которой он предложил численный метод для решения локальной проблемы Римана на границе раздела двух объемов сплошной среды. Им представлена очень важная теорема о монотонности: не существует монотонных линейных разностных схем с порядком аппроксимации по пространству выше первого. Монотонность является очень желательным свойством разностной схемы, ее наличие предотвращает образование вычислительных осцилляций в окрестности сильных градиентов решения. Лаке и Венд-рофф [4] предложили для решения уравнения переноса схему второго порядка аппроксимации, обеспечивающую хорошее предсказание резко меняющегося поведения приближенного решения, но сильно осциллирующую. Для уменьшения вычислительных осцилляций необходимо в эту схему вводить дополнительные дис-сипативные члены. В настоящее время для того, чтобы обойти жесткое ограничение теоремы Годунова, используется построение нелинейных разностных схем для решения уравнения переноса. В качестве одного из способов построения нелинейных схем является применение кусочно-полиномиальной реконструкции приближенного решения по значениям сеточной функции с введением некоторых ограничений (условий), обеспечивающих уменьшение вычислительных осцилляций при переходе на следующий момент времени. Одной из первых работ этого направления была работа Колгана [5], который предложил применять принцип минимальных значений производных, чтобы получить неосциллирующую схему второго порядка типа схемы Годунова. Другой известный подход предложил ван Лир в своей серии статей [6, 7, 9], в которой развивается теория разностных схем высокого порядка аппроксимации на основе полиномиальной реконструкции и контроля за значениями производных этих полиномов с использованием TVD-свойства (свойства невозрастания полной вариации решения при переходе на новый слой по времени) [10, 11]. На основе полиномиальной TVD-реконструкции ван Лиром, в частности, были предложены схемы MLU (Monotonic Linear Upstream) [7, 8], MUSCL (Monoto-nic Upstream centered Scheme for Conservation Laws) [9], которые используют кусочно-линейную интерполяцию и обеспечивают точность второго порядка в случае гладкого монотонного решения. Подход TVD ввел новые усовершенствования в построении схем высокого порядка точности для уравнения переноса и помог избежать ограничений для линейных разностных схем, предоставляемых теоремой Годунова. Но порядок этих схем ограничивается вторым на участке гладкого монотонного изменения решения (см. например, схему Хартена с ограничителем minmod [12] и схему с ограничителем superbee [13]). Для повышения порядка аппроксимации разностных схем для решения уравнения переноса были введены так называемые существенно неосциллирующие схемы ENO (Essentially Non-Oscillatory schemes). Схемы ENO для получения ос-редненного по ячейке значения решения содержат алгоритм построения нескольких полиномов высокого порядка на различных сеточных шаблонах, выбирающих наиболее гладкий из полиномов. Такой подход был введен в 1986 г. Harten, Osher, Engquist и Chakravarthy [14]. При численной реализации моделей атмосферного пограничного слоя и переноса примеси важную роль играет выбор качественной разностной схемы для аппроксимации конвективного члена уравнения. Целью данного исследования является разработка построенной на локальных весовых сплайнах разностной схемы для аппроксимации конвективных членов уравнения переноса. Полученная аппроксимация используется в явной разностной схеме, которая сравнивается с другими известными аппроксимациями конвективных слагаемых (MLU, MUSCL, ENO и др.) при численном решении двух одномерных нестационарных тестовых задач о распространении примеси. 1. Интерполяционный весовой кубический сплайн _ Г _я-1 ] Введем сетку сой = ■< xt■; xi+l = xi + ht■, / = 0, и -1; x0 = a; ^ ht■ = b - a >■ на отрезке [д./)]. На каждом элементарном отрезке [хг,хг+1] интерполяционный весовой кубический сплайн есть полином третьей степени ,S'( (х): Si (Х) = ai0 + ai\(Х ~Х,) + а,2 (Х ~ Х, f + ai3 (Х ~ Х, )3 > хе[хг,хг+1], коэффициенты которого надо определить из следующих условий [15]: £;(*;) = fi,'=0,n; S>(x1) = S'_1(x1)J = l,n-l; w^1(x1) = w1S»(x1)J = l,„-l. Здесь \г( > 0. /' = 0, п -1, - весовые коэффициенты. Построение весового кубического сплайна будет проводиться через наклоны (mi = S'(xi). / = 0. n ). для получения единственного решения поставленной задачи будем использовать дополнительные условия на границах: S'(x0) = /0', S\x„ ) = /„'. Из приведенных выше условий находим формулу для построения сплайнов через наклоны: h2 hf h2 2m, +m, !+1 (х-х;)2 +mi(x-xi) + f{, хг *i] ' Wi-1 /[*i>*i+l] то S'(x) > О для всех x e [a, ft], т.е. сплайн S монотонен на \a,b\. В [16] для монотонно убывающих данных был сформулирован аналог этой теоремы. Для нахождения весовых коэффициентов можно следовать алгоритму [15]. Пусть весовой коэффициент w^ известен. Полагаем \г( = м>г. Проверяем неравенства из теоремы. Если какое-то из них нарушается, то выражаем из соответствующего неравенства wi, заменяя знак неравенства на знак равенства. Алгоритм начинаем с w0=l и легко находим все коэффициенты {\г( j. обеспечивающие монотонность весового кубического сплайна для произвольных монотонных данных. Если на некотором шаге \г( < е, > - j. то полагаем \г( = е, = - j , где е - достаточно малое положительное число, позволяющее предотвратить зануление (переполнение). Рис. 1. Интерполяция таблично заданной функции алгебраическими многочленами Fig. 1. Algebraic polynomials interpolation Продемонстрируем преимущество весовых кубических сплайн-функций при интерполировании таблично заданных функций. Рассмотрим случай, когда заданы таблично четыре значения функции. Требуется восстановить значение интер-полянты для промежуточных значений х [16]. Для решения этой задачи были применены два способа: кубический интерполяционный сплайн, построенный через наклоны, и такой же сплайн, но с весами, подобранными по алгоритму, описанному выше. При построении сплайнов на концах выбранного отрезка использовались дополнительные условия равенства нулю производной интерполируемой функции. Заметим, что табличные значения функции выбирались таким образом, чтобы при интерполировании могла возникнуть проблема монотонности построенных интерполяционных многочленов. Данный набор табличных значений интерполируемой функции соответствует профилю аналитического решения задачи «бездиффузионного» переноса примеси, рассмотренной ниже (Тест 1). Из рис. 1 видно, что интерполяционный кубический сплайн осциллирует на рассматриваемом промежутке между табличными значениями, предсказывая отрицательные значения на отрезке [1,0; 2,0] и больше единицы на отрезке [3,0; 4,0]. Лучший результат демонстрирует интерполяционный весовой кубический сплайн. В данном случае осцилляций визуально не наблюдается и сглаживание в области резкого изменения значений интерполируемой функции меньше, чем в предыдущем случае. 2. Численный метод решения уравнения переноса Используем полученный выше результат для весового кубического сплайна при построении разностной схемы со сплайн-аппроксимацией конвективного члена уравнения переноса примеси. Рассмотрим одномерное нестационарное уравнение конвекции - диффузии с постоянными коэффициентами для концентрации примеси: 5Ф 5Ф д2Ф -+ u-= D-- + БФ, и > О, D >= 0; 00, [ф-+1 - 0,5 max [о, min(2©, (2 + ©) / 3,2)] (®f+1 - Ф- ), и < 0, Ф^ г+1/2 Ф^ г+1 / 2 0, [ф*+1 - 0,5(1 + Kh) min mod (ф*+1 - Ф*, Ф* - Ф*_,), и < 0 и ограничителя superbee [13]: - 0,5(1 - Kh) max mod [min mod (2 (®f+1 - Ф- ), Ф- - Ф^), min mod ( 0, Ф*+1 - 0,5(1 + Kh) max mod [min mod (2 (ф*+2 - Ф*+1), Ф*+1 - Ф*), minmod(®f+2-Ф*+1,2(ф*+1 -Ф*))], и < 0. Здесь Kh = ux/h. Тест 1. Решение уравнения переноса в случае «отсутствия» диффузии (D - 0) и Зф = 0 . Начальное условие для задачи 1 (Х= 2): ф (1,хе [0.25, 0.75], ooW~ К* £[0.25,0.75]. Граничное: Ф0 (0 = 0-Аналитическое решение задачи 1 имеет вид [11] Фехас^*) = Ф00 Расчеты производились при следующих условиях: U= 1, М = 100, М = 200, тм//г = 0.2, D = 0 . I-, маанаааань О II 0.8g 0.6- * 0 0.4 0.8 1.2 1.6 2 х, м Рис. 3. Рассчитанные (кружки) по сплайновой схеме значения концентрации примеси при t = Т = \с Fig. 3. Calculated values (circles) of the impurity concentration calculated at / = T = lc by the spline scheme a Из рис. 3 видно, что численное решение достаточно хорошо соответствует точному, однако в областях резкого изменения решения наблюдается некоторое «размытие» профиля |ф( j. Данная задача решалась и с использованием перечисленных выше методов аппроксимации конвективного члена уравнения переноса. Для сравнения степени согласования расчетов с точным решением рассматривались нормы погрешности: I х К-Ф(п*г)|2 Нх = max bf -Ф(Т,хг)|, Н2 \ --^-, ;=о,...,м1 г г I 2 V М значения которых для рассматриваемой задачи представлены в табл. 1. Таблица 1 Сравнение норм погрешности различных методов аппроксимации конвективных членов нестационарного уравнения переноса для первой тестовой задачи Метод М= 100 М= 200 Я2 Я2 1 Upwind 0,48747 0,17209 0,49108 0,14451 2 MLU 0,32540 0,05414 0,32540 0,03828 3 MUSCL 0,32989 0,05468 0,32898 0,03867 4 EN03 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 Superbee 0,37625 0,06941 0,38222 0,05024 На основании полученных данных можно сделать вывод, что для рассматриваемой тестовой задачи применение весовой сплайн-функции для аппроксимации конвективных членов имеет преимущество перед остальными схемами. Худший результат показывает противопотоковая схема. Эти результаты также согласуются с теоретическими данными о порядке аппроксимации этих схем. Сходимость схемы проверялась удвоением количества узлов. Нормы погрешности методов Н2 при этом уменьшились. Тест 2. Распространение примеси от мгновенного точечного источника. Рассмотрим следующую нестационарную задачу, моделирующую распространение примеси от мгновенного источника, расположенного внутри области [17]: ЗФ 5Ф ^ ^ -+ и-+ стФ = £> dt дх t = 0: Ф(0,х) = О, Ф(£0) = 0, ЗФ - (t,X) = 0. дх Здесь 5(z) - дельта-функция. Точное решение данной задачи имеет вид [17] д2Ф дх2 + Q5(x-x0)5(t-t0), 0 < х < X, 0

Скачать электронную версию публикации