Г. П. Прокопов Контроль энтропии в алгоритмах и расчетах газодинамических течений Ордена Ленина Институт прикладной математики имени М.В.Келдыша Российской академии наук Г.П. Прокопов Контроль энтропии в алгоритмах и расчетах газодинамических течений Москва, 2006 год УДК 519.6 Контроль энтропии в алгоритмах и расчетах газодинамических течений Прокопов Г.П. Препринт Института прикладной математики им. М.В.Келдыша РАН На примере одной недавно опубликованной схемы проведено аналитическое исследование поведения энтропии при численном интегрировании уравнений газовой динамики. Оно подтвердило опасение о возможности реализации в ходе расчета разрывов, которые аналогичны ударным волнам разрежения. Следовательно, может происходить нарушение постулата о неубывании энтропии. Предлагаются простые алгоритмы контроля энтропии при расчетах газодинамических течений и по другим известным численным методам. Работа выполнена при поддержке Российского Фонда Фундаментальных Исследований (грант № 05-01-00097). Control of entropy in algorithms and calculations of gas-dynamics flows. Prokopov G.P. Preprint of Keldysh Institute of applied mathematics, RAS Based on one recently published approach, an analytical investigation of entropy behaviour in process of numerical integration of gas-dynamics equations is implemented. It confirms an existence of discontinuous solutions that similar to rarefaction shock-wave in numerical calculations. So, the postulate of non-increasing entropy is disturbed. Some simple methods to control of entropy in gas-dynamic flows and other known numerical methods are suggested. This work is support by RFFI (grant N 05-01 -00097) Содержание стр. Введение ………………………………………………………….. 3 § 1. Схема С для расчета газодинамических течений …………. 4 § 2. Энтропийное исследование схемы С ….. ………………….. 7 § 3. О схеме С и ударной волне разрежения……………………. 10 § 4. Нужен ли энтропийный контроль схемам типа Годунова? . 12 § 5. Как быть с другими схемами?…………………….………… 13 § 6. Снова о схеме C и ее экономичности ...……………………. 15 § 7. Предложение по улучшению схемы С ………….…………. 19 § 8. О проблеме сложных уравнений состояния ……………….. 23 Заключение ………………………………………………………. 26 Литература ………………………………………………………... 28 Введение Настоящая работа является непосредственным продолжением краткой публикации автора [1]. Одной из ее целей было напоминание о том важном обстоятельстве, что в широко известном методе Годунова (см., напр., монографию [2]) успех достигается благодаря использованию дополнительного закона сохранения энтропии на гладких решениях и постулату о ее неубывании на разрывах. В связи с тем, что имеется много публикаций различных авторов, заявляющих о методах «типа Годунова», интересен вопрос, с должным ли вниманием они относятся к упомянутому обстоятельству. Естественно, что это касается и не только таких методов, но и вообще алгоритмов для расчета газодинамических течений. Проясняя для себя этот вопрос, автор обратился к изучению численных методов с такой точки зрения. Весьма полезной оказалась монография [3], в которой можно найти как содержательную информацию о работах многих авторов, так и ссылки, по которым можно разыскивать их публикации. В первую очередь автора интересовали методы расчета газодинамических течений (как нестационарных, так и стационарных), обеспечивающие выполнение законов сохранения массы, импульса и энергии, а потому основанные на использовании дивергентной формы уравнений газовой динамики. При ознакомлении с таким кругом работ удалось обнаружить слишком мало информации, касающейся поведения энтропии. Исключение составляет разве что раздел 2.10 в монографии [3], что и будет отмечено в §5. Попытавшись выполнить такие исследования аналитически, автор потерпел полную неудачу из-за весьма сложного характера известных и широко используемых алгоритмов. И вдруг в этом «темном царстве» блеснул луч надежды. На семинаре им.К.И.Бабенко в ИПМ им.М.В.Келдыша 26 октября 2006 г. был сделан доклад по результатам работы [4]. Эта очень небольшая по объему публикация производит впечатление на специалиста, знающего в деталях содержание вопроса, чрезвычайной простотой предлагаемого алгоритма для приближенного решения задачи Римана о распаде произвольного разрыва. А она является главным элементом при реализации численного метода С.К.Годунова, а также его усовершенствований и модификаций, предлагаемых различными авторами. Правда, следует отметить, что представленный в [4] разностный метод для уравнений газодинамики, для краткости именуемый далее схемой С (название взято из [4]), сразу же вызвал определенные сомнения . Именно в силу упомянутой его простоты по сравнению с другими методиками удалось провести аналитические исследования по интересующему вопросу о поведении энтропии, которые, к сожалению, подтвердили справедливость возникших опасений. Они будут представлены в настоящей работе. Эти результаты позволяют автору выступить с предложением об «энтропийном контроле» , которым следует сопровождать расчеты газодинамических течений, выполняемых по методам «типа Годунова» , схеме С и другим аналогичным методам. В § 5 автор рассматривает возможность использовать «энтропийный контроль» и при расчетах по другим методикам. Например, по представленным в монографии [3] сеточно-характеристическим схемам типа КИР (Куранта-Изаксона-Рисса), Роу, TVD-схемам и т.п. Основная часть работы для простоты изложена для идеального газа. Однако идея энтропийного контроля тем более необходима и может быть реализована и для сложных уравнений состояния . Но есть и проблемы, связанные с их заданием. Этим вопросам посвящен § 8. § 1. Схема С для расчета газодинамических течений. В настоящем параграфе будет практически полностью представлено описание разностной схемы из работы [4]. Сделать это целесообразно, во-первых, потому, что все приводимые формулы потребуются для дальнейшего изложения. Во-вторых, сочтено необходимым изменить некоторые обозначения, с тем, чтобы обеспечить единообразие с монографией [2], известной широкому кругу читателей, и которая тоже будет использоваться в дальнейшем изложении. Наконец, в-третьих, к счастью, нужное нам описание сделано в [4] с удивительной краткостью, занимая всего одну печатную страницу. Описание схемы С рассматривается на примере хорошо известных нестационарных уравнений одномерной газовой динамики в дивергентной форме: (1.1) , где U , F – векторы для законов сохранения массы, импульса и энергии: (1.2) , Здесь r - плотность газа, u – скорость, р – давление, (1.3) - полная внутренняя энергия, e - внутренняя энергия единицы массы, определяемая уравнением состояния газа. Пока, до § 8, ограничимся простейшим уравнением состояния идеального газа с показателем адиабаты g: (1.4) . Разностная схема для (1.1) записывается в виде: (1.5) Здесь, как было принято в монографии [2], - номер ячейки сетки по оси х , ограниченной узлами с номерами j -1 и j . Формула (1.5) описывает «пересчет» величин на одном шаге по времени – переход от величин с «нижнего» слоя по времени t к «верхнему». В методе С.К.Годунова потоки на границах ячеек F j вычисляются посредством численного решения задач Римана о распаде разрыва с параметрами газа в соседних ячейках сетки. Этот расчет делается либо приближенно (например, в виде «звукового» приближения), либо (если необходимо) итерационным методом до достижения результата с предписываемой точностью. Основная идея предложенной в [4] схемы С также состоит в вычислении потоков на границах ячейки. Задача Римана решается приближенно на основе соотношений на разрывах, аналогичных описанным в [2] на стр.103: (1.6) Здесь D – скорость распространения разрыва, а – массовая скорость (вместо обозначений w,m соответственно в [4]). Квадратными скобками [ ] обозначается разность значений величин, заключенных в скобки, по обе стороны от разрыва. В отличие от метода Годунова, при приближенном решении задачи Римана рассматривается упрощенная схема течения с распадом на левую волну, контактный разрыв и правую волну. Она изображена на рис.1. Обозначим индексами: 1 – параметры в левой ячейке сетки, 2 – в правой, 3 – между левой волной и контактным разрывом, 4 – между контакт-ным разрывом и правой волной, как показано на схеме течения (см. рис.1). Рис. 1. Рис. 2. Расчетные формулы имеют вид: , , (1.7) Назначение параметров газа ( R , U , P , E ) j на границах ячейки для вычисления векторов потока F j (в методе Годунова они назывались «большими» величинами) определяется значениями U , D 1 , D 2 . В случае U ³ 0: при D 1 ³ 0 эти параметры равны параметрам в зоне 1, а при D 1 < 0 – параметрам в зоне 3. В случае U < 0: при D 2 £ 0 – параметры равны их значениям в зоне 2, а при D 2 > 0 – параметрам в зоне 4. Определяющую роль играет назначение массовых скоростей а1 , а2 . Для этой цели предлагается использовать формулы: (1.8) Присутствующие в них величины с1 ,с2 представляют скорости звука, определяемые уравнением состояния. Для идеального газа (1.4) они вычисляются по формулам: (1.9) , По утверждению [4], такое назначение массовых скоростей а1 , а2 обеспечивает отсутствие осцилляций на разрывах при произвольных физических параметрах в соседних ячейках разностной сетки. Пока ограничимся этим описанием схемы С и перейдем к ее исследованию. § 2. Энтропийное исследование схемы С . Итак, одним из принципиальных элементов конструирования схемы С является замена волны разрежения на фронт разрыва . Напомним, что, как отмечалось в [1], в монографии [2] на стр.115-116 была рассмотрена в качестве примера задача о распаде разрыва с симметричными данными относительно границы раздела. Если в ней заменить возникающую волну разрежения фронтом разрыва , на котором будут выполняться соотношения (1.6) , то при этом реализуется так называемая ударная волна разрежения . А ее следует запрещать как нарушающую постулат о неубывании энтропии . В схеме С сначала «волевым» образом по формулам (1.8) назначаются массовые скорости а1 , а2 , а затем, исходя из тех же соотношений на разрыве, производится вычисление параметров, характеризующих состояние газа на разрыве. Поскольку на этом процесс их назначения заканчивается , говорить о реализации ударного фронта разрежения нет оснований. Но в каком случае можно полученный результат считать разумным и допустимым? Предлагается постулировать следующее положение : газ с вычисленными параметрами , определяющими величины потоков через разрыв, должен иметь энтропию, не меньшую, чем второй (исходный, «напарник») газ, участвующий в вычислении потоков . Представляется, что только такое требование позволит контролировать неоправданное и скрытое занижение энтропии в газодинамических расчетах. Назовем его «жестким контролем энтропии». Более просто было бы, например, ограничиться вычислением и контролем энтропии в ячейках сетки на верхнем слое по времени. Однако тогда, в соответствии с формулой (1.5), появляется возможность компенсировать «потери энтропии» в одном узле сетки за счет другого. Поэтому назовем его «мягким контролем энтропии». Нарушение его создает возможность конструирования, по крайней мере, «патологических» примеров недопустимых течений. Но отсюда недалеко и до возможности реализации таких ситуаций и в практических расчетах. Если согласиться с провозглашенным постулатом, то следующим шагом исследования должен быть вопрос: обеспечивает ли назначение массовых скоростей (1.8) его выполнение для рассматриваемой схемы С ? Это исследование будем проводить на уже упомянутом упрощенном примере распада разрыва с симметричными данными. Итак, пусть: (2.1) , , , . Поскольку два газодинамических параметра можно задавать независимо, можно было бы упростить формулы, полагая , . Мы сознательно не будем этого делать, чтобы не потерять возможности контролировать их с точки зрения размерности. Тогда в задаче остается один свободный параметр u * (если не считать показателя адиабаты g в уравнении состояния (1.4), которым мы пока ограничиваемся). Отметим, что в качестве аналитического решения этой газодинамической задачи при u * >0 (встречные потоки) реализуются ударные волны, а при u * <0 (разбегающиеся потоки) – волны разрежения. Они распространяются симметрично от границы раздела. Благодаря симметрии (2.1) расчетные формулы (1.7)-(1.9) принимают более простой вид: (2.2) , (2.3) (2.4) (2.5) В соответствии с назначением, описанным вслед за формулами (1.7), для «больших» величин в узле, которые обозначим R , E , получаем: (2.6) Вот теперь приступаем к исследованию энтропии. Обычно для этого привлекается энтропийная функция . Практически будет удобнее работать с величиной . Заметим, что именно эта величина, с точностью до множителя с v (удельная теплоемкость при постоянном удельном объеме), принималась за «настоящую» энтропию (см. [5], стр.33). Она остается постоянной вдоль адиабаты Пуассона. Поэтому, чтобы не вводить новых обозначений, полагаем, что (2.7) , Условие ее неубывания принимает вид: (2.8) Следовательно, из формул (2.4),(2.6) и (2.2) получаем: (2.9) Введя безразмерный параметр и с учетом (2.3) получаем: (2.10) Проанализируем полученную формулу сначала для случая m < 0. Тогда и формула (2.10) приобретает вид: (2.11) . Ее производная определяется формулой: (2.12) В точке m =0 имеем . Производная на интервале . Для значений функция не определена . имеет вертикальную асимптоту при . Стоит заметить, что в появлении «виновата», конечно, формула для давления (2.4): при получалось бы Р< 0. Теперь обратимся к случаю m>0. Поскольку тогда , формула (2.10) приобретает вид: (2.13) . Ее производная определяется формулой: (2.14) После несложных преобразований эта формула приводится к виду: (2.15) Функция определена для всех . . Квадратный трехчлен в числителе (2.15) положителен для g < 7. Поэтому для всех m >0, если g < 7. В случае g>7 квадратный трехчлен достигает минимума при и этот минимум равен , т.е. отрицателен. Следовательно, при g>7 есть некоторый интервал, на котором . Поскольку , получается, что и на некотором интервале , образуя «провал» ниже оси m . Эскиз графика функции - сплошная кривая на рис.2 для g < 7. При m =0 он касается оси m и при имеет асимптотику . Главным (и тревожным) фактом является поведение функции f ( m ) при m < 0. Она отрицательна на интервале , а при вообще не определена. Это означает, что в случае разбегающихся потоков в расчете распада разрыва получатся параметры потока, имеющего меньшую энтропию , т.е. будет нарушен сформули-рованный постулат . Налицо нарушение жесткого контроля энтропии. Поскольку опасный участок примыкает к m =0, такая ситуация реализуется при сколь угодно малых различиях параметров с индексами 1 и 2. В случае g>7 аналогичная ситуация наблюдается и для встречных потоков, если m попадает внутрь «провала» (см. рис.2). § 3. О схеме С и ударной волне разрежения Полученному результату не приходится удивляться. Обратимся к уже упомянутому в начале § 2 примеру на стр.115-116 монографии [2]. В нем возникала ударная волна разрежения . Получим для нее соответствующую функцию f ( m ) , аналогичную (2.10). Ограничимся только ее частью, отвечающей m < 0. Если вместо m временно принять за параметр величину , то, с учетом описания в [2] на стр.116 и определения энтропии (2.8), можно сразу выписать формулу: (3.1) Для ее производной после несколько кропотливого упражнения получается красивый результат, приведенный на стр.116 в [2]: (3.2) Нам не хватает только формулы, связывающей m и q . Вся информация для этого содержится в формулах (2.2)-(2.5): (3.3) В результате получаем: (3.4) Поскольку , для ее производной будем иметь: (3.5) . Отсюда следуют те выводы, которые уже приведены в [2] на стр.116-117. При малых значениях m получаем, что (3.6) , , где . Следовательно, в ударной волне слабой интенсивности скачок энтропии является малой величиной третьего порядка малости по сравнению со скачком давления. Для функция f ( m ) не определена. Поведение f ( m ) на интервале аналогично изображенному на рис.2. То тревожное обстоятельство, что f ( m )< 0 для , означает нарушение там упомянутого выше постулата , т.е. реализацию ударной волны разрежения , если скорость разбегания газов попадает на этот участок. В какой-то степени утешительным можно считать тот факт, что (в силу касания третьего порядка) величина остается выше некоторого условного порога - d для значений m на протяженном отрезке. Длина его оценивается, исходя из формулы (3.6), как . Аналогичная оценка для примера, рассмотренного в § 2, с использованием формулы (2.12) получается несколько хуже. Поскольку при малых m имеем оценку: , протяженность аналогичного отрезка (3.7) . Полученные результаты свидетельствуют о том, что при m <0 схема С реализует разрыв , на котором уменьшается энтропия, т.е. по существу похожий на ударную волну разрежения. § 4. Нужен ли энтропийный контроль схемам типа Годунова? Обратимся теперь к вопросу о том, как срабатывает энтропийный контроль для схемы Годунова . И сразу обнаруживаем, что в ней нет обсуждаемой проблемы , потому что «честно» отслеживаются границы области, занимаемой волной разрежения, в которой обеспечивается закон сохранения энтропии. Причем это обеспечивается не только в окончательном решении, но и в ходе всего итерационного процесса. Оборвать его можно в любой момент, который выберет исполнитель. Достигается это следующим образом. Как описано в монографии [2] на стр.114-115, получается очередное значение для давления . Далее в случае, если в конфигурации распада разрыва присутствует волна разрежения, значения плотности R и остальных величин, описывающих эту волну, досчитываются, исходя из адиабаты Пуассона , что тождественно сохранению значения энтропии. Даже при расчете исходного «звукового» приближения, получив начальное приближение для давления по формуле (13.26) на стр.113 в [2], не стоит «соблазняться» расчетом и других величин по аналогичным линеаризованным формулам, а следует воспользоваться только что упомянутым алгоритмом с учетом адиабаты Пуассона. «Копеечная» экономия от замены их ударной адиабатой Гюгонио (или линеаризованными формулами) может обернуться ударной волной разрежения, интенсивность которой и соответствующее уменьшение энтропии еще нужно оценивать. Да и нужно ли? А описанный «честный» досчет величин даже для грубого «звукового» приближения гарантирует неубывание энтропии. Изложенное, конечно же, не означает, что у метода Годунова вообще нет проблем. Ввиду первого порядка схемы к ним относятся прежде всего проблемы относительно низкого уровня точности, что особенно существенно при расчете гладких течений. Именно поэтому усилия многих авторов направлены на повышение порядка аппроксимации и другие его усовершенствования для повышения точности расчетов. Одно из направлений деятельности такого рода связано с заменой кусочно-постоянной аппроксимации текущего распределения газодинамичес-ких параметров, используемой в исходном варианте метода Годунова. Идея – подать задаче о распаде разрыва более совершенные данные, благодаря которым ожидается повышение точности. Работу такого рода алгоритмов можно условно разделить на этапы, выделив в отдельный подготовку таких данных, а в следующий этап – расчет распада разрыва. Обсуждаемый в этой работе «жесткий энтропийный контроль» предлагается применить к подготовленным входным данным для расчета распада разрыва. Цель его – убедиться, что не «нахимичили» и с энтропией все в порядке. А уж алгоритм расчета распада разрыва не подведет, если применять так, как было описано выше. Сказанным автор ни в какой мере не пытается принизить роль этапа подготовки совершенных данных. Напротив, по-видимому, именно это – наиболее перспективный путь развития алгоритмов такого типа. Теперь самое подходящее время отметить, что в качестве заключительного этапа выполняется «пересчет» значений величин в ячейке сетки с учетом потоков через ее границы. Интуитивно ясно, что это эквивалентно «смешиванию» нескольких газов, каждый из которых имеет свою энтропию. Если все новые компоненты имеют энтропию, не меньшую, чем исходный газ, неубывание энтропии для результата «смешивания» гарантировано . Во всяком случае, ограничение шага по времени по условию Куранта направлено именно на это. Сказанное означает, что «мягкий контроль энтропии» в ячейках сетки становится тогда излишним. Такой контроль однако не теряет своей актуальности, если «подмешиваются» компоненты с низкой энтропией . Очевидна и недостаточность только «мягкого» контроля . Результат контроля может быть благоприятным, хотя «кое-где у нас порой» происходило убывание энтропии , но это удалось скрыть за счет соседних узлов сетки. § 5. Как быть с другими схемами? Есть еще широкий круг алгоритмов расчета газодинамических течений, реализующих законы сохранения в форме (1.5). При этом величины потоков F j вычисляются, исходя совсем из других соображений, чем в рассмотренных выше методах «типа Годунова». Сюда, в частности, можно отнести многие схемы, рассмотренные в главах 2-3 монографии [3]. Не будем повторять и расширять их перечисление, сделанное во введении. В разделе 2.10, который называется «Энтропийная коррекция» и занимает всего лишь стр.139-144, специально рассматриваются «алгоритмы, которые позволяют избежать появления в численных результатах нефизических решений, в частности, газодинамических ударных волн разрежения. Такая коррекция вводит дополнительный механизм отбора физически приемлемого решения» (см. [3], стр.139). Ввиду того, что аналитическое исследование поведения энтропии слишком затруднительно из-за сложности алгоритмов, автор предлагает пойти по следующему пути. Будем рассматривать любой применяемый алгоритм расчета газодинамических течений посредством уравнений (1.1) как «черный ящик». Пусть (безразлично каким способом ) в некотором узле сетки получены значения потоков, т.е. вектор F из трех величин (F 1 , F 2 , F 3 ) компонент (1.2). Будем рассматривать его как результат для некоторого условного, гипотетического газа. Параметры этого газа должны восстанавливаться , исходя из уравнений для потоков, дополненных нужным уравнением состояния, в рассматриваемом пока случае – уравнением (1.4): (5.1) (5.2) (5.3) Имеем систему четырех уравнений для четырех величин . Существует ли ее единственное решение и можно ли его восстановить? Предположим, что решение существует . Тогда из (5.1)-(5.3) последовательно получаем: (5.4) , (5.5) Запишем уравнение состояния в виде: (5.6) Подставляя в него результаты (5.4)-(5.5), получаем: (5.7) Это квадратное уравнение для u приводится к виду (конечно, ): (5.8) Его дискриминант (5.9) должен быть неотрицательным. Поскольку предполагается , что решение существует , воспользуемся соотношениями (5.1)-(5.3) и подставим их в формулу (5.9): (5.10) , что и требовалось. Следовательно, уравнение (5.8) может быть решено. Остается правильно распорядиться его двумя вещественными корнями: (5.11) , Опять же, поскольку решение существует, воспользуемся соотношениями (5.1)-(5.3) и подставим их в формулы (5.11). Получим: Одно из этих равенств превращается в тождество: (5.12) , При работе с формулами (5.11) определить, какой из корней является содержательным, проще всего, если вычислить оба корня и проверить их . По значениям и (или) могут быть вычислены (5.13) , . Возможность восстановить параметры условного газа по полученным значениям потоков реализуется по (5.9),(5.11) и (5.13) при условии . Конечно, после этого прежде всего проверяется «разумность» полученного результата: (5.14) , . Если эти условия выполнены, вступает в силу постулат об энтропии. По формуле (2.8) с привлечением параметров газа-«напарника», участвующего в расчете потоков, вычисляется величина изменения энтропии . Результат такого «жесткого энтропийного контроля» признается правомочным только в случае, если энтропия неубывает: . Конечно (хотя бы с учетом ошибок округления), его стоит несколько ослабить: (5.15) , или , где . Хочется пойти и на большее ослабление: в связи с обсуждением в конце § 3 и рис.2 допустимое значение увеличить. Однако уместно заметить, что «ложка дегтя может испортить бочку меда». Наконец, главное: из двух корней (5.11) «хороший» должен быть и только один – иначе «грозит» неединственность или нефизичность решения основной задачи. Заметим также, что речь идет о контроле , который (пока) не влияет на расчет. § 6. Снова о схеме C и ее экономичности Вернемся снова к тревожным результатам «жесткого энтропийного контроля», полученным для схемы С уже на простейшей модельной задаче с симметричным распадом разрыва. Чем же можно тогда объяснить предъявление успешных результатов расчетов по схеме С ? В частности, цитируем [4]: «Численное решение различных одномерных и пространственных задач показало, что в разработанной схеме нет проблем прохождения звуковой точки и не требуется снижение числа Куранта на сильных разрывах… Предложенный подход имеет фундаментальный характер и обобщен на стационарный случай. Разработаны монотонные схемы 2-го порядка… В качестве иллюстрации представлены результаты расчетов пяти вариантов типовых задач Торо [6]». Что можно выдвинуть в качестве возможных версий такого благополучного развития событий? Начнем с того, что, по-видимому, энтропией просто не интересовались . Оценки, сделанные в конце §3 для интервала значений m , гарантирующих попадание в условный пороговый интервал d , показывают, что, как правило, это попадание наверняка происходило. В таком случае возможность компенсировать «потери энтропии» в одном узле сетки за счет другого вполне реальна. Возможно, если ввести «мягкий энтропийный контроль» только по ячейкам сетки, то он и не обнаружит нарушений энтропийного режима. А даже, если и обнаружит, вовсе не обязательно, что они должны проявляться визуально. Внешне результаты могут выглядеть вполне удовлетворительно и с «неправильной» энтропией (и вот это тревожит и вынуждает ставить вопрос о контроле). В качестве еще одной версии можно предполагать, что круг решенных задач был достаточно «благоприятным» для схемы С . Интуитивно не вызывает сомнений, что она может успешно работать при расчете течений, у которых нет резких различий параметров в соседних ячейках разностной сетки. По-видимому, так и было. (Заметим, например, что для формул (1.7) даже не оговаривается , что делать в случае или ). Во всяком случае, для одномерных задач неизмеримо выросшая производительность вычислительных машин и многопроцессорных комплексов позволяет считать с таким числом узлов сетки, что достигнуть отсутствия этих «резких различий» вполне реально. Сложнее обстоит дело с двумерными и, тем более, пространственными задачами. К тому же главную роль играет еще и то, каковы эти задачи. В связи с этим стоит коснуться еще одного вопроса. Во введении был упомянут семинар, на котором автором работы [4] докладывались результаты, полученные по схеме С . По традиции объявление о семинаре сопровождалось аннотацией, которую, естественно, предоставляет докладчик. Приведем цитату из этой аннотации: «Метод является самым экономичным работоспособным приближенным решением задачи Римана распада разрыва и сокращает вычисления в 5-6 раз по сравнению с точным». В тексте [4] формулировка более сдержанная: «При существенно меньших затратах машинного времени на контактных разрывах, в зонах разрежения и при ударных волнах слабой и «средней» интенсивности результаты очень близки к схеме Годунова , сильный разрыв размазывается на 2-3 ячейки сетки, в то время как в схеме Годунова на 1-2». Отдав уже во введении должное этому действительно очень простому и компактному алгоритму, целесообразно сделать следующее замечание. Конечно же сокращение времени в несколько раз достигается, как за счет гораздо более простых формул, так и за счет того, что предлагаемая схема С , условно говоря, делает только одну итерацию, удовлетворяясь приближенными решениями. Так может быть для тех задач, которые рассматривались, и при использовании схемы Годунова достаточно было бы при расчете распадов разрывов делать одну итерацию, а то и вовсе ограничиваться «звуковым» приближением? Такое сравнение по затратам машинного времени было бы более корректным. Безусловно, формулы (1.7)-(1.8) действительно очень просты и экономны. Но уж слишком на рафинированную ситуацию они рассчитаны: сетка с постоянными шагами Dx (равномерная), неподвижная. Да, наверно, такие задачи так и нужно считать («зачем из пушки палить по воробьям»?). Конечно, очевидно, что эти формулы легко обобщаются на неравномерные, подвижные сетки. Но это будут уже другие формулы: для них потребуются скорости движения границ, пересчет координат узлов сетки и т.д. В свое время при расчете нестационарных задач с двумерной геометрией по методикам, основанным на обобщении одномерного варианта схемы С.К.Годунова, очень долго обходились «звуковым» приближением при расчете внутренних ячеек . И это несмотря на то, что итерационный метод расчета распада разрывов присутствовал уже в первой, сильно задержавшейся с публикацией в 1959 г., работе [7]. К этому времени метод успешно работал уже несколько лет. И основан этот метод расчета распада разрыва был как раз на использовании двух массовых скоростей а1 , а2 . Они обслуживали и ударную волну и волну разрежения (по своим формулам в каждом случае), и изменялись в ходе итерационного процесса. Этот алгоритм был позже воспроизведен в монографии [2], см. стр.110. Как отмечено еще в [7] на стр.291-292, детальное исследование показало, что итерационный процесс сходится, если в результате распада не получается очень сильных волн разрежения . Чтобы сделать его сходящимся, в [7] было предложено вести итерации по видоизмененным формулам, которые приведены и в [2] на стр.110. Исследование сходимости было опущено и в [7] и в [2], так как оно проводится стандартным способом исследования сходимости итерационных процессов типа , что сводится к вычислению и исследованию громоздких выражений для производных. Именно «ужасный» вид этих формул в свое время «вдохновил» автора на разработку того алгоритма расчета распада разрыва, который приведен в монографии [2] на стр.110-115. Вторым стимулом была реализация мечты о квадратичной скорости итерационного процесса за счет использования метода касательных Ньютона. Метод был разработан задолго до публикации монографии [2] и до сих пор присутствует в производственных программах. А настоятельная потребность в использовании такого алгоритма возникла, в частности, в связи с расчетами задач, представленных в работе [8], где этот алгоритм расчета распада разрыва и был впервые опубликован. Речь идет о расчетах нестационарных процессов, возникающих при взаимодействии набегающей ударной волны с течением за отошедшей головной ударной волной перед тупым телом, движущимся со сверхзвуковой скоростью. При их расчете возникали ситуации, когда в соседних ячейках двумерной сетки из-за сложной структуры и взаимодействия идущих в разных направлениях волн иногда возникали параметры, про которые иначе, как «нарочно не придумаешь», и не скажешь. Так что на приближенные расчеты распада разрыва надеяться не имело смысла. Потребовался такой алгоритм, который бы надежно выдавал «разумный» результат в ситуациях с любыми входными данными . Он и был разработан. Естественно, при его использовании предусмотрены различные методы контроля, которые позволяют завершать итерационный процесс при достижении уровня точности, назначаемого исполнителем конкретного расчета (уровнем, может быть, заведомо «грубым» ), ограничивать (на всякий случай) максимальное число итераций и т.п. Другое дело, что на практике далеко не всегда эти возможности используются в полной мере. Конечно, возможно и такое, что алгоритм будет делать «лишнюю» работу. Важно одно: в случае необходимости имеется возможность получить точный результат (конечно, в разумных достижимых пределах). Даже, если это не спасает аварийную ситуацию, которая может быть обусловлена и физическим содержанием рассчитываемого процесса. Наконец, остановимся еще на одном вопросе. Точный расчет распада разрыва оказался универсальным инструментом для расчета движения границ физических областей, например, разделяющих области с совершенно разными параметрами. Более подробно обсуждение этого вопроса можно найти в §15 монографии [2]. Результаты, которые дал бы приближенный метод расчета распада разрыва, как правило, были бы неприемлемыми . Это можно проиллюстри-ровать уже на том элементарном примере, который рассматривался в §2. Аналитическое решение задачи о распаде разрыва с входными данными (2.1) показывает, что (см., напр., [2], стр.112) разрежение до полного вакуума наступает при скорости разлета , т.е. при значении параметра . Использование точного распада разрыва позволяет его реализовать . А приближенный метод распада разрыва (как и вариант с ударной волной разрежения) по вине формулы (2.4) для давления всегда дает полный вакуум уже при . Для воздуха g=1.4, а поэтому получаем , а . Комментарии просто излишни. Легко представить, какими были бы погрешности в определении скоростей движения границ при использовании приближенного распада. § 7. Предложение по улучшению схемы С Как уже отмечалось, определяющую роль в схеме С играет назначение массовых скоростей формулами (1.8). Ввиду получения результатов, неблагоприятных с точки зрения поведения энтропии, естественно возникает вопрос, нельзя ли исправить ситуацию изменением этих формул. Простейшим способом реализации такой возможности является введение в формулы (1.8) произвольных коэффициентов. Что мы и сделаем, введя временно сразу три таких коэффициента : (7.1) Возникшую схему для краткости назовем СП-схемой, т.е. схемой С с параметрами. Исследуем ее на том же примере симметричного распада разрыва (2.1), который рассматривался в § 2. Будем предполагать, что , , . Как легко видеть, в симметричном случае Следовательно, для массовых скоростей вместо формулы (2.3) получаем: (7.2) Остальные формулы (2.2)-(2.9) сохраняются. В зависимости от знака безразмерного параметра формула (7.2) приобретает вид: (7.3) , где введено обозначение: (7.4) Чтобы при m =0 реализовать «звуковые» волны ( , в общем случае), надо полагать . Тогда для оценки энтропии вместо (2.10) с учетом (7.3) получим: (7.5) Изучение функции f( m) начнем с вычисления ее производной: (7.6) Для m =0 имеем , . При m <0 область определения функции (7.6) ограничена тремя условиями: , , Они будут выполнены, если: (7.7) , . Наша «энтропийная мечта» - достигнуть выполнения условий (7.8) при , при . Наибольший интерес представляет окрестность m =0: удастся ли переломить тенденцию при m <0? Поэтому начнем с разложения (7.6) в окрестности m =0. Следовательно, для (7.6) получаем: , (7.9) , Конечно, заманчиво было бы назначить (7.10) . Тогда А =0, в игру вступает следующий член и получается: , . Чтобы при m< 0 получить f( m )>0 для малых m , должно быть В< 0. Подставляя в формулу (7.9), получим: , что и требовалось: B< 0, если . Пока ограничимся этим случаем. Итак, для m< 0 с учетом формулы (7.4) полагаем, согласно (7.10): (7.11) Перейдем к случаю . Поскольку для придется иметь дело со всей бесконечной полуосью m , обращаемся к формуле (7.6). Ее можно записать в виде: (7.12) , где h( m) – произведение трех знаменателей в (7.6), g( m) – многочлен второго порядка от m. . Кропотливой выкладкой получаем: Следовательно, (7.13) , , , . Из (7.12) следует, что для получения при необходимо и достаточно . В свою очередь для этого необходимо и достаточно выполнение трех требований: (7.14) , , ( , если ). Последнее связано с разрешением отрицательных действительных корней уравнения g( m)= 0 , поскольку полуось m <0 нас сейчас не интересует. Первые два из требований (7.14) порождают ограничения: , . Для учета третьего ограничения воспользуемся достаточным условием . Тогда из (7.13) получим . Следовательно, в соответствии с (7.4), должно быть: (7.15) . Чтобы одновременно выполнить (7.11) и (7.15), предлагается такая формула назначения q1 и q2 : (7.16) , . Отметим, что в формуле (7.7), поскольку , получаем . В качестве иллюстрации рассмотрим пример. Пусть g=1.4. Тогда , Для будем иметь а =1.92 , b =0 , c =2×1.082 ×0.08 , для всех . Это хорошо. При получаем: , А =0, В =-0.84 . В <0, как и требовалось. По формулам (7.13): а =0, b =-1.44 , с =-2×0.36×0.4 для всех m из интервала –5<m <0. Это хорошо. (7.17) , . По формуле (7.12) получаем, что для всех m из интервала , т.е. на всем интервале, где определена формула (2.6) для плотности. Следовательно, неубывание энтропии гарантировано. На рис.2 пунктиром изображен эскиз графика при для рассмотренного примера. Благодаря параметрам q1 и q2 в СП-схеме ситуацию в буквальном смысле удалось переломить. Интересно сравнить описанное в примере со схемой С . Отличие только в том, что в ней задаются q1 = q2 =1. Тогда для будем иметь q = q1 + q2 =2. По формулам (7.13) а =5.6 , b =11.2 , c =8. Следовательно, для всех . Это хорошо . При получаем q = q1 - q2 =0. По формулам (7.9) А=- 2.4 , В =0.96 . Это плохо , так как , т.е. при малых m . Эти факты обсуждались в §§ 2-3, и отображены на рис.2 сплошной кривой . Подведем итоги. В СП-схеме назначаем параметры по формулам (7.15) и (7.16): (7.18) , . Легко проверить, что первая из них может быть записана в таком виде: (7.19) Напомним, что изучение случая отложено, и в настоящей работе рассматриваться не будет, чтобы не загромождать изложения. Кроме того, и полученные результаты не совсем обоснованы, поскольку формулы (7.18) получены для симметричного варианта распада разрыва. Возможно, они потребуют уточнения. В частности, чтобы в асимметричном случае распада при реализовать «звуковые» волны, необходимо в (7.1) заменить скорости на , где , Полагая , массовые скорости будем вычислять так: (7.20) Тогда, как проверено для симметричного примера (2.1), дальнейшие вычисления в СП-схеме по формулам схемы С , описанным в § 1, должны позволить избежать разрывов, на которых убывает энтропия. Конечно, такое предложение должно быть тщательно проверено в численном эксперименте. Хотелось бы, чтобы описанный в настоящей работе контроль энтропии СП-схема выдерживала автоматически (как схема С.К.Годунова с соответствующим алгоритмом расчета распада разрыва). В предлагаемой СП-схеме, аналогично схеме С , речь идет только о приближенном расчете распада разрыва. Для модельного примера симметричного распада разрыва сохраняется роль вместо значения в аналитическом решении. Поэтому все, что было изложено в конце § 6 по поводу необходимости использовать для счета границ и в случаях резкого различия параметров в соседних ячейках сетки алгоритм точного распада разрыва, остается в силе. Роль алгоритма точного расчета распада разрыва в качестве универсального инструмента решения таких вопросов по-прежнему сохраняется. Что касается больших затрат машинного времени, то что же делать – «бесплатный сыр бывает только в мышеловке». Эти затраты объективно вызваны сложными формулами точного расчета распада разрыва в алгоритме, описанном в [2] на стр.110-115. Фактически это – плата за возможность в ходе итерационного процесса продвинуться от «барьера» к правильному значению . Если численный эксперимент не обнаружит скрытых недостатков предлагаемой СП-схемы, можно было бы вместо «звукового» приближения рассчитывать начальное приближение для алгоритма точного распада разрыва посредством СП-схемы. (Напомним, что можно остановиться на вычислении давления , если из-за асимметрии распада разрыва энтропия все же уменьшится. В этом случае – нарушения условия (5.15) – следует продолжить расчет, как в методе Годунова – см. § 4). Более того, в случаях, когда нет резкого различия между параметрами в соседних ячейках (это типично для сквозного счета газодинамических течений) полученное приближение может оказаться вполне достаточным, чтобы им и ограничиться . § 8. О проблеме сложных уравнений состояния. При решении практических задач зачастую приходится иметь дело с ситуацией, когда уравнение состояния среды задается не в виде идеального газа (1.4), а описывается гораздо более сложными зависимостями. Опыт работы с методом С.К.Годунова выявил особую роль так называемых двучленных уравнений состояния: (8.1) . В отличие от (1.4), оно содержит три параметра . Удобно вместо параметра r0 использовать параметр р0 , определяемый формулой: (8.2) Для дальнейшего будут нужны также соответствующие (8.1) формулы для скорости звука с и энтропии s : (8.3) , . Очевидно, что идеальный газ (1.4) будет частным случаем ( ). Нетрудно убедиться, что все изложенное в §§ 1-4 реализуется и для уравнения состояния (8.1) после внесения очевидных изменений в те формулы, которые его используют. Описание в [2] на стр.110-115 задачи о распаде разрыва это подтверждает, поскольку алгоритм был разработан сразу именно для двучленного уравнения состояния (8.1). Поэтому обратимся сразу к § 5, где это еще предстоит сделать. Замена идеального газа (1.4) на (8.1) приводит к тому, что уравнение (5.7) придется заменить на следующее: . Полученное уравнение тоже является квадратным относительно u и, с учетом (8.2), приводится к виду: (8.4) Его дискриминант (8.5) тоже должен быть неотрицательным. Формулы (5.11) заменяются на следующие: (8.6) Нужный корень или определяется так же, как в § 5, с привлечением параметров газа-«напарника». При этом условия (5.14) заменяются на следующие: , , а условие неубывания энтропии (5.15) формулируется в виде: (8.7) Присутствие в формуле (8.7) могло бы привести к изменению важной формулы (7.5), которая столь громоздко исследовалась. Но этого не произойдет : величина назначается формулой (7.2), в которой участвует лишь через . Поскольку , как и было ранее, если полагать (8.8) , то формула (7.5) сохранит свой вид . Следовательно, все изменения сводятся только к определению величины m формулой (8.8) вместо прежней в случае идеального газа (1.4). Это позволяет в случае двучленного уравнения состояния (8.1) воспользоваться всеми результатами проведенных в § 7 исследований, имея лишь в виду определение (8.8). Упомянутая выше особая роль двучленного уравнения состояния состоит в следующем. Как выяснилось при многолетней эксплуатации метода С.К.Годунова, при работе с уравнениями состояния, которые задаются сложными зависимостями , оказывается полезным следующий прием, предложенный А.В.Забродиным. Для решения задачи о распаде разрыва эти сложные зависимости аппроксимируются локально двучленными уравнениями состояния (8.1). Формулы для вычисления параметров можно найти, напр., в [3] на стр.177. Такой прием успешно используется в практике расчетов, но не является универсальным . Трудности возникают, если для параметров , вычисленных указанным способом, не выполняются требования (8.9) , , . Как указано в [3] на стр.179, «двучленная аппроксимация может быть использована для уравнения состояния общего вида при условии, что оно удовлетворяет условиям нормального газа… Для других типов уравнений состояния решение задачи Римана может быть неединственным , иметь неклассическую структуру или даже нарушить гиперболичность системы уравнений… В этих случаях использование двучленной аппроксимации может, вообще говоря, давать посторонние или физически неприемлемые решения». Есть ли выход из такой ситуации? Есть. Сейчас уместно упомянуть об универсальном подходе к решению задачи о распаде разрыва в случае сколь угодно сложных уравнений состояния, изложенном в небольшой монографии [9], изданной в 1970 году. Правда, этот подход требует организации довольно сложного «хозяйства», и его рассмотрение выходит за рамки настоящей работы. Еще одна проблема связана с тем, что уравнения состояния могут не удовлетворять термодинамическому тождеству. В таком случае для них не может быть определена энтропия , а, следовательно, теряет смысл постановка вопроса о ее контроле. В качестве конкретного примера можно привести реализацию уравнений состояния, заданных в табличной форме в виде зависимостей (8.10) , , где Т – температура. В работе [10] эти таблицы аппроксимировались независимыми билинейными функциями вида (8.11) , со своими коэффициентами в каждой ячейке таблиц. При этом, за исключением случайных ситуаций , не выполнено термодинамическое тождество (см. [11], стр.6). В работах автора [11-12] предложен громоздкий, но вполне реализуемый путь преодоления этого принципиального затруднения. Уравнение состояния в отдельной табличной ячейке конструируется как линейная комбинация частных решений, удовлетворяющих условию термодинамического согласования. Чтобы такая комбинация обеспечивала нужные табличные значения функций (8.10) в четырех углах ячейки, она должна содержать 8 свободных параметров. При удачном выборе частных решений для вычисления параметров выписываются явные формулы. Поскольку для используемых частных решений энтропия определена , она будет определена и для построенного решения. Это открывает дорогу для реализации энтропийного контроля, описанного в настоящей работе. Для уравнений (8.11), возникающих при реализации табличных данных, гипотетически возможен и другой путь. Из уравнения для энергии e выражаем температуру Т: и подставляем в уравнение для давления: (8.12) Далее полученное уравнение состояния с помощью упомянутого выше алгоритма локальной аппроксимации «превращаем» (локально) в двучленное уравнение состояния (8.1) с некоторыми параметрами (или р0 ). Если нам повезет и будут выполнены условия (8.9), то все в порядке. Можно будет вычислить энтропию по формуле (8.3) и реализовать описываемый ее контроль. А если нет? Упираемся в уже упомянутую выше проблему. Судя по изложенному в статье [10], стр.782, «трудности возникают при наличии областей, где среда находится в состоянии фазового перехода между жидкостью и газом (кипящая жидкость). В такой области некоторые термодинамические величины, например, теплоемкость при постоянном давлении, обращаются в бесконечность, что недопустимо для двучленного уравнения состояния» . Это вынудило автора [10] отказаться от такого пути. В связи с этим возникает серьезный вопрос о допустимости вообще уравнений состояния, не удовлетворяющих термодинамическому тождеству . Не собираемся ли мы бороться с «нефизическими» решениями (типа разрывов, на которых убывает энтропия) в «нефизических» условиях, которые сами же допускаем? Заключение Как следует из названия настоящей работы, главной ее целью автор считает изложение метода контроля поведения энтропии в ходе выполнения расчетов газодинамических течений. Описанные в ней практические алгоритмы такого контроля просты и могут быть использованы при выполнении расчетов для разнообразных известных и широко используемых методов. Как правило, сами эти методы настолько сложны, что вполне могут «скрытно» допускать при численной реализации нефизические решения, в частности, газодинамические ударные волны разрежения. Решающую роль сыграла возможность аналитически разобраться в вопросе о поведении энтропии на примере очень простой и изящной разностной схемы С , опубликованной в работе [4]. При этом, к сожалению, получили подтверждение опасения о возможности реализации таких разрывов, на которых уменьшается энтропия. Следовательно, они по существу аналогичны ударным волнам разрежения. Автор предвидит возможные возражения, напр., такого характера: «А зачем нам энтропия? Аэродинамические нагрузки определяют давление и плотность, а их мы считаем хорошо». Ну, что же, возможно , в таких задачах и можно не обращать внимание на энтропию. (Хотя и здесь есть подозрение: не слишком ли хорошо? Но ведь этим не исчерпывается круг задач, которые приходится решать. В качестве поддержки своей позиции автор предлагает обратиться к стр.5 работы [13]. Там сказано, что дополнительное требование неубывания энтропии «позволяет в классе всех возможных кусочно-гладких решений газодинамических уравнений выделить единственное физически реализуемое , отбросив решения с термодинамически аномальным поведением, вроде ударных волн разрежения. Таким образом, задача построения разностного алгоритма расчета газодинамических течений общего вида не сводится к выбору разностной схемы для системы дифференциальных уравнений, а является гораздо более сложной задачей» . Идентичное по содержанию мнение С.К.Годунова автор цитировал на стр.100 работы [1]. На стр.7 работы [13] отмечаются трудности правильного моделирования изменения энтропии, поскольку на разностной сетке гладкие волны сжатия формально неотличимы от размазанных ударных фронтов . В то же время на гладких участках энтропия должна сохраняться , а на фронтах разрывов – расти . Автор обращается с предложением к специалистам, конструирующим алгоритмы, и тем, кто считает газодинамические задачи. Осуществлять в процессе расчета описанный в настоящей работе контроль энтропии для специалистов, располагающих реализованными алгоритмами , не должно составлять большого труда. (Конечно, из соображений экономии возможен и выборочный контроль наиболее «ответственных» мест по пространству и времени.) При этом автоматизированный контроль энтропии совсем не обязательно предполагает остановку расчета (если этого не хочет исполнитель, даже в неблагополучной ситуации). Важно, что о таких ситуациях контроль может сообщать исполнителю и подготовить его к тому, что результаты расчета могут оказаться сомнительными (или, напротив, что с энтропией все благополучно). Одновременно будет накапливаться информация по обсуждаемому вопросу. Автор благодарен М.С.Гавреевой и А.В.Северину за помощь в оформлении работы. Литература 1. Прокопов Г.П. Загадка метода Годунова.// Вопросы атомной науки и техники. Сер.: Математ. моделир. физ.процессов. 2005, вып.4, с.98-101. 2. Численное решение многомерных задач газовой динамики. Под ред. С.К.Годунова.// М., «Наука», 1976, 400 с. 3. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений//М., ФИЗМАТЛИТ, 2001, 608 с. 4. Сафронов А.В. Разностный метод для уравнений газодинамики из соотношений на разрывах для вектора потока.// Материалы XVII Школы-семинара «Аэродинамика летательных аппаратов», ЦАГИ, 2006. 5. Зоммерфельд А. Термодинамика и статистическая физика.// М., Изд-во иностр. литер., 1955, 480 с. 6. Toro E.F. Riemann Solvers and Numerical Method for Fluid Dynamics.// Springer – Verlag. Second Edition, June 1999. 7. Годунов С.К. Разностный метод численного расчета разрывных решений уравнений гидродинамики.// Матем. сб., 1959, 47, вып.3., с.271-306. 8. Прокопов Г.П., Степанова М.В. Расчет осесимметричного взаимодействия ударной волны с затупленным телом, движущимся со сверхзвуковой скоростью.// Препринт АН СССР, 1974, №72, 26 с. 9. Алалыкин Г.Б., Годунов С.К., Киреева И.Л., Плинер Л.А. Решение одномерных задач газовой динамики в подвижных сетках.// Изд-во «Наука», Гл.ред.физ.-мат.лит., М., 1970, 112 с. 10. Чарахчьян А.А. Об алгоритмах расчета распада разрыва для схем С.К.Годунова.// ЖВМиМФ, 2000, 40, №5, с.782-796. 11. Прокопов Г.П. Аппроксимация табличных уравнений состояния для расчета газодинамических задач.// Препринт ИПМ им.М.В.Келдыша РАН, 2004, №80, 28 с. 12. Прокопов Г.П. Исследование формул аппроксимации табличных уравнений состояния в переменных температура-плотность.// Препринт ИПМ им.М.В.Келдыша РАН, 2005, №26, 28 с. 13. Забродин А.В., Софронов И.Д., Ченцов Н.Н. Адаптивные разностные методы математического моделирования нестационарных газодинамичес-ких течений. (Обзор).// Вопросы атомной науки и техники. Серия: Методики и программы численного решения задач математической физики, 1988, вып.4, с.3-22.