Содержание

                                                                                                                стр.

Введение

         Целью и содержанием работы является изложение численного метода решения задач о распаде произвольного разрыва, которые в свою очередь составляют основной структурный элемент известного метода С.К.Годунова для численного решения задач гидродинамики сжимаемых сред. Привлекаемые при этом уравнения состояния (УРС) предполагаются заданными в виде таблиц. Оказалось, что в случае линейной аппроксимации УРС на треугольных ячейках все нужные для расчета распада разрыва функции могут быть представлены в явном виде, как для двучленного уравнения состояния.

§ 1. Постановка задачи.

         1.1. В недавно изданной монографии [1] на с.75 изложено следующее. Основное уравнение термодинамики

                                                ,                                                   (1.1)

являющееся обобщенной формой обоих начал термодинамики, содержит пять термодинамических функций, определяющих состояние исследуемой среды: Т – температура, s – энтропия, e - внутренняя энергия единицы массы,  р – давление, v – удельный объем.

         Поскольку (1.1) содержит три дифференциала, можно утверждать, что две из термодинамических функций являются независимыми. Чтобы определить три других – зависимых, необходимо к (1.1) присоединить еще два уравнения. При изучении процессов, зависящих от температуры Т, в качестве таких уравнений могут выбраны, например:

                                   ,                                                    (1.2)

Их называют соответственно термическим и калорическим УРС.

         Вместо удельного объема v зачастую удобнее использовать плотность . Естественно, возможны и используются на практике и другие варианты выбора независимых термодинамических переменных, в том числе, и не входящих в число перечисленных выше.

         1.2. Система нестационарных дифференциальных уравнений газовой динамики обычно основывается на законах изменения массы, количества движения и энергии. Наряду со скоростью среды u (скалярной в одномерном случае или векторной в пространственном) она включает термодинамические величины  (или v). Поскольку, как уже сказано, из них независимыми являются две, система уравнений газовой динамики замыкается заданием дополнительного соотношения (УРС), например, вида:

                                    или                                             (1.3)

Одной из простейших форм (1.3), широко используемой на практике и особенно в методических целях, является так называемая двучленная форма УРС:

                                                                                 (1.4)

Она содержит три числовых параметра: g - показатель адиабаты,  - плотность и скорость звука для «холодной» конденсированной среды. Это позволяет рассматривать «растяжение» среды до конкретного допустимого уровня отрицательных давлений:

                    ,            ,                                     (1.5)

Еще более простым является УРС идеального или политропного газа:

                                                                                                      (1.6)

Оно получается как частный случай (1.4) при .

         Система дифференциальных уравнений газовой динамики, замкнутая заданием УРС, становится источником разнообразных задач, представляющих интерес для практики, если ее дополнить соответствующими начальными и краевыми условиями. Уравнения состояния при этом могут задаваться своими для каждой из геометрических областей, на которые «раскраивается» задача.

         1.3. Основным структурным элементом и массовой операцией широко известного метода С.К.Годунова (см., напр., [2]) является решение задач о распаде разрыва. Эта задача настолько хорошо известна, что мы ограничимся лишь кратким описанием. Речь идет о расчете конфигурации, которая возникает при взаимодействии (соприкосновении) двух сред с параметрами  слева и  справа.

         Знакомство с обширной литературой (см., напр., [2]-[3]) обнаруживает, что практический алгоритм решения этой задачи, относительно простой и пригодный для массового счета, разработан для двучленного УРС (1.4). Но и в этом случае задача сводится к некоторому нелинейному уравнению для определения давления Р на контактном разрыве, общего для двух сред:

                                                       (1.7)

Для его решения зачастую или иногда приходится привлекать итерационный процесс (например, метод касательных Ньютона). Поэтому даже в этом относительно простом случае УРС (1.4) или (1.6) обнаруживается стремление при практических расчетах решать задачу о распаде разрыва не точно, а лишь приближенно. Иногда такого рода приближенные (а потому менее затратные) алгоритмы дают основания для серьезных сомнений (см. по этому поводу, напр., [4]).

         1.4. Сложнее обстоит дело при задании УРС в форме (1.3). Тогда при конструировании упомянутого уравнения (1.7) для давления Р на контактном разрыве возникают следующие ситуации.

         В случае ударной волны (, равного  или ), кроме соответствующего УРС, связывающего «большие» величины (терминология из [2]) на контактном разрыве, привлекается ударная адиабата Гюгонио:

                                                                    (1.8)

В случае волны разрежения () вместо нее должно конструироваться гладкое решение на основе характеристик и соотношений на них:

                            вдоль () – характеристики,                        (1.9)

                           вдоль линии тока.                                            (1.10)

Скорость звука с в случае УРС (1.3) определяется формулой:

                                                                        (1.11)

или                                                                    (1.12)

Она должна давать положительный результат, чтобы обеспечивать гиперболичность системы уравнений газовой динамики с УРС (1.3). Заметим, что для обеспечения единственности решения задачи о распаде разрыва на функцию (1.3) накладываются более жесткие ограничения, известные как условия выпуклости и Бете-Вейля. Мы их сейчас обсуждать не будем.

         Проинтегрировать (1.9) в явном виде и получить формулы для соответствующих инвариантов Римана, как правило, не удается. Поэтому при работе с (1.9)  в случае сложных УРС приходится прибегать к численному интегрированию.

         С соотношением (1.10) дело обстоит так. Если из термодинамического тождества (1.1) для УРС (1.3) удается получить явную формулу для энтропии s или произвольной монотонной возрастающей функции , то (1.10) принимает форму изэнтропы:

                          ds=0 ,    т.е.  s=const  вдоль линии тока.

         Именно так будет в случае двучленного УРС (1.4), для которого

                            ,            .                     (1.13)

