Разработка универсальной техники реализации неоднородных граничных условий Дирихле и Неймана в методе спектральных элементов Текст научной статьи по специальности «Математика»

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2008

Математика и механика

№ 2(3)

УДК 532+681.3

А.М. Бубенчиков, В.С. Попонин, В.Н. Мельникова

РАЗРАБОТКА УНИВЕРСАЛЬНОЙ ТЕХНИКИ РЕАЛИЗАЦИИ НЕОДНОРОДНЫХ ГРАНИЧНЫХ УСЛОВИЙ ДИРИХЛЕ И НЕЙМАНА В МЕТОДЕ СПЕКТРАЛЬНЫХ ЭЛЕМЕНТОВ1

Разработан обобщенный метод спектральных элементов, использующий универсальную технику реализации неоднородных граничных условий Дирихле и Неймана, позволяющий повысить точность и качество решений плоских нелинейных задач динамики вязкой жидкости по сравнению с аналогами. На основе обобщенного спектрального метода разработан алгоритм решения плоских нелинейных задач динамики вязкой жидкости. Описан математический аппарат, позволяющий получать решения высокого порядка точности в областях сложной геометрии на грубых неструктурированных сетках. Предложен способ решения систем линейных алгебраических уравнений, получающихся при дискретизации уравнений Навье - Стокса обобщенным методом спектральных элементов, позволяющий существенно сократить время расчета за счет подбора предобуславливающей матрицы. В работе представлены тестовые расчеты для течений, имеющих аналитическое решение, также представлены расчеты течений в двумерной прямоугольной полости, произведен расчет стационарного потока за уступом

Ключевые слова: метод спектральных элементов, уравнения Навье - Стокса, неструктрурированные сетки, вязкая жидкость.

Как и в методе Галеркина, в спектральном методе приближенное решение исходного нелинейного дифференциального оператора

разыскивается в виде комбинации по линейно независимой системе элементов {ф;}, называемой базисом в пространстве искомых функций или координатной системой:

В качестве элементов ф; используются ортогональные функции, являющиеся собственными функциями задачи Штурма - Лиувилля на единичном отрезке. Числовые коэффициенты с0..., ем находятся из системы алгебраических уравнений:

Здесь {уу} - базис проекционной системы. В принципе, базисы {ф,} и {уу}, имеющие одинаковую размерность, различны, однако ниже мы принимаем

1. Спектральный метод

G (u) = f .

(1)

(2)

i=0

Фi =Vi (і = 0-N).

1 Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант РФФИ № 08-01-00484-а)

Итак, упомянутая выше задача Штурма - Лиувилля формулируется следующим образом:

—— | а (х)dф£- | + Ь (х)фг- = w(х)фг-, а > 0, Ь > 0,

(ix ^ dx )

Фг. (-1) = Фг (1) = 0. (4)

В общем случае решением задачи (4) являются полиномы Якоби. Так как полиномы Якоби взаимно ортогональны на интервале [-1,1], можно доказать, что

Ум € и: ||м - Ржйм|| ^ 0 при N ^ да . (5)

, N

здесь Рк и = 2 СФг- .

г=0

Более того, если и е Нт (О), иными словами, если искомое решение и является т раз непрерывно дифференцируемым, то, согласно [1], ошибка аппроксимации будет следующей:

IIи - руЙи|| - СуЫ-т ||и|| нт . (6)

II 1И2

Таким образом, используя спектральное разложение, для достаточно гладких функций можно получить экспоненциальную скорость сходимости приближенного решения к точному. В этом и состоит основное преимущество спектрального метода: очень точные приближенные решения могут быть получены при небольшом числе слагаемых, входящих в Ржйи, причём ошибка аппроксимации будет

уменьшаться экспоненциально с ростом N. Таким образом, определив коэффициенты приближенного разложения с,, можно получить значения искомой функции в любой точке области с заданным порядком точности. Указанная схема характерна для классического глобального спектрального метода.

