Г. П. Прокопов Оприближенных реализациях Ордена Ленина Институт прикладной математики имени М.В.Келдыша Российской академии наук Г.П. Прокопов О приближенных реализациях метода Годунова Москва, 2007 год УДК 519.6 О приближенных реализациях метода Годунова Прокопов Г.П. Препринт Института прикладной математики им. М.В.Келдыша РАН В целях экономии объема вычислений при реализации метода С.К.Годунова для расчета газодинамических течений используются приближенные алгоритмы. Если при этом исключается адиабата Пуассона, возможно привлечение таких конструкций, на которых энтропия может убывать. Вопрос об их допустимости является, по меньшей мере, дискуссионным. Работа выполнена при поддержке Российского Фонда Фундаментальных Исследований (грант № 05-01-00097). On approximate realizations of Godunov method. Prokopov G.P. Preprint of Keldysh Institute of Applied Mathematics, RAS To economize the computational costs of Godunov method realizations one uses approximate algorithms to calculate gas-dynamics flows. If Poisson adiabate is not considered in calculations , there is possibility entropy decrease. Admissibility of such solutions is under a discussion. Содержание стр. Введение …………………………………………………………... 3 § 1. Об упрощенной схеме течения при распаде разрыва ……... 3 § 2. Фронт разрыва – возможный источник убывания энтропии. 7 § 3. О единственности решения и выполнении уравнения состояния ……………………………………………………… 9 § 4. Алгоритм приближенного расчета распада разрыва ………. 12 § 5. Некоторые замечания о практической реализации приближенного расчета ……………………………………… 18 § 6. Снова об энтропии …………………………………………… 23 § 7. О сравнении результатов …………………………………….. 24 Заключение …………………………………………………….…. 26 Литература ……………………………………………………….... 27 Введение Настоящая работа является непосредственным продолжением [1]. Продолжением в определенной степени вынужденным, поскольку в ее § 7 не были завершены исследования предложения по улучшению схемы С , опубликованной в [2]-[3] . При проведении этих исследований возникли некоторые новые результаты, которые по мнению автора заслуживают внимания и станут содержанием настоящей работы. Они касаются не только схемы С , но и вообще приближенных алгоритмов «типа метода Годунова», предлагаемых различными авторами. Описывается один из таких алгоритмов, не допускающий убывания энтропии. § 1. Об упрощенной схеме течения при распаде разрыва Для численного решения газодинамических задач широкое применение получил известный метод Годунова. Его первая развернутая публикация была дана в [4]. Затем популярности метода способствовала монография [5], в основном посвященная именно ему и его применению для решения разнообразных задач. 1.1. Основным структурным элементом метода и главной «массовой операцией» при его практической реализации является решение задачи о распаде разрыва с параметрами газа в соседних ячейках разностной сетки. Эта задача Римана настолько хорошо известна специалистам, что мы не будем сейчас описывать ее постановку. Как правило (например, при «сквозном» счете газодинамических течений), параметры в соседних ячейках достаточно близки между собой. Это, естественно, создает благоприятные условия для того, чтобы использовать упрощенные алгоритмы для решения обсуждаемой, хотя и элементарной, но достаточно громоздкой задачи. Но так бывает, во-первых, не всегда, а, во-вторых, далеко не все методы ее приближенного решения можно считать приемлемыми для расчета произвольных газодинамических задач. Как описано, например, в [5] на стр.106-107, в точной постановке при распаде произвольного разрыва образуется контактный разрыв и две волны, которые распространяются от него. Условно их называют правой и левой волной (условно потому, что в неподвижной системе координат они могут двигаться и в одну сторону). Каждая из волн может быть либо ударной, либо волной разрежения. В случае ударной волны речь идет о движущемся фронте разрыва, по обе стороны которого параметры газа полагаются постоянными (своими для каждой из сторон) и связаны определенными соотношениями. В случае волны разрежения ситуация иная: имеется область переменного течения, в которой параметры газа постоянны вдоль прямолинейных лучей, играющих роль характеристик системы дифференциальных газодинамических уравнений, а значения параметров зависят (определенным образом) от наклона в «веере характеристик», описывающем волну разрежения. Эта зона разрежения граничит с областями постоянного течения, аналогичными описанным выше для ударной волны. В качестве примера на фиг.1 изображена одна из таких возможных конфигураций с левой ударной волной и правой волной разрежения. 1.2. В методе, описанном в [2]-[3], который для краткости называется схемой С , используется более простой алгоритм для приближенного решения задачи Римана о распаде разрыва. Его принципиальное отличие состоит, прежде всего, в том, что переменное гладкое течение в зоне разрежения (если она присутствует в возникающем конкретном варианте распада разрыва) заменяется на разрывное течение с постоянными параметрами по обе стороны фронта разрыва. Таким образом, рассматривается упрощенная схема течения, включающая контактный разрыв и левый и правый фронты разрыва (называть их волнами было бы не всегда неправильно). В узле сетки с некоторым номером j такая схема изображена на фиг.2. 1.3. Обозначим индексами: 1 – параметры в левой ячейке сетки, 2 – в правой, 3 – между левым и контактным разрывом, 4 – между контактным разрывом и правым разрывом. Параметры газа представляют соответственно: u – скорость, r - плотность, р – давление, e – внутреннюю энергию единицы массы газа. Последние три величины связаны так называемым уравнением состояния. В простейшем случае оно имеет вид, которым мы пока и ограничимся: (1.1) . Здесь g - числовой параметр (показатель адиабаты). Фиг. 1. Фиг. 2. Сразу же отметим, что скорость звука с , которая потребуется в дальнейшем, в этом случае определяется формулой: (1.2) Далее, как указано в [3] на стр.155, характерная особенность схемы С заключается в априорном задании двух параметров: массовых скоростей через поверхность левого и правого разрывов. Эти параметры, которые мы обозначим a 1 , a 2 (предполагается, что a 1 >0, a 2 >0 ) назначаются по состоянию газа в соседних ячейках сетки, т.е. в зонах 1 и 2, с помощью некоторых конкретных формул. Временно (до § 2) отложим описание этих формул. Далее, исходя из отношений, которые должны быть выполнены на разрывах, выписываются формулы для «простых» (терминология [3]) переменных P , U : (1.3) , Они задают давление Р и скорость U газа на контактном разрыве и в двух прилегающих к нему зонах 3 и 4. Затем, в силу тех же соотношений на разрывах, вычисляются другие параметры в зонах 3 и 4: (1.4) и скорости движения разрывов: (1.5) , Назначение параметров на границах ячейки для вычисления векторов потока (в методе Годунова они назывались «большими» величинами) определяется значениями U , D 1 , D 2 . Соответствующие формулы мы опустим. 1.4. В [1] было проведено «энтропийное» исследование случая распада разрыва с симметричными данными. Для дальнейшего целесообразно исследовать с «энтропийной» точки и произвольный («асимметричный») случай распада разрыва. Оказалось, что такое исследование не так уж и громоздко. Введем в рассмотрение величины: (1.6) Тогда формулы (1.3) можно записать в виде: (1.7) , , Из формул (1.4) получаем: (1.8) Энтропию S для уравнения состояния (1.1) определим формулой (1.9) Для изменения энтропии на левом и правом фронтах разрывов получаем: (1.10) Используя формулы (1.6)-(1.8), имеем: (1.11) Чтобы рассматриваемая упрощенная схема описывала реальное физическое течение, согласно постулату о неубывании энтропии должны быть одновременно выполнены два условия: (1.12) , . § 2. Фронт разрыва – возможный источник убывания энтропии. 2.1. В качестве первого примера рассмотрим схему С . В [2] для нее предлагается задавать массовые скорости формулами: (2.1) , где с1 ,с2 - скорости звука, вычисляемые с помощью (1.2). Если ввести величины (2.2) , , они записываются в виде (см. [1], с.8): (2.3) , , В случае - «разбегающихся» потоков с точки зрения системы координат, движущейся со скоростью , - будем иметь и, следовательно, (2.4) , . В [1] было показано, что уже для частного случая распада разрыва с симметричными данными в некоторых ситуациях на разрывах может уменьшаться энтропия . Этого уже достаточно для тревоги относительно физической реализуемости, а, следовательно, достоверности получаемых результатов. Заметим, что формулы (2.4) предлагались как альтернативный вариант задания массовых скоростей в работе [3] наряду с (2.1). Поскольку в случае разбегающихся потоков они совпадают, примеры убывания энтропии годятся и для варианта (2.4). В случае - «встречных» потоков – из (2.3) получаем (2.5) , . Исследование в [1] для симметричного распада разрыва показало, что эти формулы в случае встречных потоков обеспечивают возрастание энтропии, но только при дополнительном предположении . 2.2. В [1] предлагался вариант улучшения схемы С посредством введения дополнительных параметров, названный условно СП-схемой. Для «спасения» энтропии предлагалось массовые скорости задавать формулами, обобщающими (2.3): , , . Параметры предлагалось назначать некоторыми специальными формулами, полученными при рассмотрении симметричного варианта распада разрыва. Однако, нетрудно убедиться, что в случае асимметричного распада разрыва обеспечить выполнение обоих условий (1.12) неубывания энтропии для произвольных значений параметров , невозможно . Проще всего это подтверждается в случае, если . Это позволяет сразу исключить роль . Получаем, как и в случае (2.4), , , и формулы (1.6) приобретают вид: , . Если , то , . В соответствии с (1.11) рассмотрим функцию (2.6) . Если величина а не зависит от параметра q , для ее производной получаем формулу: . В точке q =0 будем иметь: , Отсюда легко получается, что, по крайней мере, на слабом разрыве , если , а следовательно . Если ограничиться случаем слабых разрывов, можно было бы воспользоваться приближенной формулой Тогда, в соответствии с формулой (2.6), получаем: (2.7) . Отсюда очевидно следует, что для достаточно слабого разрыва , если и или при любом для . Возможны и другие ситуации, когда , если «положительного запаса» первого слагаемого недостаточно для «нейтрализации» отрицательного второго. Старательные оговорки автора в [1] на стр.22-23 о необходимости тщательной проверки «улучшенной» схемы в численном эксперименте, который может обнаружить скрытые недостатки СП-схемы, к сожалению, подтверждаются и без эксперимента. 2.3. Подводя итог изложенному, можно констатировать, что замена течения в области волны разрежения фронтом разрыва может приводить к тому, что на этом разрыве убывает энтропия. Такое течение следует считать физически нереальным. Подчеркнем, что это может происходить при сколь угодно близких значениях параметров в соседних ячейках. Автору представляется уместным следующее замечание. В словаре русского языка [6] на стр.1128 читаем: «Фикция – положение, построение, которому ничто не соответствует в действительности, но которым пользуются как допущением с какой-нибудь определенной целью». Именно такой фикцией может оказаться результат расчета распада разрыва при использовании упрощенной схемы течения. § 3. О единственности решения и выполнении уравнения состояния. 3.1. Корректность математической постановки задачи, гарантирующая успешную работу вычислительного алгоритма, в качестве неотъемлемых элементов включает существование и единственность искомого решения. Для точного расчета распада произвольного разрыва эта серьезная проблема обсуждалась, и мы не будем на ней останавливаться. (Во избежание недоразумений отметим, что в ситуации, если, например, задача имеет два решения – как в случае сильного и слабого отражения волн – алгоритм должен предусматривать, какое именно решение будет искомым). Запрет ударных волн разрежения как дополнительное требование, накладываемое на обобщенное решение системы газодинамических уравнений для обеспечения его единственности , появилось в [7], опубликованной еще раньше, чем [4].Пример, представленный в [5] на стр.115-116, о распаде разрыва с симметричными данными относительно границы раздела, является иллюстрацией такой ситуации. Как показано в предыдущем § 2, схема С в качестве вспомогательного инструмента, именуемого приближенным решением задачи Римана, может привлекать разрывы, на которых энтропия убывает. Фактически это могут быть физически нереальные течения. Наряду с единственным реальным обобщенным решением (которое дает точный расчет распада разрыва) таких фикций может существовать сколько угодно. Достаточно априорно задать пару массовых скоростей – и пожалуйста. Эта фикция потом «лакируется пересчетом» и может давать внешне прекрасный результат. Не будем обсуждать достоинств формул (2.1) – они, конечно, есть. Но как все-таки быть с единственностью? Как быть, если используя другие формулы для массовых скоростей, получится другой результат для некоторой важной характеристики рассчитываемой конструкции (например, о роли ее конкретной детали)? Апеллировать к точному расчету распада разрыва? Не лучше ли действовать так сразу (причем разумно распоряжаясь закладываемыми в него возможностями экономить время конкретного расчета – позже вернемся к обсуждению этого вопроса)? 3.2. У схемы С есть еще одно уязвимое место. В [2] отмечается как ее достоинство тот факт, что величины , по которым вычисляется вектор потока на границе ячеек, получаются без привлечения уравнения состояния . Более того, в [3] на стр.155 читаем: «Отсутствие необходимости использования уравнения состояния для определения вектора потока на границе ячеек… облегчает приложение предложенной схемы к сложному уравнению состояния газа с переменными свойствами. Поэтому первоначально данный метод успешно применялся для расчета струйных течений реагирующих газов на основе решения уравнений Навье-Стокса». А вот это уже серьезное и обязывающее заявление. Поскольку из трех величин R , P , E независимыми являются только две, то ситуация далеко не столь проста. Параметры в зонах с номерами 3 и 4 должны удовлетворять уравнению состояния (1.1): (3.1) , . Если подставить в эти формулы выражения, представленные в (1.3)-(1.4), то получаются два соотношения, которые следует рассматривать как уравнения для величин массовых скоростей а1 ,а2 . Именно они и определяются в ходе итерационного процесса для расчета точного распада разрыва (естественно, если этот процесс сходится). Высказывание в [2] о том, что точное вычисление массовых скоростей не требуется в силу выполнения соотношений на разрыве, можно оправдать только тем обстоятельством, что ставится цель получения лишь приближенного решения. Никаких оснований ожидать выполнения (3.1) при априорном («волевом») назначении массовых скоростей, вообще говоря, нет. 3.3. Каковы могут быть последствия невыполнения условий (3.1)? Фактически это означает, что при последующей реализации формул «пересчета» будет происходить «смешивание» нескольких газов, удовлетворяющих не уравнению состояния (1.1), а, условно говоря, различным. А результирующему газу снова «навязывается» уравнение состояния (1.1). Вопрос о том, насколько корректна такая процедура (при ее систематическом употреблении), остается открытым. Вычислительный опыт в таких ситуациях обычно допускает конструирование неблагоприятных примеров. Возможно , иллюстрацией к обсуждаемому вопросу могут служить результаты расчета теста 2 в упомянутых работах [2]-[3]. В нем рассматривался симметричный разлет двух газов с параметрами , , , . При этом образуются сильные волны разрежения. График внутренней энергии, полученный для разностного счета, настолько сильно отличается от аналитического решения, что ставит под сомнение уравнение состояния. Гипотеза состоит в том, что выполненный расчет в зоне сильного разрежения фактически реализует, условно говоря , некоторый «эффективный» показатель адиабаты, отличный от заданного значения g , несмотря на то, что в расчетных формулах только оно и присутствует. Мы еще раз вернемся к этим результатам позже (в § 7). Заметим, что при реализации схемы С.К.Годунова с точным расчетом распада разрыва, описанным в монографии [5] на стр.110-115, подобная ситуация не допускается. Во-первых, величина Е досчитывается по Р и R именно с привлечением уравнения состояния. Во-вторых, в случае доведения итерационного процесса до сходимости, будут удовлетворены точные соотношения на ударных волнах или волнах разрежения, т.е. реализованы единственные правильные значения их массовых скоростей. В таком случае результаты вычисления внутренних энергий по формулам (3.1) и (1.4) просто совпадают. 3.4. Автор считает уместным в связи с обсуждаемым вопросом напомнить о давней работе [8], которая представляет текст доклада С.К.Годунова на IV Всесоюзном математическом съезде в 1961 году. В ней, в частности, на стр.157 приведен пример, показывающий, что последовательность решений уравнений газовой динамики в форме законов сохранения, удовлетворяющих уравнению состояния, может иметь предел, уравнению состояния не удовлетворяющий. В те времена становления прикладной математики проблема существования представлялась чрезвычайно важной и была предметом серьезных дискуссий. Приведем цитату с той же стр.157: «…Имеющийся в нашем распоряжении опыт работы с кусочно-гладкими аналитическими решениями относится только к случаю с конечным числом областей гладкости. Однако даже если в начальных данных содержится всего 3-4 разрыва, может случиться, что в дальнейшем число этих разрывов будет увеличиваться с течением времени по геометрической прогрессии и станет бесконечным уже при конечных t . Произойдет как бы сгущение особенностей. Никакого опыта в изучении таких сгущений у нас нет. Продолжимо ли решение за это сгущение?…» Не удивительно поэтому, что разностный алгоритм может выдавать «острую реакцию», отражающую реально складывающуюся в рассчитываемом течении. Так ли уж правильно заранее заботиться о сглаживании всех возникающих эффектов, исходя из соображений «монотонизации» искомого решения? Конечно, сейчас вычислительного опыта накопилось существенно больше. Однако окончательного ответа на такого рода вопросы, по-видимому, так и не появилось. § 4. Алгоритм приближенного расчета распада разрыва. Исходя из изложенных выше соображений, автор считает целесообразным описать один из возможных алгоритмов приближенного расчета распада разрыва, не допускающий убывания энтропии. Конечно, многие его детали уже описывались и в монографиях [5], [9] и в работах других исследователей. Автор считает своим долгом отметить также [10], на которую обратил внимание при подготовке настоящей работы, благодаря ссылке на стр.170 монографии [9]. Для полноты и связности изложения известные уже детали придется повторить. Но некоторые элементы являются новыми и могут представлять практический интерес. Особое внимание уделено аргументации предлагаемых решений. 4.1. Итак, в качестве исходной информации задаются: - параметры течения в соседних ячейках , ; - три термодинамических параметра связывает уравнение состояния (1.1) с заданным показателем адиабаты g ; вопрос о более сложных уравнениях состояния будет обсуждаться позже (в § 5); - задается безразмерный параметр , который будет играть роль критерия применимости «звукового приближения». 4.2. Расчет начинается с вычисления скоростей звука (4.1) , и «пробного» назначения массовых скоростей (4.2) , Заметим, что «волевое» назначение массовых скоростей a 1 , a 2 – операция достаточно ответственная и достигает цели только в случае, если обеспечивает достаточно аккуратную аппроксимацию исходных уравнений. Поэтому из возможных различных вариантов, упоминавшихся ранее, выбирается именно этот, поскольку представляет математически обоснованную линеаризацию уравнения для давления Р : (4.3) Отметим, что это предлагалось, например, в [5] на стр.113. 4.3. В соответствии с изложенным в § 1, вычисляем по формулам (1.6). Если выполняются два условия: (4.4) , , полагаем, что достаточным является «звуковой» расчет распада разрыва. В этом случае, в соответствии с (1.7), назначаем: Для оценки изменения энтропии можно воспользоваться приближенной формулой (2.7): (4.5) Очевидно, что совершенно аналогично: Следовательно, при выполнении условий (4.4) получаем: (4.6) , . Если «примириться» со столь незначительным (при достаточно малом значении параметра ) убыванием энтропии, можно вычислять величины r3 , r4 по формулам (1.8), величины е3 ,е4 – по формулам (1.4), величины D1 , D2 – по формулам (1.5). Это те формулы, которые предлагались в схеме С , однако со следующим существенным замечанием. Именно благодаря назначению a1 , a2 по формулам (4.2) достигается обращение в нуль первого слагаемого в формуле (4.5). При назначении a1 , a2 формулами (2.1) или (2.4), по крайней мере, для одной из величин этого не будет, и она будет иметь не второй , а лишь первый порядок относительно (вот где сказывается аппроксимация!). Соответственно, «энтропийный шум» при назначении (2.1) или (2.4) будет больше, чем при назначении (4.2). К вопросу о допустимости, хотя и незначительного, но все же убывания энтропии мы еще вернемся. Ситуацию (4.4) для краткости будем дальше называть «звуковым приближением» для распада разрыва. 4.4. Если хотя бы одно из условий (4.4) не выполнено, придется применять более сложный алгоритм вычисления величины Р . Как описано в [5] на стр.110 (см. также [9], стр.168), она должна удовлетворять уравнению: (4.7) , где для значений индекса k=1,2: (4.8) , Для упрощения описания алгоритма предполагается, что (см.[5], стр.111). 4.5. Начинаем (см. замечание ниже в разделе 5.3) с вычисления величин, определяющих «содержание» конфигурации распада разрыва. (4.9) Возможны следующие ситуации: 10 если , то решения уравнения (4.7)-(4.8) не существует: образуется зона вакуума, в которой следует полагать Р =0, R =0, более подробное рассмотрение этого случая можно найти в [11]; 20 если , то искомое решение находится в интервале , т.е. образуются две волны разрежения; 30 если , то ; слева имеем ударную волну, справа – волну разрежения; 40 если , то ; образуются две ударные волны. Случаи равенства значениям или сознательно не оговариваются, чтобы при программной реализации (см. ниже раздел 5.3) охватывать их при «первой возможности». 4.6. Начнем с ситуации с образованием двух волн разрежения. Тогда уравнение для давления Р имеет вид: Легко видеть, что из него явно определяется величина и поэтому сразу выписывается формула для его решения: (4.10) Эта формула приведена в [5] на стр.114 и в [10] на стр.44. Можно было бы считать это исчерпывающим решением для рассматриваемой ситуации, если бы не следующее обстоятельство. При работе со сложными уравнениями состояния приходится, например, прибегать к их локальной аппроксимации двучленными уравнения состояния. Тогда параметры (а с ними и другие параметры) оказываются различными для индексов 1 и 2. Поэтому описанный способ становится неприемлемым и придется прибегнуть к другому алгоритму. К его описанию применительно к следующей ситуации мы и переходим. К рассматриваемой ситуации вернемся после нее. 4.7. Пусть . Тогда и уравнение (4.7) будет содержать обе разновидности формул (4.8). Предлагается следующий алгоритм его приближенного решения. Аппроксимируем F(P) линейной функцией , принимающей нужные значения на концах отрезка: , . Очевидно, что имеет вид: Далее решение линейного уравнения принимаем в качестве приближенного искомого решения. В качестве полезного практического совета в [10] на стр.45 рекомендуется поменять роли переменных Р и j и сразу определять обратную зависимость Р от j : Тогда искомый приближенный корень определяется формулой: (4.11) В [10] на стр.45 говорится, что «не представляет сложности и повысить точность интерполяции, используя значения производных в этих же точках, при этом удобно интерполировать полиномом третьей степени зависимость ». К сожалению, это не совсем так. Формально действительно трудностей нет (соответствующий полином можно найти в курсах вычислительной математики). Однако он достаточно громоздок и (ввиду его единственности) должен иметь вид: (4.12) , где коэффициенты А,В определяются условиями, чтобы производные этого полинома принимали заданные значения в точках на концах отрезка: , . Естественно, значения , для функции , определяемой формулами (4.7)-(4.8), выписанные в [5] на стр.111, должны быть предварительно вычислены. Для определения А,В получается система двух линейных уравнений и, конечно, она легко решается. Но выписывать формулы мы не будем, поскольку при произвольных значениях и соответствующих им значениях гарантировать монотонность полинома (4.12) на интервале не представляется возможным. А без этого его практическая ценность весьма сомнительна. Для описанного выше варианта линейной интерполяции вопрос о монотонности очевиден. Таким образом, весь расчет в случае сводится к формуле (4.11). Напомним, что такой прием мог быть использован и в предыдущем разделе 4.6. Тогда вместо точной формулы (4.10) возникла бы приближенная, которая получается из (4.11) заменой р2 на р0 =0, т.е. формулой: (4.13) . Как ни удивительно, она должна приближенно заменить (4.10). 4.8. В [5] на стр.111 были обоснованы свойства монотонности и выпуклости каждого из двух слагаемых в левой части уравнения (4.7)-(4.8). Поэтому наряду с описанной линейной интерполяцией естественный интерес представляет также интерполяция в виде параболы, проходящей через точки , , , причем именно в виде обратной зависимости Р( j ). Тогда она определяется формулой: (4.14) Искомое приближенное значение Р , отвечающее значению , будет определяться формулой: (4.15) Эта формула может использоваться на всем интервале значений . С точки зрения точности она предпочтительнее, чем (4.11) и (4.13). Исключение составляют «вырожденные» случаи: если , то предпочтительнее (4.11), если , то предпочтительнее (4.13). Условие можно практически реализовывать, например, как , условие , как . 4.9. Обратимся теперь к ситуации . Тогда и образуются две ударные волны. Положение осложняется тем, что Р может быть сколь угодно велико (при достаточно больших значениях ). Наряду с описанием итерационного процесса решения соответствующего этому случаю уравнения (4.16) , в монографии [5] на стр.113-114 и в [9] на стр.169-170 обсуждался вопрос о замене его приближенным квадратным уравнением. Еще один интересный вариант получения квадратного уравнения предложен в [10] на стр.45. (При его описании будем иметь в виду необходимость обобщения на случай различных значений g 1 слева и g 2 справа). Запишем уравнение (4.16) так: (4.17) , введя вспомогательные обозначения: (4.18) , , Заметим, что множитель Q близок к единице, так как р1 и р2 не превышают Р . Если этот множитель опустить, (4.17) превращается в уравнение (4.19) Оно сводится к квадратному и легко решается. Заметим далее, что точно таким же приемом получается уравнение: (4.20) Ввиду уже упоминавшихся свойств монотонности и выпуклости каждого из двух слагаемых в левой части уравнения (4.16) нетрудно обосновать, что (4.19) и (4.20) дают для его искомого корня приближения снизу и сверху соответственно. Можно было бы, например, использовать это для контроля точности определения корня (и, в зависимости от этого, решать вопрос, стоит ли прибегать к итерационному процессу или можно ограничиться найденным значением). Выберем все же для вычисления приближенного корня уравнение (4.19), поскольку именно оно дает при Р=р2 точное значение , и это важно. После возведения в квадрат, несложных преобразований и дополнительных обозначений (4.21) , , получаем квадратное уравнение: В качестве решения выбирается его больший корень (4.22) . Таким образом, в рассматриваемой ситуации с двумя ударными волнами расчет приближенного значения давления Р состоит в реализации формул (4.18), (4.21) и (4.22). Кроме упомянутой выше возможности замены р1 на р2 , заслуживает внимания с точки зрения точности замена в формуле (4.22) величины р1 на «взвешенное» значение . Тогда формула (4.22) принимает вид: . Безитерационный алгоритм приближенного вычисления давления Р описан для всех ситуаций, которые могут возникнуть при распаде разрыва (если ). Как уже отмечалось, его следует рассматривать как один из примеров таких алгоритмов. В частности, еще один пример описывался в [11] в расчете на более сложную ситуацию, связанную с учетом пористости среды. § 5. Некоторые замечания о практической реализации приближенного расчета. 5.1. После получения значения Р могут быть рассчитаны остальные величины, определяющие возникающую конфигурацию распада разрыва. Как специально отмечалось в [1] на стр.12, не следует «соблазняться» расчетом этих величин по линеаризованным формулам, а следует руководствоваться адиабатой Гюгонио (если ) или адиабатой Пуассона (если ) даже в случае «звукового приближения». Соответствующие формулы приведены в [5] на стр.112-113 и в [9] на стр.171-172. Для случая ударной волны мы не будем их повторять. Отметим лишь следующее важное обстоятельство. Вычисление плотностей по адиабате Гюгонио: (k=1,2) гарантирует возрастание энтропии . По объему вычислений эти формулы не хуже, чем (1.3). Между тем эти последние в случае назначения (4.2) могут энтропию уменьшать, как следует из (4.5). Что касается волны разрежения, то при реализации итерационного процесса мы располагали массовыми скоростями, полученными на последней итерации. Теперь их нет, и предпочтительнее вычисления осуществлять в таком порядке. Прежде всего получаются значения плотностей в зонах 3 и 4: (5.1) , если (k=1,2) Затем – значения скоростей звука в зонах 3 и 4: (5.2) (k=1,2) Величина U определяется, исходя из того, что в случае точного решения должно было бы быть: Например, так (5.3) Границы «веера характеристик», задающие волну разрежения, определяются условиями: (5.4) , , 5.2. В [10] на стр.42 отмечено следующее. Хотя формулы «звукового» распада разрыва получаются для случая, когда параметры газа по обе стороны разрыва близки между собой, они приводят к устойчивому счету не только в этом случае. Несмотря на то, что вычисленные параметры P , U могут оказаться далекими от их истинных значений на соответствующих контактных поверхностях, «звуковое приближение» можно использовать гораздо чаще, чем можно было бы ожидать, исходя из предпосылок, положенных в основу его вывода. Отсюда следует важный для практики вывод, что безразмерный параметр в критерии (4.4) можно назначать не так уж малым, охватывая подавляющую часть ситуаций, возникающих при практических расчетах. Возможно, при расчете некоторых задач описанный «дорогой» алгоритм вообще не будет использоваться. 5.3. Для экономии объема вычислений при программной реализации описанного алгоритма полезно иметь в виду следующее. При получении давления Р после «отказа в звуковом приближении» (4.4), если и , можно, не вычисляя координат «опорных» точек (4.9), сразу обратиться к случаю двух ударных волн в разделе 4.9. Если и , то сначала можно вычислить только F0 и лишь убедившись, что , вычислять F 1 , F 2 и рассматривать оставшиеся случаи. Далее, целью расчета является получение «больших» величин P , U , R , E на единственном луче, отвечающем значению в случае неподвижной сетки (или задаваемому значению , если сетка подвижна). Трудность состоит в том, что необходимо определять, в какой именно из секторов, образуемых конфигурацией распада разрыва, он попадает. Очевидно, что нет необходимости рассчитывать эту конфигурацию полностью . Например, если U <0, нет необходимости рассчитывать «левые» величины и, наоборот, если U >0 – «правые». В первую очередь естественно проверять возможности, требующие меньшего объема вычислений. Самым «трудоемким» является получение величин в случае попадания внутрь сектора, занимаемого волной разрежения (см. [5], стр.120 или [9], стр.171). Заметим, что «протяженность» сектора, охватываемого левой волной разрежения, определяется формулой: Аналогично, для правой волны разрежения: Если эти величины малы, можно внутрь такого сектора и не «влезать», а ограничиться формулами типа линейной интерполяции. Однако в случае сильных волн разрежения выбор соответствующего луча становится существенным. Об этом свидетельствует и опыт расчетов и необходимость довольно громоздкой меры для обеспечения сходимости итерационного процесса, описанной в [4] (см. также [5], стр.110). 5.4. Достоинством описанного приближенного алгоритма является то, что он во всех случаях качественно правильно описывает возникающие конкретные конфигурации распада разрыва. Учитывая, что использоваться будет лишь небольшая часть информации в виде «больших» величин на одном луче, на основе которых будут вычисляться газодинамические потоки, можно надеяться, что и количественно его результаты вполне приемлемы. 5.5. В [1] на стр.18 отмечалась особая роль точного расчета распада разрыва как универсального инструмента для расчета движения границ физических областей, например, разделяющих области с совершенно разными парамет-рами. По нашему мнению, и описанный приближенный безитерационный метод, как правило , вполне может претендовать на такую роль. 5.6. Для успешного исполнения этой роли, а также для практической возможности использования описанного метода не только для уравнения состояния (1.1), но и для более сложных случаев, требуется его модификация. В § 8 работы [1] рассматривались некоторые из возникающих при этом проблем. Рассмотрим важный для практики расчетов случай, когда газодинамические параметры с индексами 1 и 2 удовлетворяют двучленным уравнениям состояния со своими параметрами: (k=1,2) Это отражается в формулах для скорости звука и энтропии: , Соответственно должны быть изменены формулы (1.9)-(1.11) и формулы (4.1), (4.4). Например, последняя должна быть заменена на следующую: , В основном уравнении (4.7)-(4.8) должна быть сделана замена на ; должны заменить c . Это приводит к столь многочисленным и громоздким изменениям в формулах § 4 и § 5, что мы ограничимся только наиболее важными расчетными формулами. Прежде всего это касается формул (4.9) для «опорных» точек, которые должны быть заменены на следующие: Что касается последней из формул (4.9), то в ней Р= 0 должно быть заменено на . Поскольку вакуума достигает первой либо левая (если ), либо правая среда (если ), формула для F0 принимает следующий вид: Применение формулы (4.10), как уже отмечалось в разделе 4.6, в случае , становится невозможным, и ее придется заменить приближенной. Поскольку , такая приближенная формула (4.13) должна быть заменена на аналогичную (4.11): Усложняется уравнение параболы (4.14), в котором добавляется еще одно слагаемое . Соответственно изменяется и формула (4.15): В разделе 4.9 основное «рабочее» уравнение (4.19) должно быть заменено на следующее: Поэтому изменения в расчетных формулах таковы. Формулы (4.18) заменяются на , , а в формуле (4.22) слагаемое заменяем на . Корректируются также формулы (5.1) и (5.2): , Подчеркнем, что перечислены изменения не всех , а только расчетных формул алгоритма для случая . 5.7. В случае сложных уравнений состояния (не двучленного вида) на практике для решения задачи о распаде разрыва широко используется прием замены их локальной аппроксимацией соответствующими двучленными уравнениями состояния (своими для каждой из зон 1 или 2). Именно поэтому такой случай стал предметом специального рассмотрения в предыдущем разделе. Однако, как отмечалось в [1] на стр.25, такой прием не является универсальным. Возможность расчета распада разрыва для сложных уравнений состояния рассматривались, например, в [12] и в [9] на с.175-176. Наконец, уместно упомянуть о достаточно универсальном подходе к решению задачи о распаде разрыва в случае сколь угодно сложных уравнений состояния, изложенном в небольшой монографии [13]. Рассмотрение упомянутых работ выходит за рамки настоящей. § 6. Снова об энтропии. 6.1. Вопрос об энтропийных эффектах и «энтропийных следах», появляющихся в численных расчетах газодинамических задач, издавна был одним из дискуссионных. В качестве примера можно привести цитату из [14]: «По нашему мнению, для расчета адиабатических течений могут применяться лишь такие разностные схемы, в которых изменение энтропии из-за погрешностей аппроксимации уравнений при малых значениях (приращение удельного объема ) хотя бы не превосходит физического изменения энтропии на слабых ударных волнах». Как известно (см., напр., [15]), на слабых ударных волнах приращение энтропии пропорционально кубу приращения удельного объема: . В [14] рассматривались разностные схемы для уравнений газовой динамики в переменных Лагранжа в адиабатическом приближении, т.е. при отсутствии вязкости, тепловых потоков и источников энергии. Проведенные в ней исследования показали, что такие разностные схемы, в которых вычислительные источники производства энтропии при расчете гладких течений по порядку величины не превосходят физических источников, существуют. В этих схемах погрешности аппроксимации суть величины порядка ( - шаг по временной координате). 6.2. Схема С.К.Годунова и в случае лагранжевых, и в случае эйлеровых переменных таковой не является, поскольку для нее погрешность аппроксимации имеет порядок . Именно поэтому и возникает желание порядок аппроксимации повысить. Эту цель ставят перед собой некоторые из авторов, предлагающих схемы «типа Годунова». Поскольку при этом обычно используется задача о распаде разрыва, естественно, чтобы она давала решение, для которого «вычислительный шум», не превосходит создаваемый реальными физическими процессами, например, слабыми ударными волнами. Точный расчет распада разрыва такими свойствами обладает. Это уже неоднократно отмечалось, и мы не будем на этом останавливаться. Что касается описанного приближенного расчета распада разрыва, то приращение энтропии имеет не третий, а второй порядок (см. выше раздел 4.3). Этого вполне достаточно для схемы первого порядка, каковой является метод Годунова, но может оказаться недостаточным, если речь пойдет о достижении более высокого порядка аппроксимации при расчете гладких течений. 6.3. В связи с обсуждаемым вопросом уместно процитировать также стр.12 работы [16]: «При расчетах сильных изэнтропических волн разрежения по моей [С.К.Годунова] схеме наблюдается значительный рост энтропии, который постоянно служил источником неприязни к результатам разностных расчетов. Существенное снижение этого эффекта было замечено только, когда приступили к расчету двумерных задач… Двумерные программы, использующие эйлеровы координаты, приводят к существенно меньшему «паразитическому» повышению энтропии в волнах разрежения. Причина этого явления состоит в том, что лагранжевы координаты (даже в простейшем одномерном случае) непригодны для формулировки понятия «обобщенные решения» уравнений гидродинамики. Эти уравнения допускают при некоторых начальных данных появление полостей, не заполненных средой, - зон вакуума, которые в лагранжевых координатах изображаются наличием у решения особенностей типа дельта-функции. Такие особенности у обобщенных решений обычно не предусматриваются и при разностной аппроксимации очень плохо моделируются…» 6.4. Отметим, что в настоящей работе, как и в [1], обращается внимание на необходимость контроля энтропии с тем, чтобы не допустить неоправданного ее занижения в газодинамических расчетах. В некотором смысле это противоречит «борьбе» за минимизацию энтропийных эффектов , порождаемых разностной схемой. На первый план выдвигаются интересы гарантированного получения физически реализуемых течений. § 7. О сравнении результатов. 7.1. В воспоминаниях [16] на стр.11-12 описывается, как «потребители вычислительной продукции – физики и инженеры – привыкли к идеальному внешнему виду получаемой информации, и это было поводом к неудовольствию, которое они стали испытывать к первым результатам расчетов по разностным схемам. Как это теперь хорошо всем известно, разностные расчеты разрывных решений уравнений газовой динамики изобилуют большим числом проявлений именно «разностной» микроструктуры, которые выглядят как погрешности, искажающие картину течения. При сравнении машинных выдач с «прилизанными» графиками из «ручных» расчетов эти разностные эффекты вызывали у потребителей неприятные впечатления, приводили к утверждениям, что в программе грубые ошибки или расчет выполнен по ошибочным исходным данным…» А вот в [2]-[3] специально подчеркивается, что схема С не имеет осцилляций на разрывах при произвольных физических состояниях газа в соседних ячейках сетки, если назначать массовые скорости по формулам (2.1). Такая позиция их автора создает впечатление, что стремление к «прилизанным» графикам превалирует над получением пусть и визуально «неприятных», но, возможно, более точных и, главное, физически реализуемых результатов. 7.2. В качестве одного из существенных доказательных аргументов в [2]-[3] приводится сравнение результатов расчетов специальных тестовых задач с расчетами по методу Годунова. При этом демонстрируется весьма хорошее их совпадение, по крайней мере, на визуальном уровне. Следует однако отметить, что сравнение производится для профилей величин по пространственной координате при некотором случайно или специально выбранном, достаточно далеком моменте времени. Поскольку речь идет о нестационарной задаче, такого сравнения безусловно недостаточно для надежных выводов о качестве вычислительного алгоритма. С этой точки зрения наиболее «ответственными» являются, например, как раз моменты прихода ударных волн или волн разрежения на объект, формирование процесса их отражения и т.п. А предельные режимы установления могут получаться правильными даже, если процесс их установления не отвечает реальному физическому процессу. (В качестве аналогии можно привести примеры по «искусственному» ускорению сходимости итерационных процессов для других уравнений математической физики). Когда запросы практики требуют точности результатов, исчисляемой процентами и даже меньше, это приобретает особое значение. 7.3. При проводившихся сравнениях результаты расчетов по методу Годунова принимались за эталонные. Между тем, как уже отмечалось в § 3, например, для теста 2 график профиля внутренней энергии весьма существенно отличается от аналитического решения. Какой же это эталон? Для схемы С и метода Годунова демонстрируются весьма близкие результаты (как для профиля внутренней энергии, так и для остальных величин – давления, плотности и скорости). Чем это можно объяснить? Осмысливая эту ситуацию, отметим следующее. Судя по списку цитируемой литературы в [2], под методом Годунова подразумевается первая публикация [4]. В ней использовался итерационный процесс, основанный на двух массовых скоростях (изменявшихся в ходе итерационного процесса). Скорость его сходимости оценивается геометрической прогрессией. В рассматриваемом случае теста 2 с сильными волнами разрежения знаменатель такой прогрессии очень близок к 1. Если число итераций было небольшим (у автора создалось впечатление, что их число было равным 3), то фактически и схема С и то, что названо методом Годунова, были весьма близки по своей реализации. Это, как можно предполагать, и обеспечило близость тех результатов, о которых шла речь выше. 7.4. Заметим, что сам метод Годунова развивался и совершенствовался. В статье [17] было описано непосредственное развитие схемы [4] для численного решения двумерных нестационарных газодинамических задач и методом установления было рассчитано обтекание сферы сверхзвуковым потоком с образованием головной ударной волны. При этом для расчета распада разрыва использовался приближенный метод («звуковое приближение») и упрощенная схема течения – та самая, которая описана выше в § 1, а массовые скорости задавались формулами Но это было почти пятьдесят лет назад. В дальнейшем, в упоминавшейся уже монографии [5], появилось изложение другого итерационного процесса для точного расчета распада разрыва. Он основан на решении уравнения для давления Р методом Ньютона, имеющем квадратичную скорость сходимости, и использовании «настоящей» схемы течения, возникающего при распаде разрыва. Более подробно это описано в [1]. Конечно, были и другие работы, посвященные непосредственно обсуждаемому кругу вопросов. Как уже упоминалось, обзор таких работ можно найти, в частности, в монографиях [5] и [9]. Заключение По признанию многих авторов метод С.К.Годунова [4]-[5] стал одним из популярных для численного решения газодинамических задач. «Массовой операцией» при его реализации является решение задач Римана о распаде разрыва. Точное решение таких задач требует значительного объема вычислительной работы. Естественно поэтому стремление ряда авторов уменьшить эти затраты. Поскольку наиболее «трудоемкая» часть связана с адиабатой Пуассона (что эквивалентно энтропии), у части авторов проявляется тенденция к ее исключению из алгоритмов. Это, в свою очередь, влечет за собой привлечение в алгоритмы таких конструкций, на которых энтропия может убывать. На примере одной из схем, опубликованной в [2]-[3], такая ситуация изучалась в [1] и в настоящей работе. Описан один из возможных безитерационных алгоритмов для приближенного расчета распада разрыва, позволяющий избежать упомянутой опасности убывания энтропии. При практических расчетах газодинамических течений в подавляющем большинстве случаев параметры в соседних ячейках отличаются мало, удовлетворяя критерию (4.4). Заметим, что это может достигаться не только за счет гладкости течений, но при достаточно большом числе ячеек расчетной сетки и «благодаря размазыванию» разрывов при сквозном счете. В этой ситуации для вычисления давления Р достаточно обходиться «звуковым приближением». Вопрос о том, можно ли при этом использовать «звуковые» формулы и для остальных расчетных величин, допуская возможность, хоть и незначительного, но все же убывания энтропии, является, по меньшей мере, дискуссионным. Не может ли систематическое «отравление» (пусть и в незначительных «дозах») убыванием энтропии приводить к нежелательным последствиям? В [1] на стр.15 такая ситуация характеризовалась пословицей: «ложка дегтя может испортить бочку меда». Заметим, что путь, как этого избежать, известен. Что касается вычислительной «стоимости» его реализации, то автор считает целесообразным отметить любопытный нюанс. Формулы метода Годунова были разработаны в 50-х годах для ЭВМ «Стрела» с быстродействием 2000 операций в секунду. А в наш век неизмеримо выросшего быстродействия вычислительных машин и многопроцессорных комплексов столь важна стала «сверхэкономичность»? Стоит ли рисковать потерять главное – физическую достоверность результатов? Не забудем, что наш классик А.С.Пушкин одну из своих известных сказок закончил словами: «Не гонялся б ты, поп, за дешевизной!» В заключение автор считает нужным подчеркнуть еще раз дискуссионный характер обсуждаемого вопроса. Отмечается лишь тревога по поводу использования в схеме С вспомогательной конструкции, которая может терять физическую реальность. Схема С обладает достоинствами, связанными и с существенным сокращением объема вычислений и с упрощением возможностей ее обобщения на случай сложных уравнений состояния. Конечно же, найдутся задачи, в которых все «обойдется» благополучно и для простейшего алгоритма. Но они не могут служить доказательством , что так можно считать всегда. В свою очередь, тревога автора не может считаться обоснованной, пока не будут предъявлены конкретные примеры расчета течений с убыванием энтропии. Их можно пытаться либо конструировать искусственно, либо искать в практических расчетах. Предложение поиска таких примеров – одна из целей настоящей публикации. Автор выражает благодарность М.С.Гавреевой и А.В.Северину за помощь в оформлении работы. Литература 1. Прокопов Г.П. Контроль энтропии в алгоритмах и расчетах газодинамических течений.//Препринт ИПМ им.М.В.Келдыша РАН, 2006, №84, 28 с. 2. Сафронов А.В. Разностный метод для уравнений газодинамики из соотношений на разрывах для вектора потока.// Материалы XVII Школы-семинара «Аэродинамика летательных аппаратов», ЦАГИ, 2006. 3. Сафронов А.В. Разностный метод решения нестационарных уравнений газодинамики на основе соотношений на разрывах.// Космонавтика и ракетостроение, 2006, вып.2(43), с.152-158. 4. Годунов С.К. Разностный метод численного расчета разрывных решений уравнений гидродинамики.// Матем. сб., 1959, 47, вып.3., с.271-306. 5. Численное решение многомерных задач газовой динамики. Под ред. С.К.Годунова.// М., «Наука», 1976, 400 с. 6. Ушаков Д.Н. Большой толковый словарь современного русского языка.// М., Альта-Принт, 2005. 7. Годунов С.К. О единственности решения уравнений газодинамики.// Матем. сб., 1956, т.40, вып.4, с.467-478. 8. Годунов С.К. Проблемы обобщенного решения в теории квазилинейных уравнений и в газовой динамике.// Успехи матем. наук, 1962, 17, вып.3(105), с.147-158. 9. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений//М., ФИЗМАТЛИТ, 2001, 608 с. 10. Матвеев С.К. Некоторые аспекты применения метода Годунова к решению задач нестационарной газовой динамики.// Газодинамика и теплообмен. Ученые записки Ленинградского Университета. Л., 1977, №393. Серия матем. наук. Вып.55, сборник 5, с.41-54. 11. Прокопов Г.П. Расчет распада разрыва для пористых и сплошных сред с двучленными уравнениями состояния.// Вопросы атомной науки и техники. Сер.: Методики и программы численного решения задач матем. фтизики. 1982, вып.2(10), с.32-40. 12. Прокопов Г.П. Задача о распаде разрыва в трехтемпературной газовой динамике. // Препринт ИПМ им.М.В.Келдыша РАН, 2004, №66, 28 с. 13. Алалыкин Г.Б., Годунов С.К., Киреева И.Л., Плинер Л.А. Решение одномерных задач газовой динамики в подвижных сетках.// Изд-во «Наука», Гл.ред.физ.-мат.лит., М., 1970, 112 с. 14. Куропатенко В.Ф. О точности вычисления энтропии в разностных схемах для уравнений газовой динамики.// Численные методы механики сплошной среды. Новосибирск, 1978, т.9. №7, с.49-59. 15. Зельдович Я.Б., Райзер Ю.П. Физика ударных волн.// М., Физматлит, 1963. 16. Годунов С.К. Воспоминания о разностных схемах. Доклад на Международном симпозиуме «Метод Годунова в газовой динамике». Мичиганский университет (США), Май, 1997// Новосибирск. Научная книга, 1997, 40с. 17. Годунов С.К., Забродин А.В., Прокопов Г.П. Разностная схема для двумерных нестационарных задач газовой динамики и расчет обтекания с отошедшей ударной волной.// ЖВМ и МФ, 1961, 1, №6, с.1020-1050.