Бабенко - Основы численного анализа (947491), страница 149
Текст из файла (страница 149)
1, 2, ..., ЛХ). Учитывая ураннение (35), замечаем, что и системе (58) отличны от нуля коэффициенты А>,>-н, Ал,— с, Лм, А> ><->, Л>,>он. Поыипяу 727 9 2. Вариационпме методы реигенггл краевых задач систему можно записать в блочном виде шгедующим образом: %,5, г 4- 11,5, -~- Ь,5,эг = п„д = 1, 2, ..., ЛХ вЂ” 1, 5 — "0; (59) причем Зг =- О.
Матрицы 55 трехдиагональные, а матрицы Зь 45, диагональные. Система (59) по виду напоминает разностнуго систему. 19.3.17) с тем лишь отличием, что теперь коэффициенты системы — матрицы, а не вещественные числа. Поэтому метод прогонки можно организовать совершенно аналогичным образом, только нужно учесть, что теперь умножение не коммутативно. Из уравнения, отвечающего г = 1, находим Х, = — 55,'сам д, =55,',, (60) 5,=Х5, дп и поэтому можем последовательно получить соотношения прямого хода про- гонки: 5 —. Х5мы —,.
дь г = 2, 3,, .. ЛХ вЂ” 1. (61) т1тобы получить формулы для матриц Х, и векторов д,, подставим протопочное соотношение в г — ' 1-е уравнение. В результате получим ю эг(Х 5 э1 + д ) + 55 ь15 -1 т ч) э15 эг = и -1, откуда — г 1,г .-г, 1,г рр — гр =р 2 ' 2 и поскольку И = ягас) Ф, то безразмерное давление равно р = —,, — —,, [(1 — ~ю ) + р„, +;,.'э 1 + р 1 1 Приведем некоторые результаты расчетов. Потенциал обтекания эллипсоидов несжимаемой жидкостью можно вычислить в явном виде.
Мы не будем приводить эти классические формулы, а отметим, что в слу.чае шара Т = = (х: х, '( 11 потенциал имеет вид Ф = гп -с хчХ'(2~хне). Расчет производился для шара и для ряда зллипсоидов вращения. Ыы приведем ре:зультаты тсггько для обтекания шара.
В табл. 1 приведены результаты расчета при Х = 25, ЛХ вЂ”.. 40. Триангуляция по всей области правая. В первом столбце таблицы указаны номера точек, расположенных на поверхности шара (нумерация, как Х,зг =- — (ЖаюХ, з-)),-г) 'б,эп дбю = (З* гХ х— )5 — 1) '(и, г — %+19.). (62) Соотношения (60), (62) позволяют рекуррентно определить матрицы Х„и векторы до т. е, сделать прямой ход прогонки, Протопочное соотношение (61) при г — -- ЛХ вЂ” 1 позволяет определить 5м —. поскольку 5м известно (равно нулю). Пользуясь теперь соотношением (61) при г = ЛХ вЂ” 2, ЛХ вЂ” 3, ..., 1, находим векторы 5,. Метод матричной прогонки требует = лргЛХ операцгпй и поэтому проигрывает по числу операции хорошим итерационным методам.
Однако он прост в реализации и хорошо распараллеливается. Как окончательный выход алгоритма нас интересует не только потегп1иал скорости, по и давление, и вектор скорости, особенно на поверхности обтекаемого тела. Давление вычисляется по формуле. вытекающей из интеграла Бернулли; 728 Глаза 10. Нскаслорые вопросы часлснного решения краееыт задач '!'аблнца 1 Номер точки р* 0,50 0,4992 — 0,023 — 0,005 — 0,056 — 0,043 0,483 0.466 0,496 — 0.116 — 0,220 -0,347 -0,127 — 0,227 — 0.350 0,462 0,446 0,419 0.,433 0,385 0,397 — 0,487 0,354 0.338 0,292 0,304 0,.250 0,241 10 0,186 0,191 0,129 0,120 12 0,060 0,065 13 — 0.00008 0,0 — 0,06196 — 1,047 — 0,989 — 1.082 — 1,009 17 19 -0,346 — 0,386 20 21 -0.420 — 0,220 — 0,116 — 0.219 22 — 0,116 — 0,045 — 0,043 —.
0,039 — 0,005 25 на рис. 6); точки отстоят друг от друга на расстоянии по дуге, равном ксс24. Во втором столбце расположены вычисленные значения потенциала Ф,, затем точное значение Ф„вычисленное давление р„и точное давление р,. При определении давления принимшюсь р = — 1с'2. Давления приведены посредине счетных интервалов. Как видим, хотя при расчете бралось 1000 узлов, точность расчета весьма скромная.
Объяснить это тем, что произошла замена бесконечной области на конечную, невозможно, поскольку окружность Гя имела радиус, равный 121,2. При этом надо учестч, что отыскивалась функпия уз = -0,1300 — 0,188 — 0,243 -0,294 --0,450 — 0,470 — 0.484 — 0,506 — 0,065 — 0,129 — 0,191 — 0,250 — 0,304 — 0,354 — 0,397 -0,462 — 0,483 — 0,496 — 0.,500 — 0,641 — 0,.636 — 0,774 — 0,778 — 0,898 — 0,905 --0,984 — 1,009 -1,073 -1,082 -1,103 -1,120 — 1,106 — 1,120 — 0,884 — 0,905 -0,766 -0,778 --0,612 — 0,636 — 0.,073 - 0,082 — 0,339 —. 0,347 9 2.
Вириационьме мет.одм реше1и л краевых задач 729 —.. — х~?? (2,х' ), регулярная вне тела. Возможно ли, что взятая триангуляция ,з была неудачной? Для проверки этого предположения были произведены расчеты с три.шгуляциями, симметричными относительно оси у. Рассматривались два варианта триангуляции: в первом смена ориентации триангуляции производилась на оси у, а во втором область разбивалась на шесть секторов равного раствора и ориентация менялась с правой па левую и затем на правую поочередно, начиная с узла Рь Расчет производился при рй —..
21, ЛХ вЂ”.. 40. Результаты расчета представлены в табл. 2. Ввиду. симметрии решения приведены значения потенциала на теле в узлах Рм ..., Рм; в таолице Ф( -- вычисленный потенциал при одной смене ориентации триангуляции, Ф, — при пяти сменах и Ф, — значение аналитического решегшя. Как видим, частая смена ориентации пользы не приносит. Результаты для симметричной и несимметричной триангуляций в данной задаче практически не различаются. К сожалению, мы пе можем привести результаты расчетов эллипсоидов вращения: отметим любопытный факт — с увеличением эксцентриситета точность расчетов увеличивается.
Таблица 2 Сделаем несколько выводов из рассмотренного примера. При проектирш вании и построении алгоритма, основанного на методе конечных элементов, прежде всего нужно абсапотно четко, не по шаблону., сформулировать вариационную зада ел Затем основная трудность, которая ожидает нас, это — построение сетки и триангуляции области. Выше об этом говорилось. К вопросу об использовании способа интерполяции нужпо подойти очень внимательно.
Отметим, что практически не реально ввиду громоздкости использовать интерполяционные многочлены высокого порядка, к тому же будет сказываться величина константы Лебега. Тем пе менее для отыскшп1я глЮГких решений вряд ли целесообразно пользоваться линейной интерполяцией. В расслютренной задаче не возникала проблема вычисления интегралов.
В общем случае это уже не так, и для вычисления 1штегравов по стандартному симплексу или прямоугольнику. (параллелепипеду) вужпо пользоватыя кубатурными формулами (см. З б гл. 6). Мы уже отмечали в гл. 9, что при решении неоднородной задачи, особенио ври не очень гладкой правой части, временная сложность алгоритма может всецело определяться сложностью вычисления правой части системы (32). 12. Итерационные методы решения линейной системы. Особую заботу при проектировании алгоритма уделяют методам решения полученной системы линейных уравнений. Итерационные методы решения являются самыми экономными, и поэтому огромпое значение 730 Глава 10. Некотории вопроси численного решения краеонк задач имеет величина числа обусловленности лзатрицы системы.
Если для триангуляции выполнено условие (48) и 06 < 6 < 6, где 6 -- максимальная сторона симплекса Й, то при линейной интерполяции получим, что матрипа А системы (32) имеет число обусловленности соЫ А < С6 Это неравенство несложно доказать, но мы тем не менее опустим доказательство. Однако сформулированные ограничения на шаги 6 очень обременительны, и, как показывает разобранный в предыдущем пункте пример, они могут не выполняться в некоторых случаях.
Так, в предьн дущей задаче при обтекании шара мы имели г, = еГ' з, 6, = е«' — е«е поэтому 6ы з,(6з .= й (при принятых параметрах, использовавших- Л! — 2 ся в расчетах, зта величина = 10г). В этом случае предыдущая оценка не верна, н в ней 6 нужно заменить на 6 кн Во всяком случае. вопрос о числе обусловленности матрицы системы требует всегда дополнительного исешедования. Матрица системы разностных уравнений или уравнений, полученных методом конечных элементов, является разреженной. Местоположение ненулевых элементов матрицы А =- (а, ) можно охарактеризовать множеством позиций, в которых расположены ненулевые элементы матрицы %(А) = ((6 «): а, ф О).
Величины саге( Я(А) и сагс( Я(,4)Х з, где Х --. порядок матрицы, характеризуют степень разреженности матрицы. Расгзогзожение ненулевых злелиентов в матрице можно охарактеризовать следующим образом: егчи даны два ненулевых элемента в позициях (ц «) и (з', «ч), то назовем пртеле множество (го,,1о)., (1ь, «ь), где зо = 1, «о = «, гь .=. г', «л — —. «', (г'„, «'„) Е %(.4), причем для двух соседних пар (ьч «„), (д„ез, «, д) либо з, =д„ез, либо «е Число 6 назовем длиной путя.
Среди путей, соединяющих элементы (г, «) и (Е, «ч), существует путь мипималыюй длины. и мы обозначим ее через Й[(1, «): (г', «')]. Если в %(А) имеются элементы (з, «6 (ю', Д, которые нельзя соединить путем, то будем полагать е( (1. «); (1'., «')] = оо. Пусть в(А) = шах еК'(6 «); (з', «о)]. 0, В 4 ч зЧшл00 Введем операцию пореслановки строк и столбцов матрицы А.