Федоренко Р.П. Введение в вычислительную физику (1185915), страница 67
Текст из файла (страница 67)
В уравнения на верхнем слое входят близкие величины трех видов. Поясним это на примере Т +,1; (то 'ке самое относится и к Ии»112, и ). Мы имеем дело с Тл»1пи Тф,112, Т<'+11))2. Величины Т", Пз (с нижнего слоя) считаются уже известными. При вычислении Т"+.,', 2 методом итераций фигурирукп уже найденные значения Т1'1»пз (1-е приближение к Т"„,.',и) и неизвестные значения Тф) . В пределе величины Т~') 1)2 — Т".,1вз.
Мы используем обозначения Т пз — — Т +„„Т лп,— — Т »112. ИменнопоотиошениюкнеизвестО) О+1) ным Т пз производится линеаризация при выполнении очередной итерации. Аппроксимация уравнений движения: ил+1 — ил 1р+д)л+1 — 1р+д)"»1 + Р Я "+гп Р ~ -' О (2) Я 326 ПЦИЕЛИЖЕИИЫЕ МЕТОДЫ ЦЫЧИ!'ЛИТЕЛЬНОЙ ФИЗИКИ (ч. и Уравнение для энергии: и+Ц2 (ц +1)2 ! (ц +1)2 (ци)2 ! (ци )2 — е"+' + Ф.~. 1 Ы + 4 ц!+Ц2 4 + ' 1(/2 + В)ц+1 иц+! (р + В)ц»1 иц»1~ Ф+1П т"+' — т"+' т"Ы т»Ф! ! ФЗП, ФФЬЯ ц+(а — !а — х "»Ф!и '( "„+! (уравнение имеет смысл прн п) = 1, 2, ..., М вЂ” 1). К уравнениям (2) — (5) следует добавить выражения для величин 2/ +,!2, р, х, не входящих в число основных счетных величин, а) Неймановская зязкосты/ » ц2 имеет вид (и) = О, 1,, М вЂ” 1) 2/ Ф! =(е/ц «ц2) ((й„,~) — й ) — )й +,— й ~)(и~~! — и ).
ПРи атем вм»ц2 = О («фиктивное» гРаничное Условие). б) Интерполяция р осуществляется по естественной формуле цт-!пцц+!и+ц»цу -т " -т+(2»«и2 Это — линейная интерполяция, имеющая формально погрешность аппроксимации О(й2), т.е. минимальную, обеспечивающую при численном дифференцировании погрешность О(й).
Онс была введена на заключительном этапе эксплуатации программы. Первоначально использовалась рассчитанная на равномерную сетку формула р = (р ц2+ р„„.цз)/2. В дальнейшем сетка стала неравномерной, а формула осталась, что привело к определенным трудностям, о которых будет сказано подробнее ниже. Заметим, что формально ошибочная формула с полусуммой рекомендуется и сейчас. Это приводит к существенным ограничениям сетки х: она должна быть «почти равномерной», т.е. й +,!2 = й цз(1 + О(й)).
в) Коэффициент теплопроводности х вычисляется по формуле, в которой совмещаются идеи линейной и «гармонической» интерполяции (см. з 21). Обсуждение разностной схемы. Приведем соображения, на основании которых были выбраны вышеприведенные формулы аппроксимации. Два типа узлов. Разделение счетных точек на целые (механические) и полуцелые (термодинамические) является очевидным следствием различного вхождения соответствующих величин в уравнения. В каждое уравнение входят производные по времени от вели- згт шхлизхпия глзиостиой схемы з гг~ чин одною сорта и производные по х от величин другого сорта. Это кажется нарушенным для уравнения энергии, но выполняется и для него, если использовать эквивалентную недивергентную форму (б) е, + (р+ ц) и„= (хТ„),.
В дальнейшем мы заменим приведенную выше аппроксимацию уравнения для энергии (5) на эквивалентную, но столь же компактную, как и очевидная аппроксимация уравнения в форме (6). Используемая сетка позволяет при минимальном шаблоне получить (на равномерной сетке) второй порядок аппроксимации по х, а при аппроксимации пространственных производных комбинацией аппроксимаций по верхнему и нижнему слоям с весами 0.55 и 0.45, например, можно получить «почти второй» порядок аппроксимации по т. При весах 0.5 и 0.5 порядок был бы вторым, однако схема стала бы 1по спектральному признаку) нейтральной, т.е.
точки спектра, соответствующие параметру р = и, оказываются на единичной окружности и высокочастотнме паразитические возмущения хотя и не нарастают катастрофически, но и не затухают. Кроме того, формальный разностный порядок производных совпадает с формальным их дифференциальным порядком. В связи с этим схема не требует дополнительных краевых условий (или «односторонних» разностных аппроксимаций в ближайших к краям счетных точках), необходимых при превышении разностным порядком схемы истинного порядка дифференциальных выражений.
Дивергентность схемы. Все разностные уравнения имеют так называемый дивергентный вид, т.е. они могут быть записаны в следующей форме (для термодинамических величин и скорости уравнение имеет ту же форму со сдвигом индекса лз на 1/2): Р«+! Р л+! л+~ +ш +из+ ". «1 "» и где Р и Ц вЂ” функции от основных счетных величин. (Такой вид имеют разностные уравнения в предположении, что уравнения на верхнем слое решены точно.
Фактически же они выполняются с точностью до погрешности итерационного процесса, обычно пренебрежимо малой.) Следствием дивергентности схемы является выполнение разностных аналогов известной интегральной формы уравнений газовой динамики (см. 3 20). Именно она является основой определения обобщенных решений, и это обстоятельство весьма существенно длв расчетов, так как программа предназначалась в первую очередь для решения задач с разрывами. згз ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ [ч. и Разностный аналог интегральных уравнений можно получить, суммируя разностные уравнения по прямоугольной (для простоты) области: и у — 1 г г г ггг Ф1 Рг 7 +1 7г-~-ггг Х г+Иг и+112 + ггпу! «й г т гг М г У ! = Ч' Р"~„„Ь „„— Ч' Р"'~„„й „„+ ~~' О"+',т — ~Х' Ц'"'т =О. (7) Значение подобных аютношеннй обычно обосновывают ссылкой на законы сохранения, разностным аналогом которых они являются.
Это, конечно, справедливо, но мы постараемся привести более четкие соображения. Как уже отмечалось, наиболее важна интегральная форма уравнений для расчета разрывных решений. Рассмотрим (в весьма упрощенной и идеализированной форме) расчет изолированной ударщгй и Л12 ЬГЗ ги Рис. 35 Рнс, 36 волны, Существенным для таких расчетов является определение правильной скорости волны и правильных скачков величин при прохождении фронта волны.
Пусть в расчетах получена следующая характерная картина: графики всех величин выглядят примерно так, как это изображено на рис. 35, где они показаны для моментов времени ги, г„, . г Отметим основные свойства численного решения: есть некоторые значения иы ны е, перед фронтом волны, они постоянны (по Л2 и л) н есть аналогичные величины иг, нг, егза фронтом волны, они тоже постоянны. Наконец, есть «размазанная» волна, точечный график которой за время от гл до 2л просто сместился на несколько узлов.
Кар- ! тина, конечно, идеализирована, но достаточно хорошо отражает свойства численного решения, которое получается по описываемой программе при расчете такого «чистого» течения, как движение первоначально покоящегося газа под действием равномерно движущегося поршня. Преобразуем соотношение (7), используя рис. 36. е22А72изкция глзностиой схвмы я 22! Запишем разностнос интегральное соотношение, учитывая, что на участках 12, 87, 18 величины Р" +и, Д' имеют постоянные значения Рг, Цг, а на участках 34, 45, 56 — значения Р„(1г ОбознаЧаЯ Т=Гн — 1н И 1ы, 1гп ... — ДЛИНЫ СООтВЕтСтВУЮЩИХ УЧаСтКОВ 2 контура, нз (7) получаем 1«трг+ХР 2+иг+ ~мр2 1мрг Х Р '«27г 1збр, + Т(Ц7 — 0г) =О.
72 гз Предполагая, что графики в зонах волн прн г„н 1, одинаковы, 2 т.е. ',7 = ~~', и учитывая, чго 1ы + 127 + 1м = 1вт + 17в + 1вз 1гз = 1зв 22 72 имеем (ИТ)(Р— Р ) — Я вЂ” (е7) = О, гДе г. = 1в7 — 17г = 1м — 1«з — РасстоЯние, пРойДенное счетной УДаР- ной волной, а 17 = 11Т вЂ” ее счетная скорость.
Полученные ссютношения суть соотношения Гюгонио, связывающие состояния до и после фронта волны и скорость волны. Подчеркнем, что этот результат является следствием не только дивергентной формы разносгных уравнений, но и правильного характера численного решения. Последний факт устанавливается экспериментально. Поэтому возникающие при расчетах по некоторым схемам явно нефизические высокочастотные осцилляцин с заметной, но не очень большой амплитудой подрывают доверие к расчету, хотя, конечно, они и не являются безусловным признаком его ошибочности.
Стремление получить решения, свободные от такого «эстетическогв> недостатка входит в цели современной практики конструирования разностных схем. Формула для кинетической энергии. При аппроксимации полной энергии (е+ иг12) бпг используется линейная интерполяция иг, хотя в принципе можно использовать величину (и + и +,)г/4. Это решение оправдывается следующими соображениями. Во-первых, квадрат полусуммы нечувствителен к возмущениям вида ( — 1)"', т.е. сеточные функции и и и + ( — 1) в этом случае имеют одну и ту же кинетическую энергию. Это дает основание ожидать появления в численном решении подобных возмущений, и эксперимент часто подтверждает зтн ожидания.
Вторая причина предпочесть именно интерполяцию квадрата скорости выяснится ниже, Полностью консервативные схемы. Кроме классических и теоретически обоснованных обязательных качеств разностной схемы (аппроксимация н устойчивость), в современной практике вырабо- ззо ЛРИБЛижеииЫЕ МетОды БычислительиОЙ Физики ~ч. и тана система дополнительных требований, улучшающих численное решение в тех нли иных ситуациях. В первую очередь это относится к ситуациям критического характера, когда искомое решение содержит нарушения гладкости, сосредоточенные в сравнительно узких зонах в плоскости (2, х), причем шаги сетки т, Ь не могут быть взяты столь малыми, чтобы упомянутые узкие зоны разрешались достаточно большим числом счетных узлов сетки. Днвергентиость схемы, роль которой мы выше разъяснили, — это одно из таких дополнительных свойств.
Перейдем к Обсуждению другого свойства, получившего название «полная консервативность». Оно связано с некоторыми деталями аппроксимации уравнения для энергии. Начнем с критики построенной выше аппроксимации (5), основанной на прямом использовании уравнения в дивергентной форме. Уравнение (6) для е аппроксимируется проще и компактнее (выбросим пока теплопроводность, не в ней дело): е" +' — е" и"~~ — и"+1 +'" + ( + )"+' +' = О. (8) й т+п2 в -Рт Это уравнение имеет определенное достоинство, хотя мы потеряли днвергентность.
В чем же оно? В различных руководствах по этому поводу можно встретить рассуждение такого типа. Дивергентное уравнение (5) «правильно» описывает эволюцию полной энергии е+ и92, однако распределение приращений отдельных ее видов (внутренней е и кинетической и2/2) может оказаться нарушенным.
Уравнение (8) «правильно» описывает приращение собственно внутренней энергии, но оно, к сожалению, неднвергентно. Из этого естественно возникает задача: построить такую аппроксимацию уравнения энергии, которая была бы и дивергентной, и «правильно» описывала бы изменение е. Такие схемы были построены, их называют полностью консервативными. Попробуем разобраться в этом вопросе. Не случайно слово «правильно» взято в кавычки. Оно не имеет сколько-нибудь определенного смысла. Требование полной консервативности уже вошло в практику конструирования разностных схем, и все более или менее однозначно понимают ее, хотя это свойство, видимо, не имеет четкого определения.
Можно дать такое объяснение. Уравнение для е получается (в дифференциальной форме) из уравнений для е+ и92 и и простой формальной выкладкой. Аналогичную выкладку можно проделать н для разностных аппроксимаций этих уравнений. Результат легко предвидеть: получается разностная аппроксимация уравнения для е типа (8), но с точностью до каких-то членов, пропорциональных т, Ь (с точностью до «аппрокснмационных источников»). Это есть очевидное следствие простою факта: разностное уравнение совпада- РеАаизАция РАзностной схемы зз! ст с дифференциальным с точностью до погрешностей аппроксимации. Схема называется полностью консервативной, если упомянутое формальное преобразование приводит к разностной аппроксимации уравнения для е, не содержащей упоминавшихся выше аппрокснмационных источников н, следовательно, «правильно» описывающей эволюцию е во времени.