Прикладная математика и механика, 2023, T. 87, № 3, стр. 475-488

Исследование волны термической детонации в смеси капель воды с расплавленным свинцом

В. И. Мелихов 1*, О. И. Мелихов 1**, Салех Башар 1***

1 НИУ “МЭИ”
Москва, Россия

* E-mail: MelikhovVI@mpei.ru
** E-mail: MelikhovOI@mpei.ru
*** E-mail: basharsaleh10@gmail.com

Поступила в редакцию 25.03.2023
После доработки 12.04.2023
Принята к публикации 24.04.2023

Полный текст (PDF)

Аннотация

Исследованы закономерности волны термического взаимодействия капель воды, находящихся в высокотемпературном расплаве свинца, или волны термической детонации. Вследствие вскипания воды на поверхности свинца обе жидкости (фазы) разделены паровой пленкой. Используется одномерная модель взаимодействующих и взаимопроникающих континуумов, которая описывает динамику каждой жидкости введением специального поля, характеризующегося своими скоростью, температурой и объемной долей. Скорость волны определяется равенством скоростей и температур фаз в плоскости Чепмена–Жуге. Параметры на скачке давления вычисляются из условий на разрыве и являются граничными условиями для интегрирования уравнений сохранения в зоне взаимодействия капель воды с расплавом. Получающаяся структура волны термической детонации характеризуется тем, что максимальное давление находится на некотором удалении от ударной волны.

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

1. Введение. Исследования распространения детонационных волн в горючих газовых смесях проводились в течение длительного времени, что позволило установить основные закономерности этого явления [1, 2]: нагрев холодного слоя газа из-за сжатия на ударной волне воспламеняет эту порцию газа, а выделяющаяся при горении энергия поддерживает распространение волны. Таким образом, детонационная волна представляет собой лидирующую ударную волну, вслед за которой движется слой горящего газа. Течения газа в такой волне могут быть неустойчивыми, исследование этого вопроса выполнено в работе [3].

В случае двухфазных горючих систем, когда горючее представляет собой твердые частицы или капли, а окислителем является газ, эффекты, возникающие из-за относительного движения фаз, усложняют картину детонационной волны, что приводит, например, к различным формам распределений давления в детонационной волне горючих аэровзвесей, что было обнаружено в работах [46], результаты которых достаточно подробно изложены в монографии [7].

В работе [8] был предложен детонационный механизм взаимодействия смеси капель высокотемпературного расплава с водой, объясняющий природу крупномасштабных паровых взрывов [9, 10]. В отличие от детонации горючих смесей в данном случае распространение волны осуществляется за счет передачи тепла от расплава к воде и ее последующего вскипания и расширения пароводяной смеси. Лидирующая ударная волна вызывает дробление капель расплава на мелкие фрагменты. Это приводит к резкому увеличению площади поверхности теплопередачи от расплава к воде и соответствующему увеличению теплового потока в воду. Таким образом, процесс дробления капель расплава и передачи его тепловой энергии в воду аналогичен экзотермическим химическим реакциям в горючих смесях. Основные результаты исследований детонации данного типа (термической детонации) опубликованы в статьях [1113].

Новый вид термической детонации расплава с водой был обнаружен во время работ над созданием энергетических реакторов со свинцовым теплоносителем [14]. Их конструкция предусматривает так называемые проектные аварии, одной из которых является разрыв теплообменной трубки парогенератора, когда струя воды высокого давления из разрыва истекает в объем, занимаемый расплавленным свинцом, находящимся под низким давлением. В этом случае в области, окружающей разрыв, может образоваться дисперсная смесь капель с расплавленным свинцом. В такой системе лидирующая ударная волна будет дробить капли воды, увеличивая площадь поверхности теплопередачи от расплава к воде и подпитывая тепловой энергией расплава распространение волны термической детонации. Такая физическая система была исследована [15] методом адиабат Гюгонио, что позволило определить скорость волны и параметры многофазной смеси в плоскости Чепмена–Жуге.

Настоящая работа является продолжением исследования [15]. На основе одномерных стационарных уравнений динамики многофазных систем [7] построена модель волны термической детонации в системе “капли воды–водяной пар–жидкий свинец”, которая позволила рассчитать структуру волны термической детонации и получить ее основные характеристики.

2. Постановка задачи. Будем рассматривать бесконечное пространство, однородно заполненное многофазной смесью, по которой распространяется стационарная плоская волна термической детонации. Выберем систему координат, связанную с распространяющейся волной. Пусть декартова ось x перпендикулярна плоскости волны и направлена в сторону движения продуктов детонации. Плоскость $x = 0$ совпадает с плоскостью лидирующей ударной волны. Исходная многофазная смесь поступает со скоростью распространения волны детонации на плоскость ударной волны из области $x < 0$. На самой плоскости $x = 0$ происходит скачок давления и некоторых других параметров смеси, после чего начинается фрагментация капель воды и взаимодействие отдельных составляющих (фаз) смеси друг с другом, в результате чего смесь продуктов детонации асимптотически достигает скоростного и термического равновесия.

Значения равновесных параметров смеси продуктов детонации лежат на адиабате Гюгонио и определяют скорость распространения волны термической детонации [9]. Значения параметров на лидирующей ударной волне вычисляются из балансовых соотношений на разрыве, дополненных рядом модельных допущений. Структура волны в зоне фрагментации определяется численным решением стационарных одномерных уравнений динамики многофазных систем.

