Содержание к диссертации
Введение
Глава 1. Многочлены Чебышева и оптимизация численных методов . 19
1.1. Многочлены Чебышева и их свойства 19
1.2. Явные разностные схемы и решение жёстких задач математической физики 23
1.2.1. Формулировка явных методов 24
1.2.2. Общий анализ явных разностных методов с переменными шагами по времени 27
1.2.3. Метод Адамса-Бэшфорта 34
1.2.4. Явные разностные схемы для мнимого спектра 38
1.3. Блочные чебышевские итерационные методы 41
1.3.1. Чебышевский набор параметров 41
1.3.2. Устойчивое упорядочивание итерационных параметров 46
1.3.3. Критерий окончания итераций 54
1.3.4. Поиск улучшенного начального приближения методом Ритца. 55
1.4. Выводы по главе 1 60
Глава 2. Вихреразрешающая модель, разностные аппроксимации и численные эксперименты . 61
2.1. Описание модели 61
2.2. Разностная аппроксимация уравнений модели 66
2.2.1. Вихревая форма записи адвективных слагаемых 69
2.3. Программная реализация модели 75
2.3.1. Задание расчётной области с помощью логических масок 76
2.3.2. Исследование операторов Ca:Da и Т 79
2.3.3. Внедрение программы DUMKA 80
2.3.4. Распараллеливание 83
2.4. Численные эксперименты 85
2.5. Выводы по главе 2 92
Глава 3. Чебышевские цифровые фильтры 93
3.1. Явный чебышевский фильтр 94
3.2. Алгоритм неявной чебышевской операторной фильтрации 96
3.3. Особенности численной реализации неявного чебышевского фильтра 99
3.4. О сравнении фильтров 101
3.5. Выводы по главе 3 105
Заключение
- Явные разностные схемы и решение жёстких задач математической физики
- Чебышевский набор параметров
- Разностная аппроксимация уравнений модели
- Алгоритм неявной чебышевской операторной фильтрации
Введение к работе
Настоящая работа посвящена изучению особенностей применения в задачах гидродинамики численных методов, общей особенностью которых является удобная представимость через многочлены перехода. Анализ свойств этих методов, таким образом, сводится к анализу свойств многочленов. Так, при применении явных разностных схем для нестационарных задач математической физики решение на каждом слое по времени в случае линейной правой части представляет собой результат действия некоторого многочлена от оператора правой части на начальное значение решения. Считается [4, 34, 41], что простые тестовые задачи такого типа в достаточной степени отражают с точки зрения устойчивости характеристики того или иного явного метода. Для жёстких задач математической физики именно ограничения на устойчивость, как правило, определяют величину шагов по времени. Поэтому, имея данные о расположении спектра задачи, естественно будет выбрать для её численного решения разностный метод, определяемый многочленом, область устойчивости которого (область комплексной плоскости внутри единичной лемнискаты) содержит в себе этот спектр при как можно большем среднем шаге по времени. В совокупности с возможностью хорошего распараллеливания явных схем это даёт основания ожидать высокой эффективности использования современных вычислительных ресурсов. В данной работе мы ограничимся исследованием именно явных схем.
Близкие идеи находят применение при построении линейных итерационных методов [17, 9], когда задача оптимизации сводится к построению многочлена перехода для серии итераций, наиболее эффективно по-
давляющего ошибку. В качестве такого многочлена обычно берётся многочлен, обладающий некоторыми экстремальными свойствами на множестве, содержащем спектр оператора решаемой задачи.
Общими проблемами в этих двух областях приложения теории многочленов являются, во-первых, организация вычислений, устойчивых по отношению к ошибкам округления внутри серии шагов метода, реализующей выбранный многочлен перехода, и во-вторых, получение количественных данных о спектре задачи. Первая из них в значительной степени решена с помощью рекурсивных алгоритмов построения устойчивых перестановок параметров [23]-[25], [17]. Вторая является гораздо более трудной, её детальное решение может оказаться сравнимо по сложности с решением исходной задачи. Однако часто можно дать общее описание спектра, исходя из физических или каких-либо иных соображений, и выбрать многочлен перехода, по возможности близко соответствующий этому описанию. При этом выборе бывает очень полезно учитывать также свойства многочленов, касающиеся удобства реализации, объёма требуемых для осуществления шага вычислительных ресурсов, а также возможности использования данных, полученных на предшествующих шагах [48, 49, 50, 72]. Перспективной темой для исследования является учёт в расчётах эффектов взаимного влияния методов различного назначения, использующих общую идею оптимизации многочленов перехода.
Идеи, аналогичные чебышевским итерационным методам, могут с успехом применяться также в цифровой фильтрации [44]. Например, в [16] для сглаживания решения или правой части системы уравнений был использован явный фильтр второго порядка, функция перехода которого обладает свойством альтернанса в области подавления. В [62] для улуч-
шения устойчивости счёта к правой части применялись псевдоразностные операторы, аналогичные исследовавшимся в [38].
Практические задачи гидродинамики, выбранные в качестве области приложения исследуемых методов, являются динамично развивающимся разделом численного моделирования [45]. Их решение осложняется многими факторами, в частности - наличием нелинейных членов. В этом свете представляет интерес вихревая форма записи таких слагаемых (форма Громеки-Лэмба), стабилизирующие свойства которой исследовались, например, в работах [30, 65], а также в [56]. Она позволяет исключить существенную часть нелинейных слагаемых из проблем построения аппроксимации и турбулентного замыкания.
Цель диссертационной работы состоит в исследовании эффективности и оптимизации применения устойчивых явных разностных схем программы DUMKA, методов с функциями перехода, основанными на че-бышевских многочленах, и вихревой формы Громеки-Лэмба в задачах гидродинамики.
Диссертация состоит из введения, трёх глав, заключения и списка литературы, включающего 74 наименования.
Работа выполнена в Институте вычислительной математики РАН. Автор выражает искреннюю благодарность своему научному руководителю В.И. Лебедеву за поддержку и постоянное внимание при выполнении работы, А. В. Глазунову за продолжительные консультации по вопросам вихреразрешающего моделирования и признателен Т.И. Шитовой, Ю.М. Нечепуренко и Р.А. Ибраеву за полезные советы по оформлению диссертации, а также сотрудникам ИВМ РАН за помощь в период выполнения работы.
Содержание работы Глава 1 диссертации посвящена теоретическому исследованию проблем
оптимизации многочленов перехода в свете решения задач гидродинамики. Для составления основы рассматриваемых методов вначале приведены базовые сведения о многочленах Чебышева, касающиеся возможных вариантов определения, формул для корней, точек экстремума и композиции многочленов, формулировки экстремальных теорем Чебышева и Маркова.
В данной работе исследование численных методов для решения нестационарных задач ограничено явными разностными схемами по времени, которые обычно легко описываются при помощи многочленов перехода для линейной задачи
du .
— = -Ait, u\t=o = Щ (1)
от начального условия к решению на выбранном слое по времени:
uN = PN{A)uQ. (2)
Для задач, приводящихся линеаризацией и отбрасыванием свободного члена в правой части к виду (1) со спектром, лежащим на или тесно вблизи мнимой оси (например, расчётов течений с большими числами Рейнольдса [14, 36, 66]), в серии работ [4, 15, 16, 17, 59, 60] предложено использовать явные разностные схемы DUMKA-7 и DUMKA-8. Этим схемам соответствует многочлен перехода 4-й степени, аппроксимирующий экспоненциальную переходную функцию решения дифференциальной задачи с 5-м порядком точности. Его корни комплексно сопряжены,
поэтому соответствующие шаги могут быть реализованы в действительной арифметике. В диссертации приведён анализ области устойчивости этих схем из указанной серии работ. Согласно [57], соответствующий многочлен перехода обладает наибольшим для своей степени отрезком области устойчивости на мнимой оси при фиксированной сумме временных шагов, а в [15] показано, что при увеличении степени многочлена перехода для задач с наличием мнимого спектра увеличение среднего шага ограничено. Поэтому для задач с большими числами Рейнольдса оправдано применение схем DUMKA-7 и DUMKA-8, реализованных в свободно доступной программе DUMKA. Следует отметить также, что, согласно [28], возможно применение и схемы DUMKA-6 [59, 4], но в этом случае требуются дополнительные данные о спектре для настройки её параметров.
Модель [7] использует для интегрирования по времени схему Адамса-Бэшфорта, описанную, например, в [41]. Для этой схемы в диссертационной работе проведён аналогичный анализ, показывающий, что её область устойчивости касается мнимой оси в нуле. При этом приближённое решение un, соответствующее мнимым собственным значениям Л оператора А, может.расти при малых |А|, как 1 + А^|А|4/4. Показано также, что в случае наличия действительного спектра использование схемы Адамса-Бэшфорта может приводить к возникновению в решении двухшагового шума. На основании этих результатов можно сделать вывод о нежелательности применения этой схемы напрямую для интегрирования по времени в моделях турбулентности.
Также в первой главе исследуется блочный чебышевский итерационный метод [17] решения уравнения Пуассона для давления. Записанный
для уравнения в "красно-чёрном" [31] виде
\-С2 D2 ) \ Р2 J ~ \ h J'
его шаг представляет собой итерацию Гаусса-Зей дел я с поправкой
Dip4+1 = ClPk2 + h,
>й>*+1/2 = Crt1 + Л,
рк2+1=р1 + аМ(рк2+1/2-ркг).
Изменение ошибки за цикл из N итераций описывается многочленом перехода Рдг для второго блока неизвестных, так что
^М^П^1-^)' 4 = PN{I-T)sl Т = D^^D^C,. (4)
г=1
Для применяемой в модели дискретной записи уравнения Пуассона в диссертации доказано с использованием результатов работы [43], что оператор I — Т обладает базисом из собственных векторов и его спектр лежит на некотором отрезке действительной оси [q, l],q > 0. Пусть с помощью дополнительных данных получена оценка спектра Sp(/ — Т) Є [га, М] (например, га = q, М = 1). Тогда в качестве P/v возьмём нормированный в соответствии с первым равенством (4) многочлен Чебышева 1-го рода, приведённый к отрезку [га, М]. Этот многочлен наименее отклоняется от нуля на данном отрезке, обеспечивая таким образом наилучшее среди многочленов вида (4) подавление ошибки. Соответствую-
щее отклонение равно
2N-
1 + сг
где а =
1- л/т/М 1 + л/т/М
и уменьшается с ростом iV, поэтому имеет смысл использовать достаточно большие итерационные циклы. Итерационные параметры а^ являются величинами, обратными к корням многочлена перехода, которые для многочленов Чебышева вычисляются по простым известным формулам. Устойчивость этого метода по отношению к ошибкам округления существенно зависит от порядка употребления параметров ск&. Для её обеспечения необходима такая перестановка параметров, чтобы величины
г? =
t[m,M]
П(і - ajt)
i=i
,N
t[m,M]
П (1 -a^)
j=i+l
(qjy = 1), а также Еі=і ?f > Еі=і ?f ff были ограничены постоянными, зависящими только от т/М [17]. Метод получения таких перестановок путём разложения многочлена перехода в произведение дробно-рациональных функций, зависящих от тригонометрического параметра, в случае N = 2Р3Г предложен в работах [23]-[25]. Его можно упростить, исключив анализ тригонометрических параметров, с помощью следующей доказанной в диссертации леммы. ЛЕММА 1. Пусть \р\ < 1, тогда
Пп{х) - (3 Тзп(в) -
= апЬпсп,
_ Тп{х) - рг _ Тп(х) - р2 _ Тп(х) - р3 М -Утп
йп - Тп(Э) -Р1'Ьп~ Г„(Є) - р2' Сп " Тп(0) - рз' М - т'
ІР1.2.3І < 1, Pi > 0, рз < О, sign(p2) = -sign(/3), рз< Р2< /01,
причём
^п < Ьп < сп, |сп| < 1, |спЬп| < 1.
Основываясь на лемме, сформулируем метод упорядочивания итерационных параметров в виде последовательного разбиения (перед приведением к отрезку [т, М]) многочленов перехода тпх^[_р степени п\, кратной 3, на множители (многочлены) степени щ/З при р < 0 в порядке спі/з> ^щ/З) апі/з- Если же в разбиваемом многочлене р > 0, то множители И СООТВеТСТВуЮЩИе ИМ КОрНИ будем браТЬ В ПОрЯДКе ЬП1/з, СП1/з, СЬщ/З-
Приведён также алгоритм аналогичного разбиения многочленов перехода чётной степени. Осуществляя последовательно эти два алгоритма, мы реализуем метод, дающий последовательность устойчивых перестановок параметров длины N = 2Р3Г, становящуюся оптимальной по Чебышеву после итераций с номерами 1, 2,4, 8,..., 2Р, 3 21\ 9 2Р,..., 3Г2Р.
Дополнительно в диссертации предложен и протестирован следующий приём удвоения длины устойчивой перестановки параметров: если Pjk - корень многочлена Тп(х), то пара ід/-2^— - корни многочлена Т2п(х). Корень с минусом определяет итерацию, многочлен перехода которой меньше единицы по модулю на отрезке [m, М]. Поэтому её будем осуществлять первой в соответствующей паре итераций.
Поскольку первое уравнение из (3) после каждой итерации выпол-
няется точно, в качестве критерия остановки итераций возьмём малость ошибки второго уравнения, разрешенного относительно р2. Соответствующая невязка определяется простой формулой, не требующей дополнительного вычисления каких-либо операторов
к к к+1/2
Г2=Р2-Р2 >
а для ошибки справедливо неравенство
\\4\\<т-т)-1\\ш<^-1ш
в некоторой норме || ||. Выбор в качестве индикатора сходимости нормы ошибки (а не невязки) является желательным, т.к. норма невязки может уменьшаться в процессе счёта не только вследствие сходимости, но и из-за изменения спектрального состава. Определение величины т является сложной задачей, она может быть приближённо решена опытным путём с помощью наблюдения скорости сходимости во вспомогательном численном эксперименте с решением заданного уравнения (3).
Для ускорения сходимости блочного чебышевского метода можно применить один из алгоритмов нахождения начального приближения по результатам предшествующих временных шагов. Например, в [72] для этого используется алгебраический подход, основанный на накоплении больших объёмов информации о векторах из подпространств Крылова в методе обобщённых сопряжённых невязок. Мы применим "физический" подход, используя полученную ранее информацию о характере поля давления, не требующую больших объёмов памяти.
Если записать уравнение Пуассона на временном новом слое для поправки к давлению
A5P(t + т) = F(t + т) - F(t) = Ф(* + г), SP(t + г) = P(t + т) - P(t),
с соответствующими краевыми условиями, то начальное приближение для SPft + r) можно взять равным a(t + r)SP(t), где коэффициент определить методом Ритца, используя знакоопределёный оператор из уравнения для второго блока неизвестных. Дополнительно учтём, что схема DUMKA-8, реализующая метод Гилла 4-ого порядка, предполагает использование нулевых и ненулевых шагов по времени. Определим для каждого из этих двух типов шагов начальное приближение по поправке с предыдущего шага того же типа. Отметим, что можно также использовать метод наименьших квадратов вместо метода Ритца, а в качестве второго базисного вектора взять величину div PU, имеющую смысл "переноса" давления потоком жидкости.
Явные разностные схемы и решение жёстких задач математической физики
Как известно, существуют явные и неявные методы решения нестационарных задач математической физики. В неявном методе для нахождения Uk+i приходится решать некоторое уравнение. Явные методы свободны от этого недостатка, однако для их устойчивой реализации могут потребоваться сильные ограничения на временные шаги.
Рассматриваемые в данном разделе явные разностные схемы с переменными шагами по времени, определяемыми как обратные величины к корням многочленов, аппроксимирующих оператор перехода на новый временной слой, были разработаны [15, 4, 16] с целью учесть в едином комплексе проблемы, возникающие при построении методов решения жёстких систем обыкновенных дифференциальных уравнений большого порядка, разностных методов решения нестационарных задач математической физики и методов распараллеливания алгоритмов для многопроцессорных вычислительных систем. Набор ссылок на литературу по этим проблемам приведён, например, в работе [60]. Примеры оптимизации численных алгоритмов, приводящие к экстремальным многочленам, можно найти в [41] и [5]. Анализируемые в разделе 1.2.4. схемы с переменными шагами по времени входят в состав программы DUMKA [15, 17].
Возникающие проблемы, связанные с реализацией явных методов, мы проанализируем на модельной задаче Коши для уравнения типа Далкви-ста [17, 41, 34], к которой придём, заменив (1.18) на линейную однородную задачу Коши du — = -Аи, u\t=to = и0, (1.22) при 0 t Т, где А - постоянная матрица (п х п). Главные проблемы, связанные с аппроксимацией и устойчивостью при применении явных разностных схем, возникают и при решении таких задач. Пусть спектр оператора А лежит в комплексной плоскости в круге Кг минимального радиуса г 0 вида Kr = {z: \z-r\ r}. (1.23) Пусть (A{, (pi) - собственные пары оператора А и { pi} образуют базис в рассматриваемом пространстве, аМ = 2г. Пусть Щ = 2а№, (1-24) і тогда решение задачи (1.22) даётся выражением u(t) = ехр(-А;)а;с , (1.25) і причём норма каждого из слагаемых в (1.25) не растёт со временем. В условиях принятых предположений назовём задачу Коши жёсткой, если Тг 1. Величина г для задач математической физики обычно имеет порядок 0(h p), где р 1, a h - характерный размер шага сетки [4]. Формула явного метода (1.21) для задачи (1.22) примет вид ик+1 = (I - тк+1А)щ. (1.26)
Учитывая вид решения (1.25) задачи (1.22), естественно потребовать, чтобы спектральный радиус оператора перехода на следующий временной слой в (1.26) при постоянном шаге Tk+i = г не превышал единицы, ибо в противном случае мы не только можем потерять свойство аппроксимации, но и столкнёмся с катастрофическим ростом нормы приближённого решения. Для оператора перехода в (1.26) это означает, что т сои = 1/г, (1.27) т.к. S{I — тА) = тах(1 — тМ\, 1). Величина сои получила в [15] название куранта. Сопоставляя (1-27) с условием Тг 1, мы видим, что при их выполнении для получения решения задачи Коши потребуется очень много временных шагов, что побуждает обратиться к явным методам с переменными шагами г/с.
Общий анализ явных разностных методов с переменными шагами по времени. Разностные схемы с переменными шагами по времени давно применяются в вычислительной математике, однако две проблемы, связанные с их использованием, исследованы сравнительно мало: это связь между условиями аппроксимации и устойчивости и устойчивая реализация вычислений внутри цикла, соответствующего набору переменных шагов. Решение второй из них предложено в [23, 24, 25] и рассмотрено в разделе 1.3. применительно к блочному чебышевскому итерационному методу, описанию первой посвящен настоящий раздел.
Отделение модуля многочлена перехода от единицы в случае ддг 1 может вводиться для обеспечения гарантированной устойчивости счёта в случае больших степеней iV, а также для лучшего приближения быстро убывающих компонент решения дифференциальной задачи. Будем предполагать также, что выполнено условие устойчивой реализации вычислений внутри цикла по отношению к ошибкам округления.
Поскольку ljy = —P N(0), для наибольшего продвижения по времени за JV шагов требуется найти решение обобщённой задачи Маркова: PN{z) = argsup(-i4(0)), (1.33) где точная верхняя грань берётся по всем многочленам вида (1.30), удовлетворяющим условиям аппроксимации (1.29) и устойчивости (1.32). Если обозначить через 2 корни экстремального многочлена (1.33), то в соответствующем оптимальном методе (1.26) множество временных шагов представляет собой множество чисел z 1.
Чебышевский набор параметров
Предположим, что операторы D\ и D2 таковы, что уравнения D{V{ = ді, і = 1, 2 легко решаются при произвольных gi. Операторы С\,С2 осуществляют связь между неизвестными и\ и и2, их будем считать ограниченными. Замечательные свойства многочленов Чебышева позволяют успешно применять их при оптимизации итерационных методов решения линейных уравнений. Итерационные методы, в которых используются свойства или параметры наименее отклоняющихся от нуля многочленов, будем в соответствии с [17, 9] называть чебышевскими. В литературе широко описаны чебышевские методы, учитывающие особенности спектра задачи. В данной работе сделан акцент на применении одношагового че-бышевского метода к вихреразрешающей модели [7] с учётом блочной структуры использованных в ней операторов.
Базовыми итерационными методами модели являются метод последовательной верхней релаксации и метод симметричной релаксации, применяемый самостоятельно или в качестве предобуславливателя для метода сопряжённых градиентов [31]. Здесь Ах, Ay, Az - шаги сетки вдоль координатных направлений, а оператор Са действует на массив величин щ а, но выдаёт значения в те ячейки, где определены иа. Справедлива основанная на результатах [43] лемма: Если оператор А в (1.48,1.49) симметричен и положительно определён, то оператор I — Т обладает базисом из собственных векторов, и его собственные значения принадлежат отрезку [q, 1], 0 q 1.
Из симметричности А следует, что матрицы D\ и D2 симметричны и С\ — С2 , а из положительной определённости А - положительная определённость симметричных матриц 1 і и 1 2.
Возьмём в качестве P/v многочлен, обладающий наименьшим отклонением от нуля на отрезке [т, М] среди многочленов вида (1.53), т.е. среди всех многочленов с действительными корнями, принимающих в точке нуль единичное значение. Согласно теореме 2, этот экстремальный многочлен можно задать формулами (1.12-1.13).
Устойчивое упорядочивание итерационных параметров. При использовании одношагового чебышевского метода важной проблемой является его устойчивость по отношению к ошибкам округления, которая существенно зависит от порядка употребления итерационных параметров [23]-[25],[35], поскольку при малых значениях отношения т/М (большом числе обусловленности А) операторы перехода итераций / — ajA с малыми j будут иметь очень большую норму.
Один из способов рекуррентного построения перестановок параметров, удовлетворяющих этому требованию, получен в [23]-[25]. Так, параметры а/г в (1.12), являясь обратными к корням величинами, формируют чебышевский многочлен перехода, приведённый к отрезку [га, М] и нормированный делением на T/v(G).
Поскольку \(3\ 1, то в силу свойств (1.5-1.6) уравнение Тз(р) = (5 имеет три различных действительных корня pi, р2 Рз, Для которых справедливо (1.59), как иллюстрирует рис. 1.5. Следовательно, для многочлена Тзп{х) — /3 справедливо разложение на множители (1.57, 1.59), откуда очевидно следует (1.58, 1.60). Двойное неравенство в (1.61) следует из того.
Многочлен перехода Ъп = г"(е\ оптимален среди многочленов степени п (обладает наименьшей нормой в С[— 1,1]), поэтому естественно будет в формируемой перестановке корней Т- п отвести корням Тп(х) первую треть позиций. Следовательно, в (1.12-1.13) параметрам а&, к = 1,...,п будут соответствовать jk Є {2, 5, 8,..., 3n — 1}. Действительно, поскольку корни Тдг(ж) определяются формулой (3j = cos ("2 -71").
Разностная аппроксимация уравнений модели
Опыт показывает, что адвективные члены должны сохранять кинетическую энергию, если мы хотим, чтобы моделирование нестационарного несжимаемого течения было устойчивым, свободным от численной диссипации и не происходило необоснованного перераспределения энергии по спектру [7, 64]. В моделях, сформулированных в терминах скорость-давление, часто используются смещённые сетки. Они позволяют получить физически корректное поле давления при экономном использовании ресурсов памяти. Используемые в данной работе октаэдральные сетки (впервые предложенные в [18, 19, 20], английское название - C-grid [46]) с постоянными шагами вдоль координатных направлений предполагают определение компонент скорости в центрах соответствующих граней ячейки-параллелепипеда, а давления - в центре ячейки.
Уравнение переноса кинетической энергии получается скалярным умножением (2.1) или (2.3) на вектор скорости. Известно [64], что косо-симметричная форма записи адвективных членов консервативна априори в уравнении переноса кинетической энергии, а дивергентная и адвективная будут консервативны в случае выпонения условия неразрывности.
Вихревая форма записи членов переноса в уравнениях гидродинамики в переменных скорость-давление (форма Громеки-Лэмба) известна давно, но в настоящее время в численных расчётах применяется относительно редко по сравнению с потоковой и дивергентной формами. Подробное исследование численных свойств вихревой формы при решении методом конечных элементов может быть найдено в [65], [30] и по ссылкам в этих работах. Целью данного раздела является исследование эффективности этой формы в конечноразностных задачах с точки зрения устойчивости счёта.
На дифференциальном уровне вихревая форма эквивалентна адвективной. Однако при численных расчётах эквивалентность будет нарушаться в силу отличий дискретизации, и запись (2.19-2.20) может обладать новыми полезными свойствами. Так, при переходе к обобщённому давлению существенная часть нелинейных слагаемых исчезает из проблем построения аппроксимации и турбулентного замыкания. В частности, в (2.20) уничтожаются диагональные члены, которые, как показывает опыт, являются одной из главных причин возникновения неустойчивости.
Выражение (2.19) представляет собой векторное произведение, ортогональное скорости U в каждой точке области. Такая запись нелинейного члена представляется наиболее приспособленной для применения центральноразностнои аппроксимации по пространству, которая не вносит дополнительной (схемной) вязкости в процесс счёта и обладает более высоким порядком аппроксимации, что является преимуществом перед односторонними схемами аппроксимации, в особенности при больших числах Рейнольдса. Ортогональность будет также обеспечивать консервативность схемы по энергии, причём в отличие от случая кососимметричной формы плотность кинетической энергии будет сохраняться схемой в каждой точке, а не только в интеграле по всей области.
Перейдём к описанию численной реализации вихревой формы. Основная проблема состоит в том, что в (2.20) скорость Uj и её производные dUi/dxj,dUj/dxi при обычной аппроксимации центральными разностями оказываются определенными в разных точках в силу смещённости сеток. В качестве её решения можно предложить два метода: 1) аппроксимация оператора ротора с помощью дополнительных интерполяций, 2) замена его на другой оператор, эквивалентный на дифференциальном уровне и легко реализуемый на смещённых сетках.
Эта величина определена в центрах ячеек сетки. Уравнение, описывающее изменение К со временем, можно получить, умножив г-е уравнение (2.3) на г-ю компоненту скорости, проинтерполировав результат в центр ячейки с помощью оператора и затем просуммировав по всем г = 1, 2, 3.
Разностный аналог свойства ортогональности вихревой формы и вектора скорости описывается следующей леммой: ЛЕММА 2. При достаточно гладком поле скорости выражение (2.23) есть малая величина порядка 0(h2), где h - максимальный из шагов сетки вдоль координатных направлений.
Будем постепенно раскладывать правую часть (2.23) по формулам Тейлора, чтобы получить выражение, не содержащее разностных операторов, а также производных от нелинейных функций.
Алгоритм неявной чебышевской операторной фильтрации
В рамках диссертационной работы построено семейство неявных чебы-шевских операторных фильтров, требующих решения системы уравнений для получения фильтрованного массива. Их переходная функция оказывается ближе к единице в области малых волновых чисел по сравнению с переходной функцией явных фильтров и обладает легко регулируемой в широких пределах длиной отрезка альтернанса в области больших волновых чисел. При этом фильтрация может производиться по собственным функциям некоторого
Заметим, что вид функции Fn(t) не зависит от типа оператора L и его собственных элементов ірі, но от них зависит разложение (3.2) и, следовательно, uF. Может возникнуть проблема с постановкой ненулевых граничных условий в уравнениях (3.3,3.4). В этом случае запишем уравнение для поправки Ф = Du — DuF: (Pn(L) + 6Ьп)Ф = 6LnDu (3.7) с однородными краевыми условиями. Например, можно положить в некоторых дополнительных для Ф граничных узлах сетки Ф = 0, ЬкФ = 0, к = 1, ..,п — 1.
Преимуществом приведённого алгоритма фильтрации является осуществимость за 0(т) действий вместо не менее 0(т log т), требуемых при использовании быстрого преобразования Фурье. При обработке серии массивов одинаковой длины часть вычислений (матрица и часть про-гоночпых коэффициентов [69] для системы (3.7)) может быть оформлена в виде отдельной подпрограммы, вызываемой один раз. Такая реализация оказывается полезной, например, при построении многомерных фильтров в виде комбинации одномерных. Переходная функция двумерного фильтра изображена на Рис. 3.2а.
С помощью фильтра (3.3) был обработан ряд тестовых задач и экспериментальных данных. Например, на Рис. 3.26 приведены графики исходной функции (суммы функции-ступеньки и двух синусоид) и результата её обработки с помощью фильтра порядка п = 3 с v, вычисленным по формуле (3.8).
Следуя этому критерию, сравним неявный чебышевский фильтр с фильтром Поттера [68] при разной ширине полосы пропускания. Из приведённых в [33] графиков следует, что чебышевский фильтр более оптимален, чем фильтр с окном Поттера во всех случаях. Так, на рис. (3.3) при фильтрации на 80 с.ш. отношение оснований трапеции, приближающей график фильтра Поттера равно 1.4, а для неявного чебышевского фильтра это отношение равно 1.6. На рис. (3.4) при фильтрации на 89 с.ш. эти же отношения равны соответственно 5.0 и 2.2.
Разработаны неявные чебышевские операторные фильтры, позволяющие осуществлять сглаживание данных, ориентированное на разложение по собственным функциям характерного для задачи оператора. Проведены тестовые расчёты и показано, что чебышевские фильтры лучше, чем фильтры с окном Поттера приближают оптимальный фильтр Валле-Пуссена.
Проведено исследование серии проблем, возникающих при применении программы DUMKA и чебышевских методов в вихреразрешающей модели, использующей вихревую форму записи адвективных слагаемых. Предложенные алгоритмы для вихреразрешающих моделей реализованы в программах на языке Фортран для многопроцессорных вычислительных систем. Построены неявные чебышеские фильтры, нашедшие применение в модели общей циркуляции океана.
Экспериментально и теоретически исследована эффективность устойчивых явных разностных схем программы DUMKA с переменными шагами по времени, ориентированных на расположенный вблизи мнимой оси спектр, в задаче вихреразрешающего моделирования и в сравнении со схемой Адамса-Бэшфорта.
Рассмотрен и реализован новый комбинированный итерационный метод для получения поправки к обобщённому давлению, учитывающий характер изменения решений на последовательности временных шагов схем DUMKA-7, DUMKA-8, а также сочетающий обобщённый метод Ритца с характерными для задачи базисными функциями и блочный чебышевский вариационный метод. Предложена новая модификация алгоритма вычисления параметров в устойчивом блочном чебышевском бесконечно продолжаемом итерационном методе с длиной цикла N = N(r:p) — ЪГ2Р.
Предложена и реализована трёхмерная разностная аппроксимация нелинейного члена уравнений Навье-Стокса на октаэдральных сетках в вихревой форме с заменой классического давления на обобщённое давление Бернулли, имеющая дивергентный вид и ортогональная с точностью не ниже 0(h2) сеточному вектору скорости. Создай и реализован в программе алгоритм неявной чебышевской фильтрации, в том числе многомерной. Проведённый анализ показал оптимальность построенных фильтров по сравнению с фильтрами Поттера.
