Содержание

1. Введение........................................................................................................................ 3

2. Разностные схемы для расчетов задач газовой динамики........................................ 4

2.1. Неконсервативные разностные схемы................................................................ 5

2.2. Консервативные и полностью консервативные разностные схемы................ 8

3. Ударные волны разрежения в численных решениях.............................................. 12

3.1. Адиабатический случай...................................................................................... 12

3.2. Изотермический случай...................................................................................... 15

4. Дифференциальное приближение............................................................................ 18

4.1. Модельный пример............................................................................................. 18

4.2. Дифференциальное приближение для изотермического случая..................... 19

5. Анализ дифференциального приближения.............................................................. 21

6. Заключение................................................................................................................. 29

Литература...................................................................................................................... 30

1. Введение

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

Дискретная и дифференциальная модели не тождественны друг другу. Например, в дискретной модели могут не выполнятся аналоги законов, характеризующих физический объект, например, законов сохранения [1,2,3]. Поэтому возможны заметные различия между решениями, полученными численно на основе дискретных моделей, и решениями дифференциальных задач. При этом поведение численного решения зависит от шагов разностной сетки, способа аппроксимации отдельных членов исходных дифференциальных уравнений, аппроксимации граничных условий, способа представления области, в которой решается задача и т. д.

Примером раздела науки, где активно используются методы математического моделирования, является газовая динамика. Разнообразие задач газовой динамики, сложность системы уравнений, ее нелинейность, - все это диктует высокие требования к методам численного решения. Примером непростых для моделирования газодинамических объектов могут служить ударные волны и волны разрежения, изучение которых представляет большой теоретический и практический интерес.

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

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

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

2. Разностные схемы для расчетов задач газовой динамики

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

Выберем, например, следующий вид системы уравнений одномерной газовой динамики в лагранжевых массовых переменных для случая плоской симметрии [2,6]:

, , .                                                              (2.1)

Систему дополняют уравнение состояния идеального газа:

,                                                                                              (2.2)

и уравнение для определения траектории частицы газа:

.                                                                                                       (2.3)

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

, где  - газовая постоянная (в описанных ниже расчетах ).

Система уравнений (2.1)-(2.3) решается в области ,  - масса газа в единичном параллелепипеде (см. [3]). В качестве модельной выберем задачу о поршне.

Рассматривается полупространство (в массовых координатах ), заполненное в начальный момент  однородным, покоящимся газом:

, , .

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

2.1. Неконсервативные разностные схемы

Напомним некоторые известные результаты исследований разностных схем газовой динамики.

В пространственной области  введем равномерную сетку:

Ниже будут использованы следующие обозначения [1]:

, , , , , ,

, , ,

где  сеточная функция.

К узлам сетки  будем относить сеточные функции скорости и эйлеровой координаты, к точкам  – давления, плотности, внутренней энергии. Аппроксимируя дифференциальные уравнения газовой динамики (2.1) – (2.3), можно получить следующую (одну из возможных) разностную схему “крест” [13, 14]:

, , , , .                                   (2.4)

Локальная погрешность аппроксимации отдельных дифференциальных уравнений системы (2.1)-(2.3) разностными соотношениями (2.4) - . Привлекательная черта схемы состоит в том, что, формально являясь неявной, она, тем не менее, разрешается явным образом. С другой стороны схема не лишена недостатков. Они проявляются, например, при расчете задачи о вдвигаемом в газ поршне.

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

.

Параметры газа за фронтом волны вычисляются по формулам:

, , , .                           (2.5)

При значениях , , ,  они составляют , , , , .

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

, ,

в качестве начальных данных возьмем исходное состояние газа

, , .

Для обеспечения возможности сквозного расчета ударных волн введем в схему линейную вязкость

с коэффициентом . Выберем равномерную разностную сетку с шагом по массе , шаг по времени будем варьировать и сравнивать с шагом Куранта:

, где  - скорость звука,

который для параметров газа за фронтом волны составляет .

На рис. 1 (на этом и всех последующих рисунках пунктиром, если не оговорено обратное, показано точное решение) приведены результаты расчетов по схеме (2.4) на момент времени  с шагами по времени  () и  (). При небольших временных шагах () решение передается хорошо, фронт ударной волны “размазан” на 3 интервала сетки. Наблюдаемый пик температуры в окрестности поршня – известный эффект, энтропийный след. Однако при увеличении  () в численном решении наблюдаются колебания за фронтом волны. Дальнейшее увеличение шага по времени приводит к сильным искажениям решения.

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