В противном  случае для реализации (1.10) придется также прибегнуть к численному интегрированию.

         1.5. Теперь перейдем к обсуждению ситуации, более общей, чем УРС в форме (1.3). Существует много вариантов параметрического задания УРС. В общем виде такой подход и его реализация для задачи о распаде разрыва был разработан, например, в [5]. Из-за громоздкости и трудоемкости он не был оценен должным образом и не получил широкого распространения. Главным практическим препятствием является то, что предполагается наличие аналитических формул для функций, задающих УРС, и их производных. Аналогично обстоит дело и с алгоритмами, предложенными для решения задачи о распаде разрыва со сложными УРС в других работах и реализованными в специальных программных комплексах (см., напр., [6]).

         При описании газодинамических течений, зависящих от температуры Т, целесообразно задавать УРС соотношениями (1.2). Фактически это неявное задание зависимости вида (1.3). Исходя из (1.2), будем иметь:

                            ,           

Следовательно,    

или           

Для обсуждаемой неявной функции  будем иметь:

                         ,                               (1.14)

Во избежание путаницы, не вводя новых обозначений для рассматриваемых функций, в формулах (1.14) указаны аргументы (v,e)   для неявной функции, а аргументы Т,v для параметрического задания (1.2) опущены.

         Полученные формулы (1.14) позволяют реализовать описанную выше технологию конструирования уравнения (1.7) для расчета распада разрыва в переменных Т,v.

         1.6. Вид (1.2) является исходным при построении широкодиапазонных УРС (см., напр., [7]). Как уже отмечено, он удобен и с математической точки зрения для замыкания системы уравнений газовой динамики при описании течений, зависящих от температуры Т.

         Важно сразу отметить, что функции (1.2) с точки зрения термодинамики не являются независимыми. Условие интегрируемости, вытекающее из тождества (1.1), налагает на них ограничение (см., напр., [8], с.139):

                                                                                     (1.15)

Оно называется условием термодинамической согласованности. В случае его невыполнения для УРС (1.2) не может быть определена энтропия со всеми вытекающими отсюда последствиями.

         Чтобы гарантировать термодинамическую согласованность, при конструировании УРС обычно привлекаются так называемые термодинамические потенциалы (см., напр., [9], с.60-67). Как правило, они имеют достаточно громоздкий вид. Конкретные модели конструирования таких потенциалов обычно рассчитаны не на весь широкий диапазон изменения независимых параметров. Поэтому встает еще практическая проблема «сшивки» аналитических выражений для таких потенциалов или самих термодинамических функций.

         Кроме того, следует принимать во внимание, что, как известно из курсов математического анализа, для существования и единственности решения системы уравнений (1.2), т.е. однозначного вычисления v,T (или r,Т) по значениям р,e , определяющим условием является необращение в нуль якобиана

                                                                                  (1.16)

Далее следует отметить еще одну важную практическую проблему (см., напр. [7], с.167). «Непосредственное применение широкодиапазонного УРС в его исходной аналитической форме затруднительно в гидродинамических расчетах из-за больших затрат машинного времени для нахождения параметров течения в каждой пространственно-временной точке. Поэтому УРС может использоваться в форме таблиц».

         Именно такой универсальный подход к заданию УРС и будет предметом исследования. Естественно, что сразу же возникает вопрос о восполнении табличных данных посредством некоторых аналитических формул, необходимых, в частности, для реализации алгоритма решения задачи о распаде разрыва.

         1.7. В самом общем виде, безотносительно к этой задаче, технология восполнения таблиц УРС рассматривалась, например, в [10], с учетом или без учета условия термодинамической согласованности (1.15).

         Для этой цели используются полиномиальные бисплайны различных типов с весьма громоздкой практической реализацией. Вопрос о невырождении якобиана (1.16) просто не обсуждается. Поэтому практическое использование такой технологии восполнения таблиц УРС для решения задач о распаде разрыва вряд ли целесообразно.

         Наиболее простой вариант восполнения таблиц УРС представляет билинейная аппроксимация, предложенная в [11]. Она не удовлетворяет, как правило, условию термодинамической согласованности и поэтому была подвергнута критике в [12]-[13]. Автором были предложены алгоритмы, позволяющие, в принципе, устранить этот недостаток. Однако эти алгоритмы оказались, во-первых, хотя и реализуемыми, но слишком громоздкими (что было очевидно и отмечено сразу). Во-вторых,  они содержат существенный элемент неопределенности, связанный с неоднозначностью выбора базисных аппроксимирующих функций. В-третьих, как отмечено в [12] на с.8, при их практической реализации можно столкнуться с «эффектом лоскутного покрытия». Суть его состоит в том, что конструирование аппроксимационных формул осуществляется в табличных ячейках независимо. Если в процессе решения задачи о распаде разрыва приходится переходить в соседнюю ячейку, это может вызвать определенные трудности.

         Учитывая эти обстоятельства, а также тот факт, что речь идет о приближенном решении задачи о распаде разрыва (а на большее и не приходится претендовать, поскольку УРС передается таблицами приближенно), сочтено целесообразным обратиться именно к простейшей билинейной аппроксимации, предложенной в [11]. Будем исходить из того, что задаваемые таблицы УРС получаются на основе функций, которые условию термодинамической согласованности удовлетворяют, и достаточно подробны, чтобы его «передать». Если они гарантируют невырождение соответствующих якобианов для табличных ячеек, то использование билинейных формул, по-видимому, можно считать вполне приемлемым.

          1.8. Мы намерены в максимальной степени использовать описание алгоритма из [11], обратив основное внимание на принципиальные вопросы, связанные с его реализацией, и предложив некоторые изменения в этом алгоритме.

         При анализе этих изменений обнаружилось следующее. Если решиться на еще более простой вариант аппроксимации УРС, а именно на линейную аппроксимацию в треугольных табличных ячейках, то удается получить алгоритм расчета распада разрыва, не требующий численного интегрирования. Его описанию будет посвящен §4.

§ 2. Табличная реализация УРС.

         2.1. Будем рассматривать случай, когда УРС задается в виде (1.2), именно в виде зависимостей от температуры Т и удельного объема v, а не плотности r. Преимущества такого выбора выяснятся чуть позже. Так же, как в [11], разделе 2, будем (1.2) задавать в табличной форме.

         Для этой цели задаются одномерные монотонно возрастающие последовательности: температуры  и плотности . При замене плотности r на удельный объем  возрастающая последовательность  превращается в монотонно убывающую . Сохраним такую ее нумерацию, что позволит в максимальной степени сохранить описание алгоритма из [11].

         В частности, будем предполагать, что заданы две двумерные таблицы  и , удовлетворяющие требованиям, изложенным в [11]. А именно, таблица  монотонно возрастает по обоим индексам, а таблица  монотонно возрастает по первому индексу. Правомерность таких требований обсудим позже в §3.

         2.2. Область на плоскости независимых переменных

                    ,         ,        ()                         (2.1)

