Производные газодинамических функций за поверхностью сильного разрыва Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 17, № 1, 2012

Производные газодинамических функций за поверхностью сильного разрыва*

Е.С. Прохоров

Институт гидродинамики им. М.А. Лаврентьева СО РАН, Новосибирск, Россия

e-mail: prokh@hydro .nsc.ru

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

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

В газовой динамике под сильным разрывом понимают двумерную поверхность в пространстве, на которой функции плотности р, массовой скорости и, давления p и других параметров сплошной среды имеют разрывы первого рода. Величины разрывов данных функций не могут быть произвольными, а удовлетворяют некоторым соотношениям — уравнениям сильного разрыва [1]. Так, для сильного разрыва типа "ударная волна" эти уравнения выражают законы сохранения массы, импульса и энергии на ударном фронте (поверхности разрыва) и могут быть представлены в следующем виде:

p*(D - и*) = po(D - ио),

p* + p*(D - и*)2 = po + po(D - ио)2,

U* + p*/p* + (D - и*)2/2 = Uo + po/po + (D - ио)2/2. (1)

Здесь D — скорость фронта, U — удельная (на единицу массы) внутренняя энергия среды; индексами 0 и * обозначены значения газодинамических величин перед фронтом (в исходном состоянии) и непосредственно за фронтом волны соответственно. Отметим, что с учетом теплового эффекта за счет химических реакций эти законы сохранения применимы и для детонационного фронта (сильного разрыва с тепловыделением) [2, 3].

* Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00634а) и Фонда Президента РФ по государственной поддержке ведущих научных школ (НШ-5770.2010.1).

Если движение среды за фронтом описывается гладким одномерным решением, а параметры перед фронтом постоянны, то можно установить однозначное соответствие между частной пространственной производной (ЧПП) функции любого газодинамического параметра (öy/ör)*, где y = {p,u,p}, и производной по времени скорости (ускорением) фронта dD/dt. Для одномерного адиабатического течения совершенного газа такие ЧПП за фронтом ударной волны приведены в [4], где они использованы для построения теории затухания ударных волн. Уместно также отметить близкую по тематике работу [5], В последней получены ЧПП газодинамических функций за искривленной стационарной ударной волной, на которую набегает равномерный сверхзвуковой поток. Однако в этой работе влияние ускорения фронта на ЧПП не учитывалось.

Область применимости формул для ЧПП, представленных в [4], в значительной мере ограничена условиями модели совершенного газа [6], в рамках которой полагается, что газ является идеальным, химически инертным и имеет постоянный показатель адиабаты y = (д lnp/д ln p)s = const, где S — энтропия вещества. Например, полученные ранее формулы не могут быть использованы для таких важных в практическом отношении случаев как движения газа за сильными ударными волнами, когда возможны 1) нарушение условий его идеальности, 2) возбуждение дополнительных степеней свободы и диссоциация молекул и 3) равновесные течения реагирующих газов за фронтом детонации, распространяющейся в химически активной среде,

В настоящей работе удалось снять эти ограничения путем использования естественного предположения к форме записи внутренней энергии U, основанного на том, что с учетом термического уравнения состояния (например, уравнений Клапейрона — Менделеева или Ван-дер-Ваальса) полную внутреннюю энергию U = Uth + Uch, включающую в себя кроме термодинамической части Uth и потенциальную химическую энергию Uch, можно представить в виде функции давления и плотности U(p,p). Это имеет место как для инертных сред, так и для продуктов реакции в состоянии химического равновесия (см, например, [7]), В частности, для совершенного газа имеем U = 1/(y — 1) ■ p/p+const. При такой форме записи внутренней энергии с учетом первого начала термодинамики равновесная скорость звука c в среде определяется уравнением

c2 = (dp/dp)s = (p/p2 — Up)/Up, Up = (dU/dp)p, Up = (dU/dp)p. (2)

При выводе искомых соотношений для ЧПП будем использовать следующие дополнительные предположения. Пусть в одномерном приближении ударная или детонационная волна с мгновенной химической реакцией за фронтом движется в сплошной (в общем случае химически активной) среде с заданными газодинамическими параметрами, Условие о мгновенности химической реакции означает, что за фронтом реализуется равновесное течение продуктов реакции. Полагаем, что U(p, p) известная, в общем виде неявная функция. Отметим, что при детонации в силу необратимости химической реакции U0 = U(p0,p0).