Модифицируя схему (2.4), увеличивая степень ее неявности и повышая запас устойчивости, можно добиться возможности расчета задачи с крупными шагами по времени ():

, , , , .                                   (2.6)

На рис. 2 представлены результаты соответствующих расчетов (в этой работе для решения неявных схем применялся метод Ньютона) на момент времени  с шагами по времени  (),  () и  (). Как видно, при малых шагах по времени (, ) численное решение достаточно хорошо соответствует точному. В случае схемы (2.6) возможно использование и больших  (); параметры течения газа за фронтом волны, полученные в результате расчетов, заметно отличаются от точного решения (для температуры при  примерно на 25%), ширина фронта волны возрастает и скорость его движения меньше единицы. Отметим, что отклонения температуры за фронтом волны от точного решения в схемах (2.4) и (2.6) – противоположны по знаку. Аналогичны отличия в скорости фронта ударной волны.

Иследуем, с чем связаны недостатки схем (2.4) и (2.6).

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

,                                                                                  (2.7)

В разностном случае для схемы “крест” (2.4) уравнение энергии приводится к виду (в этой и всех последующих формулах черта над сеточными функциями плотности, давления и внутренней энергии, означающая, что значение функции относится к полуцелой точке, опущена):

, где ,

для схемы (2.6) к виду:

, где .

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

Естественно требовать, чтобы в схеме выполнялись разностные аналоги законов сохранения, то есть разностная схема была консервативной [2,3].

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

Рассмотрим семейство консервативных схем газовой динамики [3]:

, , ,

,                                                            (2.8)

.

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

На рис. 3 представлены результаты расчетов задачи о вдвигаемом поршне по одной из схем семейства (2.8): , . Начальные данные, граничные условия и шаг сетки по пространству были взяты теми же, что и в предыдущих расчетах. Расчеты проводились с шагами по времени  () и  (). Легко видеть, что решение в обоих расчетах передается хорошо. Фронт волны размазан на 3-4 интервала сетки; с увеличением  ширина фронта возрастает.

Однако и у консервативных схем есть недостатки. Связаны они с тем, что в реальной газодинамической среде выполняются не только законы сохранения, но и другие важные соотношения. Так, уравнение неразрывности (второе уравнение в (2.1)), используя уравнение для траектории частицы газа (2.3), можно записать в виде:

.                                                                                                     (2.9)

Это соотношение называют “законом сохранения объема” [3]. Уравнение энергии помимо дивергентной может быть записано в энтропийной и в недивергентной формах:

,                                                                                      (2.10)

.                                                                                              (2.11)

Рассматривая семейство консервативных схем, легко придти к следующим результатам:

, где ,                                              (2.12)

, где ,        (2.13)

, где .                                    (2.14)

Требуя, чтобы дисбалансы , ,  тождественно обращались в нуль, получим условия [3]

, ,                                                                  (2.15)

которые выделяют из семейства консервативных схем (2.8) однопараметрическое семейство. Построенные схемы обычно называют полностью консервативными схемами газовой динамики [2].

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

Пусть в холодный газ () вдвигается поршень с нарастающей скоростью . Начальные данные и граничные условия выберем следующим образом:

, , , .

Эта задача допускает автомодельное решение. В качестве автомодельной переменной удобно взять

.

Автомодельность сохранится, если ввести вязкость следующего вида:

,

для расчетов положим . Решение будет иметь вид:

, , , , ,                          (2.16)

где автомодельные функции , , , ,  легко найти из системы обыкновенных дифференциальных уравнений [8].

Выберем равномерную разностную сетку с шагом  по массе, и  по времени. На рис. 4 в автомодельной переменной даны результаты расчетов поставленной задачи по консервативной схеме с весами ,  на три момента времени: , , . Значения  накладываются примерно на одну линию, что свидетельствует о автомодельности численного решения; однако эта линия не совпадает с точным решением. Причем, отличия от точного решения в некоторых точек составляют более 100% (см. ).

Таким образом, рассмотренная консервативная схема дает существенно искаженное решение задачи.

На рис. 5 показаны результаты расчетов этой же задачи по полностью консервативной схеме (). Легко видеть, что использование полностью консервативной схемы приводит к хорошему результату; численное решение автомодельно, значения  накладываются на точное решение.

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

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

Для завершения настоящего обзора результатов, приведем необходимые и достаточные условия устойчивости семейства схем газовой динамики (2.8) в акустическом приближении [7]:

