CC BY
18
УДК 532.5
УСТОЙЧИВОСТЬ ФРОНТА ИСПАРЕНИЯ ПРИ ИЗМЕНЕНИИ ЧИСЛА МАХА В ПОТОКЕ ИСПАРЕННОГО
ВЕЩЕСТВА
И. Н. Карташов, А. А. Самохин
На основе полученного ранее дисперсионного уравнения для возмущений фронта испарения исследуется поведение дисперсионных кривых при различных значениях числа Маха М в потоке испаренного вещества. Показано, что температурные возмущения, вызывающие неустойчивость при М = 1, оказывают стабилизирующее влияние на неустойчивость типа Даррье-Ландау при М < 1.
В работе [1] получено дисперсионное уравнение для возмущения плоского фронта испарения в диапазоне изменения числа Маха М в потоке пара от нуля до единицы. В предельном случае М <С 1 и при постоянстве температурного профиля в конденсированной среде это уравнение переходит в известное уравнение Даррье-Ландау для устойчивости фронта медленного горения, а при стремлении М к единице и условии экстремальности потоков массы дг, импульса д^ и энергии дз получается рассмотренное ранее уравнение [2], в котором газодинамические эффекты не проявляются. Случай произвольных чисел Маха 0 < М < 1 требует численного анализа, который проводится в настоящей работе для поверхностного поглощения излучения в конденсированной среде.
Пусть слева от плоского невозмущенного фронта испарения, движущегося в отрицательном направлении оси г, находится испаряемая под действием поверхностного источника нагрева конденсированная среда с постоянными температуропроводностью х, теплоемкостью с и плотностью описываемая уравнениями теплопроводности, Эйлера и непрерывности для несжимаемой жидкости. В системе отсчета, движущейся вместе с фронтом испарения с постоянной скоростью —V по неподвижному веществу, стационарный температурный профиль определяется известным выражением:
Г/ = АТ[ехр(к02) — 1] + Т3, где АТ = Т$ — Т^, Т3 и Тоо - соответственно температура на поверхности г = 0 и в глубине конденсированной среды, к0 = и/Х- Скорость V, зависящая от Г, и М, определяется балансом энергии в конденсированной фазе I = р[ь[Ьпе + сДТ], где I - поглощаемая интенсивность излучения, Ьпе - теплота испарения. Значения давления р/ и скорости протока V в этом случае постоянны, поскольку тепловое расширение конденсированного вещества не учитывается.
В области г > 0 находится испаренное вещество, которое считается идеальным одноатомным газом с показателем адиабаты 7 = ср/с„ = 5/3 и описывается системой газодинамических уравнений непрерывности, Эйлера и адиабатичности, невозмущенным решением которой являются постоянные значения плотности пара р, давления р и компонент скорости их — 0 и2 = и = д1/р — Мис, где и2с = 7р/р.
Область вблизи плоскости г = 0, называемая кнудсеновским слоем, требует кинетического описания. Размеры этой области порядка длины свободного пробега и для возмущений с достаточно большой длиной волны ее можно рассматривать как разрыв, на котором формулируются законы сохранения массы, импульса, энергии и два дополнительных соотношения, учитывающих неравновесность процессов в кнудсеновском слое.
Для возмущений физических параметров вида ехр(о;£—1кх) линеаризованная система уравнений для конденсированной среды формулируется следующим образом:
дбу
—гк8ух + —— = О, ог
с д8уг 1 дбр!
ШбУх + и-7— =---,
ог р1 ог
2 к
Ш8Ух + и-д-^ = —8р1, (1)
ог р1
д6Т'мХ дТ,м 1,4т дЧТ1 п
и при условии отсутствия волн, идущих к границе раздела из бесконечности, имеет решение:
6иг = 6ьекг, 8ух = — ¡8ьекг, 6р1 — — -¿-(ш + ку)8гекг,
к
¿Г, = коАТ6у(еЯ2 - е(к+ко)2) + (6Тв - (к0АТ)е"% (2)
ш — кг>
где q удовлетворяет уравнению
Я2 -дко-к2-и/Х = 0 (3)
с условием Лед > 0, а 6Т3 = ¿71(0) + ^дТ\\дг - модуляция температуры с учетом смещения £ возмущенной поверхности раздела относительно г — 0.
Решение линеаризованной системы уравнений движения для пара
г их, д8и* ■ д8Р п шор - гкроих + р—--^ и—— — 0,
ог ог
д8их гк шдих + и—-— = —¿>р, ог р
д8иг 1 дёр
шди2 + и—— =---—, (4)
ог р ог
2 I
ш + и— (6р - ис6р) = 0
записывается в виде:
6р = £ А,-е-*", 8р = + А3е~кзг,
i=l,2 Uc
Suz = £ , к\ + (5)
8и* = £ , гк U +
»=1^2 - ад
где волновые числа с г = 1,2 и кз определяются соотношениями:
к3 = ш/и, (и\ - и2)к2 + 2u>uki - и2 - к2и2с = 0. (6)
Устремляя во втором соотношении Reo; к бесконечности, т.е. переходя к пределу мгновенного включения источника в момент времени t — 0, можно убедиться, что при дозвуковом течении одна из волн (с амплитудой А\) соответствует распространению от границы раздела, а другая с А2 соответствует волне, идущей из бесконечности [3]. Аналогично, сносовые решения Аз и Л4 также имеют направление распространения от фронта перехода к бесконечности в потоке пара. Следует отметить, что приведенные в [4] рассуждения о том, что для неустойчивого фронта пламени сносовое решение должно быть отброшено в области, где возмущения нарастают от границы, являются не совсем корректным. Если следовать такой логике, то для затухающих возмущений должно быть оставлено нарастающее в пространстве решение [5]. Однако, как легко проверить, такой выбор для затухающих на фронте возмущений соответствует волне, идущей из бесконечности, и не соответствует постановке задачи. Выбор решений по
направлению распространения волны в общем случае производится, исходя из анализа дисперсионного уравнения, если устремить Яео; —► оо, т.е. в пределе мгновенного включения источника возмущений в начальный момент времени, что соответствует начальной постановке задачи.
Условие отсутствия волн, идущих из бесконечности, позволяет положить Аг — 0 в диапазоне М < 1. Остальные шесть постоянных определяются из граничных условий для возмущенных величин в плоскости г = О, которые могут быть записаны в виде:
8у2 = + бдг/р1, 8иг = + 8и, 8р1 = 8д2 - 2у8дх + ак2£, 8р = 8(д2 - дги),
+ ср1Х + = 0, (7)
8их — 6ух + 1к{и —
Эти условия представляют собой законы сохранения потоков массы, импульса, энерп л, а также условие непрерывности тангенциальной к поверхности возмущенного фронта компоненты скорости. Поток массы д\, как и связанная с ним величина V = д\/Р1, и поток импульса <72 зависят от Т„ и М, причем зависимость эта определяется кинетическими процессами в кнудсеновском слое. Таким образом может быть получена однородная система из четырех уравнений относительно 8у, 8Т3 и ¿М, условие разрешимост:; которой и определяет дисперсионное уравнение для возмущений фронта испарения [ 1 ]:
det\\Dij\\ =0, (8)
в котором элементы определителя определяются следующими выражениями:
п 1 п п 1 п 1 Бц = 1, £>12 = -ш, £>13 =--^г, Оы =--—
Р1о15 р1()М
к0АТ(д - к - к0)
Аи = ср1Х-:-, -032 = ср,хк0АТ{к0 -
ш — ки
д(д1Ьпе) д(д1Ьпе)
А* = СР1ХЧ + Д34 = (9)
ш2
£>41 = 1, -£>42 = ^ - ¿(м - и),
_ ш du F <9(ff2 — giu) _ О) ди F д(д2 - д\и)
43 ~ ku dT3 ри dTs ' 44 _ ки дМ ри дМ '
где F = к(ш~-к\и) • Уравнение (8) приводится к алгебраическому уравнению 16 степени относительно величины q, определяемой соотношением (3).
При решении этого уравнения использовались испарительные граничные условия, определяемые следующим образом [6]. Для частиц, вылетающих с поверхности, функция распределения бралась в виде распределения Максвелла с температурой Ts и концентрацией ns: /W = /o("s;?s), а для возвращающихся частиц в виде: = a7/0(ns, а2Т,). Потоки массы g\ и импульса д2 определяются как соответствующие моменты полной функции распределения + причем параметр а находится из соотношения
(l-a8)(l-a10) _х72М2[(7-1)М2+2]
(1 + а9)2 8 (7-1)(1+7М2)2 ' 1 ;
следующего из условия баланса энергии на кнудсеновском слое. Использование такой модели для получения испарительных граничных условий обеспечивает экстремальность потоков при М — 1 и удобно в применении. Для численного анализа использовались следующие значения параметров:
pi = Юг/см3, х = 0Лсм2/с, с = 1.5 ■ Ю6эрг/г • К, а = 220дин/см,
L = 8.8 ■ 109эрг/г, 7 = 5/3, cp = 5kB/2m, т = 3.46 • 1(Г22г, Ps = nskBTs = pbexp(rj(l - Tb/Ts)), г] = 10.9, рь = Ю6дин/см7, (11)
Too = 300А', Ть = 2023 К.
При М = 1 задача для конденсированной среды оказывается замкнутой, поскольку газодинамические эффекты в потоке пара не оказывают влияния на возмущения фронта перехода. В такой постановке проблема устойчивости исследовалась в работе [2]. На рис. 1 этому предельному случаю М —> 1 при выбранных испарительных граничных условиях и температуре поверхности Ts = 1.5Ть соответствует кривая 1, а кривые 2- i относятся к значениям чисел Маха М = 0.9; 0.7; 0.5 соответственно. Скорости фронта испарения v в этих случаях равны 36.83 (1), 36.75 (2), 35.9 (3) и 33.1 см/с (4). На рис. 2 приведены аналогичные кривые для тех же значений М и другой температуры поверхности Ts = 2Ть- Скорости фронта испарения при этом оказываются больше: 196.2 (1), 195.8 (2), 191 (3), 176.1 см/с (4).
Рис. 1. Зависимость инкремента неустойчивости от волнового числа при Т, = 1.ЪТЬ для М = 0.999 (1), 0.9 (2), 0.7 (3), 0.5 (4).
Рис. 2. Зависимость инкремента неустойчивости от волнового числа при Т= 2Ть для М = 0.999 (1), 0.9 (2), 0.7 (3), 0.5 (4).
Отметим, что при таких скоростях испарения дисперсионные кривые для случая, рассмотренного в работе [2], уже не содержат области неустойчивости. Отличия в границах области неустойчивости, а также в величине максимального инкремента в этих случаях связаны с различиями в параметрах, используемых при расчетах, и с довольно резкой зависимостью решений от этих параметров. Дальнейшее увеличение скорости испарения приводит, как и в [2], к уменьшению, а затем и исчезновению этой области неустойчивости.
Кривые 1 на рис. 1,2 с графической точностью совпадают со случаем М = 1, когда поведение конденсированной фазы можно рассматривать отдельно от газодинамических возмущений и дисперсионное уравнение существенно упрощается [2]. При уменьшении числа Маха максимальная величина инкремента и размер области исходной неустойчивости уменьшаются. На рис. 2 кривые 3 и 4 эту область уже не содержат.
В то же время в длинноволновой части появляется и растет дополнительная область неустойчивости, связанная с газодинамическими эффектами. Дисперсионные кривые, которым соответствуют значения числа Маха 0.3 (1); 0.2 (2); 0.1 (3) приведены на рис. 3 для Тя — 1.5Ть, а на рис. 4 для Т3 = 2Ть при тех же значениях М. На этих рисунках исходная область неустойчивости уже полностью отсутствует, а максимум инкремента в "газодинамической" зоне неустойчивости и ее ширина уменьшаются с уменьшением М.
Рис. 3. Зависимость инкремента неустойчивости от волнового числа при Т, = 1.5Ть для М = 0.3 (1), 0.2 (2), 0.1 (3).
Рис. 4. Зависимость инкремента неустойчивости от волнового числа при Т, = 2Ть для М = 0.3 (1), 0.2 (2), 0.1 (3).
Рис. 5. Зависимость инкремента неустойчивости от волнового числа при Т, = 1.5Ть, ДТ1 = 0 для М = 0.3 (1), 0.2 (2), 0.1 (3).
Рис. 6. Зависимость инкремента неустойчивости от волнового числа при Т, = 2Ть, АТ = О для М = 0.3 (1), 0.2 (2), 0.1 (3).
Предельный случай малых чисел Маха М С 1 без учета температурных эффектов описывается известной задачей Даррье-Ландау [4]. В рассматриваемом нами более общем случае для перехода к этому пределу необходимо исключить влияние температурных эффектов, что достигается наложением условия постоянства температуры в невозмущенной конденсированной среде, т.е. при ДГ = Т„ — Тоо = 0. На рис. 5, б приведены зависимости инкремента неустойчивости возмущений фронта перехода от волнового числа к для тех же значений М и Г3, что и на рис. 3 и 4, но при ДТ = 0.
Кривые 3 на рис. 5, 6 с графической точностью совпадают с дисперсионной кривой в задаче Даррье-Ландау с учетом стабилизационного эффекта от поверхностного натяжения на фронте перехода. В то же время исключение температурных эффектов на кривых из рис. 5, 6 приводит к увеличению максимального инкремента и области неустойчивости по сравнению с рис. 3, 4, поскольку температурные эффекты в данном случае оказывают стабилизирующее влияние.
Дисперсионное уравнение для этого случая имеет вид:
Отметим, что при произвольных значениях М условие АТ = 0 не исключает, вообще говоря, влияние температурных эффектов, поскольку модуляция числа Маха 8М вызывает соответствующую модуляцию температуры поверхности даже при А / = 0. Однако, при АТ = 0 и М = 1 влияние температурных эффектов, как и при М «С 1, исчезает. Дисперсионное уравнение в этом случае принимает вид:
Это уравнение отличается от обычного дисперсионного уравнения для капиллярных волн на поверхности неподвижной жидкости наличием дополнительного члена шку, связанного с протоком массы через границу раздела и определяющего затухание поверхностных волн в этом случае. В то же время в уравнении Даррье-Ландау аналогичны: член, связанный с затуханием, имеет вид 2шку, и может возникнуть вопрос о том, како ва же роль потока массы через границу раздела (абляции) в затухании поверхностных волн на этой границе.
Отметим, что вопрос о стабилизирующем влиянии абляции неоднократно обсуждался в ряде работ по неустойчивости Рэлея-Тэйлора в проблеме лазерного управляемого термоядерного синтеза [7-9]. При этом, в частности, приводились качественные соображения, поясняющие природу абляционного затухания, из которых следует, что в дис персионном уравнении появляется член 2шку, т.е. скорость затухания равна ку (см., например, [7], формула (3)). Смысл этих соображений сводится к тому, что за счет абляции поверхность раздела уходит вглубь конденсированной среды за время на величину vAt и попадает в ту область, где амплитуда поверхностных волн оказывает ся уменьшенной на фактор ехр(—в соответствии с законом убывания амплитуды капиллярной волны в глубь неподвижной жидкости.
Подобные соображения, однако, не могут служить основанием для количественного определения величины абляционного затухания. Прежде всего, следует иметь в виду.
ш2 + 2 шку + ак3 / р1 — к2иу — 0.
(12)
ш2 + шку -+- ак3/р1 — 0.
(13)
что уравнение (12) отличается от уравнения (13) не только дополнительным членом шку, который дает коэффициент 2, но и еще одним членом к2иу, который также связан с газодинамическими эффектами, обусловленными в том числе и протоком массы. Другими словами, для определения роли собственно абляции в уравнении (12) требуется дополнительное исследование, тогда как в уравнении (13) затухание обусловлено только абляцией, поскольку другие газодинамические эффекты в этом случае просто не проявляются из-за соответствующих граничных условий.
В коротковолновом и длинноволновом пределах уравнения (12)—(13) дают соответственно:
<^1,2 = — ку ± ¡^ак3/р1, к^>р[Уи/а,
и>1,2 = —ку ± ку/тг, к <С р^и/сг, (14)
и>1,2 = —кг/2 ± г^/стк3/р[, к^>р[У2/а,
^1,2 = 0; —ку, к < р1У2¡ст. (15)
Из этих соотношений видно, что разные предельные случаи (коротковолновый и длинноволновый пределы) дают разные значения для декремента затухания. По этой причине любые соображения о величине обусловленного абляцией затухания, в которых не учитываются конкретные значения величины к, являются, по крайней мере, недостаточными.
Что же касается "вывода" величины затухания, основанного на использовании затухающего в глубину профиля поверхностной волны на неподвижной жидкости, то в этой процедуре имеются неявные предположения, которые не согласуются с более строгим рассмотрением задачи. В частности, получение фактора ехр(—АгиД£) основано на предположении, что поле скоростей привязано к неподвижной жидкости. В то же время в подобных задачах результирующее поле соответствующих физических величин перемещается по неподвижной жидкости. Например, температурный профиль в невозмущенной конденсированной среде содержит зависимость АТ ехр(ког), которая привязана к движущемуся фронту, а не к неподвижному веществу.
С учетом этого обстоятельства, приведенные в [7] рассуждения можно модифицировать следующим образом. Очевидно, что затухание экспоненциального профиля ~ ехр(ш£ — к\г\) эквивалентно его движению в сторону границы раздела со скоростью уэф, которая отличается от скорости потока у. Если предположить, что величина уэф = у — уэф, то величина затухания соответствует коротковолновому пределу
формулы (15), а при уэф = и, когда профиль скорости привязывается к неподвижному веществу, можно получить затухание, соответствующее длинноволновому пределу в (14), или коротковолновому пределу в (15). Однако подобные предположения, в свою очередь, требуют дополнительных обоснований, по существу выходящих за рамки элементарных качественных оценок. Недостаточность таких оценок связана фактически с неадекватным учетом неодномерности рассматриваемой задачи, когда анализ проблемы затухания целиком основывается на одномерном распределении исследуемой величины по г, причем распределение возмущений вдоль поверхности явно не рассматривается.
Существенную роль неодномерных эффектов в подобных "одномерных" оценках можно продемонстрировать также на примере неустойчивости фронта сублимации в случае объемного нагрева конденсированной фазы (см., например, [2] и цитированную там литературу). Для максимального инкремента неустойчивости такая оценка дает: со = -^г • Эта формула получается при использовании следующих соотношений Дг = Дг?Д£, Ау = -^¡гАТ3, АТ3 = д| = ^Д-г. Такой же результат получается
из соответствующего дисперсионного уравнения в пределе больших к за счет температуропроводности в плоскости фронта. При исключении же этого неодномерного эффекта дисперсионное уравнение дает совершенно другое выражение для инкремента, определяющее истинно одномерную (к = 0) неустойчивость фронта в этой задаче.
Таким образом, полученные в настоящей работе результаты свидетельствуют о существенном влиянии газодинамических эффектов на устойчивость фронта испарение при М < 1, которое приводит к появлению дополнительных областей неустойчивости на дисперсионных кривых, отсутствующих при М = 1. Уменьшение М сопровождается уменьшением инкрементов и сужением области неустойчивости, существовавшей при М = 1 и обусловленной тепловыми и гидродинамическими эффектами в конденсирован ной фазе. Однако роль этих эффектов оказывается существенной и при исчезновении исходной области неустойчивости, что проявляется в их стабилизирующем влиянии на развитие неустойчивости типа Даррье-Ландау.
ЛИТЕРАТУРА
[1] К а р т а ш о в И. Н., Мажукин В. И., Перебейнос В. В., Самохин А. А. Краткие сообщения по физике ФИАН, N 9-10, 22 (1996).
[2] С а м о х и н А. А. Труды ИОФАН, 13, 3 (1988).
[3] Ландау Л. Д., Л и ф ш и ц Е. М. Физическая кинетика. М., Физматлит, 2001.
[4] Л а н д а у Л. Д., Лифшип Е. М. Гидродинамика. М., Наука, 1988.
[5] 3 е л ь д о в и ч Я. Б., Баренблатт Г. И., Л и б р о в и ч В. Б., Махвиладзе Г. М. Математическая теория горения и взрыва. М., Наука, 1980.
[6]Мажукин В. И., Прудковский П. А., Самохин А. А. Математическое моделирование, 5(6), 3 (1993).
[7] К i 1 k е n п у J. D. et al. Phys. Pias, 1(5), 1379 (1994).
[8] V e 1 i к о v i с h A. L. et al. Phys. Pias, 5(5), 1491 (1998).
[9] Mikaelian K. 0. Phys. Rev. A, 46(10), 6621 (1992).
Институт общей физики РАН Поступила в редакцию 4 марта 2002 г.
