The Design of Univrsal Algorithm Which Implements Inhomogeneous Dirichlet and Neumann Boundary Conditions in Spectral Element Method .pdf 1. Спектральный методКак и в методе Галеркина, в спектральном методе приближенное решение исходного нелинейного дифференциального оператораG (и) = / .(1) разыскивается в виде комбинации по линейно независимой системе элементов {фг-}, называемой базисом в пространстве искомых функций или координатной системой:P = I С Ф; .(2)В качестве элементов фг- используются ортогональные функции, являющиеся собственными функциями задачи Штурма - Лиувилля на единичном отрезке. Числовые коэффициенты с0cN находятся из системы алгебраических уравнений:G (I C Ф, j, V = ( f, V j), j = 0...N .(3)Здесь - базис проекционной системы. В принципе, базисы {ф,} и {у,-}, имеющие одинаковую размерность, различны, однако ниже мы принимаемф; =Х|/, (i = 0...N) .Итак, упомянутая выше задача Штурма - Лиувилля формулируется следующим образом:-- ( a (x)| + b (x)ф; =Xjw(x)ф;, a > 0, b > 0, dx \ dx JФг (-1) = фг (1) = 0. (4) В общем случае решением задачи (4) являются полиномы Якоби. Так как полиномы Якоби взаимно ортогональны на интервале [-1,1], можно доказать, чтоVu € U :||и - PNhu\ - 0 при N -^со .(5)здесь PN u = 2 СФг- .Более того, если u е Hm (Q), иными словами, если искомое решение u является m раз непрерывно дифференцируемым, то, согласно [1], ошибка аппроксимации будет следующей:и -- QN т ||UIнт (6)п iii2Таким образом, используя спектральное разложение, для достаточно гладких функций можно получить экспоненциальную скорость сходимости приближенного решения к точному. В этом и состоит основное преимущество спектрального метода: очень точные приближенные решения могут быть получены при небольшом числе слагаемых, входящих в PNhu, причём ошибка аппроксимации будетуменьшаться экспоненциально с ростом N. Таким образом, определив коэффициенты приближенного разложения с,, можно получить значения искомой функции в любой точке области с заданным порядком точности. Указанная схема характерна для классического глобального спектрального метода.Недостатком глобального спектрального метода является то, что многочлены Якоби являются ортогональными на отрезке [-1, 1] и, следовательно, утверждения об экспоненциальной скорости сходимости приближенного решения к точному имеют место только в случае, если область интегрирования представляет собойотрезок Q = [-1,1]. Для того чтобы решить задачу с произвольной областью интегрирования, необходимо найти замену координат, переводящую исходную область интегрирования в единичный отрезок.Спектральный метод обобщается на случай двух и более измерений путём использования в качестве базисных функций тензорного произведения соответствующих одномерных базисных функций.В случае двух и более измерений задача о нахождении преобразования координат становится достаточно сложной, особенно если область интегрирования имеет сложную форму. Поэтому имеет смысл разбить исходную расчетную область на конечные элементы и искать локальные представления решения через специальные функции, определенные на этих элементах. Однако при дискретизации, полученной на основе локального спектрального метода, матрица системы становится плохо обусловленной, что приводит к медленной сходимости итерационных методов. Эта проблема, как и в случае обычных локальных аппроксимаций, решается с использованием методов на основе пространств Крылова и подбором предобуславливателей.В настоящей работе в качестве базисных функций мы использовали интерполяционные многочлены, представляющие собой комбинации полиномов Лежанд-ра и их производных и получившие название полиномов Гаусса - Лежандра - Ло-батто [2]:NU (Х) = 2 UiCi (х)-1(1 - x2) L'NC(X} = N (N + 1)Ln (хг.) ~^-XT,(8)C (*/) = 8(/ ,здесь 5j - символ Кронекера, Xj - точки Лежандра - Гаусса - Лобатто, определяемые формулойxo = -1, Xj - нули LN ,1 < j < N -1,xN = 1 .(9)Веса, необходимые для численного интегрирования при использовании интерполяционных функций (7), (8), определятся соотношением21 j = од,^ N .(10)В то же время сами полиномы Лежандра определяются следующим рекуррентным соотношением:L0 (x) = 1 ,L( x) = x,2n +1 nL+i (x) =xLn (x)-Ln_x (x),n = 1,2...N .(11)n +1n + 1Полиномы (11) получаются из решения задачи (4) при1a(x) = (l - x2 )2 , b (x) = 0, w (x) = 1.2. Несжимаемые вязкие течения. Математическая постановка задачиУравнения, описывающие двумерные стационарные несжимаемые ламинарные течения, имеют вид[du w/ г \ 1 . w - + V(u u)Au = -Vp,5? 1 ; Дe (12) V-u = 0.d (13)-j = 0.dxjЗдесь j - индекс суммирования, Re = UL - число Рейнольдса, U, L - характернаяvскорость и характерный линейный размер соответственно, v - кинематическая вязкость, u = {u1u2) - векторная функция, представляющая скорость жидкости в плоском сечении, p - скалярная функция давления жидкости.3. Метод проекцийДля численного решения уравнений (12) применили метод установления совместно с широко известным методом проекций [3]. Суть этого метода состоит в следующем:1. На первом этапе находим промежуточное значение для скорости из уравнения и - И 1-AЬ = (иn -V)иn. (14)AtRe2..На втором этапе производим расчет давления по формулеn+1 nu - uV2 pn+l =-Vь .(15) At3..На третьем этапе производим расчет скорости для временного слоя n +1 следующим образом:Vpn+l .(16)AtШаги 1 - 3 выполняются до установления решения.4. Дискретизация уравнений Навье - Стокса методом спектральных элементовПусть x = (xj,x2xm)ё Rm и Qc Rm . Рассмотрим проблему построения решения уравнений (14) в слабой постановке:JvdV- - jV2wvdV = j/VdV .(17) n AtRe n nЗдесь f = u"V«" , u - одна из компонент скорости, v е L2 (Q) - некоторая произвольная функция. Ко второму слагаемому из левой части уравнения (17) применим формулу Грина:- Jv2ьvdv = - J jwьwvdr - 4 v -dMRe 0Re ^gO dn )У 'Здесь n = (n1; n2 nm) - нормаль к поверхности SQ. В результате наша задача (14) переформулирована для пространства функций Соболева:u,v€ H (П) = ju(x)\uga = 0, € L2(П)| .(19)Формулу (17) с учётом (18) можно преобразовать следующим образом:f^-^- vdV + - f VwVvdV = f fVdV + - С v- dS. (20) a AtRe ££Re /n dnРассмотрим различные типы областей ьc Rm . Наипростейший случай - это кубическая область Q = (-1,1 )m с ^m . Поскольку мы хотим получить метод, позволяющий аппроксимировать решение с высоким порядком точности, то и интегралы, входящие в формулу (20), необходимо вычислить с высоким порядком точности. Для этого используем квадратурные формулы Гаусса. В нашем случае интерполяционная формула может быть рассмотрена в виде прямого произведения одномерных базисных функций:N mVqb...,qm (x) = п Cq, (x; ) .(2)ь (x) = I ьPlPm п CPi (x) .(21) На элементе рассмотрим тестовые функции видаm( x) =Ii,-i=1Для простоты, предположим, что мы имеем граничные условия Дирихле. В таком случае нет необходимости производить расчёт интеграла d v - dS. Конкретную реализацию граничных условий мы рассмотрим более подробно ниже.Подставим (21) и (22) в уравнение (20). Рассмотрим второе слагаемое, стоящее в левой части уравнения (20). Тогда, с учётом того, что область интегрированияQ = [-1,1]m с Жт, получим•·f VwVvdV = Re О•·1 11 mNd md m•·= ■R- i - i ii "д,...,ft, T-п Cft (Xi)T-п Cq; (Xi )dx1 - dxm•·Далее, выражение (23) может быть переписано как•·f VwVvdV = Re О1 mN11 Л mЛ mR ii ьPi,...,Pm j -j f-П cPi (x пВ конечном итоге получаем(23)(24)11 d шd шI " iCp. (Xiп С,. (xi ) dx1 " dxm =-1 -1 dxk /=1 1 dxk i=1 1\ dCpk (Xk ) dCqk (xk )» 1, ч, чj-?7dxk п i Cp. (x )C,. (x,)dx,.Воспользуемся квадратурными формулами Гаусса:IC (x,) C (x,) dx. = £ ш/; Сд (x1 )Cqi (x1) .(26)Используя ортогональность интерполяционных функций и формулу (26), выражение (24) перепишем следующим образом:1 \\ ™ du dv ,- j J Lt-■^~dx1 dxm =Re _! _! k=i dxk dxk1N1 dCpk (xk ) dCqk (Xk ) , ЛЯ (27)Re д =1,...,pm=1k=1 _1 dxkdxki=U*kiSC. (xk) dCq (xk) Проинтегрируем выражение j--dxk с использованием квад--1 dxk dxkратурных интегральных формул Гаусса:i dC» (X' dCq. ' „x, = £ Bj dCp (X) SC(X). (28)-1 dxkdXkj=ldXk dXkЗдесь ю7- - веса квадратур, xJk - нули исходных интерполяционных функций. Внастоящей работе мы использовали интерполяционные функции и соответствующие им квадратурные веса и нули, определяемые формулами (7) - (11). Распишем первое слагаемое, входящее в левую часть формулы (20):-1 -1 (29)INm1N m=- 5] U P n П „ 8 „ qU P „ П „ 8 „ q .'-'Pi Pm*-*- Pi Pi Ai \t '-'Pi Pm*-*- Pi Pi ,1iАналогично распишем правую часть формулы (20), используя выражения (21), (22), а также свойство ортогональности базисных функций и квадратурные формулы Гаусса:II11 тI I f ( x ) v (x )dx1 dxm = I I f (x )ь Cq,. (xi )dx1 dxm = 1 q/=1I ... I f (x'1, x22,-, )П Щ C, (x' )= (30)NN Iii' \ mЗдесь со/. - веса, x1 - нули полинома в i-м измерении.Используя формулы (25), (28) - (30), а также меняя индексы q =qm = 1..N, получим матрицу системы и вектор правой части.Остановимся подробнее на методе аппроксимации граничных условий. В настоящей работе мы воспользовались технологией метода конечных элементов, применяемого для решения задач с ненулевыми граничными условиями. Мы полагали неизвестными все точки расчётной области, включая граничные. В случае граничных условий Дирихле нет надобности производить расчёт интегралаcf v -dS, поскольку значения искомой функции определены в точках границы.an dnФиксируя индекс q,, можно получить строку коэффициентов спектрального разложения для нахождения узловых значений зависимой функции, в том числе и в точках границы. В случае граничных условий Дирихле для граничной точки, соответствующей индексу q,, получаем строку матрицы, состоящую из одного ненулевого элемента, равного единице и расположенного на диагонали матрицы. В правую часть в таком случае мы записываем значения функции на границе. Для внутренней точки, соответствующей индексу q;, мы получаем строку коэффициентов, где число ненулевых элементов в строке равно Nm, где N - степень полинома, m - размерность задачи. Коэффициенты разложения определяются формулой (27). Если же какая-либо точка й _;Р попадает на границу, то, в случае граничных условий Дирихле, мы переносим значение функции в этой точке с соответствующим спектральным коэффициентом в правую часть. Таким образом, мы достигаем согласованности интерполяционных функций и граничных условий. В случае граничных условий Неймана, для каждой граничной точки, соответствующей индексу q,, мы получаем строку матрицы, состоящую из спектральных коэффициентов, определяемых формулой (27), к этим коэффициентам мы также должны добавить вклады от расчёта интеграла cf v -dS. Так как производная -an dn 5« является заданной, то произвести расчёт вышеупомянутого интеграла не представляет трудностей. Далее, для каждой внутренней точки, соответствующей фиксированному индексу, мы вновь получаем строку матрицы, состоящую из спектральных коэффициентов, определяемых формулой (27). Таким образом, мы опять получаем согласованность интерполяционных функций и граничных условий. Такой подход позволяет получать решения высокого порядка точности для граничных условий любого типа, кроме того, обработка граничных условий становится универсальной и не составляет никаких трудностей.Для аппроксимаций уравнений (15), (16) используем технологию, аналогичную описанной для аппроксимации уравнения (14).Необходимо отметить тот факт, что интерполяционные функции являются ортогональными на интервале Q = [-1,1]m с Rm . Следовательно, если область интегрирования отлична от единичного куба, необходимо найти преобразование координат, преобразующее исходную область к единичному кубу. Поскольку в качестве сеточных элементов мы имеем четырехугольники, то преобразовать сеточный элемент к единичному квадрату не составляет сложности.5. Результаты расчетовСопоставление результатов с известным аналитическим решениемДля того чтобы показать спектральную скорость сходимости численного алгоритма, рассмотрим течение Коважного [4]. Течение Коважного является точным решением плоских стационарных уравнений вязкой жидкости, т.е. плоских уравнений Навье - Стокса и имеет следующий вид:U (X, У ) :v (x,y ) =а - e(-Xx) cos (2ny), ^ e-Xx sin(2ny),2n(31)p (x,y ) = - 2eПараметр X, в свою очередь, определяется следующим соотношением:Х ='Re2 л 2 Re 4 2Численное решение было найдено в области прямоугольной формы [-0,5; 1] х [-0,5; 1,5]. Граничные условия задавали согласно аналитическим формулам (31). Расчеты велись для числа Рейнольдса, равного 40, на сетке, состоящей из 81 конечного элемента. Относительные погрешности для компонент скоростей, давления, а также невязка уравнения неразрывности показана на рис. 1. Из рис. 1 видно, что скорость сходимости приближенного решения к точному носит экспоненциальный характер.ю-1ю-2 ю-3Ю-4 Ю-5 Ю-6 Ю-7 Ю-8 Ю-9ю-10012Расчет течения в плоской кавернеРассмотрим двумерную полость, представляющую собой область квадратной формы с длиной грани / = 1. Нижняя и боковые грани являются твердыми стенками, верхняя грань - подвижной стенкой, перемещающейся с постоянной скоростью. Граничные условия для данной задачи задавались следующим образом: u = 0 на твердых неподвижных стенках, щ = 1, u2 = 0 на подвижной стенке,- = 0 на всех границах. dnЧисло Рейнольдса связано с физическими параметрами задачи следующей формулой: Re = Ul / v, где U - размерная скорость перемещения подвижной границы.В табл. 1 приведены результаты расчетов минимального значения функции тока для течения при числе Рейнольдса 400. Сравнение с данными других авторов показывает достоверность данного подхода к аппроксимации уравнений Навье - Стокса.На рис. 2 приведены линии тока для чисел Рейнольдса 400 и 1000. Результаты для числа Рейнольдса, равного 1000, можно сравнить, например, с данными работы [5]. На рис. 3 приведены линии тока для числа Рейнольдса, равного 1800.Стационарное течение за уступомРассмотрим область, представляющую собой канал с уступом. Пусть L - длина канала за уступом, d - диаметр входного сечения, а также высота уступа, D -диаметр канала за уступом (рис. 4).■ L + d-и2d_____Dи-L-иxРис. 4. Область расчета стационарного потока за уступомГраничные условия задаются следующим образом. На входе течение считается установившимся, поэтому первая компонента скорости имеет параболический профиль, а вторая компонента скорости тождественно равна нулю. Проинтегрировав уравнение Пуассона для продольной компоненты скорости по диаметру входного сечения, в предположении, что отсутствуют продольные изменения рассматриваемой величины, получим следующее граничное условие для первой компоненты скорости:дРЩ (О, y) = Re-дх' У-2■- y +22(32)гдеардх■ заданный и постоянный градиент давления, Re = U 2d / v - число Рей-нольдса, U - средняя скорость на входе, 1Х, /2 - пределы интегрирования входного сечения в направлении у. Граничное условие на входе для второй компоненты скорости задается следующим образом:«2 (о, У) = о.Граничные условия для давления на входе в этом случае(33)дРdx -a,(34)где a - некоторая положительная константа.На выходе течение также считается установившимся и граничные условияdnр (х y ) = 0.На твердых стенках для скорости ставится условия прилипания, а также однородные условия Неймана для давления:u (x, y) = 0;dp0.dn(36)0Часть расчетной области с сеткой, состоящей из четырехугольников, представлена на рис. 5. Отметим, что линейный размер L составлял 20d.В работе были проведены расчеты для чисел Рейнольдса 100, 200, 300, 400. Степень полинома равнялась 3, а число конечных элементов не превышало 3 000. За уступом формируется глобальный вихрь, продольный размер которого / увеличивается с ростом числа Рейнольдса. В табл. 2 представлена зависимость длины большого вихря за уступом от числа Рейнольдса. Здесь следует отметить, что экспериментальные данные для двумерных течений отсутствуют при числах Рей-нольдса больше 200.ЗаключениеТаким образом, в работе был предложен эффективный алгоритм расчета уравнений Навье - Стокса, использующий неструктурированные сетки, преобразования координат и спектральные разложения искомых функций. Также был предложен универсальный алгоритм реализации неоднородных граничных условий Дирихле и Неймана, позволяющий без потери спектральной точности применять данный метод к широкому спектру задач. Решения, полученные с помощью метода спектральных элементов, могут иметь самостоятельные значения, если нас интересуют тонкие эффекты, сопровождающие исследуемый процесс, или же они могут применяться как тестовые результаты при построении алгоритмов метода конечных разностей.

Скачать электронную версию публикации