Будем полагать, что исходная многофазная система состоит из трех фаз. Первая (несущая) фаза представляет собой сплошной расплавленный свинец, который будем идентифицировать индексом $m$ (melt). Вторую (дисперсную) фазу составляют исходные (крупные) капли воды, окруженные паровой пленкой. Паровая пленка возникает вследствие вскипания воды на поверхностях водяных капель, граничащих с высокотемпературным свинцом. $d$ (droplet) – индекс второй фазы. Третьей (тоже дисперсная) фазой являются мелкие капельки (фрагменты) воды, образующиеся при дроблении исходных капель. Эти капельки тоже окружены паровыми пленками, отделяющими их от высокотемпературного расплава. Индекс третьей фазы – $f$ (fragment). Предполагается, что вода и пар, из которых состоят частицы второй и третьей фазы, находятся в состоянии теплового и скоростного равновесия. Каждая фаза характеризуется распределением в пространстве своей объемной доли $\alpha $, и каждая фаза описывается своей плотностью, скоростью и температурой. Давления всех фаз одинаковы.

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

Уравнения неразрывности фаз:

(2.1)
$\frac{d}{{dx}}\left( {{{\alpha }_{m}}{{\rho }_{m}}{{V}_{m}}} \right) = 0$
(2.2)
$\frac{d}{{dx}}({{\alpha }_{d}}{{\rho }_{d}}{{V}_{d}}) = - {{\Gamma }_{f}}$
(2.3)
$\frac{d}{{dx}}({{\alpha }_{f}}{{\rho }_{f}}{{V}_{f}}) = {{\Gamma }_{f}},$
здесь используются следующие обозначения: $\rho $ и $V$ – плотность и скорость фазы, соответственно. Величина ${{\Gamma }_{f}}$ выражает массовую скорость фрагментации исходных капель воды.

Уравнения импульсов фаз:

(2.4)
$\frac{d}{{dx}}\left( {{{\alpha }_{m}}{{\rho }_{m}}V_{m}^{2}} \right) = - {{\alpha }_{m}}\frac{{dP}}{{dx}} + {{K}_{{dm}}}({{V}_{d}} - {{V}_{m}}) + {{K}_{{fm}}}\left( {{{V}_{f}} - {{V}_{m}}} \right)$
(2.5)
$\frac{d}{{dx}}\left( {{{\alpha }_{d}}{{\rho }_{d}}V_{d}^{2}} \right) = - {{\alpha }_{d}}\frac{{dP}}{{dx}} + {{K}_{{dm}}}({{V}_{m}} - {{V}_{d}}) - {{\Gamma }_{f}}{{V}_{d}}$
(2.6)
$\frac{d}{{dx}}\left( {{{\alpha }_{f}}{{\rho }_{f}}V_{f}^{2}} \right) = - {{\alpha }_{f}}\frac{{dP}}{{dx}} + {{K}_{{fm}}}({{V}_{m}} - {{V}_{f}}) + {{\Gamma }_{f}}{{V}_{d}},$
здесь $p$ обозначает давление. Величины ${{K}_{{dm}}}$ и ${{K}_{{fm}}}$ описывают межфазное трение дисперсных фаз с несущей фазой.

Уравнения энергий фаз:

(2.7)
$\begin{gathered} \frac{d}{{dx}}({{\alpha }_{m}}{{\rho }_{m}}{{V}_{m}}{{h}_{{sm}}}) = {{R}_{{dm}}}({{T}_{d}} - {{T}_{m}}) + {{R}_{{fm}}}({{T}_{f}} - {{T}_{m}}) + \\ + \;{{V}_{m}}{{K}_{{dm}}}({{V}_{d}} - {{V}_{m}}) + {{V}_{m}}{{K}_{{fm}}}({{V}_{f}} - {{V}_{m}}) \\ \end{gathered} $
(2.8)
$\frac{d}{{dx}}({{\alpha }_{d}}{{\rho }_{d}}{{V}_{d}}{{h}_{{sd}}}) = {{R}_{{dm}}}({{T}_{m}} - {{T}_{d}}) + {{V}_{d}}{{K}_{{dm}}}({{V}_{m}} - {{V}_{d}}) + {{K}_{{dm}}}{{({{V}_{m}} - {{V}_{d}})}^{2}} - {{\Gamma }_{f}}{{h}_{{sd}}}$
(2.9)
$\frac{d}{{dx}}({{\alpha }_{f}}{{\rho }_{f}}{{V}_{f}}{{h}_{{sf}}}) = {{R}_{{fm}}}({{T}_{m}} - {{T}_{f}}) + {{V}_{f}}{{K}_{{fm}}}({{V}_{m}} - {{V}_{f}}) + {{K}_{{fm}}}{{({{V}_{m}} - {{V}_{f}})}^{2}} + {{\Gamma }_{f}}{{h}_{{sd}}}$

Здесь ${{h}_{s}}$ и $T$ – энтальпия торможения и температура фаз. Величины ${{R}_{{dm}}}$ и ${{R}_{{fm}}}$ описывают межфазный теплообмен дисперсных фаз с несущей фазой.

Все три фазы в совокупности полностью занимают произвольный локальный объем рассматриваемой области, поэтому имеет место следующее уравнение:

(2.10)
${{\alpha }_{m}} + {{\alpha }_{d}} + {{\alpha }_{f}} = 1$

Энтальпия торможения определяется следующим образом:

(2.11)
${{h}_{s}} = e + \frac{P}{\rho } + \frac{{{{V}^{2}}}}{2} = 1,$
где $e$ – внутренняя энергия.

Плотность и внутренняя энергия жидкого свинца являются функциями его температуры, которые определялись по данным [16]. Плотности $d$ – фазы и $f$ – фазы определяются объемными долями пара в каждой фазе и давлением:

