CC BY
23
УДК 536.4
О МОДЕЛИРОВАНИИ СПИНОДАЛЬНОГО РАСПАДА ПЕРЕГРЕТОЙ ЖИДКОСТИ
С. Н. Андреев, М. М. Демин, В. И. Мажукин, А. А. Самохин
Исследуется возможность применения разностных численных методов для решения задачи спинодалъного распада, в которой используется неустойчивая ветвь уравнения состояния (Ван-дер-Ваальса). В рамках данного подхода начальный этап фазового превращения оказывается близким к изотермическому. При этом максимальная величина взрывного скачка давления слабо зависит от параметров расчетной сетки, в отличие от общей пространственно-временной картины этого процесса.
При численном моделировании процесса интенсивного испарения лазерным облуче нием веществ с достаточно большим коэффициентом поглощения а > 10' м~1, напри мер, металлов, широко используется модель поверхностного испарения. В этой модели максимум температурного профиля Ттах совпадает с температурой поверхности Т3, с которой происходит испарение.
При меньших значениях а, характерных для диэлектриков, а также для металлов при их облучении электронным пучком либо рентгеновскими импульсами или в случае перехода металл-диэлектрик [1, 2], объемный нагрев вещества и поверхностное испарение приводят к формированию максимума температуры под поверхностью вещества. В результате этого в районе температурного максимума может происходить объемное взрывное вскипание перегретой жидкости. Исследование процесса объемного вскипания оказывается весьма сложной задачей, которая пока еще далека от своего окончательного решения, хотя этому вопросу посвящено множество экспериментальных и теоретических работ [3-8].
В недавней работе [9] рассматривалась простейшая модель, описывающая эволюцию паровой полости, возникающей при взрывном объемном вскипании. В этой модели
без учета гидродинамики рассчитывался температурный профиль в веществе мишени, нагреваемой импульсным лазерным излучением. Когда находящийся под поверхностью максимум температуры достигал температуры предельного перегрева (спинода-ли), примерно равной 0.9 Тс, где Тс - критическая температура [10], вводилась паровая полость, давление в которой принималось равным давлению насыщенного пара при тем пературе предельного перегрева. Далее рассматривалась эволюция паровой полости пр1 ■ ее расширении. В предложенной модели оставался открытым вопрос о начальном этапе формирования полости - о времени ее возникновения, ее характерных размерах, скорости нарастания давления пара в полости.
Возникновение паровой полости в подобных условиях определяется эволюцией лабильного состояния вещества, которое соответствует неустойчивой ветви непрерывного уравнения состояния. Анализ поведения вещества в лабильной области возможен с использованием методов молекулярной динамики или Монте-Карло [11]. Использование более простого континуального подхода с непрерывным уравнением состояния [12] имеет определенные преимущества, однако область его применимости остается не вполне ясной.
Целью настоящей работы является рассмотрение возможности моделирования начального этапа процесса взрывного вскипания в рамках непрерывного уравнения со стояния с использованием его неустойчивой ветви.
Постановка задачи. Рассматривается задача лазерного испарения жидкости (воды) в среду с противодавлением (воздух). Конденсированная и газовая среды разделены контактной границей. Воздух полностью прозрачен для лазерного излучения, которое поглощается слоем воды. Выделение энергии лазерного импульса носит объёмный характер. По мере нагрева воды на её поверхности начинает формироваться поток испаренного вещества, который оттесняет воздух и образует новую фазу - пар. Новая область, занимаемая паром, ограничена, с одной стороны, подвижной фазовой (испаряющаяся поверхность) границей, а с другой - подвижной контактной границей пар воздух. Для описания поведения каждой из трёх сред использовалась полная система уравнений гидродинамики, дополненная уравнением переноса лазерного излучения и соответствующими уравнениями состояния. На границах раздела сред производилась сшивка решений уравнения гидродинамики с помощью соответствующих граничных условий
Э£ дг
э(рц) _ о
дх
д(ри) , д(ри2) , д_Р _ а
дг ' дх Эх ~ и'
+ = - (р т дх \
ди | дх дх
ЭИ-^ . ЭЩ. , + Эх +
^ Т «ь(/>, Т)& = о, с = + С-^г = -А(Г)Й
9С\ Эх У
Д = 1,у,д.
(1)
[Р = Р(р,Т), е = £(/»,Г)
Здесь р,и,е,Т,Р - плотность, газодинамическая скорость, внутренняя энергия, температура и давление вещества соответственно, и (7 - коэффициент поглощения и плотность потока лазерного излучения, И^х - плотность потока тепла, А - коэффициент теплопроводности. Индексы 1,у,д обозначают принадлежность величин соответственно к жидкой, парообразной и воздушной средам. В конденсированной фазе величина е/ имеет смысл энтальпии жидкой фазы //;.
Граничные условия. В качестве граничных условий на левой (неподвижной) границе используется условие равенства нулю потока массы и тепла:
х = Г, : и = 0, XV = 0.
(2)
На испаряющейся поверхности Г{„ формулируются три закона сохранения в системе координат, движущейся со скоростью движения жидкой фазы = — и/, где скорость распространения фронта испарения в неподвижной (лабораторной) системе координат:
■Т 1 ¿тп
-31 + 3ь
х = Г/„ : - рту = рь(щ - + VI,,),
Л = р1 + ЗЪЧу = РУ+ - «„ + и/„),
(«г -иь - иг„)21
(3)
(4)
н л. '*
н1 + т
. -тп
- (т! = -]У + ]1у
Нь +
где = Ж/ = —А(Г/)^-, ^ = = —А(Т„)^. С учётом этих выражений закон
ЧуЧУ.Л.^.ДЛХХУ^ХХХХ./Х ^/ХХЧ^^УХ XX XX ХУХЧ-Г./ХЧ.ХХЧУ XX X СХХЭХХ X хэ х> о .
(5)
где Ь™ = 2£(Г|) + Сру(Ть - Г,„) + ^^ ("'7")2 - неравновесная теплота испарения, Ть -температура кипения, а - постоянная закона Стефана-Больцмана.
Процесс поверхностного испарения в приближении кнудсеновского слоя описывается тремя законами сохранения и тремя дополнительными параметрами на внешней стороне кнудсеновского слоя (температурой Г„, плотностью ру и скоростью и). Два из этих параметров (обычно Т„ и р„) в общем случае определяются из решения уравнения Больцмана, а третий (обычно число Маха М = и/ис) - из решения уравнений газовой динамики. Как правило, для определения Ту и ру кинетическое уравнение Больцмана напрямую не используется, а решение его находят с помощью некоторых аппроксима-ционных соотношений.
В данном случае Ту, рь определялись из соотношений на неравновесном кнудсе-новском слое по модели Крута [13]:
Т„ = Т/ат(М), ру = рз^ар(М), (6)
где М = с, = ат(М) = ар(М) = •
■> значение т определяется из уравнения ^(М)(т2+ 0.5)2 — т2(т2-\-1.5 +а) =
о, ™ р(м) = 1+а=- -1, (= ^+
- плотность насыщенного пара.
На контактной границе пар-воздух формулировались условия непрерывности плотности, давления и температуры:
X — Г у^д . Чу — Ру — Рду Ту — Тд. ( / )
На границе раздела конденсированная фаза-газ до начала испарения используются такие же условия, кроме условия для баланса энергии:
х = Гс>5, с = 5, / :
яс +
21
Нс +
М
21
В данном случае значением теплового потока в газе можно пренебречь:
ИЛГ = _стГ4 + (8)
Численное решение рассматриваемой задачи осуществлялось с помощью конечно-разностного метода динамической адаптации [14, 15], позволяющего рассматривать задачи с разрывными решениями и автоматическим выделением подвижных межфазных
границ. Применение метода динамической адаптации в данном случае позволяет полу-
чить ™ез^льтаты качественно отличающиеся от глез,гльтатов полеченных на сетках с
фиксированными узлами и без явного выделения подвижных границ.
Система уравнений в частных производных (1) с граничными условиями (3) - (5) аппроксимировалась семейством консервативных разностных схем [16]. Полученная си стема нелинейных разностных уравнений решалась с помощью специального вычисли тельного алгоритма, построенного по аналогии с [17]. Основу алгоритма составляют внешний и два вложенных итерационных цикла, с автоматическим выбором шага интегрирования по времени.
В качестве уравнения состояния в конденсированной среде использовалось уравнение Ван-дер-Ваальса [18]:
где критическая температура для воды Тсг = 647.15 К, а критическая плотность р^ задавалась из условия Р(Т = 273 К) = 1 бар, что дает рст « 0.390 г/см3. Для коэффициентов теплоемкости и теплопроводности, а также для зависимости плотности насыщенного пара от температуры использовались данные в виде таблиц, взятые из [19]. В расчетах использовался лазерный импульс с гауссовским распределением интенсив ности по времени С = (?0 ехр( — (¿/г)2), = 5 • 107 Вт/см2, т = 2 • Ю-7 с. При таком задании импульса моменту времени 2 = 0 соответствует максимум интенсивности, а начало расчета выбирается равным —4т. Длина волны лазера 1.06 мкм, коэффициен 1 поглощения лазерного излучения /с^ = 104 .и-1, отражения от поверхности нет.
Обсуждение рзультатов. На рис. 1 (а-г) показаны пространственные профили плотности (кривая 1), температуры (2) и давления (3) вблизи поверхности жидкости в различные моменты времени после попадания части перегретой жидкости в области неустойчивости.
Кривые на рис. 1а соответствуют моменту / = 25 нс, когда параметры перегретой жидкости в одном узле расчетной сетки (узел 285) достигают спинодали при Т = 551 А', Р = 26.7 бар на глубине Н = 0.5 мкм от свободной поверхности жидкости. В результате дальнейшего нагрева тонкий слой жидкости вблизи температурного максимума оказывается в лабильном состоянии, и его поведение определяется неустой-чивой ветвью уравнения состояния.
На графике (б) показаны те же профили, спустя 1 не (50 шагов по времени). К этому моменту времени в неустойчивой части фазовой диаграммы оказывается область
Р =
рВТ 9 ТсгЯр2
(9)
1 - п/(3Рсг) 8 рст
1000 800 600
н"
400 200 0
1=25 не -
1__-
2
.......3
20018.6 20018.8 20019.0 20019.2 20019.4
X, мкм
0.05^ 0.00
20018.6 20018.8 20019.0 20019.2 20019.4 X, МКМ
1000 800 о- 600
Н"
400 200 О
1=26,3 не
-1—лЬ—
____2......1ц.
----------
0.10^ 0.05 ей' 0.00
20018.6 20018.8 20019.0 20019.2 20019.4
х, мкм
0.10"
0.05
0.00
20018.6 20018.8 20019.0 20019.2 20019.4
X, мкм
Рис. 1. Профили плотности (кривая 1), температуры (2) и давления (3) в различные моменты времени при (а) I = 25 не, (6)2 = 26 нс, (в) ¿ = 26.3 не, (г) ¿ = 26.4 не.
жидкости толщиной /г = 0.2 мкм, содержащая 5 узлов расчетной сетки с номерами 282 - 286.
Распад лабильной области начинается в узле 285, в котором плотность жидкости имеет минимум р = 505 кг/м3, а давление возросло до величины Р = 68.5 бар. Эти перепады плотности и давления фактически реализуются между соседними узлами, и в газообразную фазу переходит только один 285-й узел. Соседние узлы (284 и 286) движутся в противоположном направлении, и плотность в этих точках возрастает (см. рис. 2, на котором приводится эволюция узлов 284, 285 и 286 на фазовой диаграмме (/э, Р)). Такое поведение р и Р можно рассматривать как начало формирования паровой полости.
График (1в) соответствует моменту времени £ = 26.3 не, когда узел 285 приближается к паровой ветви спинодали, оставаясь еще в лабильной области, а узлы 282 - 284 и 286 после пересечения жидкостной ветви спинодали оказываются в области
Рис. 2. Эволюция во времени трех узлов расчетной сетки на фазовой диаграмме. Пунктирной кривой показана спинодаль, ограничивающая область неустойчивых фазовых состоянии.
метастабильных состояний. Профиль скачка давления на рисунках (б, в) расширяется (по основанию) приблизительно со скоростью звука от значения Аг = 0.07летел« до Аг = 0.13 мкм за время Д£ - 0.3 не.
Давление в паровой полости в момент времени t = 26 н с достигает своего максималь ного значения Ртах = 254 бар. Ему в этом же узле соответствует локальный минимум плотности р = 283 кг/мг, тогда как температура в узлах 282-286 мало изменяется по сравнению со значением на рис. 1а и составляет Т = 549 К с точностью до ДТ = 5 К. Такая величина АТ согласуется с аналитическими оценками изменения температуры при возникновении паровой полости в перегретой жидкости вблизи спинодали [8] и позволяет считать, что процесс образования полости является практически изотермическим. Максимальная величина давления оказывается больше давления насыщенного пара Р3 — 221 бар, соответствующего этой изотерме. В то же время, как отмечалось ранее [8], при возникновении паровой полости давление в ней не должно превышать давление насыщенного пара, соответствующее температуре жидкости непосредственно перед образованием в ней паровой полости. Поэтому превышение давления АР — 33 бар является, по-видимому, свидетельством очевидной недостаточности модели с непрерывным уравнением состояния, в которой не учитывается возникновение границ фазового раздела. Необходимость учета межфазной кинетики на возникающих границах раздела
в подобных подходах, использующих непрерывные уравнения состояния, отмечалась и в пл.бпте [121.
1 I J
Скорее всего, именно этот недостаток данного подхода приводит к резкому уменьшению плотности, давления и температуры в паровой полости при ее дальнейшем расширении, как это видно из рис. Ы, соответствующего моменту времени I = 26.4 нс. К этому моменту времени все узлы расчетной сетки уже покинули область неустойчивости. Очевидно, что при учете испарения в полость давление и плотность пара в ней не могли бы быть существенно меньше их значений на линии насыщения при данной температуре стенок и коэффициенте конденсации /3 = 1. Поведение импульса давления на рис. 1г качественно согласуется с эволюцией начального импульса давления в неограниченной сплошной среде, когда изначально заданный импульс давления разделяется на две волны, движущиеся в противоположные стороны. Однако при этом профиль плотности качественно не изменяется, т.е сохраняет один локальный минимум.
Существенным недостатком данной модели является также зависимость пространственно-временных характеристик плотности, давления и температуры жидкости от параметров расчетной сетки. Например, изначальное увеличение (уменьшение) количества пространственных узлов в расчетной сетке в четыре раза по сравнению со схемой, используемой на рис. 1, чему соответствуют новые значения пространственного шага в области неустойчивости 0.005 мкм (0.08 мкм) мало меняет максимальную величину скачка давления Ртах — 254 бар, но существенно отражается на пространственно-временном поведении профиля давления. Время нарастания давления от величины Р = 100 бар до Ртах при этом уменьшается (увеличивается) от величины АЬ = 0.25 не до At — 0.14 не (0.6 не). Кроме того, уменьшается (увеличивается) начальный размер паровой полости от величины /г = 0.04 мкм до к = 0.01 мкм (/г = 0.16 мкм), который во всех этих случаях соответствует удвоенному расстоянию между соседними узлами расчетной сетки в области неустойчивости.
Если же увеличить количество пространственных узлов расчетной сетки, например, в восемь раз по сравнению со схемой, используемой на рис. 1, в момент I = 26 не (рис. 16), т.е. уже после того как часть узлов попала в неустойчивую область, то пространственная структура профилей давления и плотности существенно изменится по сравнению с рис. 1в и г, как показано на рис. 3. Как видно из рис. За, спустя 0.15 не по еле увеличения количества узлов, профили давления и плотности остаются подобными профилям на рис. 1в, однако затем внутри области неустойчивости начинается фор мирование субструктуры рис. 36, характерный масштаб которой определяется новым
600- 1=26.15 не 1
2 Л
а400"
н" 1 1 ч 1
200- 1 1 з % % % %
0-
0.25 0.20
0.15 ГО я
0.103 0.05 0.00
20018,90 20018,95 20019,00 х, мкм
а
Е-"
20018,90 20018,95 20019,00
X, мкм
Рис. 3. Профили плотности (1), температуры (2) и давления (3) в различные моменты времени при: (а) 2 = 26.15 не, (6)1 = 26.175 не.
пространственным шагом, аналогично тому как это было при формировании первоначальной полости на рис. 1в, г.
Таким образом, полученные в настоящей работе результаты численного моделирования с использованием неустойчивой ветви непрерывного уравнения состояния показывают, что на начальной стадии спинодального распада перегретой жидкости наблюдается резкое возрастание давления и падение плотности в той части жидкости, которая пересекла границу абсолютной термодинамической неустойчивости. Такое поведение можно рассматривать как начало процесса формирования паровой полости в перегретой жидкости. При этом максимальная величина взрывного скачка давления приближенно согласуется с оценками в работе [8] и слабо зависит от параметров расчетной сетки, в отличие от пространственно-временных характеристик профилей давления и плотности, которые по этой причине нельзя рассматривать как адекватную модель реальной физической картины процесса.
В заключение отметим, что в рассматриваемой здесь простейшей постановке задачи не учитывалась зависимость коэффициента поглощения излучения от плотности, а также возможность сингулярного поведения теплопроводности жидкости вблизи спи-нодали. Оценки показывают, что уменьшение коэффициента поглощения излучения в узкой области пониженной плотности не оказывает существенного влияния на температурный профиль в рассматриваемом временном интервале. Влияние сингулярного поведения теплопроводности на процесс спинодального распада перегретой жидкости будет рассмотрено отдельно.
Работа поддержана грантами РФФИ N 04-02-16452 и N 04-01-0701.
ЛИТЕРАТУРА
[1] Y о о J. Н., J е о n g S. Н., G г е i f R., and R u s s о R. E. J. Appl. Phys., 88. 1638 (2000).
[2] А н д p e e в С. H., Мажукин В. И., Никифорова Н. М., С а м о х и н А. А. Квантовая электроника, 33, N 9, 771 (2003).
[3] С а м о х и н А. А. Труды ИОФАН, 13, 1 (1990).
[4] М i о t е 1 1 о A. and Kelly R. Appl. Phys. A: Material Science and Proceeding, 69, 67 (1999).
[5] В u 1 g a k о v a N. M. and Bulgakov A. V. Appl. Phys. A: Material Science and Proceeding, 73, 199 (2001).
[6] X u X. and S о n g K. Appl. Phys., A69, 869 (1999).
[7] К i m D. and Grigoropoulos C. Appl. Phys. A: Material Science and Proceeding, 67, 169 (1998).
[8] Андреев С. H., О р л о в С. В., С а м о х и н А. А. Труды ИОФАН, 60, 127 (2004).
[9] А н д р е е в С. Н., С а м о х и н А. А. Краткие сообщения по физике ФИАН, N 8, 26 (2004).
[10] Скрипов В. П. Метастабильная жидкость. М., Наука, 1972.
[11] Байдаков В. Г., Проценко С. П. Теплофизика высоких температур. 41, 231 (2003).
[12] V i d а 1 F., J о h n s t о n Т. W., В a r t h e 1 e m у О., et al. Phys. Rev. Lett., 86, 2573 (2001).
[13] Crout D. J. Math. Phys., 15, 1 (1936); Мажукин В. И., Прудковский П. А, Самохин А. А. Математическое моделирование, 5, N 6, 3 (1993).
[14] Д а р ь и н Н. А., М а ж у к и н В. И. Дифференциальные уравнения, 23, N 7. 1154 (1987).
[15] Бреславский П. В., Мажукин В. И. Математическое моделирование, 7, N 12, 48 (1995).
[16] Самарский А. А. Теория разностных схем. М., Наука, 1989.
[17] Бреславский П. В., Мажукин В. И. Математическое моделирование, 3, N 10, 104 (1991).
[18] Ландау Л. Д., Л и ф ш и ц Е. М.
1 (\пп 1» I и.
[19] Р и в к и н С. Л., Александров водяного пара, М., Энергия, 1980.
Институт общей физики им. А. М. Прохорова РАН
Статистическая физика, ч. 1, М., Наука, А. А. Теплофизические свойства воды и
Поступила в редакцию 12 июля 2005 г.