назовем табличной ячейкой с индексами i,j .

Внутри каждой ячейки функции (1.2) аппроксимируем билинейными функциями вида:

                                     ,                                     (2.2)

                                    ,

которые принимают в четырех угловых точках ячейки заданные значения  ,  из таблиц. Тогда коэффициенты в (2.2), свои в каждой табличной ячейке, определяются следующими формулами:

                                                                                           (2.3)

Очевидно, что коэффициенты В₀, В₁, В₂, В₃ определяются аналогичными формулами с заменой на .

         Мы не сочли целесообразным переходить к «сдвинутым» переменным, как это было сделано в [11], чтобы не вводить новых переменных и сохранить большую наглядность формул для последующего изложения. Возможно, при практической реализации «сдвинутые» переменные (причем и «нормированные» шагами ячейки  , ) заслуживают внимания.

         Воспользуемся также случаем обратить внимание, что в [12] по небрежности автора на с.5-6 в формулах (1.3)-(1.4) для коэффициентов при Т и r пропущены слагаемые, порождаемые «сдвигом» переменных. Правда, это не имело никаких последствий для последующего изложения.

         2.3. Рассмотрение алгоритма расчета распада разрыва начнем со случая ударной волны (см. [11], с.790). Уравнения состояния (1.2) в билинейной форме (2.2) будем рассматривать совместно с ударной адиабатой

                                                                    (2.4)

Их нужно решить относительно v и Т при заданном р=Р и , равных  или , характеризующих левую и правую среду в решаемой задаче о распаде разрыва. Следовательно, для конкретной табличной ячейки (со своими коэффициентами билинейной аппроксимации) речь идет об уравнениях:

Оба они линейны относительно V и Т, и из них легко получаем:

                                                (2.5)

или

                       (2.6)

Следовательно, можно из (2.5) сразу найти Т, вычислив корни квадратного уравнения, а затем досчитать V. Либо из (2.6) найти V из своего квадратного уравнения, а затем досчитать Т.

         Найденное решение будет «правомочным», если удовлетворяет условиям (2.1) «попадания» в ячейку. Это достигается для одного из двух корней квадратного уравнения и только в том случае, если табличная ячейка выбрана правильно. Более подробно об этом – ниже, в §3.

         Последовательность перебора ячеек таблицы по индексам i,j подробно описана в [11] на с.790-791, и мы не будем ее повторять. Заметим, что в описанной реализации, использующей удельный объем v вместо плотности r, алгоритм для случая ударной волны существенно упрощается. В том числе, упрощается и процедура перебора ячеек, поскольку решение (2.5) или (2.6) дает «наводящую» информацию сразу по обоим индексам i,j, даже если ячейка таблицы выбрана неправильно (если, конечно, не окажется, что корни квадратного уравнения – комплексные).

         Остановимся сразу на вопросе о функции , используемой при конструировании уравнения (1.7), из которого должна быть определена величина Р.

         В случае ударной волны () она имеет вид, описанный в [11] на с.788:

                                    ,                            (2.7)

где функция V(P) определяется посредством (2.5) или (2.6). Если для решения (1.7) будет применяться метод касательных Ньютона, то для его реализации требуется производная . Для нее легко получается формула:

Для вычисления  необходимо выписать на основании (2.6) формулу для V(P) и продифференцировать ее. Чтобы не загромождать изложение, мы не будем описывать эту элементарную процедуру.

         2.4. Теперь перейдем к случаю волны разрежения (см. [11], с.789-790).

         Из соотношения (1.10) вдоль линии тока  с привлечением УРС в форме билинейной аппроксимации (2.2) получаем:

                            (2.8)

Это дает линейное дифференциальное уравнение:

                                                ,                                        (2.9)

                  ,             .

Оно интегрируется в квадратурах (см., напр., [14]).

Его решение такое, что , имеет вид:

                              ,                        (2.10)

Если , то          .

Тогда для m(V) получаем формулу:

            ,                            (2.11)

В случае ,  получаем:

Можно выписать для этого случая и явную формулу для T(V). Не будем этого делать, чтобы не загромождать изложения.

         Наконец, в случае ,  непосредственно из (2.8) получается:

                                                                                     (2.12)

Из (2.12) видно, что есть и совсем уж вырожденный случай, если и . В таком случае на всей табличной ячейке p=const=A₀, e=const=B₀.

Это противоречит оговоренным условиям монотонности таблиц. Однако в практике расчетов необходимость работы с таблицами, содержащими такие ячейки, вполне может возникнуть, и в алгоритме требуется предусмотреть и эту возможность.

         Полученные формулы должны заменить аналогичные формулы из [11] на с.789. Описанные случаи их «вырождения», конечно же, не являются следствием замены переменной r на v. Просто соответствующие случаи «остались за кадром» в [11], по-видимому, чтобы не загромождать изложения.

         2.5. Остановимся теперь на вопросе о функции , используемой при конструировании уравнения (1.7). В случае волны разрежения

                                        ,                                  (2.13)

где  - массовая скорость. Для билинейной аппроксимации УРС (2.2) формула скорости звука (1.12) с учетом формул (1.14) приводит к следующему результату:

            .        (2.14)

Она аналогична представленной в [11] на с.790.

         Главной трудностью алгоритма, как и в [11], является реализация формул (2.10) и (2.13), требующих численного интегрирования. Соответствующий алгоритм достаточно подробно описан в [11], и мы не будем его повторять.

         Напомним только, что в ходе его реализации может происходить выход на границу соседней табличной ячейки и последующий переход в новую соответствующую ячейку. Это требует снова реализации билинейной аппроксимации (2.2)-(2.3) и работы с ней по такому же алгоритму.

         2.6. Столь громоздкий вычислительный алгоритм (требующий и соответствующих вычислительных затрат) привлекает естественный интерес к так называемому акустическому приближению. Его можно использовать, если окажется, что

                                                                                 (2.15)

