ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2007, том 47, < 9, с. 1591-1601
УДК 519.634
НЕОБХОДИМОСТЬ КОНТРОЛЯ ЭНТРОПИИ В ГАЗОДИНАМИЧЕСКИХ РАСЧЕТАХ1)
© 2007 г. Г. П. Прокопов
(125047 Москва, Миусская пл., 4, ИПМатем. РАН) Поступила в редакцию 27.12.2006 г.
На примере одной недавно опубликованной схемы проведено аналитическое исследование поведения энтропии при численном интегрировании уравнений газовой динамики. Оно подтвердило опасение о возможности реализации в ходе расчета разрывов, на которых может происходить нарушение постулата о неубывании энтропии. Предлагается простой алгоритм контроля энтропии при расчетах газодинамических течений и по другим известным численным методам. Библ. 12. Фиг. 2.
Ключевые слова: газодинамика, расчет, энтропия, контроль.
ВВЕДЕНИЕ
Как было показано в [1], для обеспечения единственности обобщенного решения системы уравнений одномерной газовой динамики в дивергентной форме, позволяющей выполнять законы сохранения массы, импульса и энергии, следует накладывать на него некоторые дополнительные требования. Они состоят в запрете ударных волн разрежения, что эквивалентно выполнению для любого замкнутого контура неравенства, из которого следует неубывание энтропии.
Опубликована [1] была еще до того, как в 1959 г. появилось описание метода (см. [2]), который теперь принято называть методом С.К. Годунова. Это было сделано с большой задержкой (см. [3, с. 10]), когда он успешно использовался уже несколько лет. Популярности метода способствовала монография [4], которая в основном посвящена именно этому методу и его применениям.
Выполнение закона сохранения энтропии на гладких решениях и постулат о ее неубывании на разрывах были заложены при конструировании разностной схемы, лежащей в основе метода Годунова. В связи с тем, что имеется много публикаций различных авторов, заявляющих о методах "типа Годунова", интересен поставленный в [5] вопрос, с должным ли вниманием они относятся к упомянутому факту. Естественно, это касается не только таких методов, но и вообще алгоритмов расчета газодинамических течений.
Попытки исследовать с такой точки зрения какой-нибудь из разнообразных численных методов, представленных, например, в монографии [6], встретили серьезные затруднения из-за сложности алгоритмов. Исследовать поведение энтропии удалось для простой и компактной разностной схемы С (название взято из [7], где она недавно опубликована). Результаты изложены в [8] и будут представлены в настоящей работе.
1. СХЕМА С ДЛЯ РАСЧЕТА ГАЗОДИНАМИЧЕСКИХ ТЕЧЕНИЙ
Опишем разностную схему из [7], поскольку все формулы потребуются для дальнейшего изложения. Кроме того, изменены некоторые обозначения, с тем чтобы обеспечить единообразие с [4]. Описание дается на примере нестационарных уравнений одномерной газовой динамики в дивергентной форме:
эиШ + д Р/дх = 0, (1.1)
^Работа выполнена при финансовой поддержке РФФИ (код проекта 05-01-00097).
1591
1592 ПРОКОПОВ
где и, F - векторы для законов сохранения массы, импульса и энергии:
и =
р р u
р u , F = 2 pu + p
J}e (pe + p)u
(1.2)
Здесь р - плотность газа, и - скорость,р - давление, £ - внутренняя энергия, е = £ + и2/2 - полная внутренняя энергия единицы массы газа.
Для простоты ограничимся совершенным газом с показателем адиабаты у:
р = (У -1 )р£. (1.3)
Разностная схема для (1.1) записывается в виде
ДГ
U1-1/2 = U
j-1/2-ЛХ( F i- F i 1 ) ■
(1.4)
Здесь, как было принято в [4], у - 1/2 - номер ячейки сетки по оси х, ограниченной узлами с номерами у - 1 и у. Формула (1.4) описывает пересчет: переход от величин с нижнего слоя по времени t к верхнему.
В методе Годунова потоки на границах ячеек F7■ вычисляются посредством решения задач Ри-мана о распаде разрыва с параметрами газа в соседних ячейках. Расчет делается либо приближенно (например, в виде "звукового" приближения), либо (если необходимо) итерационным методом до достижения результата с предписываемой точностью.
Основная идея схемы С из [7] также состоит в вычислении потоков на границах ячейки. Задача Римана решается приближенно на основе соотношений на разрывах, аналогичных описанным в [4, с. 103]:
(1.5)
D = u - a/р, a[u] + [p] = 0,
[u] - a[ 1/р] = 0, a[e] + [pu] = 0.
Здесь D - скорость распространения разрыва, a - массовая скорость (вместо обозначений w, m соответственно в [7]). Квадратными скобками [ ] обозначается разность значений величин, заключенных в скобки, по обе стороны от разрыва.
В отличие от метода Годунова, при приближенном решении задачи Римана рассматривается упрощенная схема течения с распадом на левую волну, контактный разрыв и правую волну. Она изображена на фиг. 1.
Обозначим индексами 1 - параметры в левой ячейке сетки, 2 - в правой, 3 - между левой волной и контактным разрывом, 4 - между контактным разрывом и правой волной, как показано на схеме течения. Расчетные формулы имеют вид
U = ( u2 a2 + u1 a1 - p2 + p1 ) / ( a1 + a2 ), u3 = u4 = U, P = (p2a1 + p1 a2- a1 a2(u ))/(a1 + a2), p3 = p4 = P,
1/р3 = 1/р1 + ( U - u1 )/a1, e3 = e1 - (PU - p1 u1 )/a1, (1.6)
1/р4 = 1/р2 - ( U - u2 ) / a2, e4 = e2 + ( PU - p2 u ) /a2,
D1 = u1-a/pb D2 = u2 + a 2/P2 ■
Назначение параметров газа (R, U, P, E)j на границах ячейки для вычисления векторов потока Fj (в методе Годунова они назывались "большими" величинами) определяется значениями U, Dj, D2. В случае U > 0 при D1 > 0 они равны параметрам в зоне 1, а при D1 < 0 - в зоне 3. В случае U < 0 при D2 < 0 - их значениям в зоне 2, а при D2 > 0 - в зоне 4.
Определяющую роль играет назначение массовых скоростей ab a2. Для этой цели предлагается использовать формулы
а1 = р1(и1-ш1и(и1 - с1; и2 - с2)), а2 = р2(тах(и2 + с2, и1 + с1) - и2). Здесь с1, с2 - скорости звука. Для (1.3) они вычисляются по формулам
с1 = 7У Р1Р1, с2 = 7У Р2> Р ;.
(1.7)
(1.8)
Фиг. 1.
По утверждению в [7], такое назначение массовых скоростей а1, а2 обеспечивает отсутствие осцилляций на разрывах при произвольных физических параметрах в соседних ячейках разностной сетки.
2. ЭНТРОПИЙНОЕ ИССЛЕДОВАНИЕ СХЕМЫ С
Итак, одним из принципиальных элементов конструирования схемы С является замена волны разрежения на фронт разрыва.
Напомним, что в [4, с. 115-116] была рассмотрена задача о распаде разрыва с симметричными данными относительно границы раздела. Если в ней заменить возникающую волну разрежения фронтом разрыва, на котором выполняются соотношения (1.5), то реализуется ударная волна разрежения. А ее следует запрещать, как нарушающую постулат о неубывании энтропии.
В схеме С сначала "волевым" образом по формулам (1.7) назначаются массовые скорости аъ а2, а затем, исходя из тех же соотношений на разрыве, производится вычисление параметров, характеризующих состояние газа на разрывах. Поскольку на этом процесс их назначения заканчивается, говорить о реализации ударного фронта разрежения нет оснований. Но в каком случае можно полученный результат считать разумным и допустимым?
Предлагается постулировать следующее положение: газ с вычисленными параметрами, определяющими величины потоков через разрыв, должен иметь энтропию не меньшую, чем второй (исходный, "напарник") газ, участвующий в вычислении потоков.
Представляется, что только такое требование позволит контролировать возможное скрытое занижение энтропии в газодинамических расчетах. Назовем его жестким контролем энтропии.
Проще было бы ограничиться вычислением и контролем энтропии в ячейках сетки на верхнем слое по времени. Однако тогда, в соответствии с (1.4), появляется возможность компенсировать потери энтропии в одном узле сетки за счет другого. Поэтому назовем его мягким контролем энтропии.
Если согласиться с провозглашенным постулатом, то следующим шагом исследования должен быть вопрос: обеспечивает ли назначение массовых скоростей (1.7) его выполнение для рассматриваемой схемы С?
Это исследование будем проводить на уже упомянутом упрощенном примере распада разрыва с симметричными данными. Итак, пусть
Р1 = Р2 = Р *, Р1 = Р2 = Р *, и1 = и*, и2 = -и *• (2.1)
В качестве аналитического решения этой газодинамической задачи при и* > 0 (встречные потоки) реализуются ударные волны, а при и* < 0 (разбегающиеся потоки) - волны разрежения. Они движутся симметрично от границы раздела.
Благодаря симметрии (2.1) формулы (1.6)-(1.8) принимают вид
С1 = с 2 = с* = Р*/Р*, Р * = Р * с*/у, (2.2)
а1 = а2 = а * = (и * + |и*| + с* )р *, (2.3)
1594 ПРОКОПОВ
U = u3 = U4 = 0, P = рз = p4 = p * + a* u*, (2.4)
1/R = 1/рз = 1/P4 = 1/p*- u*/a*, £ = 63 = 64 = e* + p*u*/a*. (2.5)
Приступим к исследованию энтропии. Обычно для этого привлекается энтропийная функция a(s) = p/pY. Практически удобнее работать с величиной lna(s) = lnр - уlnр.
Заметим, что именно эта величина, с точностью до множителя cv (удельная теплоемкость при постоянном удельном объеме), принималась за настоящую энтропию (см. [9, с. 33]). Она остается постоянной вдоль адиабаты Пуассона. Поэтому полагаем, что S = lnP - ylnR, s* = lnp* - ylnp*. Условие
ее неубывания принимает вид
АS = S - s * = ln (P/p*) - у ln (R/р * )> 0. (2.6)
Используя формулы (2.2)-(2.5), получаем
А S = ln (1+ a *u */ p *) + у ln (1- u * р */ a*). Введя безразмерный параметр ц = u*/c*, с учетом (2.3) получим
А S = ln [ 1+ уц( 1 + ц + |ц|)] + у ln [ 1-ц/( 1 + ц + |ц|)]. (2.7)
Рассмотрим ее для случая ц < 0. Тогда ц + |ц| = 0 и (2.7) имеет вид
AS = /_(ц) = ln (1+ уц) + у ln (1- ц). Ее производная определяется формулой
f (||) = _!___Y = ц ( Y + 1 ) Y (28)
f - (ц) 1 + уц -1 - ц -( 1 - ц) ( 1 + уц ). (2.8)
При ц = 0 имеем /-(0) = /- (0) = 0. Производная /- (ц) больше нуля на интервале -1/у < ц < 0. Для значений ц < -1/у функция не определена, /-(ц) имеет вертикальную асимптоту при ц = цв = -1/у. Стоит заметить, что в появлении цв = -1/у виновата формула (2.4): при ц < цв получалось бы P < 0.
Анализ (2.7) при ц > 0 сделан в [8]. Эскиз графика функции AS = /(ц) - сплошная кривая на фиг. 2 при у < 7. В случае у > 7 возникает фрагмент, изображенный штриховой кривой. Такое поведение /(ц) при ц < 0 означает, что в случае разбегающихся потоков в расчете распада раз
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.