(2.12)
${{\rho }_{d}} = (1 - {{\varphi }_{d}})\rho {\kern 1pt} '(P) + {{\varphi }_{d}}\rho {\kern 1pt} ''(P)$
(2.13)
${{\rho }_{f}} = (1 - {{\varphi }_{f}})\rho {\kern 1pt} '(P) + {{\varphi }_{f}}\rho {\kern 1pt} ''(P),$
где ${{\varphi }_{d}}$ и ${{\varphi }_{f}}$ – объемные доли пара в $d$-фазе и $f$-фазе соответственно, $\rho \left( P \right)$ и $\rho {\kern 1pt} '\left( P \right)$ – плотности воды и пара на линии насыщения при давлении, которые вычисляются по [17].

Внутренние энергии $d$-фазы и $f$-фазы определяются массовыми долями пара в каждой фазе и давлением:

(2.14)
${{e}_{d}} = (1 - {{\chi }_{d}})e{\kern 1pt} '(P) + {{\chi }_{d}}e{\kern 1pt} ''(P)$
(2.15)
${{e}_{f}} = (1 - {{\chi }_{f}})e{\kern 1pt} '(P) + {{\chi }_{f}}e{\kern 1pt} ''(P),$
где ${{\chi }_{d}}$ и ${{\chi }_{f}}$ – массовые доли пара в $d$-фазе и $f$-фазе соответственно, $e{\kern 1pt} '(P)$ и $e{\kern 1pt} '{\kern 1pt} '(P)$ – внутренние энергии воды и пара на линии насыщения при давлении $P$, вычисляемые по [17].

Если в результате фазовых переходов $d$-фаза и/или $f$-фаза станут чисто однофазными (пар или вода), то уравнения (2.12)(2.15) заменяются зависимостями от давления и температуры каждой фазы, которые вычисляются по [17].

3. Описание межфазного взаимодействия. Рассмотрим силу трения между крупными каплями воды ($d$-фаза) и несущей фазой жидкого свинца. Прежде всего, следует принять во внимание, что в несущей фазе еще присутствуют мелкие фрагменты воды ($f$-фаза), которые изменяют плотность среды, окружающей крупную каплю воды, в результате чего эффективная плотность этой среды ${{\rho }_{{{\text{ef}}}}}$ становится равной

(3.1)
${{\rho }_{{{\text{ef}}}}} = \frac{{{{\alpha }_{m}}{{\rho }_{m}} + {{\alpha }_{f}}{{\rho }_{f}}}}{{{{\alpha }_{m}} + {{\alpha }_{f}}}}$

С учетом (3.1) сила трения между крупными каплями воды и несущей фазой записывается в виде [7]:

(3.2)
${{F}_{d}} = \frac{1}{2}{{C}_{{D,dm}}}{{\rho }_{{{\text{ef}}}}}\pi \frac{{L_{d}^{2}}}{4}\left| {{{V}_{m}} - {{V}_{d}}} \right|({{V}_{m}} - {{V}_{d}})$

Здесь ${{C}_{{D,dm}}}$ – коэффициент сопротивления капли, ${{L}_{d}}$ – ее диаметр.

Известно [3], что при фрагментации капли происходит ее сильная деформация, она принимает форму диска, расположенного перпендикулярно набегающему потоку. Это приводит к увеличению коэффициента сопротивления по сравнению с типичной величиной 0.4 для турбулентного обтекания сферы до величины ${{C}_{{D,dm}}} = 2.5$ [18], которое и используется в настоящей работе.

Если концентрация крупных капель равна ${{n}_{d}}$, то сила трения между этими каплями и несущей средой равна ${{F}_{{d,\Sigma }}} = {{n}_{d}}{{F}_{d}}$. Поскольку объемная доля $d$-фазы связана с концентрацией капель соотношением

(3.3)
${{\alpha }_{d}} = \pi \frac{{L_{d}^{3}}}{6}{{n}_{d}},$
то силу трения ${{F}_{{d,\Sigma }}}$ можно записать следующим образом:

(3.4)
${{F}_{{d,\Sigma }}} = \frac{{6{{\alpha }_{d}}}}{{\pi L_{d}^{3}}}{{F}_{d}}$

Подставляя (3.2) в (3.4), получим:

(3.5)
${{F}_{{d,\Sigma }}} = \frac{{3{{\alpha }_{d}}}}{{4{{L}_{d}}}}{{C}_{{D,dm}}}{{\rho }_{{{\text{ef}}}}}\left| {{{V}_{m}} - {{V}_{d}}} \right|({{V}_{m}} - {{V}_{d}})$

Сила ${{F}_{{d,\Sigma }}}$ в уравнении импульса $d$-фазы была обозначена членом ${{K}_{{dm}}}({{V}_{m}} - {{V}_{d}})$. Тогда с учетом (3.5) величина ${{K}_{{dm}}}$ имеет вид:

(3.6)
${{K}_{{dm}}} = \frac{{3{{C}_{{D,dm}}}}}{{4{{L}_{d}}}}{{\rho }_{{{\text{ef}}}}}{{\alpha }_{d}}\left| {{{V}_{m}} - {{V}_{d}}} \right|$

Аналогичным образом определяется величина ${{K}_{{fm}}}$, описывающая силу трения между $f$-фазой и несущей средой

(3.7)
${{K}_{{fm}}} = \frac{{3{{C}_{{D,fm}}}}}{{4{{L}_{f}}}}{{\rho }_{m}}{{\alpha }_{f}}\left| {{{V}_{m}} - {{V}_{f}}} \right|,$
где ${{L}_{f}}$ – диаметр образующегося фрагмента воды, ${{C}_{{D,fm}}}$ – коэффициент сопротивления фрагмента. Поскольку образующиеся фрагменты воды имеют маленький размер и далее не дробятся, то для их коэффициента сопротивления принята стандартная величина ${{C}_{{D,fm}}} = 0.4$ [7].