,                                                                                               (2.17)

.                                                             (2.18)

3. Ударные волны разрежения в численных решениях

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

3.1. Адиабатический случай

Рассмотрим задачу о поршне, выдвигаемом из газа с постоянной скоростью . В этом случае по газу от поршня будет распространяться волна разрежения; решение может быть получено в автомодельном виде [8]:

,  и т.д.,

где ,  - скорость звука,

, ,

, ,

, , .                                     (3.1)

В лагранжевых массовых переменных это решение для любого  “подобно” автомодельному, только ширина волны разрежения  линейно растет со временем. При этом передний фронт волны движется со скоростью звука на “фоне” , а задний со скоростью, вычисляемой по параметрам в зоне, примыкающей к поршню:  (см. рис. 6).

Выберем для расчета следующие начальные данные:

, ,  (),

и граничные условия:

 - левое граничное условие,

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

Для численного решения воспользуемся схемой “крест” (2.4) без искусственной вязкости.

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

, где  - скорость звука.

Расчеты, проведенные с шагом по времени равным , , показывают, что схема хорошо передает решение задачи (см. рис. 7).

Однако при расчете по той же схеме, например, с шагом ,  в численном решении возникает разрыв - ударная волна разрежения (рис. 8 и рис. 9).

В дифференциальном случае возникновение ударных волн разрежения запрещает теорема Цемплена [6], однако дискретные среды обладают иными свойствами, в частности, в них могут возникать такие нефизичные эффекты как ударные волны разрежения. Пример такой среды дает рассмотренная схема “крест”. Следует отметить, что ударные волны разрежения наблюдаются при расчетах не только по схемам рассматриваемого семейства, но и по другим используемым в настоящее время схемам, например по схеме Роу [9].

Эффект возникновения ударных волн разрежения присущ не только схеме “крест”, но и другим схемам семейства (2.8). Например, этим же свойством обладает явная полностью консервативная схема (). На рис. 10 и рис. 11 приведены результаты соответствующих расчетов с шагами по времени ,  (рис. 10) и ,  (рис. 11). Обратим внимание, что в явной полностью консервативной схеме ударная волна разрежения возникает при других величинах шагов сетки по времени, чем в схеме “крест”. Так, в “кресте” она наблюдается лишь при , а в явной полностью консервативной схеме уже при  (см. рис. 11).

Однако не для всех схем семейства (2.8) в численных решениях задачи о выдвигаемом поршне возникают ударные волны разрежения. Например, для схемы (2.6) – неявного “креста” и для неявной полностью консервативной схемы () этот эффект отсутствует при любых . Соответствующие результаты расчетов приведены на рис. 12. В обоих случаях расчеты проводились при , . Легко видеть, что решение задачи заметно размазано. Это объясняется наличием в схемах сильных внутренних диссипативных и дисперсионных свойств. При увеличение  схемная диссипация усиливается, решение размазывается на большее число интервалов сетки.

3.2. Изотермический случай

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

Схема (2.8) в этом случае принимает вид:

, , , ,                                               (3.2)

 - скорость звука. Начальные данные и граничные условия выберем следующим образом:

, ,  () - начальные данные,

 - левое граничное условие,  - правое граничное условие,  - изотермическая скорость звука.

Рассмотрим аналог схемы “крест” для изотермического случая: , и проведем для нее подробный анализ численного решения задачи о выдвигаемом поршне при различных шагах по времени. Шаги по времени будем сравнивать с шагом Куранта.

При  получаем решение, хорошо согласующееся с точным. Результаты соответствующего расчета для  () представлены на рис. 13.

При , где  - экспериментально найденное значение для данного расчета, получаются решения, показанные на рис. 14 для  (). Возникает ударная волна разрежения, распространяющаяся по газу. Особенностью ударной волны разрежения в решении по схеме “крест”, как и в адиабатическом случае, является практически неразмазанный скачок параметров газа на фронте волны (см. рис. 15, где окрестность ударной волны разрежения дана в укрупненном виде).

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

При  () решение помимо ударной волны разрежения содержит ударную волну сжатия (см. рис. 17). Осцилляции за ее фронтом происходят в силу отсутствия в схеме искусственной вязкости.

Проверим, выполняются ли условия Гюгонио - соотношения, выражающие законы сохранения массы и импульса - на ударной волне разрежения. В случае лагранжевых массовых переменных они могут быть записаны в виде:

 - следствие закона сохранения массы,                           (3.3)

 - закон сохранения импульса.                                     (3.4)