Недостатком глобального спектрального метода является то, что многочлены Якоби являются ортогональными на отрезке [-1, 1] и, следовательно, утверждения об экспоненциальной скорости сходимости приближенного решения к точному имеют место только в случае, если область интегрирования представляет собой отрезок ^ = [-1,1]. Для того чтобы решить задачу с произвольной областью интегрирования, необходимо найти замену координат, переводящую исходную область интегрирования в единичный отрезок.

Спектральный метод обобщается на случай двух и более измерений путём использования в качестве базисных функций тензорного произведения соответствующих одномерных базисных функций.

В случае двух и более измерений задача о нахождении преобразования координат становится достаточно сложной, особенно если область интегрирования имеет сложную форму. Поэтому имеет смысл разбить исходную расчетную область на конечные элементы и искать локальные представления решения через специальные функции, определенные на этих элементах. Однако при дискретизации, полученной на основе локального спектрального метода, матрица системы становится плохо обусловленной, что приводит к медленной сходимости итерационных методов. Эта проблема, как и в случае обычных локальных аппроксимаций, решается с использованием методов на основе пространств Крылова и подбором предобуславливателей.

В настоящей работе в качестве базисных функций мы использовали интерполяционные многочлены, представляющие собой комбинации полиномов Лежандра и их производных и получившие название полиномов Г аусса - Лежандра - Ло-батто [2]:

N

и (х) = 2 щСг (х) , (7)

¡=0

— 1 (1 — X 2 ) -^дг

С(х} = N(+ 1)Ьд (х-) х — X- , 7)

С (X') = 8' ,

здесь 5 у - символ Кронекера, х,- - точки Лежандра - Гаусса - Лобатто, определяемые формулой

хо =-1,

Ху - нули ,1 < у < N -1,

хм = 1. (9)

Веса, необходимые для численного интегрирования при использовании интерполяционных функций (7), (8), определятся соотношением

I П ! \2, } = О’1’-’N. (Ю)

N ( ^ (ху )2

В то же время сами полиномы Лежандра определяются следующим рекуррентным соотношением:

А) (х) = 1,

А( X ) = X ,

2п +1 п

4+1(х) = —- х1п (х)-------7 4-1(х )>п = !>2-^. (11)

п +1 п + 1

Полиномы (11) получаются из решения задачи (4) при

1

а(х) = (1 - х2 )2 ,

Ь (х) = 0, w (х) = 1.

2. Несжимаемые вязкие течения. Математическая постановка задачи

Уравнения, описывающие двумерные стационарные несжимаемые ламинарные течения, имеют вид

+ У(ит и)—— Ли = -Ур, д1У’Яе (12)

У-и = 0.

Уравнения (12) можно переписать в координатном виде следующим образом:

ды, ды ,• Щ 1 д2ы, др . , „

—- +— ----------т = ——,1 = 1,2,

д? дх .• Ле дх2 дх,

' д ' (13) ды ,•

= 0.

дху

Здесь у - индекс суммирования, Ле = - число Рейнольдса, и, Ь - характерная

V

скорость и характерный линейный размер соответственно, V - кинематическая вязкость, и = (м1,и2) - векторная функция, представляющая скорость жидкости в плоском сечении, р - скалярная функция давления жидкости.

3. Метод проекций

Для численного решения уравнений (12) применили метод установления совместно с широко известным методом проекций [3]. Суть этого метода состоит в следующем:

1. На первом этапе находим промежуточное значение для скорости из уравнения

и — П 1

—— -------Ли = (ип -V)и”. (14)

Лг Яе у 1

2. На втором этапе производим расчет давления по формуле

V2 рп+1 =—V«. (15)

М

3. На третьем этапе производим расчет скорости для временного слоя п +1 следующим образом:

П+1 П

----— = -Урп+1 . (16)

м

Шаги 1 - 3 выполняются до установления решения.

4. Дискретизация уравнений Навье - Стокса методом спектральных элементов

Пусть х = (х1, х2 ,•••, хт )е Кт и Ос Жт . Рассмотрим проблему построения решения уравнений (14) в слабой постановке:

|1у2иvdV = |рсСУ . (17)

о ^ ^е о о