с некоторым заданным достаточно малым .

Тогда имеет место асимптотика

               ,            ,                   (2.16)

одинаковая для ударной волны и волны разрежения. Это не случайно хотя бы потому, что ударную адиабату (2.4) можно рассматривать как дискретизацию  .

         Для ударной волны это позволяет сразу по формулам (2.5) или (2.6) получить все необходимые величины. Заманчивой представляется возможность в случае акустического приближения (2.15) воспользоваться этими формулами и для волны разрежения. Однако автор считает необходимым обратить внимание на отмеченный в [15] факт, что это допускает возможность, хоть и незначительного, но все же убывания энтропии. Вопрос о том, не может ли систематическое «отравление» (пусть и в незначительных «дозах») убыванием энтропии приводить к нежелательным последствиям (физически нереализуемому течению), остается открытым и дискуссионным. Тем более неприемлема замена разрывом волны разрежения, не удовлетворяющей условию акустического приближения (2.15).

§ 3. О требованиях к таблицам УРС.

         Как уже отмечалось в начале §2, в [11] предполагалось, что таблица давления  монотонно возрастает по обоим индексам (т.е. по температуре Т и плотности r), а таблица  монотонно возрастает по первому индексу (по температуре Т):

                               ,          ,                                  (3.1)

если        ,       (т.е. )

         Обратимся к уже упоминавшейся монографии [1], с.82-85.

         Говорят, что система находится в термодинамическом равновесии, если ее состояние не меняется с течением времени.

         Поскольку термодинамика изучает свойства равновесных систем, то УРС могут описывать изменение состояния вещества со временем только как смену равновесных состояний. Это ограничение носит принципиальный характер ([1], с.92).

         Равновесное состояние может быть устойчивым или неустойчивым. В термодинамике доказывается, что в состоянии устойчивого равновесия свободная энергия Гиббса минимальна. Исходя из этого, получается, что (см. [1], с.84-85):

                                       ,                                            (3.2)

Далее (см. [1], с.79), практически для всех известных веществ коэффициент объемного расширения положителен:

                                                                                               (3.3)

a - относительное изменение объема тела при изменении его температуры на один градус при постоянном давлении. Такие вещества называют нормальными и для них (см. [1], с.90):

                                                                                                         (3.4)

Совокупность условий (3.2) и (3.4) и является основанием для упомянутых выше требований к таблицам (3.1).

         Отметим, что в силу условия термодинамической согласованности (1.15) знакоопределенность производной  не предполагается. Например, для идеального газа .

         3.2. Обратимся к рассмотрению отдельной табличной ячейки на плоскости переменных (р,e).

         Производные билинейных функций (2.2) по переменным Т и v являются линейными функциями:

                                  ,                                              (3.5)

                                   ,        

Легко видеть, что якобиан отображения

        (3.6)

является линейной функцией переменных Т,v.

         Для упрощения описания присвоим четырем угловым точкам ячейки более простые номера:

                           ,             ,                            (3.7)

                           ,     

Целесообразно также доопределить вспомогательные точки:

                           ,                                              (3.8)

Обозначим  значения якобиана J в четырех вершинах ячейки. Нетрудно убедиться, что величина  с точностью до коэффициента () представляет удвоенную площадь  ориентированного треугольника при вершине ячейки с номером k после проведения одной из ее диагоналей (см. рис.1а-1б).

                                                  (3.9)

Единообразная формула  достигается с помощью (3.8).

         В силу линейности якобиана J по переменным Т,v, если  для всех вершин ячейки , то это гарантирует положительность якобиана J на всей ее «территории» (2.1). Аналогично, если  для всех , то это гарантирует отрицательность J на всей ячейке. В обоих этих случаях четырехугольная ячейка на плоскости (р,e) будет выпуклой. В одном случае обход ее вершин совершается против часовой стрелки (см. рис.1а), в другом – по часовой (рис.1б).

         В остальных случаях – если не все четыре  имеют одинаковые знаки – четырехугольная ячейка становится либо невыпуклой (см. рис.1в и 1г), либо самопересекающейся (см. рис.1д и 1е). При этом якобиан J будет обращаться в нуль внутри ячейки.

         Требования (3.1) для ячейки (3.7) приобретают вид:

              ,          ,          ,                        (3.10)

         На рис.1а-е проиллюстрировано, что эти условия не препятствуют возможности реализации любой из перечисленных ситуаций.

         3.3. Заметим, что требование невырожденности произвольной регулярной двумерной сетки в описанном виде появилось впервые в [16]. Его использование открывает возможность визуального и аналитического исследования дискретного отображения, порождаемого таблицами , задающими УРС вида (1.2).

         Построим два семейства линий сетки, одно из которых отвечает постоянным значениям индекса i, т.е. заданным значениям температуры , - семейство изотерм, а второе – постоянным значениям индекса j, т.е. заданным значениям  (или плотности ), - семейство изохор.

         Если эти два семейства образуют такую сетку на плоскости (р,e), что все ячейки сетки являются выпуклыми, то с такими таблицами можно работать. В частности, квадратные уравнения (2.5) или (2.6) в случае выбора правильной ячейки дают единственное решение, удовлетворяющее условиям (2.1).

3.4. Присутствие в сетке самопересекающихся ячеек («бантиков») является неприемлемым, так как в точке самопересечения (и не только в ней) нарушается однозначность определения Т,v по заданным р,e.  Квадратные уравнения (2.5) или (2.6) даже для правильно выбранной ячейки могут иметь либо комплексные корни, либо неприемлемые с точки зрения физического смысла. Для таких ячеек якобиан  имеет один знак в двух вершинах и противоположный – в двух других.

         В тех случаях, когда якобиан  имеет один знак в трех вершинах и противоположный – в одной из вершин, возникает невыпуклая ячейка. В этом случае якобиан  будет положительным на некоторой части ячейки и отрицательным на остальной.

         Оставлять такую ситуацию без внимания нельзя, поскольку в этом случае квадратное уравнение (2.5) или (2.6) для правильно выбранной ячейки может иметь два решения, удовлетворяющих нужным условиям (2.1). Эта неоднозначность является недопустимой с точки зрения корректности задачи о распаде разрыва.

         В качестве одного из возможных подходов, позволяющих устранить такой недостаток, рассмотрим следующий. Разрежем невыпуклый четырехугольник с помощью «хорошей» (лежащей внутри ячейки) диагонали на два треугольника.  Для  рис.1в-1г  она  должна  соединять  точки  2 и 4.