где индексом 0 обозначены значения параметров газа слева от разрыва, а индексом 1 – справа. В таблице приведены полученные в расчете для  (см. рис. 14) значения параметров среды.

Для определения скорости фронта ударной волны воспользуемся соотношением: , где  - массовый путь, пройденный фронтом волны за время . Для данного расчета эта формула дает .

Подставляя далее найденное значение  в (3.3) и (3.4), убеждаемся в выполнении соотношений Гюгонио для ударной волны разрежения.

Легко также проверить, что для ударной волны разрежения выполняются неравенства для скорости фронта ударной волны: . В нашем случае: , , .

4. Дифференциальное приближение

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

4.1. Модельный пример

Рассмотрим линейное уравнение переноса:

, , , .                                                        (4.1)

Аппроксимируем это уравнение на равномерной сетке с помощью следующей разностной схемы:

, , ,                                                         (4.2)

или в индексной форме:

.

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

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

.                                                                         (4.3)

Будем предполагать, что это соотношение выполнено всюду в области . В этом случае (4.3) представляет собой так называемую Г-форму первого дифференциального приближения схемы (4.2).

Используя (4.3) легко перейти к П-форме дифференциального приближения схемы:

, .                                                                  (4.4)

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

Заметим, что при нарушении условия Куранта, то есть при

,                                                                                                         (4.5)

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

4.2. Дифференциальное приближение для изотермического случая

Рассмотрим теперь общее семейство изотермических схем газовой динамики, с учетом искусственной линейной вязкости :

,                                                                                              (4.6.a)

,                                                                                            (4.6.b)

, , .                                                                (4.6.c)

Преобразуем последние три уравнения в одно, тогда схема примет вид:

,                                                                                              (4.7.a)

,                                                                                               (4.7.b)

,                                                                                       (4.7.c)

где обозначено .

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

В итоге получим Г-форму дифференциального приближения [11] разностной схемы (4.7):

,         (4.8.a)

,                       (4.8.b)

.                                                                                       (4.8.c)

Перейдем к П-форме дифференциального приближения. Пользуясь уравнениями (4.8.a) и (4.8.b), подставим в их правые части производные по времени  и . Приводя подобные члены, получим Г-форму в ином виде:

,                                               (4.9.a)

,                                                (4.9.b)

.                                                                                       (4.9.c)

Необходимо отметить, что полученные уравнения имеют одинаковую структуру. При этом члены с , отличающие уравнение (4.8.a) от (4.8.b), возникли за счет того, что разностное уравнение (4.7.a) записано симметрично для узла , а (4.7.b) для . Поскольку разложение в ряд Тейлора было произведено в узле , в разложение функций из (4.7.a) вошли дополнительные члены с . Однако дальнейший переход к (4.9) привел к взаимному уничтожению этих членов.

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

,                                         (4.10.a)

              (4.10.b)

, где , ,                                            (4.11)

, где ,                                                (4.12)

, ,                (4.13)

.                                                                (4.14)

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

В дифференциальном приближении присутствуют также члены, содержащие третьи производные , , которые, как показано в [12], характеризуют дисперсионные свойства схемы.

5. Анализ дифференциального приближения

1. Уравнение структуры фронта ударной волны в вязкой среде

Рассмотрим первое дифференциальное приближение семейства схем (4.6):

, , ,                                                                  (5.1.a)

,                                                                                    (5.1.b)

, где , ,                                             (5.2)

.                                                                                 (5.3)

Будем трактовать уравнения (5.1) – (5.3) как математическую модель некоторой гипотетической среды. Очевидно, при  и  она близка к обычной изотермической газодинамике. При конечных шагах сетки эта среда обладает некоторыми своеобразными диссипативными свойствами, роль которых необходимо оценить. Построим для системы уравнений (5.1) – (5.3) автомодельное решение типа “бегущей волны” [8], которое используется для получения структуры фронта ударной волны в обычных газодинамических средах.

Будем искать решение в виде

, .                                                             (5.4)

Решение рассматривается в неограниченной области , при , при этом . Роль правого граничного условия при  играет “фон” ; будем считать, что при  решение “выходит” на постоянные значения , которые неизвестны. Также считаем, что

.                                                                                              (5.5)

Таким образом, мы ищем установившееся решение, перемещающееся по массовой координате с постоянной скоростью .

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

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

,

и тем самым свести систему уравнений в частных производных (5.1) к системе обыкновенных дифференциальных уравнений:

,                                                                                              (5.6.a)

.                                                                                   (5.6.b)

Равенство (5.5), играющее роль граничных условий, преобразуется к виду

.                                                                                             (5.7)

Проинтегрируем уравнения (5.6) по  от фона  до некоторого текущего значения:

,                                                                                             (5.8.a)

.                                                                                   (5.8.b)

Отсюда легко получить

,

или для случая идеального газа

,                                                                       (5.9)

где  - значение полученное из классических соотношений Гюгонио.

Заметим, что обе части этого равенства обращаются в нуль на фоне при  () в полном соответствии с (5.7). Аналогично левая, а, следовательно, и правая части должны обратиться в нуль и при . Поэтому значение плотности, которое устанавливается при , должно быть равно .

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

Вернемся к уравнению (5.9) и подставим в него формулы (5.2) и (5.3). Имеем:

,

.

Поскольку мы ограничились рассмотрением первого дифференциального приближения и считаем, что , то можно пренебречь членами второго порядка малости по  и . Тогда

и уравнение (5.10) может быть преобразовано к виду:

.                                             (5.10)

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

2. Случай

Рассмотрим решения, которые дает уравнение (5.10) при отсутствии искусственной вязкости (). В этом случае уравнение (5.10) имеет вид

                                                            (5.11)

Если  (это условие выполнено, например, для схемы “крест”), то уравнение (5.11) не является дифференциальным и с точки зрения структуры фронта необходимо учитывать члены второго порядка малости. Из (5.11) получим:

.                                                                                   (5.12)

Возможные решения данного уравнения приведены на рис. 18 (соответствует ударной волне сжатия) и рис. 19 (соответствует ударной волне разрежения).

Отметим, что решение уравнения (5.12) имеет разрыв. Это соответствует расчетам по схеме “крест” задачи о выдвигаемом поршне, в которой возникало практически разрывное решение - ударная волна разрежения (см. пункт 3).

Если , то уравнение (5.11) – дифференциальное. Оно легко интегрируется:

, где , .

Однако решение удобнее анализировать, рассматривая само уравнение (5.11), дающее знак производной плотности:

.                                                          (5.13)

Рассмотрим случай, когда . Если  (при ), то уравнение (5.13) имеет решение, представляющее собой ударную волну сжатия; оно показано на рис. 20. Для  решение есть ударная волна разрежения, как показано на рис. 21. Таким образом, мы построили два формально возможных решения задачи.

Учитывая, что, при  плотность должна “выходить” на “фоновое” значение , можно выбрать единственно допустимое решение задачи (см. [15]). В данном случае в среде существует ударная волна сжатия (рис. 20), а ударная волна разрежения существовать не может. Математически последнее означает, что при  () решения задачи типа “бегущей волны” нет, и его следует искать в ином виде. Для обычной газовой динамики такое решение есть римановская волна разрежения, являющееся автомодельным ().

Аналогично строится решение при  - см. рис. 22. В среде возникает ударная волна разрежения с , при этом , .

Таким образом, при  в среде будет существовать ударная волна сжатия, а при  - ударная волна разрежения. Отметим, что условие возникновения ударной волны сжатия совпадает с условием (2.17) устойчивости разностной схемы в акустическом приближении.

Проведем расчеты задачи о выдвигаемом поршне по схемам (4.6) со следующими весами: , . Легко видеть, что для этих схем выполнено условие , а значит, в решении может содержаться ударная волна разрежения.

Выберем следующие начальные данные для расчета:

, ,  (),

пусть скорость поршня равна ; значение изотермической скорости звука . Выберем значения шагов по пространству и времени: , .

На рис. 23 приведены результаты расчетов для трех схем , , . Действительно, в решениях присутствует ударная волна разрежения, размазанная схемной вязкостью, действие которой увеличивается с ростом . Отметим, что с увеличением  ударная волна разрежения становится все более “размазанной” и при  практически пропадает. На рис. 24 приведены результаты расчетов по схеме  с двумя шагами по времени:  и . Легко видеть, что с увеличением шага  ширина размазывания фронта ударной волны возрастает. Это связано с увеличением коэффициента схемной вязкости . Как показывают результаты расчетов (см. пункт 3), для схемы “крест” она появляется при ; для схем ,  при , причем .

3. Случай

Оценим теперь влияние искусственной вязкости на решение задачи о выдвигаемом поршне. Проведем расчеты с различными коэффициентами : без вязкости (, ) и с вязкостью (, ), при этом потребуем выполнение следующего соотношения (см. (5.10)):

