CC BY
107
УДК 519.633.6:519.612.2
МЕТОД СПЕКТРАЛЬНЫХ ЭЛЕМЕНТОВ ДЛЯ РЕШЕНИЯ ПЛОСКИХ ЗАДАЧ ДИНАМИКИ ВЯЗКОЙ ЖИДКОСТИ НА НЕРАЗНЕСЕННЫХ НЕСТРУКТУРИРОВАННЫХ СЕТКАХ
В.С. Попонин*, А.В. Кошеутов, В.П. Григорьев, В.Н. Мельникова*
Томский политехнический университет *Томский государственный университет E-mail: [email protected]
Описан алгоритм построения решения плоских задач динамики вязкой несжимаемой жидкости методом спектральных элементов. Алгоритм позволяет получать решения высокого порядка точности на грубых неструктурированных сетках. Разработан эффективный алгоритм расчета нелинейных уравнений в частных производных с различными типами ненулевых граничных условий методом спектральных элементов. Показана высокая эффективность применения метода обобщенных невязок совместно с неполным LU разложением для расчета систем линейных алгебраических уравнений, получающихся при дискретизации уравнений Навье~Стокса.
Ключевые слова:
Уравнения Навье-Стокса, спектральный метод, метод проекций, обобщенный метод минимальных невязок, неполная LU-факторизация.
Key words:
Navier-Stokes equations, spectral method, projection method, generalized minimal residual method, incomplete LU-factorization.
Два фундаментальных качества производимых вычислений составляют доктрину всей вычислительной математики - это точность и быстродействие расчетов [1]. Если быстродействие может быть увеличено, в том числе, за счет технических решений, то точность вычислений повышается, главным образом, за счет применения математических.
В настоящей работе описан математический аппарат, позволяющий получить для плоских задач динамики вязкой жидкости решения высокого порядка точности в областях сложной геометрии при решении их методом спектральных разложений искомых функций на каждом из элементов неструктурированной сетки, накрывающей исходную область. Преимущества такого метода очевидны: на достаточно грубой сетке можно получать решения с заданным порядком точности.
Высокоточные решения могут иметь самостоятельное значение, если интересуют тонкие эффекты, сопровождающие исследуемый процесс или же они могут применяться как тестовые результаты при построении алгоритмов метода конечных разностей.
Принцип, лежащий в основе всех сеточных методов, заключается в сведении исходных дифференциальных уравнений в частных производных к системе алгебраических уравнений, которые могут быть решены известными методами. Однако в спектральных методах [2-4] процедура, реализующая эти принципы, аналогична используемой в аналитических методах решения линейных дифференциальных уравнений в частных производных.
Спектральный метод
Как и в методе Галеркина, в спектральном методе приближенное решение исходного нелинейного дифференциального оператора 0(ы)=/ находится в виде комбинации по линейно независимой системе элементов {щ}, называемой базисом в простран-
стве искомых функций или координатной систе-
N
мой PN = ^ с. В качестве элементов щ использу-
= 0
ются ортогональные функции, являющиеся собственными функциями задачи Штурма-Лиувилля на единичном отрезке [5].
Таким образом, используя спектральное разложение, для достаточно гладких функций можно получить экспоненциальную скорость сходимости приближенного решения к точному. В этом и состоит основное преимущество спектрального метода: очень точные приближенные решения могут быть получены при небольшом числе слагаемых, входящих в Р/ы, причём ошибка аппроксимации будет уменьшаться экспоненциально с ростом N. Таким образом, определив коэффициенты приближенного разложения с, можно получить значения искомой функции в любой точке области с заданным порядком точности.
Известно также, что экспоненциальная скорость сходимости приближенного решения к точному достигается лишь в том случае, когда базисные функции удовлетворяют граничным условиям [2]. В данной работе применили универсальный алгоритм аппроксимации ненулевых граничных условий [6] и показали, что возможно использовать его без потери точности для решения задач динамики вязкой жидкости.
В настоящей работе в качестве базисных функций использовали интерполяционные многочлены, представляющие комбинации полиномов Лежандра и их производных.
Несжимаемые вязкие течения.
Математическая постановка задачи
Уравнения, описывающие двумерные стационарные несжимаемые ламинарные течения, имеют вид [2]:
г 1
V (и и)----------Аи = —Ур
Я-
(1)
Уи = 0
Уравнения (1) можно переписать в координатном виде:
8и]и1
8х]
8и]
дх,
_1_д\ = —др = . 2 Яе дх2 дх ’ ’’
■= 0.
А?
■ = —Урп
(4)
где и - промежуточное решение.
Шаги 1-3 выполняются до установления решения.
Дискретизация уравнений Навье-Стокса
методом спектральных элементов
Пусть х=(хьх2,...,хт)еК
ти ОсКт. Рассмотрим проблему построения решения ур. (2) методом взвешенных невязок [4]:
Ги_и_ —- [V 2тс1У= Г^У.
J А? Я- J J
Яе О Яе чо
1 ¡VИСУ = -1-1 \VvVvdy — ОV 88и-СБ
Яе дО дп
(6)
и,V е Н0(О) = \и(х) | и8О = 0, -дие ¿2(О)1.
I дх )
Формулу (5) с учётом (6) можно преобразовать следующим образом:
Ги_и_ vdУ + — ГVйVvdУ -'А? Я- О
Г У + — 4 V 8и СБ.
(7)
Здесь} - индекс суммирования; Ке= ПЬ/у - число Рейнольдса; П, Ь - характерная скорость и характерный линейный размер; у - кинематическая вязкость; и=(и1,и2) - векторная функция, представляющая скорость жидкости в плоском сечении; р - скалярная функция давления жидкости.
Метод проекций
Для численного решения ур. (1) применяли метод установления совместно с широко известным методом проекций [2]. Суть метода состоит в следующем:
1. Находим промежуточное значение для скорости из уравнения
и — ип 1 , _ , п„. п
-----------Аи = (и V)u . (2)
Д Я-
2. Производим расчет давления по формуле
V2 рп+ = — Vй. (3)
А,?
3. Вычисляем значение скорости для временного слоя и+1:
о А? Я-о
1 г ди
— О V—!
о Я- да дп
Рассмотрим различные типы областей ОсКт. Н-аипростейший случай - это кубическая область О=(-1;1)тсКт. Поскольку необходимо получить метод, позволяющий аппроксимировать решение с высоким порядком точности, то и интегралы, входящие в формулу (7), необходимо вычислить с высоким порядком точности. Для этого используем квадратурные формулы Гаусса. В нашем случае интерполяционная формула может быть рассмотрена в виде прямого произведения одномерных базисных функций:
и( х) = X и‘р1,-,рт ПСр, (х-). (8)
Р1 =1.-.рт =1 ¡=1
На элементе рассмотрим тестовые функции вида:
VÍ1...,Í.(х) = Йс»,(х). (9)
¡=1
Для простоты, предположим, что заданы граничные условия Дирихле. В таком случае нет необходимости производить расчёт интеграла
О V СБ. Конкретную реализацию граничных
да дп
условий рассмотрим более подробно ниже.
Подставим (8) и (9) в ур. (7). Рассмотрим второе слагаемое, стоящее в левой части ур. (7). Тогда, с учётом того, что область интегрирования О=[-1;1]тсКт, получим:
— \vUVvdУ=—\ ■■■[ X X и —х
Я- яе J . р1--р» дх
Я-
1 —1 к=1 р1=1,..., рт =1
8х,г
д
(5) -А.
ХП Ср, (х ) ^П СЧ, (х С1 ■Схт . (10)
¡=1 дхк ¡=1
Далее, выражение (10) может быть переписано как:
1 _ 1 т N 1
=— X X и„,..р, /•••
Здесьf=unVun; и - одна из компонент скорости; уеЬ2(О) - некоторая произвольная функция. Ко второму слагаемому из левой части ур. (5) применим формулу Грина:
Здесь и=(иьи2,...,ит) нормаль к поверхности 8О. В результате задача (2) переформулирована для пространства функций Соболева:
1
Я- О Я- *=1 р1=1,...,рт=1 ,
¡дгП с р. (х) д- Г! с (х к --а,. ш)
—1°Хк ¡=1 °Хк ¡=1
В итоге получаем:
} ■ Г1Ср,( х ) 8т П^ х>Сх1-‘Схт =
— 1 —1иХк ¡=1 °Хк ¡=1
V 8Срк (х.) 8Ск (х,) т 1
1 8хк 8хк ¡=1)Фк — 1
Воспользуемся квадратурными формулами Гаусса:
)с, (х,)С (х,)Сх, = £щеР1 (х1 С(х1). (12)
-1 ‘ ‘ I,=1
Используя ортогональность интерполяционных функций и формулу (12), выражение (11) примет вид:
± [-[£“- .^х,-* =
Яе -1 -1 = дхк дхк
V дС (хк) дС (хк)
дх,,
дх
х П ®рА,я •
/=1,/^к
Проинтегрируем выражение
V дСРк (хк) С (хк) -1 дхк дхк
(13)
Сх,
с использованием квадратурных интегральных формул Гаусса:
V дС (хк) дС (хк)
дх.
дх.
Схк = £
дСРк (х!) —с (4)
здесь а - веса квадратур, хк - нули исходных интерполяционных функций.
Распишем первое слагаемое, входящее в левую часть формулы (7):
1 1
| | АуСх1 — Схт =
At
_ 1 At
-1 -1
N
£ и
/_! Р1,-, Р,
Р1 = 1.-.Рт =1 '=1
£
Р1»-. Р,
Рх =1,..., Рт =1 1 = 1
Па 8 _
1 1 Р, Р, Я,
1=1
П а 8 .
1 1 Р, Р, я,
Аналогичным образом аппроксимируется правая часть формулы (7).
Используя полученные аппроксимации ур. (7) и меняя индексы ^=1,...Д; дт=1,...,И получим матрицу системы и вектор правой части.
Остановимся подробнее на методе аппроксимации граничных условий. В настоящей работе была использована технология метода конечных элементов, применяемого для решения задач с ненулевыми граничными условиями. Полагались неизвестными все точки расчётной области, включая граничные. В случае граничных условий Дирихле нет надобности производить расчёт интеграла
4 V ^ СБ, поскольку значения искомой функции
дп дп
определены в точках границы. Фиксируя индекс #;, можно получить строку коэффициентов спектрального разложения для нахождения узловых значений зависимой функции, в том числе и в точках границы. В случае граничных условий Дирихле для граничной точки, соответствующей индексу #;, получаем строку матрицы, состоящую из одного
ненулевого элемента, равного единице и расположенного на диагонали матрицы. В правую часть записываем значения функции на границе.
Для внутренней точки, соответствующей индексу #;, получаем строку коэффициентов, где число ненулевых элементов в строке равно Ыт, где N - степень полинома; т - размерность задачи. Коэффициенты разложения определяются формулой (13). Если какая-либо точка и{, попадает на границу, то, в случае граничных усл1 овт ий Дирихле, переносим значение функции в этой точке с соответствующим спектральным коэффициентом в правую часть. Таким образом, достигаем согласованности интерполяционных функций и граничных условий.
В случае граничных условий Неймана, для каждой граничной точки, соответствующей индексу #;, получаем строку матрицы, состоящую из спектральных коэффициентов, определяемых формулой (13), к этим коэффициентам необходимо добавить вклады от расчёта интеграла 4 V СБ. Так
дп дп
как производная — является заданной, то произ-д п
вести расчёт вышеупомянутого интеграла не представляет трудностей. Далее, для каждой внутренней точки, соответствующей фиксированному индексу, вновь получаем строку матрицы, состоящую из спектральных коэффициентов, определяемых формулой (13).
Таким образом, получаем согласованность интерполяционных функций и граничных условий. Данный подход позволяет получать решения высокого порядка точности для граничных условий любого типа, кроме того, обработка граничных условий становится универсальной.
Для аппроксимаций ур. (3), (4) используем технологию, аналогичную описанной для аппроксимации ур. (2).
Необходимо отметить тот факт, что интерполяционные фу-нкции являются ортогональными на интервале П=[-1;1]тсКт. Следовательно, если область интегрирования отлична от единичного куба, необходимо найти преобразование координат, преобразующее исходную область к единичному кубу. Поскольку в качестве сеточных элементов имеем четырехугольники, то преобразовать сеточный элемент к единичному квадрату не составляет сложности.
Методы решения систем линейных
алгебраических уравнений
В результате дискретизации ур. (2)-(4) получим системы линейных алгебраических уравнений. Известно, что матрицы в таких системах являются разреженными. Запишем систему линейных алгебраических уравнений в виде Лх=Ь, где Л - разреженная матрица; х - вектор неизвестных; Ь - вектор правой части.
п
Для расчета таких систем применили методы сопряженных градиентов (BICGStab) и обобщенных невязок (GMR.ES) совместно с неполной LU-факторизацией [7]. Для разреженных матриц используется неполное разложение ^и с исключением элементов, принадлежащих определенному множеству. В методе ^и(0) в качестве такого множества берется множество пар индексов, соответствующих нулевым элементам матрицы А. Полученное разложение является предобуславливаю-щей матрицей к исходной матрице системы: и-х1-хАх=и-х1-хЪ.
Данную систему можно быстро решить методом BICGStab [7]. Применение метода BICGStab с пре-добуславливателем, полученным разложением ^и(0), позволяет существенно сократить число итераций в расчете систем линейных уравнений.
Наряду с методом BICGStab для решения систем линейных уравнений применили также метод GMRES с неполной LU-факторизацией [7].
Результаты сравнений эффективности различных итерационных методов приведены ниже.
Результаты расчетов
Сопоставление результатов
с известным аналитическим решением
Для того, чтобы показать спектральную скорость сходимости численного алгоритма, рассмотрим течение Коважного [8]. Численное решение было найдено в области прямоугольной формы [—0,5;1]х[—0,5;1,5]. Граничные условия задавали согласно аналитическим формулам. Расчеты велись для Re=40 на сетке, состоящей из 81 конечного элемента. Невязка уравнения неразрывности и относительные погрешности для компонент скоростей и давления приведены на рис. 1.
В табл. 1 приведены результаты сравнения эффективности применения различных итерационных методов для расчета систем линейных алгебраических уравнений, получающихся при дискретизации ур. (2)-(4). Следует отметить, что временные затраты на одну итерацию методом GMRES превосходили на 30 % временные затраты на одну итерацию ме-
тодом сопряженных градиентов. Одна итерация методом Гаусса-Зейделя примерно в 4 раза быстрее итерации метода сопряженных градиентов. Из таблицы 1 видно, что наиболее эффективен метод GMRES совместно с неполной LU-факторизацией.
Таблица 1. Оценка эффективности применения различных итерационных методов для расчета систем линейных уравнений, получающихся при дискретизации ур. (2) методом спектральных элементов
Степень полинома Число итераций, необходимое для сходимости метода
Гаусса- Зейделя В^1аЬ В^1аЬ+ + !Ш(0) GMRES GMRES+ +11_11(0)
3 650 420 120 250 48
5 2100 1150 350 440 110
7 6700 1620 510 580 156
9 21000 2350 740 890 241
11 37400 4300 1050 1100 357
Расчет течения в плоской каверне
Рассмотрим двумерную полость, представляющую собой область квадратной формы с длиной грани /=1. Нижняя и боковые грани являются твердыми стенками, верхняя грань является подвижной стенкой, перемещающейся с постоянной скоростью. Граничные условия для данной задачи задавались следующим образом: и=0 на твердых неподвижных стенках; «1=1, и2=0 на подвижной стенке; др/дп=0 на всех границах.
Рассмотрим результаты расчетов минимального значения функции тока у/^ для течения при Re=400. В настоящей работе было получено Утт=_0,111; значение работы [9] у/т^-0,11068; для [10] утп=-0,114; для [11] у„,т=-0,108. Сравнение с данными других авторов показывает достоверность данного подхода к аппроксимации уравнений Навье-Стокса.
На рис. 2 приведены линии тока для Re=400 и Re=2000.
На рис. 3 представлен профиль продольной компоненты скорости при Re=100.
Степень полинома
Рис. 1. Относительные погрешности расчетов для компонент скоростей, давления и невязка для уравнения неразрывности
-0,4 -0,2 В 0,2 0,4 -0,4 -0,2 0 0,2 0,4
Рис. 2. Линии тока в каверне при Не=400 (слева) и Не=2000 (справа)
Рис. 3. Двумерное распределение продольной компоненты скорости при Не=Ю0. ♦ - метод спектральных элементов; а - результаты [12]
Стационарное течение за уступом Рассмотрим область, представляющую собой канал с уступом. Пусть Ь и Б - длина и диаметр канала за уступом, й - диаметр входного сечения, а также высота уступа, рис. 4.
Граничные условия задаются следующим образом. На входе течение считается установившимся, поэтому первая компонента скорости имеет параболический профиль, а вторая компонента скорости тождественно равна нулю. Проинтегрировав
Рис. 4. Область расчета стационарного потока за уступом
уравнение Пуассона для продольной компоненты скорости по диаметру входного сечения, в предположении, что отсутствуют продольные изменения рассматриваемой величины, получим граничное условие для первой компоненты скорости:
Таблица 2. Зависимость длины большого вихря I за уступом от числа Не
dP
ul{0,y) = Re —
ox
(.
У_
2
где дР/дх - заданный и постоянный градиент давления; Re=2(Ш/v) - число Рейнольдса; и - средняя скорость течения на входе; /ь/2 - пределы интегрирования входного сечения в направлении у. Граничное условие на входе для второй компоненты скорости: и2(0,у)=0. Граничные условия для давления на входе:
дР
— = -а, дх
где а - некоторая положительная константа.
На выходе течение также считается установившимся, граничные условия задаются по формулам:
ди
— = 0, р(х у) = 0.
дп
На твердых стенках для скорости ставится условие прилипания, а также однородные условия Неймана для давления:
др
и (х, у) = 0; — = 0. дп
В работе были проведены расчеты течения за уступом для различных чисел Рейнольдса, табл. 2. Степень полинома равнялась 3, а число конечных элементов не превышало 3000.
Источник информации Re
100 200 300 400
Настоящая работа 5,0 7,9 9,85 10,9
[9] 5,1 7,2 9,5 10,49
[10] - - 9,08 10,04
[13] 5,0 8,5 - -
[14] 5,0 8,2 10,1 14,8
За уступом формируется глобальный вихрь, продольный размер которого увеличивается с ростом числа Re.
Выводы
Предложен алгоритм, использующий неструктурированные сетки, преобразования координат и спектральные разложения искомых функций, который позволяет получать высокоточные решения нелинейных эллиптических задач с нелинейностью типа конвективного члена, характерной для уравнений Навье-Стокса. Для решения системы линейных алгебраических уравнений рекомендуется использовать метод обобщенных невязок совместно с неполной LU-факторизацией. Метод позволяет достичь наивысшей скорости расчета по сравнению с общепринятым методом сопряженных градиентов и методом Гаусса-Зейделя. Описан также универсальный алгоритм аппроксимации неоднородных граничных условий Дирихле и Неймана методом спектральных элементов, не приводящий к потере точности.
Работа выполнена при финансовой поддержке гранта РФФИ № 08-01-00484-а.
СПИСОК ЛИТЕРАТУРЫ
1. Самарский А.А. Введение в теорию разностных схем. - М.: Наука, 1971. - 553 с.
2. Флетчер К. Вычислительные методы в динамике жидкостей. -М.: Мир, 1991. - 552 с.
3. Boyd J.P. Chebyshev and Fourier Spectral Methods. Second Edition. - N.Y.: Dover Publications, 2001 - 688 p.
4. Van de Vosse F.N., Minev P.D. Spectral Element Methods: theory and application. - Eindhoven: Eindhoven University of Technology, 1999. - 68 p.
5. Patera A.T A spectral element method for fluid dynamics: laminar flow in a channel expansion // Journal of Computational Physics. -1984. - V. 54. - Iss. 3. - P. 468-488.
6. Бубенчиков А.М., Попонин В.С., Фирсов Д.К. Спектральный метод решения плоских краевых задач на неструктурированной сетке // Математическое моделирование. - 2007. - Т. 19. -№ 10. - С. 3-14.
7. Saad Y. Iterative Methods for Sparse Linear Systems. - Philadelphia: Society for Industrial and Applied Mathematics, 2000. -447 p.
8. Helenbrook B.T A two-fluid spectral element method // Computer Methods in Applied Mechanics and Engineering. - 1999. -V. 191. - № 4. - P. 273-294.
9. Бубенчиков А.М., Фирсов Д.К., Котовщикова М.А. Численное решение плоских задач динамики вязкой жидкости методом контрольных объемов на треугольных сетках // Математическое моделирование. - 2007. - Т 19. - № 6. - С. 71-85.
10. Ozawa S. Numerical studies of steady flow in a two-dimensional square cavity at high Reynolds numbers // Journal of the Physical Society of Japan. - 1975. - V. 38. - № 3. - P. 889-895.
11. Wan D.C., Patnaik B.S.V., Wei G.W Discrete Singular Convolution-Finite Subdomain Method for the Solution of Incompressible Viscous Flows // Journal of Computational Physics. - 2002. - V. 180. -Iss. 1.- P. 229-255.
12. Junk M., Rao S. A new discrete velocity method for Navier-Stokes equations // Journal of computational physics. - 1999. - V. 155. -№ 1. - P. 178-198.
13. Armaly B.F., Durst F., Pereira J.C.F., Schonung B. Experimental and theoretical investigation of backward-facing step flow // Journal of Fluid Mechanics. - 1983. - V. 127. - P. 473-496.
14. Елизарова Т.Г., Серегин В.В. Аппроксимация квазидинамиче-ских уравнений на треугольных сетках // Вестник Московского университета. Сер. 3: Физика. Астрономия. - 2005. - № 4. -С. 15-18.
Поступила 13.04.2010 г.
