Математическая постановка и решение пространственных краевых задач методом спектральных элементов

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

2008 Математика и механика № 3

УДК 532+681.3

А.М. Бубенчиков, В.С. Попонин, В.Н. Мельникова МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА И РЕШЕНИЕ ПРОСТРАНСТВЕННЫХ КРАЕВЫХ ЗАДАЧ МЕТОДОМ СПЕКТРАЛЬНЫХ ЭЛЕМЕНТОВ1

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

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

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

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

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

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

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

С

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

г. К

Р = IСФ; .

/=0

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

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

Фг = V; (г ).

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

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

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

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

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

и

/=0

—1

С

С (X') = 5„ ,

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

хо =-1 , х] - нули , 1 < у < N -1 ,

хк = 1.

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

! \2 ,•/ = О’1’-’N.

^ ^ ^ (х7. )2

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

А) ( Х ) = 1 ,

А( X ) = X ,

А+1

п +1 п + 1

2. Уравнение Пуассона

Рассмотрим уравнение Пуассона:

-У2и = /(х), х е О,

«50 =Ф (х ).

Пусть X = (х1;х2,х3)е^3 и . Рассмотрим проблему построения решения

(9) в слабой постановке:

-1 У2им^ = | /У

п п

К левой части уравнения (10) применим формулу Грина:

ди

-1Vгиус1¥ = -1УиУус1¥- 4 V = |№У .

дО дП О

Здесь п = (щ, п2, п3) - нормаль к поверхности дО. В результате наша задача (9) переформулирована для пространства функций Соболева [3, 4]:

и, V €

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

ди

-1ЧиЧъИУ = |/ус1У + 4 .

О О дО дп

Рассмотрим кубическую область 0 = (-1,1 )3 с ^3. В этом случае интерполяционная формула может быть рассмотрена в виде прямого произведения одномерных базисных функций:

N 3

и

г1=1,г'2=1,гз =1 г=1

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

3

^ЪЧ2т

г=1 '

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

4 V

50

С учетом (14), (15), левую часть формулы (13) можно переписать следующим образом:

/ ди 5у ди 5у ди Зу! , , ,

11--------1-------1-----I ахауаг =

0 ^йх дх ду ду & & )

111 ( д ( N N N

Л

= Ш ТТ Т С( х)Сч (У)Сз (*) Т- С (Х)СР2 (У)Ср 3 (*)) +

I к=1 к =1 і =1 )СХ

-1-1-1