Здесь / = ипУи” , и - одна из компонент скорости, V е I2 (О) - некоторая произвольная функция. Ко второму слагаемому из левой части уравнения (17) применим формулу Грина:

-— ГУ2тс}¥ = — [ Г4 V• —¿Б I. (18)

Яе ^ Яе ^ £ дп ) У '

Здесь п = (щ,п2,•••,пт) - нормаль к поверхности дО. В результате наша задача

(14) переформулирована для пространства функций Соболева:

и, V € Н (О) = |«(X) | И5П = 0, ^ € 1} (П)| . (19)

Формулу (17) с учётом (18) можно преобразовать следующим образом:

Г й-й_ vdV + — Г VйVvdV = Г + — 4 V ■ — ¿Б. (20)

о & Ке о о Ке до д

Рассмотрим различные типы областей О с Шт . Наипростейший случай - это кубическая область О = (-1,1 )т с Кт . Поскольку мы хотим получить метод, позволяющий аппроксимировать решение с высоким порядком точности, то и интегралы, входящие в формулу (20), необходимо вычислить с высоким порядком точности. Для этого используем квадратурные формулы Гаусса. В нашем случае интерполяционная формула может быть рассмотрена в виде прямого произведения одномерных базисных функций:

N т

й (X) = I йР1Рт П СР1 (х1). (21)

Р =1,-,Рт =1 1=1

На элементе рассмотрим тестовые функции вида

т

%,...,дт (х) = п с„ (х). ()

г=1

Для простоты, предположим, что мы имеем граничные условия Дирихле. В таком случае нет необходимости производить расчёт интеграла 4V■ —dS. Кон-

50 д

кретную реализацию граничных условий мы рассмотрим более подробно ниже.

Подставим (21) и (22) в уравнение (20). Рассмотрим второе слагаемое, стоящее в левой части уравнения (20). Тогда, с учётом того, что область интегрирования & = [—1,1]т с , получим

— Г =

Яе о

1 11 т N д т д т

= — I • ” I Е Е и — П с (х1)— П Сч. (х1 Ук1 ■■■ •

Яе-1 -1 к=1 р =1 ,..,рт=1 Л т дхк;=1 ' дхк1=1 '

Далее, выражение (23) может быть переписано как

(23)

— Г =

о

1 т N 11 д т д т

— 2 2 «р^ I -1т-П ср, (х;)-—псч, (х ))х1 -ск„

к=1 р =1,...,рт =1 -1 -1 дхк ¡=1 дхк ¡=1

В конечном итоге получаем

(24)

1 1 д _т д _т_

I •" I т—П Ср. (х)—П С (х) dx1 ■■■dxm =

-1 -1 дхк I=1 1 дхк 1=1 1

1 дСРк (хк ) дСдк (хк ), т 1 -1 дхк дхк 1=1,1 * к-1

