Федоренко Р.П. Введение в вычислительную физику (1185915), страница 63
Текст из файла (страница 63)
[1 н пгивлижзииыз мвтоды вычислитвльиоя»ивики Расчет ударных волн основан иа предложенном фон-Нейман м и Рихтмайером приеме — на введении в уравнения искусстве)1ной вязкости. Она вводится таким образом, чтобы уравнения мало йскажались вне зоны ударной волны. Сама же ударная волна при.этом «размазывается» на пренебрежимо малую ширину. Таким образом мы заменяем уравнения газовой динамики на слабо возмущенные, но уже не имеющие разрывных решений уравнения.
Конкретно, новые уравнения имеют вид «,+(р+я)„=О...— «.=а, (21) где д — малая искусственная вязкость. Нейман и Рихтмайер предложили очень удобную конструкцию (она наиболее популярна и чаще других используется в расчетах): « = Й («. — 1«.!)«; (22) Здесь в — малый параметр, выбор которою мы в дальнейшем уточним. Очевидно, о вь О, если и„с О.
Это условие выделяет те участки течения, на которых происходит сжатие вещества (плотность растет: из и„< 0 следует и, с О, т.е. р, > 0). Наоборот, там, где течение сопровождается «разрежением» (и„> О, я = 0), вязкость «выключается». Ударная волна — зто как раз участок течения, на котором происходит ударное, скачкообразное, сжатие вещества.
Достоинства неймановской искусственной вязкости проще всего продемонстрировать, сравнив точные решения простых модельных задач для исходной системы уравнений (14) и для уравнений с вязкостью (21). Эти простые решения относятся к важному классу автомодельных решений, в которых характерные особенности реальных решений (в данном случае, ударные волны) проявляются, так сказать, «в чистом виде», без взаимодействия с какнми-то гладкими течениями. Забегая вперед, опишем результат. Оказывается, решением исходной системы будет ударная волна — «ступенька», движущаяся с постоянной скоростью, а решением системы уравнений с искусственной вязкостью (21) будет движущаяся с той же скоростью «размазанная ступенькав, причем ширина зоны размазывания фиксирована, она зависит, разумеется, от величины в.
Вие зоны «размазанной ударной волны» все функции и, и, р, е для обеих систем уравнений совпадают. Итак, рассмотрим следующую ситуацию. Пусть имеется решение уравнений (14) типа чисюй ударной волны, т.е. ~«п ип е„х> й1, одномввныв гелвнвния глзовой диялмякн зог й зо1 где Р— скорость ударной волны, и выполнены соотношения Гюгонио (в массовых лагранжевых переменных): а) — Ри, + р, =С,= — Ри + рч, (23) б) Ри~+ и~ Сг Риг+ иг г1 «г) в) — Р е, + —.'! + Р,и| = Сз = — Р ~ег+ г'! + Р,и,. Построим аналогичное решение для уравнений с искусственной вязкостью.
Это решение будем искать в классе автомодельных решений типа бегущей волны», т.е. когда все функции и, е, е зависят от одного аргумента $ = х — Рк Уравнения в частных производных (14) превращаются в обыкновенные после замены операторов д/д/ = — Р ЫИ$, д/дх = г/И$. Выписывая эту систему и интегрируя ее один раз, мы получаем систему «первых интегралов»: — Ри + (р+ д) =О ~ — Ри+ р+ е= С,, (23*) — Ре — и =О~Ри+и=С г' „гг / г'~ — Р е+ г~ + ~(р+ ч)и1 =О ~ -Р~е+ — ) + и(р+ ч) = Сз Так как мы не решаем какую-то определенную задачу, а просто конструируем нужное нам решение, постоянными интегрирования можем распоряжаться так, как нам будет удобно, В частности, здесь Сы Сг, Сз — те же самые, что и в соотношениях Гюгонио (23).
Дальнейшие выкладки будем проводить для идеального газа: е = ри/(у — 1). Руководящей идеей последующих преобразований является стремление получить уравнение только для и(г,), остальные переменные будем исключать через зл Следующие ниже выкладки оправданы лишь при и„( О, т.е. при »1 > О. В этом случае е = «Рг(е )г/в. Из полученных «первых интегралов» имеем и= Сг — Ри, (р+ 4) = С4 Р~г, р= С« — Рги — д. Точные значения постоянных С н других (А, В; см.
ниже) нас пока не интересуют. Из «интеграла энергии» для идеального газа — Р— ', + —, + (Р+ ч) и = Сз исключим р, и, (р+ е) по уже найденным формулам. После простых преобразований получаем уравнение для сч дг = — т Рггг+ Ае+ В. (24) 1ч1 и пгиелнженные методы вычислительной ьизики эоз Утверждение. Многочлен в правой части (24) обращается в нуль при в = и, и и = иг; поэтому он может быть записан в виде г Р (" и~)(" "г)' (25) Доказательство. Соотношение (24) является следствием соотношений (23«). Покажем, что полагая в этих соотношениях в = ип можно получить в качестве следствий равенства р= р„и = и„у = О.
В самом деле, из Рн + и = С = Ри, + и, прн г = и, следует, что и = и,. Из -Ри + р+ ~у С, = — Ри, + р, при и = и, имеем р+ а= р,. Из 21 — Р ~ ~— + — ~ + ( р + 4) и = С = — Р ~ ~— + — ~ + р ре и ~ в1«~ — г~ — з= )т — ~ г~ при и= но и = ио р+ а= р, получаем р= р, и, следовательно, я = О. Таким образом, при и = и, левая, а следовательно, и правая части (24) обращаются в нуль. Точно так же рассматривается и случай и =вг.
Используя выражение для д, получаем дифференциальное ураве(иг)г = 0.5(у + 1) (и — э,)(иг — в) (26) (и условие гг> 0). После замены переменных г =Ле/('у+Т~т1 и и = (и, + иг)/2+ г(и, — иг)/2 уравнение (26) принимает вид (г)г 1 г Решение угадывается: г(т)) = »1 нли е(п) = з(п и. Последним решением можно пользоваться (в силу условия иг>0, т.е.
г„> 0) лишь при Ч ~ ( — н/2, х/2). Возвращаясь к прежним переменным, получаем решение (продолжая его постоянным за пределами выделенного интервала й): ь=х — Рус — т —— 2е к у+1 2' зш 1/~(х — Рг), )$~ с Ъ— 2е т+1 г' ~/ 2е и т+1 г' ь,+ь ь,— е + г г сп Такие же формулы легко получить и для функций и(з) и р(е) + уД) с заменой ом иг на ип иг или р,, ргаютветственно. Итак, получено решение уравнений с вязкостью, совпадающее с решением типа «чистая ударная волна» всюду, за исключением узкой полосы вдоль фронта ударной волны ~х — Щ и(я/2)т/2«7(у+ 1). Ширина «размазанной ударной волны» есть щ/2«/((т + 1) .
Опыт показал; что хорошие результаты дает выбор е, при котором волна разре- 1ш1 одномгшьи гглвнвния глзовоя динАмихя шается четырьмя-пятью счетными точками. Например, 5Ь = пЛ з7(у + Т~, т е. з = 256з(у + 1)/(2 из) ш 26з. ' Отметим, что можно выписать формулы для р(х) в зоне волны. Детали нам не очень нужны, отметим лишь, что фактически область плавного перехода от р, к рь примерно в два раза уже, чем для остальных функций и, и, р+ д.
При расчетах влагранжевых координатах по хдифференцируются только функции и и р+ я, каждая из которых имеет стандартную ширину зоны размазывания. Функции р и Э отдельно по х не дифференцируются. Поэтому сокращение фактической ширины волны для р не играет роли.
Правда, мы должны еще обеспечить должное размазывание фронта волны на четыре-пять точек сетки по времени. Практика расчетов, проводившихся в пятидесятых годах привела к такому рецепту. Ударная волна за один шаг времени должна проходить примерно половину интервала по массовой координате, т,е. по времени зона размазанной волны захватывает примерно в два раза больше счетных интервалов. Обратим внимание на то, что по г дифференцируютсз все функции. Таким образом, и более крутой график р также «разрешаетсяь четырьмя-пятью точками сетки по 1. Этот рецепт (половнна шага по пространству за шаг по времени) не был связан с условием устойчивости, так как.мы применяли неявные абсолютно устойчивые схемы.
Он был связан с тем, что для необходимой точности разностной аппроксимации нужно разметить на крутых профилях функции четыре-пять счетных точек сетки (как по времени, так и по пространству). Попытки расчетов с ббльшими шагами по г (расчет, например, со скоростью один шаг по пространству за шаг по времени) приводили к ухудшению результатов: на графиках появлялись осцилляции явно счетного происхождения. Заметим, что все проведенные выкладки можно повторить и в эйлеровых координатах.
Однако ситуация осложняется тем, что в эйлеровой форме уравнений по х дифференцируются все функции. Поэтому приходится брать в два раза более широкую зону размазывания (по х), чтобы фактическая ширина зоны размазывания р была покрыта четырьмя-пятью интервалами сетки. Это уже не очень приятно, так остальные величины при этом размазываются на десять точек. В 1962 г.
автором проводился расчет ударной волны в эйлеровых координатах. Чтобы избежать слишком широкого счетного фронта волны, была вмбрана довольно экстравагантная форма записи уравнений: в качестве основных функций были взяты и, р, с где с — скорость звука. Уравнения (в недивергентной форме) оказались такими, что по х дифференцировались только функции, профили которых имели в зоне волны стандартный синусоидальный вид. Это позволило вести расчеты с шириной размазывания порядка 4Ь. (Подробнее об этом см. з 22 в связи с применением гибридной схемы для решения уравнений газовой динамики в эйлеровых координатах,) 31О неивлижвнныв методы вычислительной»заики ~ч. и й 21. Нелинейное уравнение теппопроводности Рассмотрим некоторые вопросы, возникающие прн численном решении. нелинейного уравнения теплопроводностн, ад ах н(х,г,и)ах +Д(х,г,и), (1) которое решается в простой области 0 и а и Т, 0 и х и Х с начальнымн данными и(0, х) = ив(х) и краевыми условиями, для простоты, первого рода: и(Ф, 0) и(а, Х) = О.
Нужно иметь в виду, что нелинейность уравнений приводит к сложностям не сама по себе, а лишь в тех случаях, когда она порождает сложные, описываемые негладкими функциями, явления. Чтобы сделать это более конкретным, рассмотрим важный в приложениях случай, когда коэффициент теплопроводности и зависит от и степенным образом: и, = [и«и„[„. (2) Уравнения такого типа встречаются при описании процессов в высокотемпературном веществе (лучистая теплопроводность), например в звездах. Аналогичные уравнения описывают н процессы фильтрации.
дз Рассмотрим характерное и очень важное в при« ложеннях явление, описываемое зтнм уравнеии1 ем, — так называемую тепловую волну, илн, иначе, тепловой фронт. На рис. 33 изображены графих кн функций и(аи х) для трех моментов времени г, < гз < га. В этом случае мы имеем процесс расРис.
33 пространения высокой температуры по «нулевому фону» (перед фронтом тепловой волны и = 0; в действительности, конечно, перед фронтом температура не нулевая, но очень маленькая по сравнению с температурой за фронтом). Тепловой фронт, Будем искать «автомодельное» решение уравнения, т.е. решение, зависшцее не от т и х, а от нх комбинации, в данном случае от $ = х — Вг, где аа — некоторая постоянная, смысл которой потом станет ясным. Тогда ди Ии д$ «и ди Ии — = — — = —  —, да Н$ дд а3' дх Ы$ Нам следует найти решение обыкновенного дифференциального уравнения — ааи'= [и«и'['.
Интегрируя, получаем -ааи = и«и'. (Так как мы ищем какое-нибудь решение, постоянную интегрирования положим равной нулю.) Уравнение и" 'и' = — Р, нли (и«)' = — хР, интегрируется и дает и«(С) = -йЩ. 311 нелинейное лАвненне теплопговодностн дгЦ Итак, мы получили решение и($) = ФО(-$) и». Но нас интересует вещественное и положительное и.
Поэтому это решение имеет смысл только при к ж О, т.е. при х <1М. Положим иД) = О при р ) О. Легко видеть, что функция и($) — О является решением. Вопрос только в том, можно ли эти два решения склеить, точнее говоря, можно ли н в каком смысле говорить, что функция (3) есть решение уравнения теплопроводносги. Естественно обратиться к понятию обобщенного решения, так как в любой точке (1, х), кроме линии х = Юг (фронт тепловой волны), эта функция удовлетворяет уравнению теплопроводности в классическом смысле. Обобщенное же решение вводится как функция, удовлетворяющая интегральному тождеству (закону сохранения).