Fedorenko-RP-Vvedenie-v-vychislitelnuyu-fiziku (810773), страница 66
Текст из файла (страница 66)
Конечно, значение и„ от и' зависит сильно, но физически существенная величина — тепловой поток ки„ вЂ” прн росте и' быстро выходит на предельное, асимптотическое, значение. й 22. Реализация разностной схемы дпя уравнений газовой динамики с теппопроводностью При создании алгоритма численного интегрирования уравнений газовой динамики возникает необходимость решения большого числа относительно «мелких» непринципиальных вопросов, относящихся, так сказать, к вычислительной технологии. Однако квалифицированное их решение существенным образом влияет па успех дела. В этом параграфе на примере одной конкретной схемы мы постараемся выделить зти вопросы и покажем, на каком уровне онн решаются.
Это, в основном, — уровень качесгвенных соображений, теоретических исследований упрощенных моделей и, конечно же, проверка принятых решений математическим экспериментом. Удобным примером представляется схема, разработанная в 1953— 1954 гг. авторским коллективом под руководством И. М. Гельфанда (это, видимо, была одна из первых схем подобного рода). Выбор этой схемы оправдан еще и тем, что ее реализация затрагивает достаточно полный набор наиболее важных моментов. Математическая постановка задачи. Область расчета.
Решение ищется в прямоугольной области 0 < х < Х, 0 < [ < р, где х— массовая лагранжева переменная, [ — время, т.е. рассчитываются события, происходящие в вьщеленном объеме вещества. Интервал ]О, Х] разбит на части точками 0 = Хв < Х, « ... Х, = Х, причем каждый из интервалов ]Х,, Х,.+,] заполнен газом того или иного сорта. Другими словами, уже геометрия задачи определяет наличие некоторых контактных разрывов (в процессе решения могут появиться и другие) Искомые функции. Расчет состоит в определении функций, описывающих состояние газа: и, г, р, э, е, Т (они имеют физический смысл скорости, эйлеровой координаты, давления, удельного объема, удельной внутренней энергии и температуры).
Из этих функций основными являются и, г, э, Т. Функции р, е связаны с э, Т уравнением состояния, которое имеет свою форму в каждом веществе (т.е. на каждом из интервалов [Х„Х,+,]. Для простоты и определенности можно иметь в виду идеальный газ, параметры которого различны для разных газов, хотя программы в таких ситуациях РЯАяи»Ация РАзностноя схемы згз ~ гг1 обычно пишутся в терминах заданных функций Р,(Т, в), Е,(Т, в), где г — номер вещества. Начальные данные при 1 = 0 задаются значениями и(0, х), г(0, х), в(0, х), Т(0, х). Уравнения. Схема строится на основе уравнений газовой динамики в массовых лагранжевых координатах с добавлением теплопроводности и искусственной вязкости: и, + (р+ ц)„=0, г,= и, в„— и„=0, (е + и92), + [(р+ а) и1„= [ ц(Т, «) Т„[„, (1) где в = (е/в)и„(и„— [и„[) — вязкость Неймана, х(Т, в) — заданный коэффициент теплопроводности.
Для дальнейшего сущесгвенна его следующая форма, явно выделяющая степенную зависимость х от Т: х(Т, в) = Т а(Т, в), где а > 1, а(Т, в) — гладкая функция, ограниченная неравенствами 0 < а «а(Т, в) к а+, Разумеется, параметры и внд функции а зависят от вещества. Ради простоты мы ограничимся «плоским» вариантом задачи, когда г не входит явно в уравнения. Краевые условия. Они могут иметь различную форму. Ради определенности ограничимся такими: при х= 0 заданы скорость и(1, 0) = й(1) и поток энергии кт„= ЩР); при х = х заданы температура т(г, х) и давление р(1, х).
Основные особенности решений. Сложность приближенного решения дифференциальных уравнений определяется прежде всего свойствами гладкости искомых функций, Ниже имеются в виду задачи, решения которых были кусочно-гладкими функциями. Точнее, область счета некоторыми линиями разбивалась на большое число подобластей, в каждой из которых решение было достаточно гладким. Число этих линий и их форма не задавались заранее, они определялись в процессе решения. Линии, на которых нарушалась гладкость решения, являются хорошо известными особенностями решений уравнений газовой динамики и нелинейной теплопроводности, Это ударные волны, границы волн разрежения (линии разрыва производных), фронты тепловых волн и фиксированные в лагранжевых координатах линии разрыва плотности и формул уравнений состояния.
Особенно сложный характер имеет течение в окрестностях точек пересечения линий нарушения гладкости (т.е., например, прохождение ударных и тепловых волн через контактные границы Хи сопровождающиеся «рождением» новых линий нарушения гладкости). Определенные трудности возникают тогда, когда разные подобласти состоят из существенно разных веществ, например если подобласти из очень тяжелых веществ разделяются значительной по эйлеровым размерам областью из очень легкого вещества, имеющей ничтожный 324 ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ [Ч. !! в массовых координатах размер (такую область условно назовем «вакуумом») .
Сложным является, например, прохождение ударной волны через вакуум. При выходе ударной волны на внутреннюю границу вакуума, она исчезает, сменяясь волной разрежения. Одновременно начинается быстрое движение этой границы вакуума в сторону другой его границы, В какой-то момент эти границы встречаются, происходит «удар», снова рождающий ударную волну. Интересный класс течений создается в следующей ситуации. На границу холодного покоящегося газа подается мощный поток энергии (задается либо большой поток иа одной нз границ, либо высокая температура).
Возникает характерная картина — прогрев газа в режиме тепловой волны, фронт которой движется с конечной скоростью. Расчет такого режима затруднен тем, что температура существенно негладкая около точки фронта (см. $21). Градиент температуры порождает градиент давления, и в дальнейшем возможно образование ударной волны, причем могут осуществиться два разных предельных режима: либо ударная волна обгоняет тепловую, либо ударная волна отстает от тепловой и распространяется как изотермическая по сильно нагретому веществу.
Такого рода процессы протекают при облучении сферических мишеней мощным потоком лазерного излучения. Разностная аппроксимация задачи. Введем основные обьекты, появляющиеся при конструировании метода приближенного решения. Сетка и счетные величины. Интервал (О, Х] покрывается сеткой (х )"', сетка (Б„) формируется в процессе решения, так как шаг по времени т„+и выбирается в зависимости от полученного на л-м временном шаге решения. Узлы сетки с координатами (л, т) образуют множество «целых» счетных точек.
В них определены «механические» величины и«н Р". Кроме того, вводятся «полуцелые» счетные точки с координатами х „=(х„+х +,)/2, Р„. В этих «полуцелых» точках определены «термодинамические» величины (т"„„з, н"„»ц2) (Рл О, 1, ..., м). Границы подобластей х! совпадают с какими-то из «целых» точек х . Особое положение занимают граничные точки. В точках (О, н) (левая граница), кроме механических величин, могут быть определены некоторые термодинамические, используемыедля реализации краевого условия н для нестандартной аппроксимации некотормх уравнений.
В точках (лз + 1/2„н) (правая граница) могут быть определены механические величиим и«м«!~, Рм+цм использУемые в тех же целЯх. Мы таких нестандартных счетных точек использовать не будем (нужда в них появляется, например, при иных краевых условиях). РеАлизАция РАзиостиой схн«ы й 221 (стандартная аппроксимация применима при т = 1, 2, ..., М вЂ” 1), „л»1 л »+1+ л з 2 (формула применима при всех т = О, 1,...,М). Уравнения для удельного обьема: Р"+1 — Рл л+1»1 "+1П "»1П "+1 — 0 " +мз (4) (формула имеет смысл при т = О, 1, ..., М вЂ” 1). Размостмая аппроксимация в стамдартмых точкак. Сначала опишем стандартные формулы разностной аппроксимации, т.е.
те, в которых не используются «термодинамические» величины в граничных узлах. Введем следующие обозначения: Ь +, —— х +, — х, й = х,)2 — х,)2 — шаги сетки. Численное интегрирование проводим по стандартной для эволюционных задач схеме счета «по слоям». Шаг интегрирования состоит в том, что значения на и-м слое (и, г, Т, в)" уже известны и надо вычислить (л + 1)-й слой (и, г, Т, в) "+', решая систему уравнений на верхнем слое. Ради простоты рассмотрим полностью неявную схему, хотя можно использовать и схему, в которой пространственные производные аппроксимируются взвешенной (обычно с весами 0.55 и 0.45) суммой аппроксимаций на верхнем и нижнем слоях. В любом случае приходим к системе нелинейных уравнений относительно неизвестных величин на верхнем слое, которая решается специальным итерационным алгоритмом.
Итерации строятся на основе неполного метода Ньютона. В уравнения на верхнем слое входят близкие величины трех видов. Поясним это на примере Т +,1; (то 'ке самое относится и к Ии»112, и ). Мы имеем дело с Тл»1пи Тф,112, Т<'+11))2. Величины Т", Пз (с нижнего слоя) считаются уже известными. При вычислении Т"+.,', 2 методом итераций фигурирукп уже найденные значения Т1'1»пз (1-е приближение к Т"„,.',и) и неизвестные значения Тф) .
В пределе величины Т~') 1)2 — Т".,1вз. Мы используем обозначения Т пз — — Т +„„Т лп,— — Т »112. ИменнопоотиошениюкнеизвестО) О+1) ным Т пз производится линеаризация при выполнении очередной итерации. Аппроксимация уравнений движения: ил+1 — ил 1р+д)л+1 — 1р+д)"»1 + Р Я "+гп Р ~ -' О (2) Я 326 ПЦИЕЛИЖЕИИЫЕ МЕТОДЫ ЦЫЧИ!'ЛИТЕЛЬНОЙ ФИЗИКИ (ч.