книги / Метод крупных частиц в газовой динамике

единице) корни уравнения

det — S G>J<*kCjk(x\ t*)— XI = 0

/,*= 1

имеют неположительные действительные части

i.k=i

где С — матрица аппроксимационной вязкости.

Если система дифференциальных уравнений (4.29) является гиперболи­ ческой, то для нее целесообразно ввести в рассмотрение характеристики и область зависимости Q точки (х, t). Для гиперболической формы п. д. п. (4.30)’ можем также ввести понятие области зависимости Q' точки (л:, t). Оказалось, что условия устойчивости разностных схем определяются из рассмотрения параболической формы их первых дифференциальных приближений и характера зависимости областей Q и Q'.

Ниже приводится ряд теорем, доказательства которых содержатся в-

Если матрица С^О , то первое дифференциальное приближение (III) является

неполной параболической системой уравнений (что и определяет устойчивость

простой разностной схемы (I), аппроксимирующей систему уравнений (4.32), необходимо, чтобы ее первое дифференциальное приближение было неполной

ференциального приближения является необходимым и достаточным условием устойчивости простой разностной схемы (I).

Т е о р е м а 3. Если А — симметрическая матрица, А +, А " — Липшицнепрерывные матрицы, то для устойчивости мажорантной схемы (II) необ­ ходимо и достаточно, чтобы ее первое дифференциальное приближение было неполной параболической системой уравнений.

метрические матрицы и первое дифференциальное приближение простой схемы расщепления является неполной параболической системой уравнений, то схема устойчива.

дифференциальное приближение мажорантной схемы расщепления является неполной параболической системой, то схема (IV') устойчива.

Аппроксимируем систему уравнений (4.29) следующей мажорантной

Т е о р е м а 6. Если A j — симметрические матрицы, A t, A j — липшиц-

Метод первого дифференциального приближения был также применен к исследованию устойчивости и аппроксимационной вязкости разностных схем,, аппроксимирующих нелинейные уравнения (в частности, уравнения газовой динамики), хотя строгого обоснования этому пока не имеется. В этом случае

и позволяет провести качественное сравнение разностных схем в смысле ха­ рактера их устойчивости [25, 62, 218—221].

Возвращаясь к ранее изложенному, заметим, что система уравнений (4.30) — г. ф. п.д. п.— получается из разностной схемы (4.28) разложением; в ней функций. ип+1(х), ип(х+х%а) в ряды Тейлора и удержанием членов по­ рядка 0(т2) включительно, а система уравнений (4.31) — п. ф. п. д. п.— получается из системы уравнений (4.30) заменой d2uldt2 через пространствен­ ные производные с использованием исходной системы уравнений (4.29). Удер­ живая в разложении ип+1{х) члены более высокого порядка, чем 0(т2), и используя при замене d2uldt2 в N-u дифференциальном приближении (N —1)-е- дифференциальное приближение, получим в матрицах «вязких членов» (Cjky дополнительные члены, зависящие от градиентов решений. Эти добавочные- члены объясняют неустойчивость некоторых разностных схем, наблюдаемую- в расчетах, которая не улавливается при анализе локальным методом Фурье,, ибо он не учитывает, как уже отмечалось, влияния градиентов.

Аналогичный подход применил К. Хирт [25] при исследовании разност­ ных схем для нелинейного уравнения газовой динамики.

4. Вернемся теперь к исследованию разностной схемы метода крупных: частиц [1, 23, 29].

Рассмотрим, например, какой вклад в неустойчивость дают разные формьг записи уравнения неразрывности, считая, что уравнения импульса и энергии' устойчивы. Здесь и далее мы будем рассматривать каждое уравнение газоди­ намической системы в отдельности, учитывая тем самым не все члены матриц. Практика показывает, что такое упрощение в ряде случаев правомочно.

Если АМп считалось по формулам второго порядка точности (3.11), то* разлагая приведенные разностные уравнения в ряд Тейлора и удерживая

Здесь (в (4.33) и в (4.34)) Аи Д£ — члены разложения порядка не выше поряд­

При практических расчетах, когда возникают ударные волны, контактныеразрывы и волны разряжения:

Отсюда видно, что в (4.34) коэффициент при d2pIdx2, положителен, а в (4,33) — отрицателен, т. е. схема (4.33) допускает быстрорастущие решения и является неустойчивой при счете, в то время как схема (4.34) устойчива.

Практический счет подтверждает это, что следует, в частности, из рис. 4.2, где показаны профили р на линиях r=const перед цилиндрическим торцом

Рис. 4.4. Профили плотности р за ударной волной для различных значений М при использо­ вании АМп второго порядка точности (цилиндрический торец, v = l) .

для этих двух случаев. Видно, что для (4.33) (сплошная линия) сразу за удар­ ной волной возникает существенная неустойчивость. Эта неустойчивость ил­ люстрируется также рис. 4.4, на котором приведены вдоль линии r = R /3, параллельной оси симметрии, профили плотности для схемы (4.33) при раз­ личных числах М я (1,1<Д1о& 1 4 |5). Устойчивость вычислений при исполь­ зовании (3.11), (3.11') достигается введением в схему явного члена вязкост­ ного давления q.

1 0 6 ИССЛЕДОВАНИЕ ЧИСЛЕННЫХ СХЕМ МЕТОДА КРУПНЫХ ЧАСТИЦ [ГЛ. IV

Из (4.39) и (4.40) видно, что при отсутствии вязкостного давления схема неустойчива, и лишь введение q делает ее устойчивой (Re со ^ 0). Тем не менее из-за простоты записи и удобства счета можно использовать заведомо не­ устойчивую на эйлеровом этапе схему с ^ = 0 (на последующих этапах проис­ ходит компенсация этой неустойчивости). На конечный результат неустой­ чивость первого этапа практически не влияет, так как важна устойчивость разностной схемы на всех трех этапах в целом.

Все вышесказанное иллюстрирует и рис. 4.5, на котором приведен ха­ рактер установления по времени плотности р0 в точке торможения ( ^ = 2 ) . Штрихпунктирная и сплошная линии соответствуют использованию для вы­ числения АМ п формул первого порядка точности (3.12) с q и без q(q—Q)y со­ ответственно, пунктирная и штриховая — формул второго порядка точности {3.11) также с q и без q. Мы видим, что во втором случае прогрессирует не­ устойчивость, которую можно устранить лишь введением q. Вычисление же АМ п по формулам первого порядка точности (3.12) делает счет устойчивым в любом случае (введение q лишь незначительно ускоряет сходимость). При

Рис. 4.5, Характер установления по времени плотности в точке торможения при сверх­ звуковом обтекании цилиндрического торца М со=2,0; v = l.

этом установление по всем параметрам практически достигалось с точностью до 0,0001%.

Следует отметить также периодический, автоколебательный характер неустойчивости, присущий данной разностной схеме с использованием формул второго порядка точности (3.11). По мере приближения к искомому (устой­ чивому) решению градиенты параметров потока в областях гладкого течения •становятся малы. Вместе с ними становятся малыми вязкостное давление и вязкостные члены первого дифференциального приближения, в результате чего в этих областях вновь начинает прогрессировать неустойчивость. Когда •она разовьется настолько, что градиенты станут велики, вязкостные члены •опять начнут ее демпфировать и т. д., т. е. возникают автоколебания. Такая

§4] МАТРИЦЫ АППРОКСИМАЦИОННОЙ ВЯЗКОСТИ 107

периодическая автоколебательная структура неустойчивости для схем вто­ рого порядка точности отчетливо видна на рис. 4.5.

Как видно из построения, везде используются явные разностные схемы,, что требует выполнения определенных условий устойчивости. Применениездесь абсолютно устойчивых неявных разностных схем кажется менее целе­ сообразным, так как в этом случае сильно увеличивается область зависимости, разностного решения (для задачи Коши она становится бесконечной). В прак­ тическом счете это обстоятельство выглядит как усиленное сглаживание раз­ ностных профилей, что может привести к потере характерных особенностей: решения [27].

Таким образом, последовательное рассмотрение нулевого, первого и вто­ рого дифференциальных приближений позволило исследовать вопросы аппрок­ симации, образования вязкостных эффектов и устойчивости вычислений, т. е. все основные свойства полученных разностных схем метода крупных частиц.

5. Итак, была поставлена задача получить схему сквозного счета, позво­ ляющую непрерывным образом проводить расчет и поверхностей сильного разрыва. Здесь мы пришли к тому, что выполнение необходимых условий устойчивости для широкого класса газодинамических задач обеспечивается (при использовании формул (3.12)) механизмом диссипации, определяемым,, как правило, лишь схемной вязкостью (диссипативно-устойчивые схемы). Как показали многочисленные расчеты двумерных (по пространству) газодинами­ ческих задач, эти схемы обеспечивают устойчивый сквозной счет с постоянными сеточными параметрами практически во всей области интегрирования (при этом, например, ударная волна размазывается на 2—3 счетные ячейки и т. п.).

В [1, 20—23] проводились аналитические и численные оценки и обосно­ вания постановки краевых условий задачи, а также имел место контроль процесса вычислений и полученных данных («склейка» с асимптотикой, расчет на разных сетках аппроксимации, «продолжение» полей течений, сравнение с расчетами по другим методикам, экспериментом, и т. п.). Некоторые данные методического характера будут приведены ниже.

§ 4. Матрицы аппроксимационной вязкости

Получим теперь конкретный вид матриц аппроксимационной вязкости

для схем метода крупных частиц. Как уже отмечалось, рассмотрение таких матриц представляет большой интерес для определения диссипативных ха­ рактеристик разностной схемы, оценок условий устойчивости, возможности оптимизации схем по диссипативным свойствам и т. д. Исследования в этой области были проведены Ю. М. Давыдовым и др. [218, 346, 348—349, 437].

Изучим в одномерном случае структуры матриц аппроксимационной вяз­ кости для различных форм записи исходной системы уравнений.

1. Рассмотрим систему уравнений газовой динамики, записанную в ди­ вергентной форме *),

где а — матрица аппроксимационной вязкости, Aw — члены, пропорциональ­ ные Дх, Дх2, At> не содержащие wxx или wxxx.

*) Для сокращения записи будем здесь использовать такие обозначения: р*=др/д/,.

(ир)х =д (ир)1дх и т. д.

Будем считать, что поток течет слева направо. Для противоположного направления потока достаточно поменять Ах на (—Ах) в дифференциальном приближении.

Выпишем элементы матрицы а разностной схемы метода крупных частиц с АМ п первого порядка точности (3 .12) (учитываются члены порядка Ал:, Ад:2, At) [437]

аи = у иА х— у ихАх2+ ± ( р р— и2) At,

Условие положительности диагональных элементов матрицы а исполь­ зовалось в качестве критерия устойчивости схем метода крупных частиц.

2. Если систему (4.41) обозначить через L 2(w), где до={р, ры, р£}г (см., например, метод «потоков» [55]), то П-форму второго дифференциального приближения по Ад: и первого по At можно записать в виде

Следует заметить, что матрица р использовалась для исследования устой­ чивости разностных схем в работах Н. Н. Яненко, Ю. И. Шокина [223], Ле­ ра, Пейро [224] и др.

Представляет большой интерес получить матрицу (А и для метода круп­ ных частиц. Выведем формулы преобразования для перехода от элементов матрицы а к элементам матрицы р [346, 348, 437].

•3) i-к столбец матрицы р. Так как дифференциальные приближения (4.42)

Таким образом можно получить все элементы матрицы р разностной схемы метода крупных частиц. Для краткости изложения ограничимся только диаго­

Если учитывать только первое дифференциальное приближение (п. д. п.) разностной схемы, то можно установить взаимно однозначное соответствие между элементами матриц а и р, в которых учитываются только члены первого

Отсюда видно, что матрицы b и ц были введены для преобразования тех членов элементов матриц а и р , которые пропорциональны Дх2.

Если рассматривать разностную схему метода крупных частиц без расщеп­ ления (т. е. на лагранжевом этапе при расчете потоков массы, импульса, энер­

гии через границы ячеек использовать а, £ , а не а, Ё, полученные на эйлеро-