CC BY
16
УДК 532.5
ВЛИЯНИЕ КОЭФФИЦИЕНТА КОНДЕНСАЦИИ НА ИСПАРИТЕЛЬНЫЕ ГРАНИЧНЫЕ УСЛОВИЯ И УСТОЙЧИВОСТЬ ФРОНТА ИСПАРЕНИЯ
И. Н. Карташов, С. В. Орлов, А. А. Самохин
Исследуются изменения газодинамических граничных условий на фронте испарения в зависимости от величин коэффициентов конденсации и аккомодации. На примере воды показано, что от величины коэффициента конденсации существенно зависит инкремент неустойчивости для малых возмущений фронта испарения.
В работах [1, 2] исследовалась проблема устойчивости фронта лазерного испарения конденсированных сред в случаях поверхностного и объемного нагрева мишени и с учетом газодинамических возмущений в потоке испаренного вещества. При этом предполагалось, что пар является одноатомным идеальным газом, взаимодействие которого с поверхностью характеризуется коэффициентом конденсации /3=1. Величину /3 также называют коэффициентом поглощения, прилипания, испарения и аккомодации, а величину г — I — в - коэффициентом отражения частиц [3 - 7]. Ниже под коэффициентом аккомодации, в отличие от упомянутого выше, будет пониматься параметр, характеризующий степень термализации частиц при их отражении от поверхности. Для металлов приближение /3 = 1 считается вполне удовлетворительным. Однако в общем случае величина /3 < 1, то есть не все частицы, сталкивающиеся с поверхностью, на ней конденсируются. Например, для воды этот коэффициент, по-видимому, значительно меньше единицы, хотя имеющиеся экспериментальные данные о величине /3 характеризуются значительным разбросом [8 - 10].
Испарительные граничные условия зависят от величины /3, которая характеризует уменьшение потока частиц, испаряющихся с поверхности по сравнению с формулой Герца-Кнудсена. В настоящей работе исследуется эта зависимость и влияние измененных газодинамических испарительных граничных условий на поведение возмущении фронта испарения.
В процессе испарения между поверхностью конденсированного вещества и газодинамическим потоком пара существует область неравновесного течения (слой Кнудсена). В рамках газодинамического подхода эта область рассматривается как поверхность раз рыва, на которой формулируются граничные условия, связывающие температуру поверхности Ts и концентрацию насыщенного пара п3 с параметрами пара при заданном числе Маха М < 1 в потоке испаренного вещества. Вопрос об определении испарительных граничных условий и структуре кнудсеновского слоя в рамках различных подходов рассматривался во многих работах (см., например, [1, 11 - 16] и цитированную там литературу). Сравнительный анализ этих подходов проводился, в частности, в [14, 16].
В работах [2] по анализу устойчивости фронта испарения использовались модельные граничные условия [14], которые легко обобщаются на случай ¡3 < 1 с любым значением коэффициента аккомодации. В этой модели при /9 = 1 предполагается, что неравновес ная функция распределения fne для одноатомного газа на поверхности раздела внутри кнудсеновского слоя имеет следующий вид:
fne = /(+)К, Т.) = Дп„ Т,, 0) = f(ns, Ts), Vz >0, (1)
fne = f(-\nsle2Ta) = e7f(na, e2Ts), к < o, (2)
где / - равновесная максвелловская функция распределения для потока частиц с кон центрацией п, массой т и температурой Т', движущихся со средней скоростью и в положительном направлении оси z от границы раздела (кв ~ постоянная Больцмана):
« Т Л ( m \ъ/2 + + (Ц - «О'Л „л
/(n'T'u)=4wü ехр(—w—)■ (3)
Положительный параметр 0 < 1 характеризует интенсивность процесса испарения. При отсутствии испарения 9 = 1 и fne совпадает с равновесной функцией распределения /(ns,Ts, 0). Выражения для потоков массы <71, импульса д2 и энергии #3, вычисленные с помощью /пе, определяются следующими формулами:
(4)
/. ™ \ 1/2
gi=mj Vzfned3V = (1 - 08)mns f ^ j ,
9 з
Интегрирование по каждой из переменных Vx, Vy, Vz проводится по всей числовой о. . На внешней стороне кнудсеновского слоя течение пара становится равновесным с функ цией распределения (3) и выражения для потоков <7, принимают обычный газодинами ческий вид
gl — m J Vzf(n,T,u)d3V = mrnx, д2 = т J V2f(n, Г, u)d3V = nkBT + mnu2, (5)
9з = jJ V2Vzf(n,T,u)d3V = nu (^f- + с.квТ^ ,
где cp - это теплоемкость при постоянном давлении в расчете на одну частицу в единицах кв■ Приравнивание соответствующих потоков из формул (4), (5) дает систему уравнений, которая связывает газодинамические переменные р и Т с температурой поверхности Ts, концентрацией насыщеного пара ns и степенью неравновесности процесса испарения, определяемой параметром в или числом Маха в потоке пара M:
92 , т тп-уМ2 (д2\2
P = TTW *вГ=(1 + 7м>у UJ ' ( )
где выражения для потоков <71, д2 даются формулами (4). Число Маха определяется обычным образом M = u/us, где квадрат скорости звука u2s = 7квТ/m, а 7 = cp/cv есть отношение теплоемкостей при постоянном давлении и постоянном объеме. Связь между 9 я M получается из равенства двух следующих выражений Ge и Gne:
с _ 9з91 _ 72М2 (7 - 1 )М2 + 2
е gl (1 +7М2)2 2(7-1) '
дгдг .... 4 (1 - 0*)( 1 - g") 9l ~ * (1+П2 ' Напомним, что рассматриваемая модель, в отличие, например, от [12], обеспечивает экстремальность потоков по их зависимости от M при M = 1, а соотношение (7), у которого экстремум также локализован при M = 1, имеет общий характер и от вида модели не зависит.
Формулы (4), (8) справедливы только для одноатомного газа с 7 = 5/3. В случае многоатомных газов, кроме поступательного движения, необходимо учитывать также дополнительные степени свободы. В интересующем нас диапазоне температур можно пренебречь колебательными возбуждениями и ограничиться только вращением молекул. Полная равновесная функция распределения / вместо (3) теперь имеет вид произведения
/(п,Г,и) = /(п,Г,и)/г(Т), (9)
где ее вращательная часть /г зависит только от угловых переменных. В классическом пределе, когда справедлив закон равнораспределения, среднее значение вращательной энергии (Ег), вычисленное с помощью (9), дается выражением
(Ег) = (10)
Замена / на / в левой части формул (5) не меняет вида их правых частей, где
величина с„ = ——— = ——— определяется обычным образом с учетом полного числа 2 7 — 1
степеней свободы отдельной частицы 5, которое равно 5 и 6, соответственно для двух и многоатомных (нелинейных) молекул.
Для полной неравновесной функции распределения будем использовать вместо (1) п (2) следующие выражения:
/пе = = /(п„Т,)/г(Т,), К > 0, (11)
/пе = ^-\пя,в2Т3) = ^/К,02Г8)Л(02Г5), К < 0. (12)
Формулы (11), (12) основаны на предположениях, что температура вращательных степеней свободы равна "температуре" соответствующих частей неравновесной функции распределения. При этом потоки <71, д2, дз из (4) модифицируются таким образом, что вместо множителей (1 — в8), (1 + 09) и 2(1 — О10) в соответствующих формулах возникают
5 + 1
выражения (1 — 0Ж+1), (1 + вх+2) и ———(1 — 01+3), а для (7пе вместо (8) получаем
£
Г 9391 * + 1(1-0*+1)(1-0*+3) ,
91 ' Т (1+^)2 • ^
Как и в исходной модели, величина х определяется из условия равенства отношения неравновесных потоков д3 и д\ соответствующему отношению газодинамических потоков в равновесном пределе, когда 9 = 1 и М = 0:
= ср-. (14)
Л/—>0 т
9з 9\
Из этого условия следует, что х = 2з + 1. Таким образом учет неодноатомности пара в рассматриваемой модели в итоге сводится к изменению показателя степени у пара.чп гра 9 в поступательной части неравновесной функции распределения.
До сих пор мы считали, что все падающие на поверхность частицы конденсируются на ней. При неполной конденсации, когда ¡3 < 1, доля 1—/3 частиц отражается. Отраженные частицы могут полностью сохранять свою энергию (упругий удар) или, в другом предельном случае, отражаться с измененной энергией, соответствующей температуре поверхности (полная аккомодация). Для простоты мы будем считать, что распределение отраженных частиц может быть представлено в виде суперпозиции этих двух предельных случаев. В рамках данного подхода положительная часть неравновесной функции распределения принимает такую форму:
/(+) = (31(п„Т3) + (1-/3) [ец/{пя,Та) + (1 -е)вх/(п„е2Та)} , У2 > 0, (15)
где параметр £ означает долю потока отраженных частиц, которые при соударении с поверхностью полностью термализуются. Остальная часть потока 1 — £ отражается упруго. В последнем случае фактически отсутствует различие между зерка л, ным и диффузным упругим отражением из-за используемой нами конкретной формы /(~\пц,92Т3), которая не содержит явно выделенной скорости.
Параметр ц < 1 в правой части формулы (15) возникает, как и аналогичный коэффициент в работах [15, 16], из-за различия вида функций распределения и Я 5 92Т3) при необходимости выполнения условия равенства прямого и обратного потоков массы для той доли е отраженных частиц, которые изменяют свою энергию при отражении (полностью аккомодируются):
(XI УЛ(щ,Т3)<13У = | (16)
где интегрирование по У2 ведется по полуоси, а по остальным компонентам скорое ! и - по всей числовой оси. Из этого условия следует ц = 9Х+1. Выражения для потоков массы, импульса и энергии в этом случае принимают вид:
9\
= д[+) + д[~\ g. = gi+) + gt\ gs = gi+) + gt\ (17)
а выражения для потоков от поверхности и на поверхность д\ ' определяются со
отношениями:
д[+) = {р + (1-р)ех+1}91к, = {/? + (!- Р)[е<Г+1 + (1 - (18)
= + (1 - /?)[£Г+1 + (1 - е)9х+3]}дзК,
= = ^ = (19)
/квТЛ1/2 . + 1 . _ ДвГЛ1/2
где о1Л- = -- и оЗА: — —-—п3Кв1* ~- представляют собой потоки
\27rm) 2 \2тгт)
массы и энергии, определяемые распределением Максвелла в полупространстве скоростей К > 0 (кнудсеновские потоки). Вместо (13) с учетом (18), (19) теперь для Спе имеем
_ дз91 _ s + 1 {/? + (1 - + (1 - е)9^) - Г+3}{/?( 1 -
gl -к {/? + (1 — (3)[e9x+l + (1 — е)0х+2] + 0Х+2}2 '
Подставляя далее соответствующие потоки из (17)—(19) в соотношение (14), получаем значение показателя степени х в коэффициенте 9х у функциии распределения, описывающей обратный поток частиц
Ж =-0-' (21)
Из условия равенства выражений (20) и (7) получается уравнение Gne{9) = Ge(M), связывающее число Маха с параметром 9. Как и в предыдущих случаях, когда для Gne использовались выражения (8) и (13), это уравнение обеспечивает экстремальность потоков <7, по их зависимостям от М при М — 1. Несмотря на некоторое усложнение выражения (20) по сравнению (8) и (13), уравнение Gne = Ge по-прежнему элементарно разрешается относительно М.
На рис. 1 показано поведение нормированного потока массы дх/дщ в зависимости от числа Маха при различных значениях показателя адиабаты 7 и коэффициента конденсации /3. Наиболее существенно на это поведение влияет изменение параметра 3, уменьшение которого от единицы до значений 0.5 и 0.1 приводит к примерно такому же уменьшению величины потока (кривые 1,1'; 2,2'; 3,3'). Зависимость величины потока от показателя адиабаты 7 оказывается значительно слабее, что видно по небольшому различию между штрихованными и не штрихованными кривыми, которым соответствуют значения 7 = 4/3 и 7 = 5/3. Малое отличие между собой кривых с различными 7 обусловлено также относительно малым диапазоном изменения этого параметра.
Влияние коэффициента аккомодации е оказывается весьма слабым, так что на рис. 1 и на всех последующих рисунках кривые с различными е в диапазоне между нулем и единицей с графической точностью совпадают (различие не превосходит нескольких процентов).
Рис. 1. Зависимость нормированного потока массы д\/д\к ОТп числа Маха М при /3 = 1 (1,1'), /3 = 0.5 (2,2'), /3 = 0.1 (3,3') для трехатомного (со штрихом) и одноатомного (без штриха) пара.
Рис. 2. Зависимость нормированного потока импульса д-г!р& от числа Маха М при /3 -1 (1,1'), /3 = 0.5 (2,2'), /3 = 0.1 (3,3') для трехатомного (со штрихом) и одноатомного (без штриха) пара.
Поведение нормированного потока энергии дз/дж при тех же значениях параметров 7 и /3 описывается примерно такими же кривыми, как и на рис. 1. Отношение нормированных потоков дз/дзк и дг/дгк, равное (5 + 2)/(з + 1) при М —* 0, в соответствии с
условием (14), при увеличении числа Маха уменьшается приблизительно на 10 — '20%.
Зависимости от М для полного импульса отдачи д2 и давления пара на внешней стороне кнудсеновского слоя, нормированных на давление насыщения р3 = п3квГя, показаны на рис. 2 и 3 для того же набора параметров 7 = 5/3 (кривые без штриха) и 7 = 4/3 (кривые со штрихом), /3=1 (кривые 1,1'), 0.5 (кривые 2,2') и 0.1 (кривые 3,3'). Как и в случае с д\ и д3, на поведении этих кривых изменение /3 сказывается гораздо сильнее, чем изменение 7.
Рис. 3. Зависимость нормированного давления р/р, на внешней стороне кнудсеновского слоя от числа Маха М при /3=1 (1,1'), /3 = 0.5 (2,2'), /3 = 0.1 (3,3') для трехатомного (со штрихом) и одноатомного (без штриха) пара.
Рис. 4. Зависимость нормированной температуры Т/Тв на внешней стороне кнудсеновского слоя от числа Маха М при /3 = 1 (1,1'), /3 = 0.5 (2,2'), /3 = 0.1 (3,3') для трехатомного (со штрихом) и одноатомного (без штриха) пара.
На рис. 4 приведена зависимость от числа Маха для температуры пара на внешней стороне кнудсеновского слоя, отнесенная к величине 7!,, при тех же обозначениях и величинах параметров 7 и /3, что и на рис. 1-3. В отличие от рис. 1-3, в данном случае наиболее существенно на поведение кривых влияет величина 7. В связи с этим напомним, что изменение 7 сопровождается значительным изменением теплоемкости пара. Из рис. 4 видно, что для многоатомных газов скачок температуры на кнудсенов-ском слое Т — Т3 оказывается меньше, чем для одноатомного. Уменьшение /3 приводит к увеличению этого скачка, то есть к увеличению эффекта неравновесности фазового перехода. Аналогичным образом изменение /3 влияет и на поведение потока д2, который
уменьшается вместе с /3. Однако при уменьшении /3 происходит не увеличение, а уменьшение потоков д\ и <7з, что соответствует уменьшению эффекта неравновесности. В то же время, с этой точки зрения изменение 7 оказывает одинаковое влияние на поведег re всех потоков и температуры Т: при уменьшении 7 (переход от одноатомного пара к многоатомному) эффект неравновесности испарительного процесса уменьшается, что аналогично качественно влиянию от уменьшения числа Маха.
Рассмотрим теперь влияние коэффициента конденсации на устойчивость стационарного плоского фронта испарения относительно малых возмущений ~ exp(uit — ikx), с комплексной частотой ш = ги\ + u>2 и волновым числом к. Координата х направлена вдоль невозмущенного фронта испарения, который в системе координат, движущейся вместе с фронтом, совпадает с плоскостью z = 0. Дисперсионное уравнение для таких возмущений ш(к) исследовалось в работах [2] при /3 = 1 и 7 = 5/3. Конденсированная фаза при этом рассматривалась как невязкая несжимаемая жидкость с постоянными величинами теплоемкости с, температуропроводности х, плотности pi и коэффициента поглощения излучения a, a nap считался идеальным одноатомным газом. Дисперсионное уравнение D(ui, к) = 0 получалось как условие разрешимости линеаризованной однородной системы уравнений непрерывности, Эйлера и теплопроводности в конденсированной фазе и уравнений идеальной газодинамики в паре для возмущений стационарного режима испарения с испарительными граничными условиями на фронте перехода.
Обобщение этого анализа на случай /3 < 1 и 7 = 4/3 осуществляется с помощью изложенных выше испарительных граничных условий, которые заменяют используемте в [2] соответствующие выражения для одноатомного пара с коэффициентом конденсации /3 = 1. В качестве конкретного вещества будем рассматривать воду, для теплофизик ских параметров которой используются следующие значения:
pi = 1г/см3, х = 0.001 см2/с, с = 5.2 х 107эрг/г • К, L = 15 х 109эрг/г, m — 3 х 10~23г, <т = Пдин/см, а = 1л«кл<_1,
Ps = nskBTs = рь exp (т/(1 - Tb/Ts)), 77 = 13, Рь = 106дин/см2, (22)
Too = 300Л', Ть = 373 К,
Ш-> мкс-1
к мкм"1
W2max, МКС"1
2000
1000
1500
500
0
—г"---1— и -1 -.-1-1-
100 к, мкм"1 0 0.2 0.4 0.6 0.8 р
' 2500 - 2000
■ 1000
• 1500
■ 500 О
О
50
Рис. 5. Поведение инкремента неустойчивости фронта испарения воды в зависимости
от волнового числа при различных значениях /3 = 1 (1), 0.8 (2), 0.6 (3), 0.4 (4), 0.2 (5) с температурой поверхности Ts = 1.5Ть и числом Маха M = 0.815. Пунктирная кривая Г соответствует условиям 1 для одноатомного пара.
Рис. 6. Поведение величины максимума инкремента <*>2тах (пунктирная кривая) и его положения kmax (сплошная кривая) при изменении коэффициента конденсации (3.
где Ть - температура кипения воды при нормальном давлении рь, Г«, - температура в глубине конденсированной фазы, m - масса молекулы воды, а а - коэффициент поверхностного натяжения, соответствующий рассматриваемой температуре поверхности Ts = 1.5X1,. Величина а примерно соответствует коэффициенту поглощения лазерного излучения с длиной волны Л ~ 3мкм. Отметим, что при таких значениях Ts и а величина максимума температурного профиля в рамках данной модели с постоянными теплофизическими коэффициентами достигает критического значения Ткр = 647 К уже при /3 ~ 0.05.
На рис. 5 показано поведение вещественной части частоты lj2 (инкремента) от волнового числа к при M = 0.815 и различных значениях /9 = 1; 0.8; 0.6; 0.4; 0.2 (кривые 1-5). При /3 = 1 выбранному числу Маха соответствует дисперсионная кривая с максимальным (в зависимости от M) значением инкремента = 2567л<кс-1. При M — 1 величина u>2 = 2088jw гсс-1. Пунктирная кривая на рис. 5 соответствует тем же параметрам, что и кривая 1, но с измененным значением 7: вместо используемой для всех остальных кривых величины 7 = 4/3, у пунктирной кривой 7 = 5/3. Как следует из рис. 5, это изменение относительно слабо меняет дисперсионную кривую, причем это изменение
качественно противоположно тому, что можно было бы ожидать, исходя из описанного выше изменения выражений для потоков.
Уменьшение коэффициента конденсации при неизменной температуре поверхности и числе Маха приводит к уменьшению инкремента и области неустойчивости, то есть к смещению в область меньших значений к точки локализации максимума инкремента и правой границы области неустойчивости.
На рис. 6 представлены зависимости положения кщах (сплошная кривая и левая ось) и величины максимального (по к и М) инкремента <^2тах (пунктирная кривая и правая ось) от значения параметра /3. При этом для каждого значения /3 выбиралось наибольшее из максимальных значений инкремента на дисперсионной кривой, которому соответствует определенное значение М(/3). Как видно из рис. 6 зависимости ктах и ^2тах от /3 сравнительного слабо отличаются от линейной. Практически такие же результаты получаются и при фиксированном значении числа Маха М & 0.8, поскольку величина М(/9) в диапазоне изменения /3 от 0.1 до 1 меняется примерно на 5%.
При /3 = 1 величина инкремента и>2 ~ 2.5нс"1 получается весьма большой, то есть время развития неустойчивости, вообще говоря, оказывается меньше времени установления стационарного режима испарения. При этом длина волны наиболее неустойчивых возмущений может стать меньше толщины кнудсеновского слоя, который по этой причине уже нельзя считать поверхностью разрыва. Кроме того, при таких малых длинах волн роль вязкого затухания перестает быть пренебрежимо малой, поскольку мри коэффициенте вязкости и = 10~3сж2/с и к — 100ж?сл{-1 величина декремента 2ик~ достигает сравнимого с инкрементом ш2 значения 2нс~1. Однако при уменьшении /3 или температуры поверхности значения инкремента ^2тах и ^тах также уменьшаются, так что описанные выше ограничения становятся менее существенными.
ЛИТЕРАТУРА
[1] С а м о X и н A.A. Труды ИОФАН, 13, 3 (1988).
[2] К а р т а ш о в И. Н., Самохин А. А. Краткие сообщения по физике ФИАН, N 3, 33 (2002); N 10, 24 (2002); Квантовая электроника, 33, N 5, (2003).
[3] К н а к е О., Странский И. Н. УФН, 68, вып. 2, 261 (1959).
[4] X и р с Д., П а у н д Г. Испарения и конденсация. Изд. Металлургия, 1966.
[5] Ландау Л. Д., Лифшиц Е. М. Статистическая физика, часть 1. М., Наука, 1976.
[6] М а р т ы н е н к о О. Г., X р а м ц о в П. П. Инженерно-физический журнал, 74, N 3, 3 (2001).
[7]Nagayama G., Tsuruta Т. Journal of Chemical Physics, 118, N 3, 1392 (2003).
[8] А л e к с e e в В. H., E г e p e в С. В., Наугольных К. А. и др. Акустический журнал, 33, вып. 6, 961 (1987).
[9] Кумзерова Е. Ю., Шмидт А. А. ЖТФ, 72, вып. 7, 36 (2002).
[10] Кузнецова И. А., Юшканов А. А., Яламов Ю. И. ЖТФ, 67, N 10, 21 (1997).
[11] С г о u t D. J. Math. Phys., 15, 1 (1936).
[12] Knight С. J. AIAA J., 17, N 5, 81 (1979).
[13] Коган M. H., Макашев H. К. Изв. АН СССР, МЖГ, N 6, 3 (1971).
[14] Мажукин В. И., Прудковский П. А., Самохин А. А. Математическое моделирование, 5, N 6, 3 (1993).
[15] S о n е Y., Т a k a t a S., Golse F. Phys. Fluids, 13, N 1, 324 (2001).
[16] G u s а г о v A. V., S m u г о v I. Phys. Fluids, 14, N 12, 4242 (2002).
Институт общей физики Поступила в редакцию 12 мая 2003 г.
им. А. М. Прохорова РАН