Теплообмен между расплавленным свинцом и водой описывается величинами ${{R}_{{dm}}}$ и ${{R}_{{fm}}}$, входящими в уравнения энергий фаз (2.7)–(2.9). Действуя таким же образом, как с величинами ${{K}_{{dm}}}$ и ${{K}_{{fm}}}$ в случае описания межфазного трения, можно получить связь величин ${{R}_{{dm}}}$ и ${{R}_{{fm}}}$ с коэффициентами теплоотдачи дисперсных фаз:

(3.8)
${{R}_{{dm}}} = 6{{\alpha }_{d}}\frac{{{{h}_{{dm}}}}}{{{{L}_{d}}}}$
(3.9)
${{R}_{{fm}}} = 6{{\alpha }_{f}}\frac{{{{h}_{{fm}}}}}{{{{L}_{f}}}},$
где ${{R}_{{dm}}}$ и ${{R}_{{fm}}}$ – коэффициенты теплоотдачи от капель $d$-фазы и от капель $f$-фазы соответственно.

Теплообмен между расплавом и водой представляет собой очень сложное и не до конца изученное явление [9, 10]. При высоких температурах расплава, что, как правило, всегда имеет место в приложениях, на границе раздела расплава и воды образуется паровая пленка вследствие вскипания воды. Теплопередача от расплава к воде осуществляется через эту паровую пленку двумя механизмами: излучением и конвекцией пара в пленке [19]. Проведенные исследования позволяют оценить коэффициент теплоотдачи в этих условиях в диапазоне 100–1000 Вт/(м2·К) [2022]. Следует отметить, что нами рассматриваются процессы непосредственно в зоне фрагментации сразу же за ударной волной. Как правило, на ударной волне, когда давление резко возрастает относительно уровня в исходной многофазной смеси, равновесное двухфазное пароводяное образование (капля воды, окруженная паровой пленкой) трансформируется в каплю воды, температура которой меньше температуры кипения, соответствующей температуре кипения при давлении на ударной волне. В этих условиях паровая пленка будет коллапсировать. Но поскольку прямой контакт воды с высокотемпературным расплавом невозможен, то водяную каплю будет окружать очень тонкий слой пара, что означает интенсификацию теплообмена между расплавом и водой. Поэтому далее в расчетах было использовано верхнее значение диапазона коэффициента теплоотдачи, ${{h}_{{dm}}} = {{h}_{{fm}}}$ = 1000 Вт/(м2·К).

Оказавшаяся сразу за ударной волной крупная капля воды будет тормозиться относительно окружающего ее расплавленного свинца, поскольку плотность воды на порядок меньше плотности свинца. Жидкий свинец, как более тяжелая жидкость будет меньше тормозиться на фронте ударной волны по сравнению с водяными каплями. Возникающая разность скоростей между каплей и окружающей несущей фазой будет приводить к возникновению возмущений на поверхности капли и срыву отдельных фрагментов воды. Для типичных параметров исследуемой системы фрагментация капли будет осуществляться по механизму обдирки поверхностного слоя [7], когда с поверхности капли срываются капельки размером, который существенно меньше размера “материнской” капли. Для описания фрагментации капли используется модель, предложенная в [23], согласно которой скорость уноса массы воды с поверхности отдельной капли определяется следующей формулой:

(3.10)
$\frac{{{{d}_{m}}}}{{dt}} = {{c}_{{{\text{frag}}}}}\left| {{{V}_{d}} - {{V}_{m}}} \right|\pi L_{d}^{2}\sqrt {{{\rho }_{d}}{{\rho }_{{{\text{ef}}}}}} ,$
где ${{c}_{{{\text{frag}}}}}$ – эмпирическая константа.

Умножая (3.10) на концентрацию капель воды ${{n}_{d}}$, получим выражение для массовой скорости фрагментации капель воды в единице объема:

(3.11)
${{\Gamma }_{f}} = {{n}_{d}}\frac{{{{d}_{m}}}}{{dt}} = \frac{{{{\alpha }_{d}}\left| {{{V}_{d}} - {{V}_{m}}} \right|\sqrt {{{\rho }_{d}}{{\rho }_{{{\text{ef}}}}}} }}{{{{L}_{d}}}}$

При выводе (3.11) было использовано значение ${{c}_{{{\text{frag}}}}} = 1{\text{/}}6$.

В ходе дробления исходной крупной капли воды ее размер ${{L}_{d}}$ уменьшается. Когда он достигнет значения размера фрагментов $({{L}_{d}} = {{L}_{f}})$, то полагается, что капля больше не дробится, а переходит в $f$-фазу.

Распределение концентрации дробящихся капель воды вплоть до момента достижения размера фрагментов описывается уравнением [24]:

(3.12)
$\frac{{d({{n}_{d}}{{V}_{d}})}}{{dx}} = 0$

Комбинируя уравнение (3.12) с уравнением неразрывности $d$-фазы (2.2), можно получить уравнение, описывающее изменение размера дробящихся капель воды в зоне фрагментации:

(3.13)
${{\alpha }_{d}}\frac{{d\left( {{{\rho }_{d}}L_{d}^{3}} \right)}}{{dx}} = - {{\Gamma }_{f}}L_{d}^{3}$

Полагалось, что фрагментация имеет место, когда число Вебера капли воды превосходит критическое значение:

(3.14)
$\operatorname{We} = \frac{{{{\rho }_{{{\text{ef}}}}}{{{\left| {{{V}_{d}} - {{V}_{m}}} \right|}}^{2}}{{L}_{d}}}}{\sigma } > {{\operatorname{We} }_{{\operatorname{cr} }}} = 12,$
где $\sigma $ – коэффициент поверхностного натяжения.

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