С учетом требований (3.10) это – типичная ситуация.

         Будем аппроксимировать УРС по табличным данным отдельно для каждого из двух треугольников. Результаты такого подхода излагаются ниже.

§ 4. Реализация табличных УРС в треугольных ячейках.

         4.1. Рассмотрим треугольную ячейку, в вершинах которой, имеющих на плоскости независимых переменных Т,v координаты

                                  ,         ,                                           (4.1)

заданы значения давления и внутренней энергии

                                  ,        ,         .                                 (4.2)

Они отбираются из величин (3.7)-(3.8) по номеру , назначаемому алгоритмом, который будет описан позже.

         Представляется естественным аппроксимировать УРС независимыми линейными формулами:

                                                                                      (4.3)

Их коэффициенты однозначно определяются информацией (4.1)-(4.2). Расчетные формулы для них будут выписаны позже. После их вычисления (4.3) можно рассматривать как частный случай (2.2) при .

         4.2. В соответствии с изложенным в §2, в случае ударной волны формулы (2.6) принимают вид:

                       (4.4)

Эти уравнения позволяют однозначно определить V, а затем и Т:

                                                      (4.5)

         В случае необходимости может быть вычислена и внутренняя энергия по уравнению состояния (4.3). Для функции  сохраняется прежний вид (2.7).

         4.3. В случае волны разрежения линейное дифференциальное уравнение (2.9) при  приобретает вид:

                                         .                                  

Его решение можно было бы выписать не в квадратурах, как в §2, а в явном виде. Однако еще проще и эффективнее следующий путь.

         Исключаем температуру Т с помощью формулы

и получаем уравнение:

                                                                           (4.6)

Его решение такое, что  , дается формулой:            

                                     ,                      (4.7)

если ввести вспомогательную величину

                                                                (4.8)

Можно обратить (4.7) и получить явное выражение V через Р:

                                        .                                                (4.9)

         4.4. И наконец, самое замечательное состоит в том, что формула (2.14) для массовой скорости при  приобретает вид:

                              .                                   (4.10)

Это позволяет и в случае волны разрежения из формулы (2.13) получить явное выражение для функции

       (4.11)

         Полученные явные формулы для функций  и  позволяют выписать уравнение для определения давления Р на контактном разрыве, аналогичное (1.7) для двучленного УРС (1.4).

         4.5. Обратим внимание на одно серьезное возможное осложнение. Полученные результаты можно считать «правомочными» только в том случае, если рассчитанная точка со значениями термодинамических параметров Р,Т,V не выходит за пределы той же треугольной табличной ячейки, для которой строилась аппроксимация УРС (4.3).

         А если это не так? Начнем с более сложного случая волны разрежения. Тогда, аналогично описанному ранее, в §2, процессу интегрирования уравнения (2.9), реализуемого с помощью квадратурной формулы (2.10), нужно «пройти по маршруту», предписываемому уравнением (4.6). Ситуация облегчается тем, что его решение выписывается явной формулой (4.7). Тем не менее, в случае, если нужная табличная ячейка окажется удаленной от исходной, могут потребоваться переходы в соседние табличные ячейки с перерасчетом коэффициентов аппроксимации (4.3). При этом вопрос о месте «выхода» линии тока в соседнюю ячейку решается в ходе процесса. Если такой «выход» происходит через ребро ячейки v=const, то можно вычислить соответствующее значение Р, используя формулу (4.7). В других ситуациях (например, при «выходе» через диагональ ячейки) ситуация сложнее, так как придется решать итерациями некоторое уравнение, но с явными функциями. Чтобы не загромождать изложение, соответствующие формулы не выписываются. Одновременно в ходе процесса «накапливается» и функция (4.11). Учитывая, что в [11] уже описывались детали процесса численного интегрирования для более сложной ситуации, ограничимся этим описанием.

         В случае ударной волны ситуация несколько проще. Полученные по явным формулам (4.4)-(4.5) значения P,V,T позволяют сразу получить «наводку» (т.е. номера i,j)  на правильную табличную ячейку. Пересчитав для нее аппроксимационные коэффициенты (4.3), повторяем расчет по формулам (4.4)-(4.5). Если «наводка» оказалась неточной, процедуру придется повторять, пока не придем в правильную ячейку. После этого определяется  по формуле (2.7).

         4.6. Для полноты изложения алгоритма вернемся теперь к вопросу о формулах, по которым должны вычисляться коэффициенты аппроксимации (4.3). Для давления  определяются системой уравнений:

                                                                                         (4.12)

Ее решение можно выписать наиболее компактно, по-видимому, так. Исключаем из (4.12) коэффициент :

                                   .

Отсюда получаем определитель

и формулы для коэффициентов  и :

                     .

