Fletcher-1-rus (1185917), страница 67
Текст из файла (страница 67)
Использование разложений в ряды Фурье Рассмотрение ошибок аппроксимации, приведенных в табл. 9.1 для разностной схемы со сдвигом против потока, схемы Лакса — Вендроффа и Кранка — Николсона в применении к примерам, рассмотренным в конце $9.1, очень упрощается, ЗТТ 4 9.2, Численная диссипация и дисперсия (9.39)' Однако более точное представление дТ/дг на и-м временнбм.
слое можно записать в виде дТ !л Та ы за дз7. (а дзз да7, (и — — 0.551 — ~ — — — ~ — о (л7 ). д~ ~ ЛЗ д7з ~ 6 д~з Из (9.2) имеем д'Т дзТ вЂ” =и' —, дР дхз ' дзТ дз дТ цз дЗз = дхз д1 ' Следовательно, более точная (чем (9.39) ) полудискретизация уравнения (9.2) дается выражением — "...1 Л1 д 1 Т"'~~ — Т дТ е дТ 6 д ' ~ М +и дх О.би Ж дхТ=О.
(9.40) Использование сетки с постоянным шагом и метода Галер- кина с конечными элементами и линейной интерполяцией ($ 5.4) приводит уравнение (9.40) к виду (̄— — Лх 1.„х ~(Т7'+' — Тз") + СЛх1 „Т7~ — 0.5С Ьх ЕзхТ; = О, (9.41) если обратиться к фурье-компонентам, входящим в решение. Результаты, указанные на рис. 9.2 — 9.4 и в табл. 9.1, показывают, что значительная диссипация численной схемы эффективно компенсирует малую дисперсию для разностной схемы против потока и ослабляет часть дисперсионной ошибки в решении Лакса — Вендроффа. В схеме Кранка — Николсона преобладают дисперсионные ошибки, которые приводят к значительному замедлению распространения каждой из фурье-гармоник.
Хотя все схемы Кранка — Николсона устойчивы, с увеличением числа Куранта ошибка растет как Сз (табл. 9.1). Для малых С и, следовательно, малых Лг дисперсионная ошибка схемы Кранка — Николсона в конечных элементах (табл. 9.1) значительно меньше, чем у других схем, рассмотренных в 2 9.1. Использование разложений в ряды Фурье и модифицированных уравнений позволяет построить более точные численные схемы, уменьшающие ошибки диссипации и дисперсии.
Это можно проиллюстрировать на методе Тейлора — Галеркина с конечными элементами [Г)опеа, 1984]. Согласно этому методу, сначала проводятся дискретизация по времени, которая для уравнения (9.2) дает Т" е~ — Т дТ дз + и — =О. дх З78 Гл. 9. Линейние задачи с преобладающим влиянием конвекпии которое является обобщением схемы Лакса — Вендроффа (9.16).
Вид массового оператора М„= Я, з, 1 ~ вытекает из формулировки метода конечных элементов. Член (С'/6) Ь (Т! — Т";) возникает из-за сохранения дополнительных слагаемых в разложении производной дТ/дх в ряд Тейлора. Уравнение (9.41) содержит ошибку аппроксимации порядка 0(Ага, Ахз), что можно сравнить со вторым порядком точности схемы Лакса — Вендроффа (9.16). В работе [Ропеа, 1984] с помощью разложения в ряды Фурье исследованы дисперсия и диссипация схемы; найдено, что уравнение (9,41) превосходит схему Лакса — Вендроффа для всех длин волн и чисел Куранта в интервале 0 < С < 1.
Приведенная выше схема может быть обобщена введением .дополнительных членов в рядах Тейлора для управления (анти-)диссипацней и (анти-)дисперсией. В работе ]Вакег, Кнп, !987] выведена такая общая схема и названа формулировкой Тейлора — Вика (ФТВ) метода Галеркина с конечными элементами. ФТВ-алгоритм содержит четыре произвольные постоянные. Специальный выбор этих постоянных дает 17 хорошо известных численных схем, сформулированных в методах конечных разностей и конечных элементов. Оптимальный выбор может быть сделан в принципе и для линейного уравнения конвекции (9.!О). Однако для более общих гиперболических уравнений сохранения, например уравнений Эйлера (10.40), оптимальным будет решение, которое обладает хотя и желательной диссипацией, но и значительной дисперсией. Мы можем сделать общий вывод о том, что следует избегать аппроксимационных формул первого порядка для производных.
Аппроксимация первого порядка для п и производной в исходном уравнении будет порождать пространственные про. взводные порядка (и+ 1) и выше в модифицированном уравнении (которое является эквивалентным решаемому исходному уравнению). Это особенно важно для конвективных членов (и = 1), так как введение фиктивных производных второго и третьего порядков может существенно изменить характер ре;шения.
9 9.3. Стационарное уравнение с конвекцией и днффузней Для многих задач обтекания диссипативные механизмы яв.ляются существенными только в узком слое, обычно вблизи границы обтекания. Численные решения, построенные на сетках, приспособленных к основному потоку, часто начинают осциллировать в пограничном слое, где точное решение претер- 379 Е 9.3. Стационарное уравнение певает резкое изменение. Стационарное уравнение конвенции — диффузии является полезным модельным уравнением, с помощью которого можно проиллюстрировать это явление. Это уравнение в безразмерных переменных имеет вид Д стяг и — „— а — „, =О.
(9.42) тГх тГхт Оно выражает стационарный баланс между конвекцией и диффузией. Время не играет никакой роли. Если к (9.42) добавить граничные условия Т(0)=0 и Т(1)=-1.0, (9.43) то в интервале 0 < х ( 1 точное решение можно представить в виде Т = [е '/ — 1.01/(енсе — 1.01. (9.44) ° Лх Е.г, Я 1,=4 о ах=в.1, Я =г се11 х ах=в.ез,я 1 ' се11 — - точное оешенне "«-гв ьо 0.5 -о. 0.5 1.О Ряс. 9.3.
Влияние числа Рейяольдса ячейки ка решение уравнения кояаекцяя — дяффуэяя. иым распределением Т почти на всем интервале, кроме тонкого пограничного слоя в окрестности х = 1. 9.3.!. Влияние сеточного числа Рейнолодса Если использовать трехточечные центральные разности для производных в (9.42), то получим следующий алгоритм: (1 + 0.5!чссв) Тг 1 + 27/ (1 0.51очсса) 71ч.1 = 0 (9 45) где тт'„и = иЛх/а, т. е. 1с„п — число Рейнольдса (строго говоря, Пекле), построенное на характерной длине ячейки Лх.
Это решение, приведенное на рис. 9.5, характеризуется постоян- 380 Гл. 9. Линейные задачи е иреобладакяиим влиянием конвекнии Разложение уравнения (9.45) в ряд Тейлора около 1-го узла показывает, что оно аппраксимирует уравнение (9.42) с точностью до 0(йхо). Разностная схема (9.45) является неявной, но если выписать уравнение для всех узлов, то получается система уравнений с трехдиагональной матрицей, для решения которой могут быть использованы экономичные способы (п.
6.2.2). Для относительно простого случая точное решение (9.45) можно представить в виде 1 + 0.5к Т; = Ао+ Во ~ 1 — 0.бо„ц 1 ' (9.46) (9,47а) АТ = В, То Тз — аТ, Ь с аЬс (9.47Ы Т,, аЬс а Ь вЂ” сТ, где а = — (1+ 0 5еоеец), Ь = 2 и с = — (1 — О 5)геоц) Собственные числа вышеприведенной трехдиагональной матрицы связаны с а, Ь и с следующими соотношениями: Лу —— Ь+2(ас)" сов( у, ), 1=1, 2, ..., У вЂ” 2.
(9.48) где Ао и Во выбираются таким образом, чтобы удовлетворялись граничные условия (9.43). Для случая и/а = 20 и Лх = 0.05, 0.1 и 0.02 решение (9.46) приведено на рис. 9.5. При Лх = 0.05, что эквивалентно Й„ц = 1, решение довольно хорошо согласуется с точным. При Веец = 2 (Лх = 0.1) решение совпадает с точным везде, кроме пограничного слоя в окрестности х = 1. Для более грубой сетки при Я„ц = 4 решение не только является неточным, но и носит колебательный характер. Из (9.46) видно, что колебательных решений не будет, если )7„ц ( 2.
Полезно также связать поведение решения, показан. ное на рис. 9.5, с собственными числами матрицы системы уравнений (9.45), записанных с граничными условиями (9.43). Эта система имеет вид 4 9.3. Стационарное уравнение 38! Видно, что для того, чтобы собственные числа были вещественными, необходимо удовлетворить условию ас ) 0; (1 + О.Исса) (1 — 0 асан) ~ )О нлн )ссса ~( 2. (9.49) Таким образом, появление осциллирующих решений связано с появлением комплексных собственных чисел.
Для матрицы А в (9.47а) обще~о вида ограничение на собственные числа Л дается теоремой круга Гершгорина (см. [)епи!ппз, 1977[) [Л вЂ” ам[( ~ [ап [. /ФЕ Таким образом, собственные числа матрицы А лежат в кругах, связанных с каждым рядом матрицы А. Центр каждого круга определяется диагональю с элементами ап, а радиус — суммой внедиагональных элементов каждого ряда. В работе [Сагеу, 5ерепгпоог1, 1980) применялась теорема круга Гершгорнна для анализа уравнений вида (9.47а) при различных значениях В„н. В интервале 0 ( Р,м|( 2 центр круга Гершгорина лежит в точке ( — 2а/Аха, 0) комплексной плоскости, имеет радиус 2а/Аха и все собственные числа расположены на отрицательной части вещественной оси.
С увеличением Вссн, т. е. с уменьшением а, радиус круга уменьшается, а центр перемещается к началу координат. При Я„ц = 2 центр круга расположен в точке ( — и/Ах, 0), радиус равен и/Ах = = 2а/Ахе. При Исси ) 2 радиус остается постоянным, но центр движется к началу координат и собственные значения являются комплексными. При Я„н = ео центр круга совпадает с началом координат, а собственные значения становятся чисто мнимыми величинами. Ограничение В„н ~ В для получения неосциллирующих решений применимо и к другим методам, в частности к методу конечных элементов. В работе [Г!п!ауооп, 1989] получено значение В для метода конечных элементов и метода ортогональной коллокации.
Однако обычно В меньше 10. Следует подчеркнуть, что осциллируюшие решения могут появиться, если )с„н ) В, а появятся ли они на самом деле, зависит от геометрии потока и граничных условий. Например, замедление потока перед препятствием часто порождает осциллирующее решение при Я„н ) В. Однако для равномерного или ускоряющегося потока можно построить неосциллирующее решение и при )есся ) В. Как и следовало ожидать, если для представления дТ/дх в (9.42) использовать не центральные разности, а разности против потока (Т; — Тг,)/Ьх, то осциллирующие решения не 382 Гл.