Складывая уравнения (2.1)(2.3) и интегрируя суммарное уравнение, можно получить уравнение сохранения потока массы многофазной смеси:

(4.1)
${{\alpha }_{m}}{{\rho }_{m}}{{V}_{m}} + {{\alpha }_{d}}{{\rho }_{d}}{{V}_{d}} + {{\alpha }_{f}}{{\rho }_{f}}{{V}_{f}} = {{\operatorname{const} }_{1}}$

Аналогично из уравнений импульсов каждой фазы (2.4)–(2.6) можно получить уравнение сохранения потока импульса многофазной смеси:

(4.2)
${{\alpha }_{m}}{{\rho }_{m}}V_{m}^{2} + {{\alpha }_{d}}{{\rho }_{d}}V_{d}^{2} + {{\alpha }_{f}}{{\rho }_{f}}V_{f}^{2} + P = {{\operatorname{const} }_{2}}$

Таким же образом из уравнений энергий фаз (2.7)–(2.9) получается уравнение сохранения потока энергии многофазной смеси:

(4.3)
${{\alpha }_{m}}{{\rho }_{m}}{{V}_{m}}{{h}_{{sm}}} + {{\alpha }_{d}}{{\rho }_{d}}{{V}_{d}}{{h}_{{sd}}} + {{\alpha }_{f}}{{\rho }_{f}}{{V}_{f}}{{h}_{{sf}}} = {{\operatorname{const} }_{3}}$

Обозначим параметры исходной многофазной смеси нижним индексом “0”, а нижним индексом “2” будем обозначать параметры равновесной многофазной смеси после завершения межфазного взаимодействия. Тогда с помощью первых интегралов (4.1)–(4.3) можно получить уравнения, связывающие параметры многофазного потока в сечениях “0” и “2”:

(4.4)
${{\rho }_{2}}{{V}_{2}} + {{\rho }_{0}}{{V}_{0}}$
(4.5)
${{P}_{2}} + {{\rho }_{2}}V_{2}^{2} = {{P}_{0}} + {{\rho }_{0}}V_{0}^{2}$
(4.6)
${{h}_{2}} + \frac{1}{2}V_{2}^{2} = {{h}_{0}} + \frac{1}{2}V_{0}^{2}$

В уравнениях (4.4)(4.6) были использованы следующие обозначения:

(4.7)
${{\rho }_{0}} = {{\alpha }_{{m0}}}{{\rho }_{{m0}}} + {{\alpha }_{{d0}}}{{\rho }_{{d0}}}$
(4.8)
${{\rho }_{2}} = {{\alpha }_{{m2}}}{{\rho }_{{m2}}} + {{\alpha }_{{f2}}}{{\rho }_{{f2}}}$
(4.9)
${{h}_{0}} = \frac{{{{\alpha }_{{m0}}}{{\rho }_{{m0}}}{{h}_{{m0}}} + {{\alpha }_{{d0}}}{{\rho }_{{d0}}}{{h}_{{d0}}}}}{{{{\alpha }_{{m0}}}{{\rho }_{{m0}}} + {{\alpha }_{{d0}}}{{\rho }_{{d0}}}}}$
(4.10)
${{h}_{2}} = \frac{{{{\alpha }_{{m2}}}{{\rho }_{{m2}}}{{h}_{{m2}}} + {{\alpha }_{{f2}}}{{\rho }_{{f2}}}{{h}_{{f2}}}}}{{{{\alpha }_{{m2}}}{{\rho }_{{m2}}} + {{\alpha }_{{f2}}}{{\rho }_{{f2}}}}}$

После алгебраических преобразований уравнений (4.4)–(4.6) получается уравнение адиабаты Гюгонио:

(4.11)
$({{P}_{2}} - {{P}_{0}})({{\nu }_{2}} + {{\nu }_{0}}) - 2({{h}_{2}} - {{h}_{0}}) = 0,$
где ${{\nu }_{0}} = 1{\text{/}}{{\rho }_{0}}$ и ${{\nu }_{2}} = 1{\text{/}}{{\rho }_{2}}$ – удельные объемы многофазной смеси в сечениях потока “0” и “2”.

В уравнении адиабаты Гюгонио (4.11) удельный объем ${{\nu }_{2}}$ и энтальпия равновесной многофазной смеси в сечении 2 определяются давлением ${{P}_{2}}$ и температурой ${{T}_{2}}$ в этом сечении. Задавая температуру ${{T}_{2}}$ и определяя из уравнения (4.11) давление ${{P}_{2}}$, можно вычислить значение удельного объема ${{\nu }_{2}}$ и на плоскости $({{\nu }_{2}},{{P}_{2}})$ построить кривую (адиабату Гюгонио), которая будет описывать все возможные состояния равновесной смеси для данных параметров исходной многофазной смеси.

Определение значений параметров ${{\nu }_{2}}$ и ${{P}_{2}}$ на адиабате Гюгонио позволяет вычислить из уравнений (4.4) и (4.5) скорость волны термической детонации, которая равна скорости исходной смеси ${{V}_{0}}$,

(4.12)
${{V}_{0}} = {{\nu }_{0}}\sqrt {\frac{{{{P}_{2}} - {{P}_{0}}}}{{{{\nu }_{0}} - {{\nu }_{2}}}}} $

5. Расчет параметров на лидирующей ударной волне. Теперь определим значения параметров многофазной смеси на лидирующей ударной волне, которая распространяется со скоростью ${{V}_{0}}$. Обозначим нижним индексом “1” сечение потока, в котором находится плоскость ударной волны. В этом месте происходит скачок давления, который приводит к разрыву некоторых параметров многофазного потока относительно их исходных значений.

