CC BY
50
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2008 Математика и механика № 3(4)
УДК 532+681.3
А.М. Бубенчиков, В.С. Попонин, В.Н. Мельникова МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА И РЕШЕНИЕ ПРОСТРАНСТВЕННЫХ КРАЕВЫХ ЗАДАЧ МЕТОДОМ СПЕКТРАЛЬНЫХ ЭЛЕМЕНТОВ1
В работе описана математическая постановка и решение пространственных краевых задач методом спектральных элементов. Основной упор сделан на дискретизацию трехмерного уравнения Пуассона. В работе приведен подробный вывод преобразований координат, необходимых для решения задачи в областях сложной геометрии.
Ключевые слова: метод спектральных элементов, уравнение Пуассона, неструктрурированные сетки, вязкая жидкость.
1. Спектральный метод
Принцип, лежащий в основе всех сеточных методов, заключается в сведении исходных дифференциальных уравнений в частных производных к системе алгебраических уравнений, которые могут быть решены известными методами. В глобальном спектральном методе решение ищется путём разложения в ряд по некоторой системе ортогональных функций, называемых базисными. Имея представления искомых функций в спектральном пространстве, то есть определив их в виде разложения по базисным функциям, в глобальном спектральном методе строится система интегральных соотношений, получающихся умножением исходных или преобразованных дифференциальных уравнений на тестовую функцию, и далее проводится интегрирование по всей области. Использование глобального спектрального метода ограничено областями простой геометрической формы, что существенно сужает его применимость к реальным физическим процессам.
Метод спектральных элементов основан на тех же самых принципах, что и глобальный спектральный метод. Основное отличие метода спектральных элементов состоит в том, что интегрирование ведётся по части пространства независимых переменных, которую отождествляют с конечным элементом. С учетом свойств базисных функций удаётся выразить все интегралы через нули и веса наивысших аппроксимирующих полиномов. Таким образом, мы приходим к системе алгебраических уравнений для определения значений искомой функции в узлах сетки, определённой способом построения конечно-элементного разбиения и положением нулей наивысших полиномов на каждом из элементов разбиения. В таком случае для достижения необходимой точности расчётов и плотности сетки, накрывающей расчётную область, нет необходимости использовать излишне большое число базисных функций на каждом из элементов. Весь комплекс этих мер приводит к существенной экономии вычислительных ресурсов без потери спектральной точности и даёт возможность проводить вычисления в геометрически сложных областях реалистичной формы.
1 Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант РФФИ № 08-01-00484-а).
Как и в методе Галеркина, в спектральном методе приближенное решение исходного нелинейного дифференциального оператора
С («) = / (1)
разыскивается в виде комбинации по линейно независимой системе элементов {ф/}, называемой базисом в пространстве искомых функций, или координатной системой:
г. К
Р = IСФ; . (2)
/=0
В качестве элементов фг- используются ортогональные функции, являющиеся собственными функциями задачи Штурма - Лиувилля на единичном отрезке. Числовые коэффициенты с0находятся из системы алгебраических уравнений:
(в (20 С Фг ^ > V ^ = ( /, V ]), У = 0--^ . ()
Здесь {уу} - базис проекционной системы. В принципе, базисы {ф/} и {уу}, имеющие одинаковую размерность, различны, однако ниже мы принимаем
Фг = V; (г ).
Используя спектральное разложение, для достаточно гладких функций можно получить экспоненциальную скорость сходимости приближенного решения к точному [1]. В этом и состоит основное преимущество спектрального метода: очень точные приближенные решения могут быть получены при небольшом числе слагаемых, входящих в Ркны, причём ошибка аппроксимации будет уменьшаться экспоненциально с ростом N. Таким образом, определив коэффициенты приближенного разложения сй можно получить значения искомой функции в любой точке области с заданным порядком точности. Указанная схема характерна для классического глобального спектрального метода.
Недостатком глобального спектрального метода является то, что многочлены Якоби являются ортогональными на отрезке [-1, 1] и, следовательно, утверждения об экспоненциальной скорости сходимости приближенного решения к точному имеют место только в случае, если область интегрирования представляет собой отрезок ^ = [-1,1]. Для того чтобы решить задачу с произвольной областью интегрирования, необходимо найти замену координат, переводящую исходную область интегрирования в единичный отрезок.
Спектральный метод обобщается на случай двух и более измерений путём использования в качестве базисных функций тензорного произведения соответствующих одномерных базисных функций [2].
В случае двух и более измерений задача о нахождении преобразования координат становится достаточно сложной, особенно если область интегрирования имеет сложную форму. Поэтому имеет смысл разбить исходную расчетную область на конечные элементы и искать локальные представления решения через специальные функции, определенные на этих элементах. Однако при дискретизации, полученной на основе локального спектрального метода, матрица системы становится плохо обусловленной, что приводит к медленной сходимости итерационных методов. Эта проблема, как и в случае обычных локальных аппроксимаций, решается с использованием методов на основе пространств Крылова и подбором пред-обуславливателей.
В настоящей работе в качестве базисных функций мы использовали интерполяционные многочлены, представляющие собой комбинации полиномов Лежандра и их производных и получившие название полиномов Гаусса - Лежандра - Лобатто [2]:
и (х) = X и£1 (х); (4)
/=0
—1 (1 — х )
С(Х)" N(Н + 1)^ (х) х — х ' ()
С (X') = 5„ ,
Здесь ду - символ Кронекера, - точки Лежандра - Гаусса - Лобатто, определяемые формулой
хо =-1 , х] - нули , 1 < у < N -1 ,
хк = 1. (6)
Веса, необходимые для численного интегрирования при использовании интерполяционных функций (4), (5), определятся соотношением
! \2 ,•/ = О’1’-’N. (7)
^ ^ ^ (х7. )2
В то же время сами полиномы Лежандра определяются следующим рекуррентным соотношением:
А) ( Х ) = 1 ,
А( X ) = X ,
А+1 (х) = х£я(х)--П-^(х),п = 1,2..^ . (8)
п +1 п + 1
2. Уравнение Пуассона
Рассмотрим уравнение Пуассона:
-У2и = /(х), х е О,
«50 =Ф (х ).
Пусть X = (х1;х2,х3)е^3 и . Рассмотрим проблему построения решения
(9) в слабой постановке:
-1 У2им^ = | /У(/Г . (10)
п п
К левой части уравнения (10) применим формулу Грина:
ди
-1Vгиус1¥ = -1УиУус1¥- 4 V = |№У . (11)
(9)
дО дП О
Здесь п = (щ, п2, п3) - нормаль к поверхности дО. В результате наша задача (9) переформулирована для пространства функций Соболева [3, 4]:
и, V € (О) = |и(х) | иж = 0, |и е Ь2(О) |. (12)
Формулу (10) с учётом (11) можно преобразовать следующим образом:
ди
-1ЧиЧъИУ = |/ус1У + 4 . (13)
О О дО дп
Рассмотрим кубическую область 0 = (-1,1 )3 с ^3. В этом случае интерполяционная формула может быть рассмотрена в виде прямого произведения одномерных базисных функций:
N 3
и(X) = I иг1,2^ П СР1 (х). (14)
г1=1,г'2=1,гз =1 г=1
На элементе рассмотрим тестовые функции вида
3
^ЪЧ2т (х) = ПСт (Х). (15)
г=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 (21 ^ср3(21 ^,р\_Щх\,р2Ю;21 • (16)
1=1 дг дг )
Аналогичным образом преобразуется и правая часть (13). Таким образом, используя формулы (16), а также меняя индексы рх,р2,р3, получим матрицу системы и вектор правой части.
3. Преобразования координат
Поскольку многочлены (5) являются ортогональными лишь на единичном кубе, то для областей более сложной геометрии необходимо осуществить замену координат, переводящую исходную область интегрирования к единичному кубу. Пусть х = х(£) - некоторое преобразование координат, где ^ є [-1,1] для всех
і = 1,2,3 . Тогда
ди ди д£ ди дп ди дх ;
дх д£ дх дп дх дх дх ’
ди ди дЪ ди дп ди дх
— =-------- +----------------------------- +- ’ (17)
ду дЪ, ду дц ду дх ду
ди _ ди д£ ди дп ди дх
дг д£ дг дп дг дх дг
Формулы (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 г.