.

И, следовательно,

.                                                                                 (5.14)

Выберем , ,  (при этом значении шага  в численных решениях присутствуют ударные волны разрежения). На рис. 25 представлены результаты соответствующих расчетов. Для сравнения на том же рисунке нанесены точное решение (пунктир) и численное решение с  без вязкости. Легко видеть, что введение в схему с весом  искусственной вязкости (5.14) приводит к изменению амплитуды ударной волны и ширины ее фронта (штрихпунктирная линия), в результате, численное решение задачи с  становится близко к решению задачи по схеме с весом . Отметим, что коэффициент вязкости (5.14) нелинеен. Если его задать постоянным: , то различие между численными решениями немного увеличивается (см. рис. 25).

3. Случай нелинейной искусственной вязкости

Пусть в схему введена искусственная вязкость () следующего специального вида:

.                                                                               (5.15)

Ограничимся рассмотрением случая . Преобразуем уравнение (5.10)

.                                                (5.16)

Анализ возможных интегральных кривых приводит к следующим результатам.

Если , то при  в среде возникает ударная волна сжатия, а при  ударная волна разрежения.

Если , то при  в среде возникает ударная волна разрежения, а при  ударная волна сжатия.

Как и в предыдущем пункте рассмотрим схемы семейства (4.6) с весами , , для которых выполнено условие .

Выберем следующие начальные данные для расчета:

, ,  (),

пусть скорость поршня равна ; значение изотермической скорости звука .

В схемы введена линейная вязкость

,

коэффициент вязкости положим равным

.                                                     (5.17)

Отметим, что переход от нелинейного коэффициента вязкости (5.15) к линейному коэффициенту вида (5.17) оправдан в силу того, что рассматриваемая ниже ударная волна разрежения имеет малую амплитуду; на волне будет выполнено условие

На рис. 26 показано влияние данной вязкости на решение, полученное по схеме с  и шагами ,  (для выбранных параметров ).

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

Приведенные расчеты показывают влияние вязкости с коэффициентом (5.17) на решение задачи о выдвигаемом поршне. Они подтверждают сделанные аналитически выводы о том, что при  и  возникает ударная волна разрежения, а при  волна сжатия.

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

6. Заключение

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

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

Литература

1. Самарский А.А. Теория разностных схем. – М.: Наука, 1989.

2. Попов Ю.П., Самарский А.А. Полностью консервативные разностные схемы - ЖВМ и МФ, 1969, Т. 9, № 4, С. 953-958.

3. Самарский А.А., Попов Ю.П. Разностные методы решения задач газовой динамики. – М.: Наука, 1992.

4. Тихонов А.Н., Самарский А.А. Об однородных разностных схемах - ЖВМ и МФ, 1961, Т. 1, № 1, С. 5-63.

5. Давыдов Ю.М. Применение метода дифференциального приближения для исследования и построения нелинейных разностных схем. – ЧМ МСС, Новосибирск: 1980, Т. 11, № 4, С. 41-59.

6. Ландау Л.Д., Лифшиц Е.М. Гидродинамика. – М.: Наука, 1988.

7. Арделян Н.В., Гулин А.В. К обоснованию устойчивости разностных схем для уравнений акустики. – Препринт/ИПМ АН СССР. – М., 1978. - №96.

8. Седов Л.И. Методы подобия и размерности в механике. – М.: Физматгиз, 1962.

9. Кузнецов О.А. Численное исследование схемы Роу с модификацией Эйнфельдта для уравнений газовой динамики. – Препринт/ИПМ АН СССР. – М., 1998.

10. Трощиев В.Е. О дивергентности схемы “крест” численного решения уравнений газовой динамики // Численные методы механики сплошной среды. – Новосибирск: Наука. – 1970.

11. Шокин Ю.И. Метод дифференциального приближения, Новосибирск: Наука, 1979.

12. Мухин С.И., Попов С.Б., Попов Ю.П. О разностных схемах с искусственной дисперсией // ЖВМ и МФ. - Т. 23, №6.

13. Рождественский Б. Л., Яненко Н. Н. Системы квазилинейных уравнений. – М.: Наука, 1978.

14. Von Neumann J., Richtmyer R. D. A method for numerical calculation of hydrodynamics shocks // J. Appl. Phys. – 1949. V. 21.

15. Гельфанд И.М. Некоторые задачи теории квазилинейных уравнений. – УМН, 1959, Т. XIV, вып.2(86), С. 87-158.