Fedorenko-RP-Vvedenie-v-vychislitelnuyu-fiziku (810773), страница 70
Текст из файла (страница 70)
Общая картина в таких алгоритмах такова, что метод оказывается тем быстрее сходящимся, чем больше в нем доля метода Ньютона. Однако итерация метода векторной прогонки требует больших вычислений. Вообще, следует отметить, что основное время работы ЭВМ связано не с прогонкой, а с вычислением коэффициентов системы (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У +/„,!/!,г'„— У На рис.
37 показаны фрагменты !ислениых решений, полученных по всем трем схемам. Онн соответствуют моменту ! = 30 (при Ь = 1), т.е. ударная волна прошла! 17 счетных точек, контактный разрыв— 88 точек. ПоложениЯ точных гРаниЦ хп хт, хм х4 изобРажены иа рис. 37. Обсудим результаты.
Схема Е. Дефекты численного решения очевидны: сильно размазанный контактный разрыв, скорость ударной волны занижена примерно на 15 % (3.3 вместо 3.9), заметно размыты и слабые разрывы (границы волны разрежения). Схема П. Повышение формального порядка аппроксимации привело к существенному ухудшению результатов: графики функций искажены сильными осцилляцнями явно нефизического характера. Схема 777 (гибридная)': Существенное улучшение качества решения очевидно, хотя и ие все дефекты численного решения ликвидированы.
В частности, контактный разрыв размыт больше, чем хо- ПРИБЛИЖЕННЫЕ МЕТОДЫ НЫЧИГЛИТЕЛЬИОЙ ФИЗИКИ 1ч, и зФЕ телось бы. ПРи хзс х с х давление Р,„,нз= 11.35й0.05 (точное значение 11.40), и = 2.93Ф0.01 (точное значение 2.92). В последние годы в вычислительной газовой динамике ведется активная работа по конструированию схем с улучшенными свойствами решений. Целью этой работы является получение таких схем, в которых контактные разрывы размываются как можно меньше и не проявляются нефизические осцилляции.
Отличительной чертой таких конструкций является использование различных анализаторов локальной гладкости решения в каждой точке (л, Рп). В зависимости от показателя гладкости решения используется либо схема первого порядка аппроксимации, либо второго, либо некоторая промежуточная («гибрид» схем разного порядка аппроксимации).
О качестве схемы судят по качеству решения задачи о распаде разрыва в начальных данных и других задач-тестов, й 23. Приближенное решение двумерных задач газовой динамики Прикладные задачи газовой динамики, кзк правило, не допускают явных решений, поэтому важное значение имеют методы приближенного решения. В настоящее время ведется интенсивная разработка таких методов. Их создано уже достаточно много, тем не менее работа продолжается. Это объясняется тем, что одни и те же уравнения газовой динамики описывают (в зависимости от тех или иных краевых условий, значений входящих в уравнения физических постоянных) качественно разные явления, Они часто очень сложны, и эффективный метод решения должен учитывать характерные особенности подлежащего расчету явления.
Именно стремлением учесть специфику явления при конструировании расчетной схемы определяется содержание научной работы в области численных методов газовой динамики. Если полагаться на простейшие разностные схемы, мощность существующих ЭВМ окажется явно недостаточной для решения задач, которые выдвигаются современной техникой и достаточно успешно решаются специализированными методами. При изучении различных схем решения уравнений газовой динамики нужно прежде всего четко представлять себе, каков класс задач, в которых эффективен именно тот, а не другой из многих известных методов.
Эту сторону вопроса мы постараемся разъяснить в процессе изложения. Формулировка задачи газовой динамики. В дальнейшем мы будем иметь дело с так называемыми двумерными задачами, т,е. с задачами, в которых искомые функции зависят от времени Г и двух пространственных координат х, у. Конечно, реальные задачи газовой гашении лиумигных злдич гизовой динимики заз й гз1 динамики трехмерны; мы ограничимся двумерными ради простоты изложения. Основные идеи построения методов можно объяснить уже в двумерном случае.
Переход к трехмерному случаю вносит осложнения, в основном, технического характера (в то же время переход от одномерного случая к двумерному вносит ряд принципиальных осложнений). Другая причина состоит в том, что большая часть современных расчетов в газовой динамике — пока чтО двумерные; освоение массовых трехмерных расчетов по существу только начинается. Итак, мы имеем дело с некоторой областью пространства .Р, разные части которой заполнены разными газами. Заметим, что термин «газ» не следует понимать слишком узко. В определенных условиях (при высоких температурах) металлы ведут себя, как газы, и описываются теми же уравнениями газовой динамики.
Короче, мы имеем дело со сплошной средой, состояние которой описывается следующими функциями: компоненты скорости и(/, х, у) и и(/, х, у), плотность р(ц х, у), давление у(/, х, у), удельная внутренняя энергия е(д, х, у), Они удовлетворяют уравнениям газовой динамики. Существуют разные формы этих уравнений, удобные в тех или иных ситуациях. Начнем с уравнений в форме Эйлера: ди ди ди ! др — +и — +и — + — — =О, да ах ду р дх аи аи ди 1 ар — + и —. + и — + — — О, д! дх ду р ду дЕ ае да р /аи аиз — +и — +и — + — ~ — + — ~ =О. аа ах ау р(ах ду) Система (1) замыкается конечным соотношением — уравнением состояния, связывающим термодинамические характеристики среды р, р, е в каждой точке (/, х, у). Уравнение состояния используется в виде е= Е(р, р) или р= Р(е, р), где Е, Р— известные функции.