= I—РРк----------Т-dxk П | Ср. (хг )С?; (хг Мхг •

Воспользуемся квадратурными формулами Гаусса:

1 N

ІСр (хі ) Сч. (х) (ІХі = 2 ю/; Ср, (х/ )С (х/). (26)

-і і і 1=1

Используя ортогональность интерполяционных функций и формулу (26), выражение (24) перепишем следующим образом:

1 1 1 т дй ду , ,

— \ - 1 -т-^Х -¿хм =

Яг . к=і дХк дХк

1 $ - гЛ дс,к(хк )дс?к (Хк), п 8 (27)

— У и р р У І-------------------(Іхк П ю р 8 „а.

П„ р1 >->Рт -1 Я„ Ят к р РіР

Яг р =1,...,Рт=1 к=1 _і дхк дхк і=1,і^к

і дСр (хк) дС_ (хк)

Проинтегрируем выражение Г------------ ---------- -----dxk с использованием квад-

-і дхк дХк

ратурных интегральных формул Гаусса:

1дс„ (х,) гс (х,) , х дсп (х,) гс„ (х,)

I ах, = X ю;----------------------------

-1 дх, дх, '=1 ' дх, дх,

Здесь Юу - веса квадратур, хк - нули исходных интерполяционных функций. В

настоящей работе мы использовали интерполяционные функции и соответствующие им квадратурные веса и нули, определяемые формулами (7) - (11).

Распишем первое слагаемое, входящее в левую часть формулы (20):

Г... Г й••• dxm =

-1 -1 (29)

I N т 1 N т

=----- У й Р п П Ю п 8 п п-------У М Р п П Ю п 8 п п ■

А, ^ Р1 >■■■> Рт^ Р/ Р/ >п/ Л > Р\ Рт ^ Р1 Р/ >п/

& Р! =1 ,.■■,Рт =1 1=1 & Р = !■■■,Рт = 1 /=1 1

Аналогично распишем правую часть формулы (20), используя выражения (21), (22), а также свойство ортогональности базисных функций и квадратурные формулы Гаусса:

II 11 т

|••• | у (х)• V(х1 •••dxm = |••• | у (х)Пс,. (х)^1 •••dxm =

-1 -1 -1 -1 1=1

= £••• £ У(X>х22хтт)Йс,. (Х')= (30)

^ =1 /„ =1 4 ' !=1 4 '

N N , , , т

= £ ••• £ У Х >Х22 >•> хтт ) ,1 •

7 1 7 1 ^ ' .-_1

Здесь Ю/. - веса, х1 - нули полинома в г-м измерении.

Используя формулы (25), (28) - (30), а также меняя индексы д1 = 1..^,..., цт = \ ..Ы, получим матрицу системы и вектор правой части.

Остановимся подробнее на методе аппроксимации граничных условий. В настоящей работе мы воспользовались технологией метода конечных элементов, применяемого для решения задач с ненулевыми граничными условиями. Мы полагали неизвестными все точки расчётной области, включая граничные. В случае

граничных условий Дирихле нет надобности производить расчёт интеграла 4 V •—dS, поскольку значения искомой функции определены в точках границы.

50

Фиксируя индекс д,, можно получить строку коэффициентов спектрального разложения для нахождения узловых значений зависимой функции, в том числе и в точках границы. В случае граничных условий Дирихле для граничной точки, соответствующей индексу д, получаем строку матрицы, состоящую из одного ненулевого элемента, равного единице и расположенного на диагонали матрицы. В правую часть в таком случае мы записываем значения функции на границе. Для внутренней точки, соответствующей индексу дг, мы получаем строку коэффициентов, где число ненулевых элементов в строке равно Ыт, где N - степень полинома, т - размерность задачи. Коэффициенты разложения определяются формулой (27). Если же какая-либо точка и попадает на границу, то, в случае гранич-

ных условий Дирихле, мы переносим значение функции в этой точке с соответствующим спектральным коэффициентом в правую часть. Таким образом, мы достигаем согласованности интерполяционных функций и граничных условий. В случае граничных условий Неймана, для каждой граничной точки, соответствующей индексу д;, мы получаем строку матрицы, состоящую из спектральных коэффициентов, определяемых формулой (27), к этим коэффициентам мы также должны добавить вклады от расчёта интеграла 4 V • —dS. Так как производная —

50 д

является заданной, то произвести расчёт вышеупомянутого интеграла не представляет трудностей. Далее, для каждой внутренней точки, соответствующей фиксированному индексу, мы вновь получаем строку матрицы, состоящую из спектральных коэффициентов, определяемых формулой (27). Таким образом, мы опять получаем согласованность интерполяционных функций и граничных условий. Такой подход позволяет получать решения высокого порядка точности для граничных условий любого типа, кроме того, обработка граничных условий становится универсальной и не составляет никаких трудностей.

Для аппроксимаций уравнений (15), (16) используем технологию, аналогичную описанной для аппроксимации уравнения (14).

Необходимо отметить тот факт, что интерполяционные функции являются ортогональными на интервале 0 = [—1,1]т с Жт . Следовательно, если область интегрирования отлична от единичного куба, необходимо найти преобразование координат, преобразующее исходную область к единичному кубу. Поскольку в качестве сеточных элементов мы имеем четырехугольники, то преобразовать сеточный элемент к единичному квадрату не составляет сложности.

5. Результаты расчетов

Сопоставление результатов с известным аналитическим решением

Для того чтобы показать спектральную скорость сходимости численного алгоритма, рассмотрим течение Коважного [4]. Течение Коважного является точным решением плоских стационарных уравнений вязкой жидкости, т.е. плоских уравнений Навье - Стокса и имеет следующий вид:

и (х, у) = 1 - є( ) 008 (2лу),

X _

V (х, у ) =---Є 8ІП (2пу),

(31)

Параметр X, в свою очередь, определяется следующим соотношением:

2

Численное решение было найдено в области прямоугольной формы [-0,5; 1]х[-0,5; 1,5]. Граничные условия задавали согласно аналитическим формулам (31). Расчеты велись для числа Рейнольдса, равного 40, на сетке, состоящей из 81 конечного элемента. Относительные погрешности для компонент скоростей, давления, а также невязка уравнения неразрывности показана на рис. 1. Из рис. 1 видно, что скорость сходимости приближенного решения к точному носит экспоненциальный характер.

Рассмотрим двумерную полость, представляющую собой область квадратной формы с длиной грани I = 1. Нижняя и боковые грани являются твердыми стенками, верхняя грань - подвижной стенкой, перемещающейся с постоянной скоростью. Граничные условия для данной задачи задавались следующим образом:

дп

Число Рейнольдса связано с физическими параметрами задачи следующей формулой: Яе = Ш / V, где и - размерная скорость перемещения подвижной границы.

ц_1

и_2

Р

ёгу(и)

0

2

4 6 8

Єтепень полинома

10

12

Рис. 1. Относительные погрешности расчетов для компонент скоростей, давления и невязка для уравнения неразрывности

Расчет течения в плоской каверне

и = 0 на твердых неподвижных стенках, и1 = 1, и2 = 0 на подвижной стенке, дР п

— = 0 на всех границах.

В табл. 1 приведены результаты расчетов минимального значения функции тока для течения при числе Рейнольдса 400. Сравнение с данными других авторов показывает достоверность данного подхода к аппроксимации уравнений Навье - Стокса.

Таблица 1

Сравнение минимального значения функции тока в двумерной полости при Ее = 400

Эксперимент Минимум функции тока

225 конечных элементов, степень полинома 3 -0,111

Данные работы [5] -0,11068

Данные работы [6] -0,114

Данные работы [7] -0,108

На рис. 2 приведены линии тока для чисел Рейнольдса 400 и 1000. Результаты для числа Рейнольдса, равного 1000, можно сравнить, например, с данными работы [5]. На рис. 3 приведены линии тока для числа Рейнольдса, равного 1800.

0,1 0,3 0,5 0,7 0,9 -0,4 -0,2 о 0,2 0,4

Ле=400 Ле=1000

Рис. 2. Линии тока в каверне при Ке = 400 и Ке = 1000

-0,4 -0,2 0 0,2 0,4

Рис. 3. Линии тока в каверне при Ке = 1800

Стационарное течение за уступом

Рассмотрим область, представляющую собой канал с уступом. Пусть Ь - длина канала за уступом, С - диаметр входного сечения, а также высота уступа, Б -диаметр канала за уступом (рис. 4).

Рис. 4. Область расчета стационарного потока за уступом

Граничные условия задаются следующим образом. На входе течение считается установившимся, поэтому первая компонента скорости имеет параболический профиль, а вторая компонента скорости тождественно равна нулю. Проинтегрировав уравнение Пуассона для продольной компоненты скорости по диаметру входного сечения, в предположении, что отсутствуют продольные изменения рассматриваемой величины, получим следующее граничное условие для первой компоненты скорости:

дР

«і (О, У) = Яе—

дх

( 2

У

/і + /2 /і /9

-----------1-------2 У + —

2 2 2

Л

(32)

дР

где — - заданный и постоянный градиент давления, Яе = и 2^ / V - число Рей-

дх

нольдса, и - средняя скорость на входе, 11,12 - пределы интегрирования входного сечения в направлении у. Граничное условие на входе для второй компоненты скорости задается следующим образом:

и2 (0 У) = 0-

Граничные условия для давления на входе в этом случае

дР

дх

= -а,

(33)

(34)

где а - некоторая положительная константа.

На выходе течение также считается установившимся и граничные условия

— (X, у ) = 0,

дп УУ (35)

р (х у ) = °.

На твердых стенках для скорости ставится условия прилипания, а также однородные условия Неймана для давления:

и (х, у) = 0;

^ = 0.

дп

(36)

Часть расчетной области с сеткой, состоящей из четырехугольников, представлена на рис. 5. Отметим, что линейный размер Ь составлял 20^.

11

0,8

0,6

0,4 -

0,2 -

0,2

0,4

0,8

1,4

1,6

Рис. 5. Область расчета стационарного потока за уступом, покрытая расчетной сеткой

В работе были проведены расчеты для чисел Рейнольдса 100, 200, 300, 400. Степень полинома равнялась 3, а число конечных элементов не превышало 3 000. За уступом формируется глобальный вихрь, продольный размер которого I увеличивается с ростом числа Рейнольдса. В табл. 2 представлена зависимость длины большого вихря за уступом от числа Рейнольдса. Здесь следует отметить, что экспериментальные данные для двумерных течений отсутствуют при числах Рейнольдса больше 200.

Таблица 2

Зависимость длины большого вихря за уступом от числа Яе

0

Эксперимент (1 / <!) Число Яе

100 200 300 400

Данные нашей работы 5,0 7,9 9,85 10,9

Данные работы [5] 5,1 7,2 9,5 10,49

Данные работы [7] - - 9,08 10,04

Экспериментальные данные работы [8] 5,0 8,5 - -

Данные работы [9] 5,0 8,2 10,1 14,8

Заключение

Таким образом, в работе был предложен эффективный алгоритм расчета уравнений Навье - Стокса, использующий неструктурированные сетки, преобразования координат и спектральные разложения искомых функций. Также был предложен универсальный алгоритм реализации неоднородных граничных условий Дирихле и Неймана, позволяющий без потери спектральной точности применять данный метод к широкому спектру задач. Решения, полученные с помощью метода спектральных элементов, могут иметь самостоятельные значения, если нас интересуют тонкие эффекты, сопровождающие исследуемый процесс, или же они могут применяться как тестовые результаты при построении алгоритмов метода конечных разностей.

ЛИТЕРАТУРА

1. Van de Vosse F.N. Spectral Element Methods: Theory and Application, 1999.

2. Boyd John P. Chebyshev and Fourier Spectral Methods. Second Edition. University of Michigan, 2000.

3. Флетчер К. Вычислительные методы в динамике жидкостей. М.: Мир, 1991.

4. HelenbrookB.T. A Two-Fluid Spectral Element Method. Department of Mechanical and aeronautical engineering, 1999.

5. Бубенчиков А.М., Фирсов Д.К., Котовщикова М.А. Численное решение плоских задач динамики вязкой жидкости методом контрольных объемов на треугольных сетках // Математическое моделирование. 2007. Т. 19. № 4.

6. Ozawa S. Numerical studies of steady flow in a two-dimensional square cavity at high Reynolds numbers // J. Phys. Soc. Jpn. 1975. V 38. P. 889.

7. Wan D.C., Patnaik B.S.V., and Wei G.W. Discrete Singular Convolution - Finite Subdomain Method for the Solution of Incompressible Viscous Flows // J. Computat. Phys. 2002. V. 180. P. 229 - 255.

8. Armaly B.F., Durst F., Pereira J.C.F., and Schonung B. Experimental and theoretical investigation of backward-facing step flow // J. Fluid Mech. 1983. V. 127. No. 473.

9. Елизарова Т.Г., Серегин В.В. Аппроксимация квазидинамических уравнений на треугольных сетках // Вестник Моск. ун-та. Серия 3. Физика. Астрономия. 2005. № 4. С. 15

- 18.

Статья принята в печать 19.07.2008 г.