Будем полагать, что температура расплавленного свинца в сечении “1” не меняется $({{T}_{{m1}}} = {{T}_{{m0}}})$, а фрагментация капель воды начинается после их прохождения сечения “1”.

Потоки масс расплава и $d$-фазы на ударной волне сохраняются:

(5.1)
${{\alpha }_{{m1}}}{{\rho }_{{m1}}}{{V}_{{m1}}} = {{\alpha }_{{m0}}}{{\rho }_{{m0}}}{{V}_{0}}$
(5.2)
${{\alpha }_{{d1}}}{{\rho }_{{d1}}}{{V}_{{d1}}} = {{\alpha }_{{d0}}}{{\rho }_{{d0}}}{{V}_{0}}$

Так как температура расплава на ударной волне не меняется, то ${{\rho }_{{m1}}} = {{\rho }_{{m0}}}$. (Отметим, что скорости расплава и $d$-фазы на ударной волне не равны между собой.)

Поток импульса многофазной смеси на ударной волне также сохраняется.

(5.3)
${{P}_{1}} + {{\alpha }_{{m1}}}{{\rho }_{{m1}}}V_{{m1}}^{2} + {{\alpha }_{{d1}}}{{\rho }_{{d1}}}V_{{d1}}^{2} = {{P}_{0}} + ({{\alpha }_{{m0}}}{{\rho }_{{m0}}} + {{\alpha }_{{d0}}}{{\rho }_{{d0}}})V_{0}^{2}$

Кроме того, на ударной волне выполняется закон сохранения энергии для каждой фазы [21]:

(5.4)
${{h}_{{m1}}} + \frac{1}{2}V_{{m1}}^{2} = {{h}_{{m0}}} + \frac{1}{2}V_{0}^{2}$
(5.5)
${{h}_{{d1}}} + \frac{1}{2}V_{{d1}}^{2} = {{h}_{{d0}}} + \frac{1}{2}V_{0}^{2}$

Многофазная смесь на ударной волне состоит из расплава и дисперсной $d$-фазы, поэтому имеет место соотношение ${{\alpha }_{{m1}}} + {{\alpha }_{{d1}}} = 1$.

Плотность и энтальпия $d$-фазы на ударной волне определяются давлением и объемной долей пара в капле (${{\rho }_{{d1}}} = {{\rho }_{{d1}}}({{P}_{1}},{{\varphi }_{1}})$ и ${{h}_{{d1}}} = {{h}_{{d1}}}({{P}_{1}},{{\varphi }_{1}})$) или, если паровая пленка сконденсировалась, давлением и температурой воды в капле (${{\rho }_{{d1}}} = {{\rho }_{{d1}}}({{P}_{1}},{{T}_{{d1}}})$ и ${{h}_{{d1}}} = {{h}_{{d1}}}({{P}_{1}},{{T}_{{d1}}})$). Энтальпия расплава на ударной волне определяется давлением на ударной волне и температурой исходного расплава, поскольку его температура на ударной волне не меняется, ${{h}_{{m1}}} = {{h}_{{m1}}}({{P}_{1}},{{T}_{{m0}}})$.

Уравнения (5.1)(5.5) позволяют определить неизвестные значения параметров на ударной волне (в сечении “1”): ${{\alpha }_{{m1}}}$, ${{\alpha }_{{d1}}}$, ${{V}_{{m1}}}$, ${{V}_{{d1}}}$, ${{P}_{1}}$, ${{\varphi }_{1}}$ или ${{T}_{{d1}}}$. Эти значения используются в качестве краевых условий для решения системы дифференциальных уравнений (2.1)–(2.9), описывающих структуру волны термической детонации.

6. Метод решения. Описанная выше задача решалась численно в следующей последовательности: 1) построение адиабаты Гюгонио и определение скорости волны термической детонации, 2) расчет параметров на лидирующей ударной волне, 3) решение системы обыкновенных дифференциальных уравнений, описывающих процессы в зоне фрагментации после ударной волны.

Построение адиабаты Гюгонио выполнялось с использованием метода бисекций. Расчеты параметров на ударной волне производились методом итераций Ньютона. Дифференциальные уравнения решались со вторым порядком точности по шагу интегрирования, при этом члены в правых частях вычислялись в центральной точке, а их нелинейность снималась итерациями.

7. Результаты расчета стационарной волны термической детонации. С помощью разработанной модели был выполнен расчет волны термической детонации в смеси капель воды с расплавленным свинцом. Использовались следующие параметры исходной системы: давление ${{P}_{0}} = 0.8$ МПа, температура расплавленного свинца ${{T}_{{m0}}} = 800$ К, объемная доля пара в каплях ${{\varphi }_{0}} = 0.7$, диаметр капель ${{L}_{{d0}}} = 10$ мм. Полагалось, что вода и пар в исходных каплях находятся в состоянии насыщения, соответствующего системному давлению 0.8 МПа. Размер образующихся при дроблении капель фрагментов воды принимался равным ${{L}_{f}} = 1$ мм.

По этим параметрам была построена адиабата Гюгонио и касательная к ней, проведенная из точки начального состояния многофазной смеси, то есть прямая Рэлея–Михельсона [1, 2]. Точка касания этой прямой адиабаты Гюгонио (точка Чепмена–Жуге) определяет параметры равновесной смеси продуктов детонации, а наклон прямой определяет скорость волны детонации, которая будет реализована [1, 2, 9]. Для данных системных параметров скорость детонации равна ${{V}_{{\det }}} = {{V}_{0}}$ = 222 м/с. Давление в точке (плоскости) Чепмена–Жуге составляет 46.61 МПа.