После этого вычисляем :

                                            .

         Коэффициенты  вычисляются по таким же формулам с заменой  соответственно на .

         Эти формулы могли бы быть несколько упрощены с учетом того, что координаты (4.1) суть углы прямоугольной ячейки. Однако это не сделано сознательно. Во-первых, такая общая форма упрощает программную реализацию алгоритма. Во-вторых, общая форма допускает возможность реализации табличного задания УРС не только на прямоугольных сетках, но и на произвольных регулярных (и даже нерегулярных) сетках. Такая возможность заслуживает особого внимания.

         4.7. Как было отмечено в конце предыдущего §3, описанная реализация табличного УРС на основе линейной аппроксимации (4.3) была задумана для невыпуклых ячеек посредством их разрезания по «хорошей» (лежащей внутри ячейки) диагонали на две треугольные ячейки. Полученный результат, реализуемый в столь законченном виде, естественно ставит вопрос: а почему бы не воспользоваться такой аппроксимацией УРС и для выпуклых ячеек?

         Для этого, прежде всего, нужно принять решение по следующему вопросу. Выпуклую ячейку можно разрезать на треугольные ячейки двумя способами. После отбора «правильной» (или «подозрительной») табличной ячейки какой из треугольников выбрать для конструирования приближенного УРС?

         Чтобы исключить возможность получить после разрезания выпуклой ячейки треугольник, «почти» вырождающийся в отрезок, что нежелательно с точки зрения вычисления коэффициентов (4.3), а также с учетом (3.10), разрезаем ячейку по диагонали, соединяющей вершины 2 и 4. (В случае необходимости разрезание по диагонали, соединяющей вершины 1 и 3, рассматривается аналогично).

         Уравнение диагонали для табличной ячейки на плоскости p,v имеет вид:

Определим функцию

                                               (4.13)

Точка (P,V) попадает в треугольник с номером , если результат подстановки ее координат в (4.13) имеет такой же знак, как результат подстановки координат вершины (), т.е.

если           ,    то  ;

в противном случае полагаем .

         В случае реализации в описанной процедуре знаков равенства при назначении  имеет смысл учитывать его «текущее» значение, чтобы не делать лишних пересчетов коэффициентов (4.3).

         Треугольник с назначенным номером  и будет использоваться для отбора величин координат (4.2) из набора (3.7)-(3.8) и соответствующих им значений (4.1).

         Обратим внимание, что сформулированные правила назначения  автоматически срабатывают и для невыпуклых (но не самопересекающихся) ячеек.

         Для самопересекающейся ячейки «хорошей» диагонали не существует, так как обе они расположены вне ячейки. Поэтому присутствие таких ячеек и считается неприемлемым.

§ 5. Алгоритм расчета распада разрыва

для «треугольной»  аппроксимации УРС.

         5.1. Рассматривается одномерная задача о распаде разрыва. Состояния среды (или двух различных сред) слева и справа от разрыва определяются скоростью u и двумя термодинамическими переменными. В качестве одной из этих переменных выбирается плотность r или удельный объем . В качестве другой переменной выбирается давление р. Такой выбор обусловлен тем, что таблица давления  предполагается монотонно возрастающей по обоим индексам i,j, если последовательности  и  являются возрастающими. Это существенно упрощает и делает однозначным выбор табличной ячейки, которая содержит задаваемую точку () или ().

         Таблица энергии  в тех же узлах предполагается монотонно возрастающей только по первому индексу i.

         Естественно, что для каждой из двух сред могут быть заданы свои таблицы  и . Более того, сейчас мы рассматриваем только такой простейший случай. В практике расчетов приходится реализовывать более сложные ситуации, связанные с наличием в реально решаемых задачах областей, характеризующихся как среды с разными УРС, учетом фазовых переходов и т.п.

         5.2. Итак, в качестве входной информации для задачи о распаде разрыва задаются величины:

                        ,        ,      

Индексом 1 отмечаются величины слева, индексом 2 – справа от разрыва.

         Предполагается, что существует решение этой задачи в виде двух волн и контактного разрыва (КР) между ними. Каждая из волн может быть либо ударной волной (УВ), либо центрированной волной разрежения (ВР).

         Предполагается заданным параметр  (), который будет использоваться в качестве критерия «акустического приближения» и, возможно, в некоторых других аналогичных целях.

         5.3. Реализация алгоритма начинается с того, что для каждой из двух сред проводится следующий расчет.

1. Для заданных значений  по своей таблице  отыскивается табличная ячейка i,j, в которую попадает эта точка, а затем и номер треугольника , как описано в конце §4.

2. Вычисляются коэффициенты аппроксимации (4.3) для УРС в найденной треугольной ячейке согласно описанию в §4.

3. С помощью уравнений (4.3) могут быть вычислены остальные параметры  , .

4. По формуле (4.10) вычисляются массовая скорость  и скорость звука .

5. Аналогичный расчет выполняется для среды 2.

6. Из уравнения для «акустического приближения»

вычисляется исходное значение Р=Р_ак.

Следует принимать во внимание, что это значение может выйти за пределы таблицы или из области допустимых значений (например, при значениях , меньших некоторого «критического» уровня, может возникать зона вакуума). Поэтому возможна корректировка, т.е. замена Р_ак  на некоторое значение Р_а.

7. Независимо для каждой из сред 1 и 2 формулы (4.5)-(4.4) для УВ или (4.9)-(4.8) для ВР дают возможность определить величины , . Величина  досчитывается с помощью УРС (4.3).

8. Проверяем, расположены ли полученные точки  и  в рассматривавшихся («своих») треугольных ячейках, привлекая алгоритм, описанный в конце §4.

9а. Пусть это так. В частности, с большой вероятностью можно ожидать, что так будет в случае выполнения условия (2.15) «акустического приближения». Тогда расчет величин (P,V,T) заканчивается. Если условие (2.15) не выполняется, переходим к п.10.

9б. Если  в п.8 оказалось, что полученные точки  и (или)  вышли за пределы своих треугольных ячеек, необходимо отыскать по соответствующим таблицам  и перейти в табличные ячейки, в которые они попали, а затем определить номера  нужных треугольников.

         После этого, так же как в п.2, вычисляются аппроксимационные коэффициенты УРС (4.3) для этих новых ячеек.

10. Как описано в §4, конструируется уравнение для давления Р на контактном разрыве. Затем оно решается приближенно, как будет описано в следующем §6. Полученное решение обозначим Р_в. Возможно, оно потребует корректировки, описанной в п.6.

11. Так же, как в п.7, независимо для каждой из сред 1 и 2 по формулам (4.5)-(4.4) для УВ или (4.9)-(4.8) для ВР определяются величины , , а затем .

12. Как в п.8, проверяем, расположены ли полученные точки  и  в «своих» треугольных ячейках. Если это не так, то реализуется п.9б.