D

относительно производных dy*/dD, где y* = {p*,u*,p*}. Решая эту систему, находим:

dp* _ р* - ро р*^ du* _и*-щ ^ dp* _ р* -р0 ^ | ^ dD D — u0 p0 ' dD D — u0 ' dD D — u0 '

_ 2 + (1/p* — 1/po)/(Up)* (

(1/M,)2-1 ' W

где М* = (Д — и*)/с* — число Маха относительного потока за фронтом волны. Соотношения (3) можно записать и в более компактной форме

¿1п(р» -р0) _ р» - Ц0) ¿Лп(р» ~Ро) _ л о

(¿1пМ0 ~ Ро ' сЛпМ0 ~ ' сЛпМ0 ~ '

где М0 = (Д — ио)/со — число Маха относительного потока перед фронтом,

В работах [1, 3] показано, что если сплошная среда обладает "нормальными" [8] термодинамическими свойствами, то для ударных и детонационных волн выполняются естественные неравенства

Р* > Ро, Р* > Ро, и* > «о, Мо > 1, М* < 1. (4)

Отметим также, что безразмерная величина А в формулах (3) всегда положительна (А > 0), поскольку только в этом случае плотность р* будет возрастать с увеличением скорости фронта > 0). В частности, для совершенного газа А = 2/(М^ — 1),

Далее при выводе соотношений для ЧПП воспользуемся методом, описанным в [4]. Так, для одномерных движений любой газодинамический параметр среды зависит только от одной пространственной координаты г (г > 0) и времени ¿, т. е. имеем у (г, ¿). Если определять текущую координату фронта волны г* в виде монотонно возрастающей функции от времени (волна распространяется от центра) г* (¿), то для параметров фронта у* = {р*,и*,р*} получим зависимости у* = у (г* (£),£) = у*(£). При этом полные производные параметров по времени связаны с первыми частными производными аналогичных параметров за фронтом по £ (ду/д£)* и коордипате г (ду/дг)* следующим образом:

(¿у* = (¿У* , п(!?У\

(Й (Ш ' <И \дг Д \дг у; 1 )

Здесь учтено, что скорость фронта Д =

Одномерные уравнения нестационарной газовой динамики для адиабатических плоских движений (М = 0), а также движений с осевой (М = 1) и центральной (М = 2) симметрией могут быть приведены к виду [9]

др др ди ди ди 1 др

дг дг дг дг дг р дг

(6)

где т = риМ/г, я = 0 е = 0. Поскольку приведенные уравнения справедливы и для течения газа, примыкающего к фронту волны, то, исключая в (6) с помощью (5) частные производные (ду/д£)*, получим систему линейных уравнений относительно (ду/дг)*. Решая эту систему с учетом (3), найдем искомые ЧПП функций газодинамических параметров за фронтом волны (поверхностью сильного разрыва):

1 f Ро /су л 0\и* - и0 (10 т* е* ]

+ 6) 3. 17 + — + — Н

М2 — 11 р* с* ^ р* с* р*с*

(9Р\ 1 f Po г Л , О , / Л , 1\/Л/г21Р*~Р0 dD , Л/[ , • , М*е*1

(äJ, = +2 + +■ -ж + м*с*т*++ —)' (7)

где m* = p*u*N/r*, i* = 0, e* = 0, По-видимому, (7) является наиболее компактной формой записи для (dy/dr)*. При попытке преобразовать с помощью уравнений (1) данные уравнения становятся более громоздкими и теряют наглядность. Косвенным подтверждением правильности вывода этих соотношений может служить следующее: при подстановке в (7) известных формул, соответствующих уравнениям (1) в совершенном газе [1]

р* = Р-щ = (7+1)М§ р* = _ Jjzl

Po D — u* (7 — l)Mg + 2' po 7+I ° 7+1'

M2- 2 + (T~l)Mg 2 P* 2_ Po ( )

_ 27M0 — (7 — 1)' C°-7po'

рассматриваемые соотношения сводятся к уравнениям, полученным раннее в [4] для случая, когда газ перед фронтом ударной волны находится в состоянии покоя (uo = 0), Наиболее просто соотношения для ЧПП выглядят при изотермических движениях газа (p/p = const). Для этого достаточно в формулах (8) положить y =1- Тогда из (3) находим

dp*/dD = 2po Mo/co, dp*/dD = 2poCoMo, du*/dD = 1 + 1/M2 (9)

и соответственно из (7) с учетом dD/dt = codMo/dt получим

dp\ _ po Mo(Mo + 3) dMo ..2

ft'ÜR df~poM°

3M2 + 1 dMo ( Mo uAN

'Mo(Mg-l)' dt CoM°{l + Ml-l' cojr*

(10)

Обратим внимание на то, что в (10) производные (ду/дг)* ^^то выражаются через М0 и при N = 0 эти формулы еще более упрощаются.

Выражения (7) могут быть обобщены и на другие виды движения среды. Так, в случае квазиодномерных течений газа в трубе с переменной площадью поперечного сечения о = а (г) достаточно заменить в выражении для т системы (6), а значит и в решении (7) отношение ^г на производную ^(1п о)/dг. Заметим, что в канальном приближении нетрудно учесть также потери на трение и теплоотвод в стенки трубы [10]. В этом случае правые части г и е уравнений системы (6) уже не равны нулю, а являются некоторыми функциями, зависящими от параметров потока р,и,р и гидравлического диаметра канала, который для трубы круглого сечения совпадает с обычным. Однако такие изменения не влияют на характерную структуру соотношений (7),

Как показывает анализ, во всех рассмотренных случаях эти соотношения представляют собой формулы типа (ду/дг)* = g* • dD/dt + Ь*, которые устанавливают линейную связь между первыми ЧПП параметров за фронтом волны (ду/дг)* и ускорением самого фронта dD/dt . В общем виде коэффициенты пропорциональноети g* легко получить из (7), Они обладают универсальностью и не зависят от внешних условий движения среды (т.е. вида выражений для т,г,е). При заданных начальных параметрах р0,и0,р0,и0 конкретные значения коэффициентов определяются функциональной зависимостью для внутренней энергии и (р, р) и скоростью фропта D. Для большинства

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

Тогда из (7), в частности, следует, что для адиабатических плоских, цилиндрически и сферически симметричных видов движения при ускорении фронта волны (¿В/¿к > 0) плотность, давление и массовая скорость газа уменьшаются при приближении к фронту, т. е. (ду/дг)* < 0.

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

В работе [11] для расчетов в качестве химически реагирующего газа использована стехиометрическая смесь ацетилена с кислородом, для которой внутренняя энергия и(р, р) определена в соответствии с [7]. На рис. 1 представлены характерные профили плотности (а) и массовой скорости (б') за фронтом ДВ, где значения газодинамических величин отнесены к соответствующим параметрам рс7 и исз-> которые реализуются при выполнении условия Чепмена — Жуге (М* = 1) [3] на фронте установившегося и самоподдерживающегося режима детонации. Координаты Ь0 и выделяют область конического сужения трубы. Время I отсчитывается от момента инициирования ДВ у закрытого конца широкой трубы (г = 0). Графики на данном рисунке позволяют также проследить за положением детонационного фронта (вертикальные линии) в различные моменты времени.

Соответствующая пространственно-временная диаграмма г* (£), описывающая положение фронта при распространении ДВ в узкой трубе (г* > изображена на рис. 2.

а б

Рис. 1. Профили плотности (а) и массовой скорости (б) за фронтом затухающей пересжатой ДВ в различные моменты времени сплошные линии — расчеты [11], штриховые — восстановленное решение за фронтом по зависимости г* (¿)

Эта диаграмма может служить аналогом получаемых в эксперименте фоторазверток самосвечения (см., например, [12]) фронта ДВ через щель, параллельную оси трубы. Поэтому, если на рис. 2 провести оцифровку точек кривой г* (Ъ) и использовать методы численного дифференцирования табличных функций [13], как при обработке экспериментальных данных, то можно для конкретного момента времени с достаточно высокой точностью (погрешность около 1 %) рассчитать скорость фронта В = ¿г*/¿Ъ и его ускорение ¿В/А; = А2 г* /¿Ъ2.

Далее, используя В, определим значения газодинамических параметров на фронте У* = {р*, и*,р*} и с помощью (7) найдем соответствующие ЧПП (ду/дг)*. Это дает возможность решить обратную задачу — восстановить решение за фронтом ДВ: у = у* + (г — г*) х (ду/дг)* (штриховые линии на рис. 1). Видно, что градиенты газодинамических функций вблизи фронта для исходного [11] и восстановленного решений практически совпадают. Таким образом, представленный вычислительный эксперимент может служить основанием для применения соотношений (7) при обработке реальных экспериментальных данных.

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

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

Сеточный шаблон, обладающий достаточной универсальностью для реализации различных разностных схем, показан на рис. 3. Координаты точек (узлов) подвижной сетки "нижнего" временного слоя Ъ0 обозначим через гп—/2 = 0,1, 2,...), а координаты узлов "верхнего" временного слоя Ъ = Ъ0 + т (г — шаг по времени) — через гп—/2. Для вычисления координат точек гп—/2 промежуточного слоя Ъ = Ъ0 + т/2 применяется естественная формула гп—/2 = (гп—/2 + гп—/2)/2. Для удобства последующего изложения полагаем шаг по координате постоянным, т. е. гп—у2 — гп = гп—/2 — гп = — ^/2. Понятно, что в случае каждой конкретной схемы будут задействованы лишь некоторые узлы для перехода со слоя Ъ0 на елой Ъ = Ъ0 + т.

0.8

0.7

0.6

0.5

0.4

г

150 200 250 300

Рис. 2. Зависимость координаты фронта Рис. 3. Сеточный шаблон разностной схемы ДВ от времени

Пусть правая граница — плоская (М = 0 ) ударная волна, распространяющаяся по неподвижному (й0 = 0) изотермическому газу, для которого выполняются уравнения (9) и (10), Задаем при ¿0 произвольно положение фронта волны г*, его скорость Б = с0М0 (а следовательно, и фронтовые параметры у*) и ускорение дБ/дЬ = с0дМ0/дЬ. Полагаем, что на "нижнем" слое координата граничного узла гп = г*. Тогда, используя соотношения для производных (9) и (10), можно построить сеточное решение для газодинамических функций вблизи границы на "нижнем" уп_у2 и "верхнем" уп_^/2 слоях с точностью до членов второго порядка малости по т и к:

где Дг = г* + Вт + (дБ/д£)т2/2 — гп + О(т3) — погрешность в расчете по разностной схеме координаты фронта на "верхнем" слое.

Предположим, что для численных расчетов течения газа используется схема Годунова [14]. Приведем для этой схемы разностные уравнения (которые аппроксимируют дифференциальные уравнения системы (6), записанные в дивергентной форме) с учетом движения узлов сетки (см. рис. 3) для перехода с "нижнего" слоя на "верхний":

2 _ рк_./2) ^{\рк(йк - Ьп)] - [рк-1(йк_1 - Ё>п)}} = 0, к = п, п — 1, п — 2,...,

\рк~з/2 йк~]/2 - рк-з/2 ик-ц2) + \{[Рк + Ркйк{йк - £>п)]~

— [рк_ 1 + Рк_1йк_1(йк_1 — Бп)]} = 0, (12)

где Бп = (гп — гп)/т — скорость движения граничного узла. Для выбранной схемы решения разностных уравнений ук_1/2 (т. е. величины рк_1/2, йк_1/2 и рк_1/2 = е^рк_1/2) определены только в узлах с полуцелым номером, а функции ук на промежуточном слое — в узлах с целым номером. Значения ук вычисляются через известные величины Ук_1/2 и Ук+1/2 с "нижнего" слоя при решении задачи о распаде разрыва.

Ограничиваясь точностью порядка О (к) при решении задач о распаде разрыва, можно оценить погрешность решения ук_1/2 разностных уравнений (12) за один шаг по сравнению с аналитическим решением (11) на "верхнем" слое. Так, для ближайшего (к = и) к границе узла получим

Р-1'2 - Р-1'2 = -А>(Мо - ' § + 0(т2 + Л2),

Г+

т.е. погрешность за один шаг имеет порядок О(т), В то же время аналогичная погрешность для внутренних (к = и — 1, и — 2,...) узлов составляет величину О(т2 + к2). Таким образом, вблизи подвижной границы типа "ударной волны" наблюдается снижение порядка аппроксимации для рассмотренной разностной схемы.

Полученные соотношения для ¿у*/^В и (ду/дг)* позволяют находить ЧПП газодинамических функций более высоких порядков. Так, из (5) и (7) следует (дифференцирование по = д/д£ + В • д/дг)

V д£2 / ¿В2 \ ¿П ¿В ^^ V д^ / V дгд£ / V дг2 /

\ / * \ / \ / * \ / * \ /*

=<г/ау\ /зу, \drdt), л\ат), \дт*); 1 '

где ¿2у*/^В2 вычисляются с помощью (3), Далее, дифференцируя (6) по времени (д/д£) и учитывая (13), получим систему линейных алгебраических уравнений относительно вторых частных производных по координате. Решая данную систему, находим вторые

(д2у/дг2)* (ду/дг)*

второй производной скорости фронта по времени ¿2В/^2 , В общем виде это решение выглядит достаточно громоздким, поэтому здесь приведем его только для частного случая, когда плоская (М = 0) ударная волна распространяется в изотермическом газе и справедливы соотношения (9) и (10):

/д2р\ 1 /д2р\ 4роМ3(М2 + 1) ¿2Мо роМ2(М0 - 10М0 - 27Мо - 12) /¿Мо\2

А2/* со\дг2/* с2(М0 - 1)2 с2(М2 - 1)3 V ^ ) '

_ (М2 + I)2 + 4М2 сРМ0 Ж0(М40 + 10Мо + 5) ((МЛ2 , ~ ср(Мр - I)2 Ж ср(Мр - 1)з •

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

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

Список литературы

[1] Ландау Л.Д., Лифшиц Е.М. Теоретическая физика: Гидродинамика. М.: Наука, 1986.

[2] Зельдович Я.Б. Теория горения и детонации газов. Москва, Ленинград: Изд-во АН СССР, 1944.

[3] Митрофанов В.В. Детонация гомогенных и гетерогенных систем. Новосибирск: Изд-во НГиЛ СО РАН, 2003.

[4] Седов Л.И. Методы подобия и размерности в механике. М.: Наука, 1963.

[5] Русанов В.В. Производные газодинамических функций за искривленной ударной волной. \!.. 1973. (Препр. АН СССР. Ин-т прикладной математики. № 18).

[6] 15vмер Ю.Б., Рыбкин М.Ш. Термодинамика, статистическая физика и кинетика. М.: Наука, 1977.

[7] Николаев Ю.А., Фомин П.А. О расчете равновесных течений химически реагирующих газов // Физика горения и взрыва. 1982. Т. 18, № 1. С. 66-72.

[8] Овсянников Л.В. Лекции по основам газовой динамики. М.: Наука, 1981.

[9] Станюкович К.П. Неустановившиеся движения сплошной среды. М.: Наука, 1971.

[10] Гинзбург И.П. Прикладная гидрогазодинамика. Л.: Изд-во ЛГУ, 1958.

[11] Ждан С.А., Прохоров Е.С. Квазиодномерный расчет детонации в канале переменного сечения // Физика горения и взрыва. 1984. Т. 20, № 5. С. 96-100.

[12] Дувовик A.C. Фотографическая регистрация быстропротекающих процессов. М.: Наука, 1964.

[13] Калиткин H.H. Численные методы. М.: Наука, 1978.

[14] Решение одномерных задач газовой динамики в подвижных сетках / Г.Б. Алалыкин, С.К. Годунов, И.Л. Киреева, Л.А. Плинер. М.: Наука, 1970.

Поступила в редакцию 29 июня 2011 г.