Рассчитанные параметры смеси на лидирующей ударной волне имеют следующие значения: объемная доля капель воды ${{\alpha }_{{d1}}} = 0.29$, объемная доля расплавленного свинца ${{\alpha }_{{m1}}} = 0.71$, скорость капель воды ${{V}_{{d1}}} = 70$ м/с, скорость расплава свинца ${{V}_{{m1}}} = 218.66$ м/с, давление ${{P}_{1}} = 8.52$ МПа. Также расчет параметров смеси на ударной волне показал, что на ударной волне пар в паровой пленке, окружающей каплю воды, сконденсировался, и капля стала состоять из одной воды, температура которой ${{T}_{{d1}}} = 452.7$ К меньше температуры насыщения (кипения) воды, которая соответствует давлению на ударной волне ${{P}_{1}}$. Эти значения использовались, как граничные условия для решения дифференциальных уравнений (2.1)–(2.9), описывающих процессы в зоне фрагментации капель за ударной волной.

На рис. 1 представлено распределение давления за ударной волной в зоне фрагментации. Видно, что давление исходной смеси после ударной волны возрастает до величины 4.2 МПа. После этого давление очень быстро увеличивается до значения 77.7 МПа на расстоянии около 3 см от плоскости лидирующей ударной волны. Этот рост вызван интенсивной фрагментацией капель воды в этой области из-за возникшей большой разности скоростей капель и расплавленного свинца. Образовавшиеся из-за этого мелкие фрагменты воды тормозятся силой трения между их поверхностью и жидким свинцом, вызывая рост локального давления в системе. Далее давление медленно уменьшается, асимптотически стремясь достигнуть давления в плоскости Чепмена–Жуге для равновесной смеси продуктов детонации.

Рис. 1.

Распределение давления в зоне фрагментации.

Рис. 2 показывает распределения скоростей фаз в зоне фрагментации. Тяжелая фаза (расплавленный свинец) входит в зону фрагментации практически со скоростью распространения волны детонации, что означает практически полное отсутствие торможения на ударной волне. Легкая фаза (капли воды) существенно замедляет свою скорость на ударной волне (с 222 м/с до 70 м/с). В результате после прохождения ударной волны из-за возникающей значительной разности скоростей капель воды и расплавленного свинца начинается интенсивная фрагментация капель, приводящая к усиленному тепловыделению в зоне фрагментации, что поддерживает распространение волны. Из-за сильного трения между несущей фазой и дисперсными фазами их скорости быстро сравниваются на расстоянии примерно 8 см от плоскости ударной волны. Распределения скоростей дисперсных фаз практически совпадают. Это объясняется тем, что образующиеся фрагменты воды имеют скорость материнской капли, а разница времен релаксации скоростей дисперсных фаз к скорости несущей фазы незначительна. При этом надо иметь ввиду, что размер материнской капли постоянно уменьшается в ходе ее фрагментации, что означает уменьшение времени ее скоростной релаксации.

Рис. 2.

Распределения скоростей фаз в зоне фрагментации. 1 – расплав свинца, 2$f$-фаза и $d$-фаза.

На рис. 3 представлены распределения температур фаз в зоне фрагментации капель. Видно, что массивный свинец охлаждается незначительно, в то время, как температуры легких фаз (капли и фрагменты воды) быстро растут, причем рост температуры $f$-фазы происходит быстрее, чем в случае $d$-фазы, что объясняется меньшим размером фрагментов воды.

Рис. 3.

Распределения температур фаз в зоне фрагментации. 1 – расплав свинца, 2$f$-фаза, 3$d$-фаза.

Сравнивая распределения скоростей фаз (рис. 2) и их температур (рис. 3), можно заключить, что скоростная релаксация фаз происходит значительно быстрее термической релаксации. Если скорости фаз выравниваются на расстоянии около 8 см от плоскости ударной волны, то температуры фаз на расстоянии метра от плоскости ударной волны еще значительно различаются.

Рис. 4 дает наглядное представление об исследуемом процессе в целом. На нем на плоскости “удельный объем многофазной смеси–давление смеси” построены адиабата Гюгонио (кривая 1) и ударная адиабата (кривая 2). Излом ударной волны на 90° означает полную конденсацию пара в паровой пленке на ударной волне и переход воды в состояние недогретой до температуры насыщения жидкости.

Рис. 4.

Адиабата Гюгонио (1) и ударная адиабата (2). Точки: A – исходное состояние, B – на ударной волне, C – точка Чепмена–Жуге, D – достижение максимального давления.

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

Точка A на рис. 4 показывает исходную смесь. Точка B представляет состояние смеси на ударной волне. Линия BD отображает состояние многофазной смеси при ее течении от плоскости ударной волны к плоскости Чепмена–Жуге. Видно, что в некоторый момент удельный объем смеси и ее давление совпадают с параметрами точки Чепмена–Жуге. Однако это не означает, что система пришла в равновесное состояние – в этот момент скорости фаз значительно отличаются между собой, то же самое имеет место и для температур фаз.

Далее система начинает двигаться от точки C к точке D, в которой она достигает максимального давления и максимального сжатия – ее удельный объем в этой точке минимален. После прохождения точки D начинается расширение многофазной системы, которое сопровождается снижением давления. Система движется в обратном направлении по линии DC, асимптотически приближаясь к точке Чепмена–Жуге.