13. Далее возможны два варианта действий. Либо расчет распада разрыва завершается получением приближенного решения . Либо предполагается продолжать итерационный процесс до тех пор, пока два последовательных приближения для Р не окажутся достаточно близкими.

         Такая необходимость может возникнуть, например, в случае, если задача о распаде разрыва решается на границе двух сред, характеризующихся сильно различающимися значениями параметров.

         5.4. Несмотря на сложность описанного алгоритма, следует отметить, что алгоритм, описанный в [11], еще сложнее. Ему приходится оперировать не с явно выписанными формулами для функций, а с функциями, получаемыми посредством численного интегрирования.

         В этом отношении некоторые преимущества перед алгоритмом, описанным в [11], по крайней мере, для случая ударной волны, имеет и алгоритм, описанный выше в §2. Однако и для его реализации в случае волны разрежения приходится прибегать к численному интегрированию. Поэтому алгоритм из §2 (назовем его для краткости «билинейным») имеет более сложную реализацию, чем описанный в настоящем параграфе (условно назовем его «линейным»).

         Стоит однако отметить одно обстоятельство. «Билинейный» алгоритм из §2 оперирует с прямоугольными табличными ячейками УРС, а «линейный» - с треугольными, получаемыми проведением в них одной из двух диагоналей. Поэтому вероятность выхода за пределы основной (исходной) табличной ячейки для «линейного» алгоритма выше, чем для «билинейного». А это влечет за собой процедуру смены аппроксимационных формул и повторения расчета. Ответ на вопрос, «съедает» ли это преимущества «линейного» алгоритма, может дать только достаточно серьезный практический опыт.

         Стоит также отметить, что благодаря локальной аппроксимации УРС в отдельных ячейках билинейными (2.2) или линейными функциями (4.3) на их общих ребрах достигается непрерывность реализуемых аппроксимирующих УРС’ов. Но не более того. В этом отношении билинейная аппроксимация (2.2) имеет некоторое преимущество перед линейной (4.3), которая непрерывную «сшивку» (с возможным разрывом первых производных УРС) реализует не только на ребрах четырехугольной табличной ячейки, но и на одной из диагоналей такой ячейки.

         Высказанные соображения и стали причиной того, что представлены оба варианта алгоритмов. В какой степени они являются новыми? Обнаружив возможности, описанные выше, для линейной аппроксимации (4.3), автор предпринял некоторые поиски. Однако отыскал лишь в [7] на с.167 краткое изложение использования линейной интерполяции табличного УРС  на треугольной ячейке, лежащей «выше диагонали» или «ниже диагонали» прямоугольной ячейки. С констатацией факта: «Ясно, что указанный способ интерполяции с сохранением непрерывности не является единственным».

§ 6. Уточнение акустического приближения.

         6.1. Как уже отмечалось, задача о распаде разрыва сводится к решению уравнения для давления Р на контактном разрыве:

                                                      (6.1)

         В случае линейной аппроксимации (4.3) УРС в треугольных ячейках удалось получить явные формулы для функций и . Согласно (2.7)

                                ,                                (6.2)

где V(P) определяется формулой (4.5). Согласно (4.11)

                    ,            (6.3)

где                                    .                                              (6.4)

         Решение уравнения (6.1), аналогично описанному, например, в [2]-[3], можно осуществлять с привлечением метода касательных Ньютона. В [17] предложено применять для этой цели метод обратной параболической интерполяции, который (по утверждению ее авторов) имеет кубическую скорость сходимости, экономичен и надежен.

         6.2. С точки зрения практики расчетов большой интерес представляют также безитерационные методы, если они дают достаточно удовлетворительное приближенное решение. Один из таких методов изложен в [15] применительно к УРС двучленного вида (1.4). Есть возможность использовать этот метод и в рассматриваемой ситуации.

         Для его реализации необходимо вычислить три опорных значения функции F(P), определяющих конфигурацию распада разрыва. Для сокращения числа вариантов предполагаем, как это делалось уже в [2]-[3], что .

         Тогда возможны следующие ситуации.

1⁰. Если , то  и реализуются две ударные волны.

2⁰. Если , то  и реализуется конфигурация из волны разрежения и ударной волны.

3⁰. Если  , то реализуются две волны разрежения.

4⁰. Если , то возникает область вакуума, в которой следует полагать r=0_, с=0.

         6.3. С формул для вычисления  и начнем.

                          ,

где  вычисляется по формуле (4.5) при  с использованием коэффициентов (4.3) для левой среды (с индексом 1) в нужной («правильной») треугольной ячейке.

Вычисляется по формуле (6.3) при  с использованием коэффициентов (4.3) для правой среды (с индексом 2) в своей («правильной») треугольной ячейке.

         Далее, формулы (6.3) теряют работоспособность при  для левой среды и  для правой среды. Эти значения вычисляются по формулам (6.4) со своими коэффициентами для левой и правой сред, вычисленными в соответствующих треугольных ячейках. Чтобы сохранить работоспособность формул для обеих сред, полагаем

Далее, аналогично описанному в [15] на с.21, для двучленных УРС, полагаем:

         6.4. Будем предполагать, что .

Как описано в [15] на с.16, по полученным трем опорным точкам , ,  строится парабола в виде обратной зависимости Р(j). Она определяется формулой

      (6.5)

Искомое приближенное значение Р, отвечающее , определяется как .

         Обратим внимание на особые случаи, когда оказывается, что  или .

Можно практически реализовать их, например, используя заданный параметр , как  или .

         Тогда вместо трех опорных точек фактически остается только две. Такие ситуации в [15] рассматривались как вырожденные случаи и разрешались посредством замены параболы отрезком прямой. Мы не будем загромождать изложение ввиду полной аналогии.

         Особое внимание в [15] было уделено и случаю двух ударных волн, поскольку при достаточно больших значениях  и искомая величина корня  может быть сколь угодно большой. Этот случай рассматривался на основе замены (6.1) приближенным квадратным уравнением (см. [15], с.17-18). Ввиду громоздкости соответствующего описания мы его также опускаем. Вид формул (6.2)-(6.3) показывает, что такой подход вполне возможен.

         Заметим, что в случае небольших значений  можно этого и не делать, «рискнув» экстраполировать параболу (6.5), т.е. полагая и в этом случае .

         На этом описание уточнения «акустического приближения» закончим.