( N N N

ду

Т Т Т СК (Х)С2 (у)Сз (2) I —(СР1 (Х)СР2(у)Ср3 (2)) +

^ ?! — 1 ?2 — 1?3 — 1

д( II I С (х)ЄІ2 (у)Єіз (2)^((Х)СР2(у)Ср3 (2))Л

д2

У Іі — 1 і'2 — 1 І3 — 1

N N N

<ІХ(Іу<І2 =

= I II ї І І {-д( (х)С!2 (У)СН (7))(Ср1 (х)Ср2(У)СР3 (2))-

І-1=112=113 =1 -і-і-ічдх дх

+ду((х)С,2 (у)С(1 (г))(С„ (()С,2(у)Срз(г)) +

+д^( (х)С2 ((г))((, ((Ср2 (у)Ср3 (г))'] =

дг

N N N 111 / о д

= I I I Ukhh I | 11 C!2 (y)C!3 (z)C^ (y)C^ (z)—C4(x) — Cpi (x) +

ii=1 i2=1/3 =1 -i-i-iV ox ox

+Ck (x)Cii3 (z)Cp2 (y)Cp3 (z)ddyCi2 (У)ddyCp2 (У) +

C (x)Ci2 (y)Cpi(x)Cp2 (y) d; C (z) -dZ CP3 (z)j dxiydz =

N N N f 111 f d d d

= Ё Ё Ё Uv'2i31 III I Ci2 (y)CP2(y)Ci3 (z)Cp3 (z)T~ Ci (x)T~Cpi(x) mdydz +

4 =1 i 2 = 1 ц =1 V-1 -1 -1V dx dx J

+ I I I^C(x)Cpi(x)C3 (z)Cp3(z)^дУс,(j)-dyCP2(y"\dxdydz +

+ f I I (C(x)Cpi(x)C2 (.у)Ср2 (yд^С!з (z)dZcp3 (z) jdxdJdzj =

N N N ( 1 Я Я 1 1

= HZ 1 I — Ci! (x)—Cpi( x)dx | Ci2 (y)Cp 2 (y)dy JC (z )Cp3dz + i-1 =1i2 = 1 i3 =1 V-1 dX dX -1 -1

1 1 д d 1 + J C (x)cpi(x)dx Idy Ci2 (ydy CP(У I C (z)Cp3 (z)dz +

+ J C (x)Cp1 (x)dx J Ci2 (J)Cp2 Cy)dV J” Ci3 (z)Cp3 (z)dz| •

-1 -1 -1 dz dz J

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

N N N ( N д д N N

Z Z Z UW3 (Z 0,-j Сг1 (xJ )— Cpl (xJ) X % Ci2 (yk )Cp2 (/ )Z ю, Ci3 (z1 )Cp3(z/)+

i1=1i2=1i3 =1 \J=1 OX OX k=1 1=1

+ T C (xJ )Cp1 (xJ) T c21 (Уk ^ Cp 2 (Уk )T щ C3 (zl )Cp33 zl) +

J=1 k=1 dy dy i=1

+ 1 ЮjCh (XJ )Cpl (xJ) jt «kC2 (yk )Cp2 (yk )jt Chl (zl )-^r Cp3 (zl) 1 =

j=1 k=1 i=1 dz dz J

N N N ( N N N д d

= 111 Ui¥3 ( XII «J % Щ—Ск (xJ )—Cpl (xJ )Ch (yk )Cp2 (yk )Ch (z1 )Cp3 (zl)+ /j=i/2=1/3 =1 v j=ik=i/=i dx dx

N N N д d

+Z 11Юj«kЮ/C4 (xJ)Cpl (xJ) — C21 (yk) — Cp2 (/ )Ci3 (z1 )Cp3 (z1) + i=1 k=11=1 ду dy

N N N д d ^

+1 IIЮjЮ,C4(xJ)Cpl(xJ)C!2(yk)Cp2(yk)-jC!31 (z1)—CP3(z1) 1 =

j=1k=11=1 dz dz У

МММ МММ /О о

= 111 ІII«у«кю,I — Сг1 (хУ)—Єр1 (хУ)СІ2 (/)Ср2(/)С (г1 )Срз(г1)+

г1 =1г'2 =1г3 =1 У=1к=11=1 \ОХ ОХ

+Сг1 (х7)Ср1 (х7)дд_с2 (/)-д- ср2 (/ )С3 (г1 )Ср3 (г1) +

+ С (хУ )Ср1 (хУ )С2 (ук )Ср2 (ук) А С3 (71) Срз (г1) 1 =

& 02 )

N N N ( N д д

= III ик¥з І I «/—уСг1 (хуСр1 (ху)5г-2,р2Юг-2З/,,рз3 + і1=і /2 =і і, =1 \у=1 дх дх

N д к д к +Е сі2(у )ттср2(у )8г1,рі%8і3,рзі3 +

к=і 5ук 2 5ук 1 3

+ Тсі3

1=1 дг дг )

Аналогичным образом преобразуется и правая часть (13). Таким образом, используя формулы (16), а также меняя индексы рх,р2,р3, получим матрицу системы и вектор правой части.

3. Преобразования координат

Поскольку многочлены (5) являются ортогональными лишь на единичном кубе, то для областей более сложной геометрии необходимо осуществить замену координат, переводящую исходную область интегрирования к единичному кубу. Пусть х = х(£) - некоторое преобразование координат, где ^ є [-1,1] для всех

і = 1,2,3 . Тогда

ди ди д£ ди дп ди дх ;

дх д£ дх дп дх дх дх ’

ди ди дЪ ди дп ди дх

— =-------- +----------------------------- +- ’

ду дЪ, ду дц ду дх ду

ди _ ди д£ ди дп ди дх

дг д£ дг дп дг дх дг

Формулы (17) можно переписать в следующем виде:

ди ди 1 ( ду дг дг ду Л + ди ( 1 Л ( ду дг дг ду Л + ди 1 ( ду дг ду дг

дх дЪ У ^дп дх дп дх) дпч УДд^дх д^дх) дх У^д^дп дп дЪ,, ди (1V , _т , ч ди (_ 1V ч ди (1V , _, .ч

_д^1 У У( 3 У32У 23 ;+дп1 У У ( 3 У3у 23 ^ + дх^ УУ( 32 У гг'131

ди ди ( 1 ^( дх дг дг дх \ди ^ 1 Vдх дг дг дх ^+дм ( 1 Л( дх дг дх дг

ду д^ч J)\дцдх дцдх) д^чУ)\д^дх д^дх) д%\ J)удхдц дцд£ ди ( 1V ч ди ( 1V ч ди ( 1V ч

_"д£ч У Г2 У 33 У'3 Уз2 ^^"дпІУ^^Уп У33 У1Ъ Узі ^+"дхІ У^^У'3 Уз2 УігУз1

du du f 1 дх dy dy dxdu f 1 дх dy dy dxdu f 1 Y дх dy dx dy

dz 5^v J дцд%) 5^1 J Ад^дх д^д%) 5x1 J Ад^дп дцд^)

du f 1V . 5u f_ 1V 5u f 1V _ .

_д^1 J J( 23 Jl3 J 22 ^+дп1 J J( 23 J 21 ^+дх1 J J( 22 Jl3 J 23

Здесь J - определить матрицы Якоби, Jy - компоненты матрицы Якоби.

Заключение

Таким образом, в работе приведена математическая постановка численного

решения пространственных краевых задач методом спектральных элементов.

Приведен вывод преобразований координат, необходимых для решения задач в

областях сложной геометрии.

ЛИТЕРАТУРА

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. Helenbrook B.T. A Two-Fluid Spectral Element Method. Department of Mechanical and Aeronautical Engineering, 1999.

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