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

Тип работы:
Реферат
Предмет:
Физика


Узнать стоимость

Детальная информация о работе

Выдержка из работы

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2008 Математика и механика № 3(4)
УДК 532+681. 3
А. М. Бубенчиков, В. С. Попонин, В. Н. Мельникова МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА И РЕШЕНИЕ ПРОСТРАНСТВЕННЫХ КРАЕВЫХ ЗАДАЧ МЕТОДОМ СПЕКТРАЛЬНЫХ ЭЛЕМЕНТОВ1
В работе описана математическая постановка и решение пространственных краевых задач методом спектральных элементов. Основной упор сделан на дискретизацию трехмерного уравнения Пуассона. В работе приведен подробный вывод преобразований координат, необходимых для решения задачи в областях сложной геометрии.
Ключевые слова: метод спектральных элементов, уравнение Пуассона, неструктрурированные сетки, вязкая жидкость.
1. Спектральный метод
Принцип, лежащий в основе всех сеточных методов, заключается в сведении исходных дифференциальных уравнений в частных производных к системе алгебраических уравнений, которые могут быть решены известными методами. В глобальном спектральном методе решение ищется путём разложения в ряд по некоторой системе ортогональных функций, называемых базисными. Имея представления искомых функций в спектральном пространстве, то есть определив их в виде разложения по базисным функциям, в глобальном спектральном методе строится система интегральных соотношений, получающихся умножением исходных или преобразованных дифференциальных уравнений на тестовую функцию, и далее проводится интегрирование по всей области. Использование глобального спектрального метода ограничено областями простой геометрической формы, что существенно сужает его применимость к реальным физическим процессам.
Метод спектральных элементов основан на тех же самых принципах, что и глобальный спектральный метод. Основное отличие метода спектральных элементов состоит в том, что интегрирование ведётся по части пространства независимых переменных, которую отождествляют с конечным элементом. С учетом свойств базисных функций удаётся выразить все интегралы через нули и веса наивысших аппроксимирующих полиномов. Таким образом, мы приходим к системе алгебраических уравнений для определения значений искомой функции в узлах сетки, определённой способом построения конечно-элементного разбиения и положением нулей наивысших полиномов на каждом из элементов разбиения. В таком случае для достижения необходимой точности расчётов и плотности сетки, накрывающей расчётную область, нет необходимости использовать излишне большое число базисных функций на каждом из элементов. Весь комплекс этих мер приводит к существенной экономии вычислительных ресурсов без потери спектральной точности и даёт возможность проводить вычисления в геометрически сложных областях реалистичной формы.
1 Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант РФФИ № 08−01−484-а).
Как и в методе Галеркина, в спектральном методе приближенное решение исходного нелинейного дифференциального оператора
С («) = / (1)
разыскивается в виде комбинации по линейно независимой системе элементов {ф/}, называемой базисом в пространстве искомых функций, или координатной системой:
г. К
Р = IСФ-. (2)
/=0
В качестве элементов фг- используются ортогональные функции, являющиеся собственными функциями задачи Штурма — Лиувилля на единичном отрезке. Числовые коэффициенты с0находятся из системы алгебраических уравнений:
(в (20 С Фг ^ & gt- V ^ = (/, V ]), У = 0--^. ()
Здесь {уу} - базис проекционной системы. В принципе, базисы {ф/} и {уу}, имеющие одинаковую размерность, различны, однако ниже мы принимаем
Фг = V- (г).
Используя спектральное разложение, для достаточно гладких функций можно получить экспоненциальную скорость сходимости приближенного решения к точному [1]. В этом и состоит основное преимущество спектрального метода: очень точные приближенные решения могут быть получены при небольшом числе слагаемых, входящих в Ркны, причём ошибка аппроксимации будет уменьшаться экспоненциально с ростом N. Таким образом, определив коэффициенты приближенного разложения сй можно получить значения искомой функции в любой точке области с заданным порядком точности. Указанная схема характерна для классического глобального спектрального метода.
Недостатком глобального спектрального метода является то, что многочлены Якоби являются ортогональными на отрезке [-1, 1] и, следовательно, утверждения об экспоненциальной скорости сходимости приближенного решения к точному имеют место только в случае, если область интегрирования представляет собой отрезок ^ = [-1,1]. Для того чтобы решить задачу с произвольной областью интегрирования, необходимо найти замену координат, переводящую исходную область интегрирования в единичный отрезок.
Спектральный метод обобщается на случай двух и более измерений путём использования в качестве базисных функций тензорного произведения соответствующих одномерных базисных функций [2].
В случае двух и более измерений задача о нахождении преобразования координат становится достаточно сложной, особенно если область интегрирования имеет сложную форму. Поэтому имеет смысл разбить исходную расчетную область на конечные элементы и искать локальные представления решения через специальные функции, определенные на этих элементах. Однако при дискретизации, полученной на основе локального спектрального метода, матрица системы становится плохо обусловленной, что приводит к медленной сходимости итерационных методов. Эта проблема, как и в случае обычных локальных аппроксимаций, решается с использованием методов на основе пространств Крылова и подбором пред-обуславливателей.
В настоящей работе в качестве базисных функций мы использовали интерполяционные многочлены, представляющие собой комбинации полиномов Лежандра и их производных и получившие название полиномов Гаусса — Лежандра — Лобатто [2]:
и (х) = X и?1 (х) — (4)
/=0
-1 (1 — х)
С (Х)& quot- N (Н + 1)^ (х) х — х '- ()
С (X'-) = 5″ ,
Здесь ду — символ Кронекера, — точки Лежандра — Гаусса — Лобатто, определяемые формулой
хо =-1, х] - нули, 1 & lt- у & lt- 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йх дх ду ду & amp- & amp-)
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
& lt-ІХ (Іу<-І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& quot-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 =
& amp- 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 ч
_& quot-д?ч У Г2 У 33 У'-3 Уз2 ^^& quot-дпІУ^^Уп У33 У1Ъ Узі ^+& quot-дхІ У^^У'-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 дцд%) 51 J Ад^дх д^д%) 5×1 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 г.

ПоказатьСвернуть
Заполнить форму текущей работой