Заключение

         Табличную форму задания УРС можно считать самой универсальной. Однако для практического использования в газодинамических расчетах эти табличные данные необходимо восполнять. Простейшей формой такого восполнения является независимая билинейная аппроксимация затабулированных функций в табличных ячейках. Она обеспечивает непрерывность конструируемых УРС на общих ребрах ячеек.

         Такая аппроксимация использовалась в [11] для решения достаточно сложных задач. При этом оказалось, что при решении задач о распаде разрыва для получения некоторых функций приходится прибегать к численному интегрированию.

         Если вместо плотности r использовать в качестве независимой переменной удельный объем , то алгоритм упрощается для случая ударной волны.

         А если использовать вместо билинейной независимую линейную аппроксимацию табличных функций ,  в треугольных ячейках, то удается получить явные формулы для всех функций, используемых при решении задачи о распаде разрыва. Это несомненно представляет интерес для практики расчетов.

         Наиболее «хлопотная» часть алгоритмов связана с возможной необходимостью переходов из исходной в соседние табличные ячейки УРС. Здравый смысл подсказывает, что можно ослабить требование такого перехода, заменив «строгое» условие выхода за пределы табличной ячейки на выход за «расширенные» пределы (например, с использованием того же параметра ). Кроме того, естественно ожидать, что переходов можно  избежать (или, по крайней мере, существенно сократить) при достаточно подробной геометрической сетке, на которой считается задача. Резко возросшие возможности вычислительной техники позволяют всерьез относиться к такому фактору.

         С табличными УРС ситуация противоположная. Для практики желательно по возможности «укрупнять» табличные ячейки. К этому вынуждает и широкодиапазонный характер УРС, охватывающего несколько порядков по каждой из переменных. Однако, как уже отмечалось в §1, нужно заботиться и о передаче табличными данными условий  термодинамической согласованности и невырождения якобиана.

         В §3 было обращено внимание на возможность визуального и аналитического исследования дискретного отображения, порождаемого таблицами, задающими УРС. Присутствие самопересекающихся ячеек («бантиков») следует считать неприемлемым, поскольку нарушается однозначность отображения. На практике это может приводить к «отказу» рабочих формул алгоритма или получению неприемлемых, лишенных смысла результатов.

         В связи с этим возникает ряд вопросов о контроле и возможном исправлении таблиц, задающих УРС. Они должны стать предметом специальных исследований.

         Автор выражает признательность Е.А.Забродиной. При ее помощи была получена информация о конкретных таблицах УРС, используемых в расчетах, которая стимулировала выполнение представленных исследований.

         Автор благодарен М.С.Гавреевой и А.В.Северину за помощь в оформлении работы.

Литература

1.           Куропатенко В.Ф. Модели механики сплошных сред.// Челябинск, Челябинский гос.университет, 2007, 302 с.

2.           Численное решение многомерных задач газовой динамики. Под ред. С.К.Годунова.// М., «Наука», 1976, 400 с.

3.           Куликовский А.Г., Погорелов А.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений.// М., ФИЗМАТЛИТ, 2001, 608 с.

4.           Прокопов Г.П. Необходимость контроля энтропии в газодинамических расчетах.// Ж. вычисл. матем. и матем. физики, 2007, 47, №9, с.1591-1601.

5.           Алалыкин Г.Б., Годунов С.К., Киреева И.Л., Плинер Л.А. Решение одномерных задач газовой динамики в подвижных сетках.// Изд-во «Наука», Гл.ред.физ.-мат.лит., М., 1970, 112 с.

6.           Куропатенко В.Ф., Коваленко Г.В. и др. Комплекс программ «ВОЛНА» и неоднородный разностный метод для расчета неустановившихся движений сжимаемых сплошных сред. Часть 1. Неявный разностный метод. Часть 2. Архитектура. Общая организация.// Вопросы атомной науки и техники, сер. «Математическое моделирование физических процессов», М., 1989, вып.2, с.9-25.

7.           Бушман А.В., Канель Г.И., Ни А.Л., Фортов В.Е. Теплофизика и динамика интенсивных импульсных воздействий.// Черноголовка, 1988.

8.           Рождественский Б.Л., Яненко Н.Н. Системы квазилинейных уравнений и их приложения к газовой динамике.// М., «Наука», Главная ред. физ.-мат. лит., 1968, 592 с.

9.           Зоммерфельд А. Термодинамика и статистическая физика.// М., Изд-во Иностр. литер., 1955, 480 с.

10.      Елисеев Г.М. Термодинамическая согласованность и сплайн-уравнение состояния.// Вопросы атомной науки и техники, сер. «Теоретическая и прикладная физика», 1993, вып.1, с.35-45.

11.      Чарахчьян А.А. Об алгоритмах расчета распада разрыва для схемы С.К.Годунова// Ж. вычисл. матем. и матем. физики, 2000, 40, №5, с.782-796.

12.      Прокопов Г.П. Аппроксимация табличных уравнений состояния для расчета газодинамических задач.// Препринт ИПМ им.М.В.Келдыша РАН, 2004, №80, 28 с.

13.      Прокопов Г.П. Исследование формул аппроксимации уравнений состояния в переменных температура-плотность.// Препринт ИПМ им.М.В.Келдыша РАН, 2005, №26, 28 с.

14.      Бронштейн И.Н., Семендяев К.А. Справочник по математике.// Лейпциг, Тойбнер, 1979; М., Наука, 1980.

15.      Прокопов Г.П. О приближенных реализациях метода Годунова.//  Препринт ИПМ им.М.В.Келдыша РАН, 2007, №15, 28 с.

16.      Иваненко С.А., Чарахчьян А.А. Криволинейные сетки из выпуклых четырехугольников.// Ж. вычисл. матем. и матем. физики, 1988, т.28, №4, с.503-514.

17.      Кобзева Т.А., Моисеев Н.Я. Метод неопределенных коэффициентов для решения задачи о распаде произвольного разрыва//Вопросы атомной науки и техники. Сер.: Математ. моделир. физ. процессов, 2003, вып.1, с.3-9.