Заключение. Для расчета внутренней структуры детонационной волны в смеси капель воды с расплавленным свинцом разработана одномерная стационарная модель для многофазной смеси, состоящей из трех фаз: жидкий свинец (сплошная фаза), капли воды (дисперсная фаза) и фрагменты воды, образующиеся при дроблении капель (вторая дисперсная фаза). Учитываются: сила межфазного трения, теплообмен и фрагментация капель воды. Результаты исследований показали, что в волне термической детонации в данной системе максимальное давление достигается в плоскости, расположенной между плоскостью лидирующей ударной волны и плоскостью Чепмена–Жуге.

Работа выполнена при финансовой поддержке Российского научного фонда в рамках проекта № 21-19-00709.

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

  1. Зельдович Я.Б., Компанеец А.С. Теория детонации. М.: Гостехиздат, 1955. 268 с.

  2. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика: Учебное пособие. В 10 тт. Т. VI. Гидродинамика. М.: Наука, 1986. 736 с.

  3. Крайко А.Н. Неустойчивость стационарных течений в каналах переменной площади поперечного сечения с детонационной волной Чемпена–Жуге // ПММ. 2019. Т. 83. № 3. С. 354–369.

  4. Вайнштейн П.Б., Нигматулин Р.И., Попов В.В. Переход конвективного горения аэровзвесей унитарного топлива в детонацию // Физика горения и взрыва. 1980. № 5. С. 102–106.

  5. Нигматулин Р.И., Вайнштейн П.Б., Ахатов И.Ш. Структура стационарных детонационных волн в смесях газа с частицами унитарного топлива // в сб.: Химическая физика процессов горения и взрыва. Детонация. Черноголовка: Ин-т хим. физ. АН СССР, 1980. 128 с. С. 96–99.

  6. Ахатов И.Ш., Вайнштейн П.Б., Нигматулин Р.И. Структура детонационных волн в газовзвесях унитарного топлива // Изв. АН СССР. МЖГ. 1981. № 5. С. 47–53.

  7. Нигматулин Р.И. Динамика многофазных сред. Ч. I. М.: Наука, 1987. 464 с.

  8. Board S.J., Hall R.W., Hall R.S. Detonation of a fuel coolant explosion // Nature. 1975. V. 254. P. 319–321.

  9. Мелихов В.И., Мелихов О.И., Якуш С.Е. Гидродинамика и теплофизика паровых взрывов. М.: ИПМех РАН, 2020. 276 с.

  10. Мелихов В.И., Мелихов О.И., Якуш С.Е. Термическое взаимодействие высокотемпературных расплавов с жидкостями // ТВТ. 2022. Т. 60. № 2. С. 280–318.

  11. Sharon A., Bankoff S.G. On the existence of steady supercritical plane thermal detonations // Int. J. Heat Mass Trans. 1981. V. 24. P. 1561–1572.

  12. Frost D.L., Lee J.H.S., Ciccarelli G. The use of Hugoniot analysis for the propagation of vapor explosion waves // Shock Waves. 1991. V. 1. P. 99–110.

  13. Iskhakov A.S., Melikhov V.I., Melikhov O.I., Yakush S.E., Le Tran Chung. Hugoniot analysis of experimental data on steam explosion in stratified melt-coolant configuration // Nucl. Engng.&Design. 2019. V. 347. P. 151–157.

  14. Dinh T.N. Multiphase flow phenomena of steam generator tube rupture in a lead-cooled reactor system: a scoping analysis // Proc. ICAPP 2007. Paper No. 7497. May 13–18, 2007. Nice, France.

  15. Iskhakov A.S., Melikhov V.I., Melikhov O.I. Hugoniot analysis of energetic molten lead water interaction // Annals of Nucl. Energy. 2019. V. 129. P. 437–449.

  16. Sobolev V. Database of thermophysical properties of liquid metal coolants for GEN-IV. Sodium, lead, lead-bismuth eutectic (and bismuth) // in: Belgian Nuclear Res. Centre. Sci. Rep. SCK CEN-BLG-1069. Boeretang, Belgium. 2010. P. 175.

  17. IAPWS (The Int. Assoc. for the Properties of Water&Steam). http://www.iapws.org

  18. Pilch M., Erdman C.A. Use of breakup time data and velocity history data to predict the maximum size of stable fragments for acceleration-induced breakup of a liquid drop // Int. J. Multiphase Flow. 1987. V. 13. P. 741–757.

  19. Meignen R., Picchi S., Lamome J. et al. The challenge of modeling fuel-coolant interaction: Pt. I – Premixing // Nucl. Engng.&Design. 2014. V. 280. P. 511–527.

  20. Fletcher D.F., Anderson R.P. A review of pressure-induced propagation models of the vapour explosion process // Prog. Nucl. Energy. 1990. V. 23. P. 137–179.

  21. Fletcher D.F. An improved mathematical model of melt/water detonations – I. Model formulation and example results // Int. J. Heat Mass Transfer. 1991. V. 34. № 10. P. 2435–2448.

  22. Безносов А.В., Пинаев С.С., Давыдов Д.В. и др. Экспериментальные исследования характеристик контактного теплообмена свинцовый теплоноситель–рабочее тело // Атомная энергия. 2005. Т. 98 (3). С. 182–187.

  23. Carachalios C., Burger M., Unger H. A transient two-phase model to describe thermal detonations based on hydrodynamic fragmentation // in: Proc. Int. Meeting on LWR Severe Accident Evaluation, Cambridge, Massachusetts, 28 Aug.–1 Sep. 1983.

  24. Ishii M., Hibiki T. Thermo-Fluid Dynamics of Two-Phase Flow. New York: Springer, 2011. 518 p.

  25. Kolev N.I. Multiphase Flow Dynamics. V. 1. Fundamentals. New York: Springer, 2015. 840 p.

Дополнительные материалы отсутствуют.