Федоренко Р.П. Введение в вычислительную физику (1185915), страница 61
Текст из файла (страница 61)
Неявно опираясь на единственность ее решения, мы получаем тождество и, р, е(аг, ах) = и, р, е(г, х). Зто возможно лишь в случае, когда решение зависит не от двух переменных г, х, а лишь от одной автомодельной переменной $ = х/г. Итак, и, р, е(г, х) = //, я, е(х/г).
Уравнение газовой динамики становятся обыкновенными дифференциальными уравнениями, которые допускают достаточно обозримый анализ. Зти уравнения выводятся после замены операторов а д аВ ~ д а д а~ аа «В ВГ Г дС' ах д~ ах Г Иц' одноме~ ныо ггявнония гюозой динямнки 297 Использум уравнения в форме (14), после простых преобразований получаем систему уравнений для автомодельиого решеним: Р'=~и', и'= — ВУ', — М+Ри'=О (штрих — символ производной по ~). Для уравнения состояния Е = РУ!(у — 1) они легко интегрируются: ч-щ .~п Р(с) "о чзтд +О о т 21о -~ ! (1б) и($) = ~ — — о ц~~ Од~+~>. В результате мы имеем общее решение с двумя произвольными постоянными. Отметим важное соотношение Сз(г) = т РД)/УД) = ~2.
Здесь С есть скорость звука (отлвчие от формулы со= у2о/р связано с тем, что С вЂ” это «массовая» скорость звука, так как мы используем уравнение с массовой лагранжевой переменной), Таким образом, линия $= С, т.е. х = С1, является характеристикой. Заметим, что при с = О решение имеет особенность, которой можно избежать, используя это решение только в интервале Ц', $" 1 при $' < С" с О нли О с г,'с г,", Построим еще одно точное решение типа «центрированной волны разрежения» (для определенности, нддтшей вправо). Пусть и, р„е, произвольны. Вычислим с" =~/~р~р~ (это будет правая граница волны).
Речь идет о непрерывном решении, поэтому нам известим значения У($") = и = 1/р, н У($') ин что позволяет без труда вычислить посюявные Ув, ио в (16). При г, < ф" решение опасывается формулами (16). Левую границу волны О< о,'< $" можно назначить произвольно. Вычислим из (2($'), вз= У(с'), р = РД'). Эта значения определяют константное решение прн г, с г,'. Заметим, что и ударная волна, н контактный разрыв входят в семейство автомодельных обобщенных решений — это константные решения, рвущиеся при некотором значении $. Теперь можно вернуться к вопросу о распаде произвольного разрыва в начальных данных.
Решение является автомодельным и состоит (в общем случае) из контактного разрыва, справа и слева от которого расположена ударная волна нли центрированная волна разрежения, причем возможаы четыре сочетааия: все определяется расположением точек в„р„е, и в, рз, ез. Основаые трудности, возникающие при численном решении уравневай овзовой динамики, связаны с наличиеы разрывов в искомых решениях.
При конструировании численных методов обычно выделяют характерные особенности решений, т.е. строят задачи, в ПРИБЛИЖЕИИЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬИОЯ ФИЗИКИ [ч. н 298 которых наиболее трудные особенности существуют в чистом виде, без взаимодействия с не очень трудными для расчетов гладкими течениями. В таких задачах известно точное решение н качество расчетной схемы оценивается по тому, как оно справляется с решением «модельной задачи». Метод Годунова.
Для расчета разрывных решений широко используется метод, основанный на решении задачи о распаде разрыва. Пусть начальные данные являются кусочно-постоянными на некоторой сетке (х ) н о, т.е. и, р, е(0, х) = (ио + гг, ро+и, ео+пг) при х Я (х, х +,). Оказывается, эта задача имеет точное явное решение. Оно строится так: в ка:кдой точке х нужно решить задачу о распаде разрыва независимо от всех остальных разрывов.
Такое решение можно использовать до того момента времени, когда прн каком-то значении Рл правая волна, образовавшаяся от разрыва в точке х встретится с левой волной, идущей от распада в точке х + г Перейдем к описанию схемы. Основными счетными величинами ЯЕЛЯютсЯ н !гг~ Р~~Ф!ле е» Фгп тле л номеР ЕРемеиного слОЯ Ве личина и" +,, например, представляет собой приближенное значение функции и(г„, (х + х +!)/2). Пусть начальные данные и +пг, р +,и, е Фн превращены в кусочно-постоянные на сетке о о о (х ) функции.
Построим точное решение уравнений газовой динамики и определим малое (порядка шгп (х„+, — х )) время т, в течение которого это явное решение существует. Один шаг численного интегрирования соответствует продвижению по ! на г. Для того чтобы сделать второй шаг, нужно и при ! = г иметь кусочно-постоянное решение. Точное решение, конечно, таковым не является.
Поэтому оно аппроксимируется кусочно-постоянным с сохранением в пределах каждого интервала основных физических величин: массы, импульса н энергии. Например, (х„+, — х ) р' +н — — ~ р(т, х) с(х. Аналогично вычисляются ! ! ! ! г Р„+нг ив+!гг рт+!гг(от+!!г+ (и, +г!г) /2!. Однако и, р, е при ! = г являются сложными функцнямн х, вычислять интегралы от них трудно. Это препятствие обходится следующим образом.
Рассмотрим ячейку (х, х + ) х (О, г) и запишем уравнения в интегральной форме (3), взяв ячейку в качестве й. В результате одномзгныг. гк«внзния гАзозой диилмики 299 получим соотношение х +$ ~ Я[0, х[ !/х — ~ Яг, х +!! сй — ~ Я[«, х[ ах+ ~ Ц[/, х [ г/г = О, «„ о (17) где Я[!, к[ ш Я(и(/, х), р(б х), е(/, х)). В (17) третий интеграл есть то, что нам требуется. Первый интеграл вычисляется элементарно, так как Я[0, к[ = Я(и~ +!н» р«+!и» е«+нз) в силу постоянства начальных данных на (х, х +,). Так же легко вычисляются второй и четвертый интегралы. В силу атомодельности решения задачи о распаде произвольного разрыва на линии х = х,„(т.е.
$ = сопв1) все функции постоянны. Обозначим их ич, р, е . Тогда второй интеграл есть Д(и~, р, е ) т. Формально схему Годунова можно записать в виде л"+' — л" К+!и™»+ 2+ а ! а 0 ю+1 « Следует иметь в виду, что «поток» Ц вычисляется решением задачи о распаде разрыва. Она сводится к решению системы нелинейных уравнений. Это относительно «дорогая» операция (ведь она проводится при всех т, л).
Значительные усилия прилагаются к тому, чтобы снизить ее трудоемкость. В частности, используется то, что в большинстве узлов (и, т) величины слева и справа от к мало отличаются друг от друга, Разработанная в середине пятидесятых годов схема до сих пор применяется в расчетах; при этом она, разумеется, обобщается и совершенствуется. Расчет контактного разрыва. Проблемы расчета течения, содержащего контактный разрыв, рассмотрим, используя известное нам точное решение уравнений газовой динамики типа «чистый контактный разрыв».
В этом случае из всех уравнений газовой динамики нетривиально только одно: р, + ир, = 0 (здесь и = сопз!). Будем решать его методом сеток. Построим равномерную сетку с шагом /! по х и т по и Узлы сетки х = л!/!, г„= пт. Приближенное решение ищем в виде сеточной функции р" . Используем простейшую явную схему (предполагая и~ 0): (р" ' — р" )/т + и(рч — р" !)//! = О.
(18) Как известно, эта схема устойчива при условии Куранта ит//! « 1 (см. 3 12). (Заметим, что такая схема весьма популярна при расчете задач в эйлеровых координатах: во все уравнения в этом саучае пгнвлнжвнныв методы вычислительной еизикн !ч. и входит характерный оператор «субстанциальной производной» з/д~ + и д/ах.) При и < О используется аппроксимация с шаблоном, ориентированным в противоположную сторону (против потока): (р"+' — р" )/т + и(р"„», — р" )//с = О. Когда решается полная система уравнений газовой динамики, функция и(д х) может менять знак.
Соответственно, и разностные формулы строятся в точках (л, сн) в зависимости от знака и". Что же получается оооо при расчете контактного разрыва по такой о ~решение д схеме? Происходит неприятное явление. В ! расчете контактный разрыв размывается, !о его «ширина» растет с ростом времени. На рис. 31 показана характерная картина расо о чета: точное решение («ступенька») и прито зз зо з«4о еч т ближенное.
Это явление особенно неприятно в тех Рис. 31 задачах, в которых контактный разрыв раз- деляет среды с разными уравнениями состояния. Нужно знать достаточно точно границу между разными газами, чтобы пользоваться в данной точке нужнмм уравнением состояния. На рис. 31 представлены результаты, полученные в вычислительном эксперименте. Но для того чтобы бороться с «размазыванием контактного разрыва», нужно иметь хоть какую-то теорию. Этим мы сейчас н займемся. Заметим, что излагаемый ниже аппарат имеет достаточно общее значение, не ограничивающееся только задачей о контактном разрыве. Он может применяться при построении разностных схем, лучших в том нли ином отношении (смотря по тому, что нам нужно в задаче). Исследование дисперсионного соотношения для разностной схемы.
Подчеркнем, что этот аппарат, строго говоря, работает только в случае линейных уравнений с посгояннымн коэффициентами. Он применяется для линеаризованных моделей реальных уравнений, однако полученные в линейной модели рекомендации затем используются и в реальных задачах. Хорошо известно, что такое дисперсионное соотношение для линейного дифференциального уравнения с постоянными коэффициентами.
В нашем случае для уравнения р, + ир, = О дисперсионное соотношение появляется, когда мы ищем решение вида ех<ер+ы', где цк) = -и/сг. Оказывается, разностные уравнения (линейные, однородные, с постояннымн коэффициентами) имеют решения того же вида, но, конечно, с другой функцией Х(/с), зависящей от вигов т, /е и вида разностной схемы. Найдем дисперснонное соотношение для напи- зо~ одномвгныв жевания гАзовой динлмвки санной выше схемы, полагая р" = е"~«М'+и'»". Подставляя эту функцию в разностное уравнение, после очевидных преобразований получаем соотношение (е"' — 1)/т + и(1 — е на)/Ь = О.
Из него легко вычислить дисперсионное соотношение для схемы «против потока»: Х(/с, Й, с)= — !и 1 — и — „' +и — „'е '"" Естественно оценивать качество разностной схемы по степени совпадения дисперсионных функций для дифференциального и разностного уравнений. Идеальным было бы их совпадение.
Оно обеспечивается условием ит/Л = 1. Этот идеальный случай, к сожалению, практически не интересен. Соотношение их/л = 1 в реальной задаче, когда значение и не является постоянным, а меняется во времени и в пространстве, во всех точках сетки выдержать нельзя. Поэтому связанными с ним преимуществами воспользоваться в практической работе не удается. Разумеется, функции Х(к, л, т) должна аппроксимировать А(к).
Параметр Ь есть малый параметр, и аппроксимация, естественно, тем лучше, чем меньше волновое число к (чем глаже по х рассматриваемое частное решение); на сетке с шагом Ь волновые числа й > 2п/л уже не реализуются. Суждения о качестве разностной схемы можно делать, сравнивая графики Х(й) и Цк, 6, т). Некоторые выводы можно получить, считая кл«1 (в смысле кл ц 0.5, например) и разлагая в ряд Тейлора хотя бы с точностью до второго члена: Ц/с л, т) ж — /ик — — ~1 — и — ~. ил« / т~ 2 ~ ь~' Сравним частные решения дифференциального и разностного уравнений: е м(ть-ило е!Н(ть-«по е — и«(ь-«оп'«2 Сделаем некоторые качественные выводы из полученной формулы. 1. Как отмечалось, малоинтересен специальный случай ит = /с 2. При и < 0 схема непригодна: решения разностного уравнения отличаются от решений дифференциального множителем порядка «*ь е~" ~» ""', который при к ж 1/л катастрофически (при л, т-»0) растет как с~и~с!л 3.