Федоренко Р.П. Введение в вычислительную физику (1185915), страница 69
Текст из файла (страница 69)
В данном слУчае, посколькУ а есть величина О(лг), Условие ом«нг — — 0 достаточно естественно. Таким образом, уравнение (2) для им можно представить в виде 12'м(и „и, Р пг, Р пг) = О. УчитываЯ, что Рм«пг задано кРаевым Условием, запишем его в окончательной форме: «Гм(им „им, Тм-пг, Тм+пг) = О (Тм«222 добавлено «для общности»).
Теперь все значения и обеспечены «своими» уравнениями. Перейдем к уравнениям для Т: е — е" «2П «2!2, и и ю+! ю + (Р+ Ч)~«222 А " +2а »«!/2 ! ! Тм+3!2 т«иг т««2 т-!/2 м а„„~ +! л„ 'и Эти уравнения можно использовать при лг= 1, 2, „М вЂ” 1, т.е. в системе пока не хватает двух уравнений. Исключая е пг через т„,+22, о «пг (УРавнение состоЯниЯ) и 22 «пг чеРез и, и +„пРидадим уравнениям форму Е«!+2!2(Т!«2пи Т«,+!22 Т„,+2!2, и,„, и„,+!) = О. ПРИБЛИЖЕННЫЕ МЕТОДЫ НЫЧИСЛИТЕЛЬНОй ФИЗИКИ 1Ч.
И 336 Уравнение для Тцг получаем, используя краевое условие (задан поток НТ„= Ц при х = 0): е — е" и — и Г Т вЂ” Т +ш»+цг,,-, -ь з о ! - зп иг,з1 Это уравнение можно представить в виде Е Гг(1 цг 1 зди ио и~) = О (ио включено тоже Фдля общности»), Включив в уравнение для е Фцг аппроксимационный источник, г +цг (вычисляемый по значениям и, а не й ), получим уравнение Е» цг = 0 точно такой же структуры. Уравнение для Т „иг в рассматриваемом случае тривиально— эта величина просто задана. Запишем это уравнение в общей форме, имея в виду и более сложные краевые условия: ем~из(тм цг, тмРиз, им) =О.
Подведем итог, выписав все уравнения, которые предстоит решатке 1го(ио из Тиг) =О (11) йI (и и ищ, ит»и Тт-цг Тт+цг) =О, им(им — з им Тм-цг Тм+цг) = О' Еиг(1'|ы Туг ио и~) =О, Еттцг(Тщ цг, ТщР цг Тщ» зтг ищ, ищ+,) = О, (12) Еи+цг(Тм-иг, Тч»цг, им) = О, где т = 1, 2, ..., М вЂ” 1. Перейдем к алгоритму их решения. Метод раздельной прогонки В атом методе сначала величины Т фиксируются как Т (т.е.
как уже найденные приближения к Т" +'), затем решаются уравнения (11) относительно и (линеаризацией по Ньютону). В результате получается линейная система уравнений относртельно и, имеющая ту же структуру, т.е. система с трехдиагональной матрнцей. Она легко решается методом прогонки (см.
5 10). Фиксируя ио»'>, можно линеаризовать вторую группу уравнений относительно Т. Линейная система с такой же трех- диагональной матрицей решается прогонкой, Далее эти процедуры повторяются до достижения требуемой точности. 337 РехлизАиии РАзностной схемы 5 221 Метод векторной нрогонки. В методике, которая описывается в этом параграфе, система уравнений на верхнем слое (11), (12) решалась методом векторной прогонки. (Раздельная прогонка была предложена позднее.) В методе векторной прогонки одновременно линеаризуются обе системы уравнений.
Эта операция приводит к следуюшим линейным уравнениям: 4оио + -4ои1 + АО Тцг = .4о 4~2Т112+ Вц,Туг+ Вцг ио+ В,12 и, = Во; 1 -ц2 1гг А и +А и+Аи++,+ о 1 А-цгТ -1- АцгТ = .4 Вт+цгТт-т+ В ьцгТ„,ицг+ В,'„+цгТ,згг+ -цг Ц2 +В +ц и +В цги +1=В +ц ', А 1им-1+ Аомим+ АмцгТм- г+ А~м~гТм. = Ам -1 по и-цг В +ц*т -т+"м+.Т + г+В + и =В +цг где нг = 1, 2, ..., М вЂ” 1. Вводя вектор г = (и„, Т +1~), запишем эти уравнения в матричной форме (т = 1, 2, ..., М вЂ” 1)1 .Вого+ Во21 ~о (13) Ф г~ 1+Я 2 +$' 2 „1=Я Л'мгм — 1 + ~мгм Рм.
Здесь использованы обозначения Формулы вычисления элементов этих матриц через значения функций У, В и их производных очевидны, но громоздки. Нет неооходим1кти их воспроизводить. Система уравнений (13) имеет «трехднагональную» форму и решается несложным обобщением алгоритма прогонки. Вывод формул алгоритма отличается от вывода, изложенного в 3 10, только тем, что теперь мы работаем с матрицами (некоммутативная алгебра) и надо аккуратно следить за порядком множителей. Решение ищется в форме ( в в'„„)' А1 0 цг 40 Ац2 -цг В ьцг в, ьцг (в „). ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ 1ч. И зза г, = Х г + У, где Х вЂ” матрица 2х 2, г — вектор.
Опуская простые выкладки, приведем результат (т = 1, 2, ..., М вЂ” 1): Х1 ~0 йо' Х +Б ('® +Л~ Х ) ~о ~о Теперь последнее уравнение (13) и прогоночное соотношение гм, = Хмгм + Ум можно РазРешить относительно гы. Эта величина вычисляется н позволяет начать «обратную прогонку» — вычисление справа-налево искомых величин г„. Мы не будем здесь обсуждать проблем разрешимости всех встречающихся в алгоритме задач (существования обратных матриц) и сходимости итерационного процесса, Укажем лишь, что легко угадать тривиальный результат: прн достаточно малом т все обстоит благополучно.
Это естественное следствие вырождения в пределе т- О всех уравнений в тривиальные. Теоретические оценки того малого т, начиная с которого герантнруется успех вычислений, в практике расчетов не используются — зто привело бы к неоправданно заниженному шагу по времени. Однако сам факт зависимости, например, скорости сходимостн итераций от т (она тем выше, чем меньше т) используется в режиме обратной связи. Считается, что требуемая точность должна достигаться за три-пять итераций.
Если итераций потребовалось больше, следующий шаг интегрирования уравнений выполняется с уменьшенным шагом т. Если точность достигается за меньшее число итераций, т увеличивается в пределах, определяемых другими критериями выбора шага. Что касается сопоставления скорости сходимости методов раздельной и векторной прогонок, то преимущество имеет последняя. Это естественно: оба метода являются комбинацией метода лннеаризации (Ньютона) и метода простой итерации (Пикара). Общая картина в таких алгоритмах такова, что метод оказывается тем быстрее сходящимся, чем больше в нем доля метода Ньютона. Однако итерация метода векторной прогонки требует больших вычислений.
Вообще, следует отметить, что основное время работы ЭВМ связано не с прогонкой, а с вычислением коэффициентов системы (13). Одним нз ресурсов экономии вычислительной работы является алгоритм с однократным вычиглением этих коэффициентов; при этом в процессе итераций пересчитываются только ЕР в (13). Поясним это на промере решения нелинейного уравнения у(х) = О. Упрощенный вариант метода Ньютона с однократным вычислением у, имеет форму ' — у.( ') у( '). 339 РЕАЛИЗАПИЯ РАЗИОСТИОЙ СХЕМЫ й 22! При достаточно хорошем начальном приближении х' он сходится «линейно», т.е.
11у(х'Н) - 9', где д ж ()Š— у„(х) /„'(хв) 11 (х — решение системы), Эту оценку предоставим вывести читателю. И наконец, подчеркнем, что приведенные выше формы организации решения уравнений на верхнем слое образуют некоторую общую схему„в рамках которой возможны различные варианты. Они появляются при различных способах отнесения тех или иных неизвестных к 1-й итерации (по таким неизвестным проводится линеаризация) или к (1 — 1)-й, а в иных случаях и к л-му слою. Выбор того или иного варианта диктуется особенностями решаемой задачи и здесь, естественно, не конкретизируется. Расчет ударной волны по неднвергентной гибридной схеме.
Сказанное выше может создать у читателя впечатление, что для правильного расчета ударных волн дивергентная форма разностных уравнений является существенным фактором. В общем это верно. Не следует только возводить это положение в ранг абсолютного, безусловного требования к используемым в расчетах схемам, Обсудим этот вопрос подробнее, опираясь на результаты вычислительного эксперимента, проведенного автором в 1962 г, Рассмотрим численное решение задачи о распаде произвольного разрыва в начальных данных. Кроме трудностей расчета ударной волны, мы имеем проблему расчета контактного разрыва, так как задача решается в переменных Эйлера.
Итак, уравнения имеют вид и,+ив„+с~ ~ =О, (с'+ 91 ч;+ ив„— ви„=О, г С+вс+ 2 и 9 †с +д В качестве основных термодинамнческих величин берутся удельный объем Р н величина с = ~~ (которая с точностью до множителя совпадает с адиабатической скоростью звука; рассматривается идеальный газ с у= 5/3), Вязкость берется в форме фон-Неймана: 4=(х,й) и„(и„— (и,1), где Ь 4 —:5. Выбор «экзотической» переменной с объясняется просто. В зоне размазанной волны (см. й 20) переменная с ведет себя так же, как функции и и с, тогда как р и е фактически размазываются на вдвое меньшую длину. Поскольку профиль волны должен быть разрешен четырьмя-пятью счетными точками, при счете в терминах е илн р пришлось бы увеличить |.
раза в два, что приводит к слишком большому размазыванию и и ж При решении задач в лагранжевых координатах этой проблемы нет, так как переменные е и р по х не дифференцируются, а шаг по времени по разным причинам таков, пгивлиженные методы вычислительноа «изики 340 что временное размазывание волны заметно больше пространствен- ного (содержит больше шагов т). Начальные данные имеют вид [2.0, 0.25, 1.688, ль< 0, и", цв+, и с" +из= ]О Можно найти точное решение задачи.
Оно состоит из: а) волны разрежения, левой и правой границей которой являют- ся линии х,(1) = — 0.891, х (г) = 0.390 б) контактного разрыва иа линии х (г) = 2.920 в) ударной волны на линии х4 — — 3.91. При хз(Г) < х < х4(Г) значения и(1, х) = 2.92, р(С, х) = 11.40. В этой области ц рвется на контактном разрыве: ц(1, х) =0.352 при х (К) < х < х (Г), и(Г, х) = 0.250 при х (~) < х < х4(Г). Мы имеем дело с так называемой «сильной ударной волнойь, идущей по «холодному газуы В этом случае скачок плотности при переходе через волну максимален (сжатие в (у + 1)/(у — 1) ж 4 раза). Используем явную схему: + и — „— ", '"+', =О, с + и — +ц« ц.
Ь Поясним некоторые, обозначения: общие индексы вынесены за квадратные скобки; и, из=0.5(и + и„„,); о =0.5(ц цз+ и, цз); и «цз есть значение (сз+ д)/ц, вычисленное по очевидной разностной аппроксимации л, причем значения ц и с берутся с (л + 1)-го слоя (сначала эти величины находятся из двух первых уравнений, затем считается и"«'); [Ьс[ «цз — — с цз — с„цз при и цз>0 и [Лс[ «цз= с «зп — с «цз пРи и «из< О.
Таким же обРазом («против потока«) берутся и разности Ьц, Ьи. Назовем вышеприведенную схему схемой 1. Ее основной дефект— первый порядок аппроксимации конвективной производной У, + иУ„. Эта величина (при и ) 0) аппроксимируется разностью типа т '[Д~„+ т, х ) — Д~„, х — ит)[. РеАлиздция Рлзностиой схемы 0 гг! 341 Значение х„— ит не попадает в узел сетки, поэтому в эту точку значение у интерполируется линейно по ближайшим узлам (л, ги — 1) и (и, гл).
Можно заранее предвидеть (см. 3 20), что схема ! приводит к размазыванию контактного разрыва. уточним схему в этом месте, вычисляя у(г„, х — ит) квадратичной интерполяцией значений У" „У", 7' ~Р Это будет схема П. Можно и здесь предвидеть неприятности, связанные с нефизическнми осцилляциями. Наконец, рассмотрим гибридную схему (схему 0,35 0,35 0,30 0,30 0,35 035 Рис. Зт И1), в которой используется линейная нли квадратичная интерполяция в зависимости от дифференциальных свойств решения в данной точке (см. 8 20). Они характеризуются отношением второй и первой разностей, например !У, — 2У +/„,!/!,г'„— У На рис.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
