Журнал вычислительной математики и математической физики, 2023, T. 63, № 4, стр. 678-693

ВВЕДЕНИЕ

Неустойчивость границы раздела сред применительно к задачам инерциального термоядерного синтеза исследовалась во многих работах (см. [1–6]). Обзор численных методов, применяемых для моделирования таких неустойчивостей, приведен в [7].

В настоящей работе рассматривается мишень в виде цилиндра с термоядерным горючим, окруженного тяжелой оболочкой. Как впервые было показано в [8], при наличии небольшого участка с достаточно высокой температурой в такой мишени возникает близкая к стационарной незатухающая детонационная волна вдоль оси цилиндра. В [9–11] аналогичная задача рассматривалась применительно к торцевому зажиганию предварительно сжатой мишени пучком тяжелых ионов в соответствии с концепцией быстрого зажигания, которому предшествует сравнительно медленная стадия сжатия (см. [12], [13]).

Та же задача, но с другим зажигающим драйвером рассматривалась в нашей работе [14]. Вместо пучка тяжелых ионов был выбран пучок протонов с примерно в 10 раз меньшей глубиной прогрева горючего. Было обнаружено, что при зажигании мишени вначале образуется нестационарная детонационная волна с растущей по времени амплитудой. Превращение такой волны в близкую к стационарной сопровождается изменением механизма ее распространения. Вместо детонационной волны с ударно-волновым механизмом распространения возникает волна горения, основным механизмом распространения которой является нелокальный нагрев $\alpha $-частицами, возникающими в реакции синтеза ядер дейтерия и трития. Далее будем называть такую волну горения безударной. В отличие от волны медленного горения (см. [15]), механизмом распространения которой является обычная теплопроводность, рассматриваемая в нашей работе волна является быстрой (т.е. распространяется со скоростью порядка скорости детонационной волны). Возможность существования таких волн и некоторые приближенные формулы, основанные на анализе энергетического баланса, обсуждаются в [16] и цитируемых там работах.

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

В [17] рассматривалась задача с условием симметрии на некоторой плоскости, перпендикулярной оси цилиндра, что соответствует столкновению на этой плоскости двух одинаковых волн горения. Для выбранного небольшого размера мишени детонационная волна не успевает превратиться в безударную.

Применительно к очень малым возмущениям, которые остаются малыми для рассматриваемого интервала времени, изучена роль двух возможных механизмов развития неустойчивости границы раздела: ее импульсного ускорения детонационной волной (неустойчивость Рихтмайера–Мешкова) и высокоскоростного скольжения горючего вдоль границы раздела (неустойчивость Кельвина–Гельмгольца). Использовались известные формулы для скорости роста амплитуды малых синусоидальных возмущений плоской границы (см. [18], [19]), коэффициенты которых полагались зависящими от времени и брались из расчета рассматриваемой задачи. Как оказалось, для выбранного размера мишени основную роль играет неустойчивость Рихтмайера–Мешкова, так как в волне разрежения, примыкающей к волне горения, касательная к границе раздела скорость горючего быстро уменьшается по мере удаления от волны. Следует отметить, что по мере движения волны влияние этого эффекта на неустойчивость границы вблизи волны уменьшается, так как уменьшение касательной составляющей скорости в волне разрежения замедляется.

В [17] проанализирована также эволюция небольшого локального возмущения границы раздела, возникающего при формировании детонационной волны. Для определения амплитуды возмущения использовалась некоторая обработка расчетного профиля границы раздела. Для выбранного размера мишени оказалось, что скорость роста амплитуды возмущения до отражения детонационной волны определяется неустойчивостью Рихтмайера–Мешкова, а значительное увеличение скорости роста амплитуды после отражения вызвано неустойчивостью Кельвина–Гельмгольца.

В настоящей работе неустойчивость границы раздела исследуется прямым численным моделированием задачи из работы [14], включая стадию образования и распространения безударной волны. Исследование ограничено двумерными возмущениями границы. Рассматривается эволюция как задаваемых начальных возмущений, так и спонтанных возмущений, которые самопроизвольно возникают при расчете неустойчивых течений (см., например, [20], где спонтанные возмущения используются для анализа неустойчивости детонационной волны). Важный вопрос об устойчивости волны горения выходит за рамки настоящей работы.

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

Используется подвижная регулярная сетка, одна из линий которой является границей раздела. Такая вычислительная технология не позволяет моделировать турбулентное перемешивание веществ на поздней стадии развития неустойчивости. Тем не менее изучение неустойчивости на ее ранней стадии может давать результаты, применимые и к стадии турбулентного перемешивания. В качестве примера приведем формулу для изменения скорости роста амплитуды немалого (отношение начальной амплитуды к длине волны возмущения лежит примерно между 0.5 и 1) возмущения при повторном импульсном ускорении границы раздела (см. [21], [22]), полученную по результатам расчета двумерных возмущений с явным выделением границы раздела в виде линии сетки, которая неплохо описывает изменение скорости роста ширины зоны турбулентного перемешивания сразу после повторного ударного сжатия (см. [23]).

1. ПОСТАНОВКА ЗАДАЧИ И ЧИСЛЕННЫЙ МЕТОД

Расчетная область в начальный момент времени в цилиндрических координатах $(z,r)$ показана на фиг. 1. Горючее имеет форму цилиндра радиусом $R$ = 64 мкм ($z \geqslant 0$, $r \leqslant R$), а оболочка является цилиндрическим слоем размером $\Delta R$ = 60 мкм ($z \geqslant 0$, $R \leqslant r \leqslant R + \Delta R$). Перед цилиндром располагается дополнительная область в виде усеченного конуса, предназначенная для расчета горючего, вылетающего из мишени. Первоначально эта область заполнена DT-газом небольшой плотности и атмосферного давления.

В качестве зажигающего драйвера рассматривается моноэнергетический пучок протонов с энергией 1 МэВ. Интенсивность ${{J}_{{\text{b}}}}$ пучка не зависит от $r$ при $r \leqslant R$ и равна нулю при $r > R$. Функция ${{J}_{{\text{b}}}}(t)$ определяется двумя параметрами: максимальной интенсивностью ${{J}_{0}} = 4.5 \times $ × 10¹⁹ Вт/см² и временем действия пучка $\Delta {{t}_{{{\text{pr}}}}} = 16$ пс, определяющим, в свою очередь, начальный интервал времени $\Delta {{t}_{{{\text{0pr}}}}} = 0.02\Delta {{t}_{{{\text{pr}}}}}$, в течение которого ${{J}_{{\text{b}}}}(t)$ увеличивается от 0 до ${{J}_{0}}$.

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

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

Границы AB и BC лежат на неподвижной прямой линии, но точки A и B могут двигаться вдоль этой линии в соответствии с движением боковой границы оболочки и границы раздела между горючим и оболочкой.

Математическая модель основана на уравнениях одножидкостной двухтемпературной гидродинамики, которые применительно к горючему имеют вид

где $d{\text{/}}dt = \partial {\text{/}}\partial t + {\mathbf{u}}\nabla $ – лагранжева производная по времени; $\rho $ – плотность; ${\mathbf{u}} = ({{u}_{z}},{{u}_{r}})$ – скорость; ${{p}_{{\text{e}}}}$ и ${{p}_{{\text{i}}}}$ – электронное и ионное давления, из которых складывается полное давление $p = {{p}_{{\text{e}}}} + {{p}_{{\text{i}}}}$; ${{\varepsilon }_{{\text{e}}}}$ и ${{\varepsilon }_{{\text{i}}}}$ – внутренняя энергия электронов и ионов; ${{{\mathbf{q}}}_{{\text{e}}}}$ и ${{{\mathbf{q}}}_{{\text{i}}}}$ – тепловые потоки электронов и ионов; ${{Q}_{{{\text{ei}}}}}$ определяет обмен энергией между электронами и ионами. Последние члены в уравнениях (1) и (2) определяют нагрев электронов и ионов пучком протонов (${{D}_{{\text{e}}}}$ и ${{D}_{{\text{i}}}}$) и $\alpha $-частицами (${{W}_{{\text{e}}}}$ и ${{W}_{{\text{i}}}}$), а также обмен энергией между электронами и излучением ($S$).

При описании DT-смеси мы опускаем стадию ударной ионизации (см. [24]) атомов смеси протонами зажигающего пучка, полагая степень ионизации смеси равной единице.

Для DT-смеси используются уравнения состояния электронной и ионной компонент, полученные из широкодиапазонного полуэмпирического однотемпературного (${{T}_{e}} = {{T}_{i}} = T$) уравнения состояния водорода $p = {{p}_{H}}(\rho ,T)$, $\varepsilon = {{\varepsilon }_{H}}(\rho ,T)$, основанного на модели [25]. Такое уравнение состояния удовлетворительно описывает сильно сжатый водород как при низких, так и при высоких температурах. Однотемпературное уравнение состояния для DT-смеси берется в виде

где $A = 2.5$ – атомный вес смеси. Двухтемпературные уравнения состояния смеси определяются формулами где $\xi = 1$ – степень ионизации. Для холодной смеси с ${{T}_{e}} = {{T}_{i}} = T$ получаем а для высокотемпературной плазмы, когда $p(\rho ,T)$ и $\varepsilon (\rho ,T)$ являются уравнениями состояния совершенного однотемпературного газа при $T = {{T}_{e}}$ и $T = {{T}_{i}}$, формулы (3), (4) дают двухтемпературные уравнения состояния совершенного газа с заданной степенью ионизации.

Ранее мы провели контрольные расчеты одномерной задачи с двухтемпературными уравнениями состояния (3), (4) и с некоторыми разумными формулами для уравнения состояния электронной компоненты (см. [26]). Разница в результатах была несущественная.

Оболочка описывается теми же уравнениями, но без правых частей уравнений (1) и (2). Как и в работах [9–11], в качестве материала оболочки выбрано золото. Уравнения состояния получены по формулам (3), (4) из однотемпературных уравнений состояния по модели Томаса–Ферми с квантовыми и обменными поправками (см. [27], [28]), которая позволяет определять и степень ионизации. Так как модель Томаса–Ферми не может правильно описать сильно сжатое холодное вещество, начальные значения давления и плотности оболочки вычисляются по уравнению состояния [29] (см. ниже), а начальное значение температуры определяется из уравнения состояния Томаса–Ферми.

Перенос излучения описывается многогрупповым по частоте $\nu $ и диффузионным по телесному углу приближением стационарного уравнения переноса. Многогрупповое приближение строится делением интервала $0 < \nu < \infty $ на $N$ групп ${{\nu }_{{l - 1}}} < \nu < {{\nu }_{l}}$, $l = 1, \ldots ,N$, ${{\nu }_{0}} = 0$, ${{\nu }_{N}} = \infty $.

По сравнению с работой [14] изменена модель переноса $\alpha $-частиц. Вместо стационарного кинетического уравнения в приближении Фоккера–Планка используется его диффузионное приближение из [30], что позволило значительно уменьшить объем вычислений. Чтобы получить решение с временем выхода на квазистационарный режим, близким к работе [14], была несколько уменьшена интенсивность ${{J}_{0}}$ зажигающего пучка протонов.

Заметим, что в силу условия Чепмена–Жуге возмущения из области зажигания не могут догнать квазистационарную волну горения, и поэтому параметр ${{J}_{0}}$, как и время действия зажигающего пучка $\Delta {{t}_{{{\text{pr}}}}}$, влияют только на сам факт зажигания и на время появления квазистационарной волны. Параметры квазистационарной волны зависят только от начальной плотности горючего, радиуса цилиндра, толщины и материала оболочки.

Горючее полагается предварительно сжатым изоэнтропически из состояния при температуре 4 К, давлении 0.1 МПа и плотности ${{\rho }_{a}} \approx 0.22$ г/см³ до плотности ${{\rho }_{0}} = 500{{\rho }_{a}}$, что дает давление ${{p}_{0}} \approx 400$ ТПа. Оболочка также предполагается предварительно сжатой изоэнтропически (использовалось уравнение состояния золота из [29]) до давления ${{p}_{0}}$, что дает плотность сжатой оболочки примерно 300 г/см³.

Более подробное описание постановки задачи приведено в [14].

Для уравнений гидродинамики используется двухэтапная схема типа [31], [32]. На лагранжевом этапе используется схема второго порядка точности [33]. На этапе пересчета искомых функций на эйлерову сетку используется схема третьего порядка точности [34].

Для уравнений диффузии и теплопроводности используется схема второго порядка точности по пространственным переменным на 9-точечном шаблоне из [35]. Система линейных уравнений решается прямым методом разложения на треугольные множители для ленточной матрицы (см. [36]).

Методика расчета движения сеточных узлов, расположенных на границе раздела сред и свободной границе, взята из книги [37]. На одном шаге по времени узлы границы вначале движутся в соответствии с полем скорости. Затем вычисляются крайние узлы как точки пересечения границы и прямых линий на фиг. 1. На последнем этапе граница аппроксимируется дугами окружностей, проходящими через узлы, после чего находятся новые положения узлов, равномерно расставленных по длине границы.

Для построения регулярных криволинейных сеток с заданными граничными узлами используется метод [38], основанный на методе [39]. В настоящей работе метод [38] модифицирован в соответствии с работой [40]. Так как внутренние узлы начальной сетки для метода [38] берутся из предыдущего временнóго слоя, начальная сетка может оказаться невыпуклой (т.е. состоящей не только из выпуклых четырехугольников), что снижает возможность метода [38] генерировать выпуклые сетки. Указанная выше модификация позволила повысить надежность метода для невыпуклых начальных сеток.

Параметры разностной сетки, если это не оговорено особо, следующие. Шаг начальной сетки в основной расчетной области вдоль оси $z$ ${{h}_{z}} = 0.25$ мкм, число узлов сетки вдоль оси $r$ ${{N}_{r}} = 40$ в горючем и 45 в оболочке. Проводились расчеты и на более подробных сетках с ${{h}_{z}} = 0.125$ мкм и ${{N}_{r}} = 80$.

2. ТЕСТОВАЯ ЗАДАЧА

Имеются два плоских слоя: один из DT-смеси, другой из золота, которые касаются друг друга, как показано на фиг. 2. Рассматривается плоское течение, зависящее от двух пространственных переменных $(x,y)$, которое описывается уравнениями бездиссипативной однотемпературной гидродинамики. В начальный момент искомые функции постоянны в каждом слое. Термодинамические функции в слое золота совпадают с начальными данными основной задачи, а в слое DT-смеси примерно соответствуют значениям в безударной волне основной задачи, ${{p}_{{{\text{DT}}}}} = $ = 600 ППа, ${{\rho }_{{{\text{DT}}}}} \approx {{\rho }_{0}}$. Слой золота неподвижен, а слой DT-смеси имеет скорость ${{U}_{n}}$ вдоль оси $y$ и скорость ${{U}_{\tau }}$ вдоль оси $x$. Скорость ${{U}_{\tau }}$ задает разрыв касательной составляющей скорости на контактной поверхности, а скорость ${{U}_{n}} \approx 880$ км/c равна скорости за фронтом ударной волны из решения задачи о распаде разрыва для неподвижного слоя DT-смеси плотностью $1.5{{\rho }_{0}}$. Эта задача определяет также значение ${{\rho }_{{{\text{DT}}}}}$. При таком выборе ${{U}_{n}}$ и ${{\rho }_{{{\text{DT}}}}}$ решение задачи о распаде разрыва для подвижного слоя DT-смеси содержит только ударную волну в слое золота, а значения искомых функций в слое DT-смеси не меняются.

Рассмотрим начальное синусоидальное возмущение границы раздела

где ${{d}_{0}}$ – амплитуда возмущения, $\lambda $ – длина волны. Возмущение предполагается малым (${{d}_{0}}{\text{/}}\lambda \ll 1$). Зависимость амплитуды возмущения $d(t)$ определяется формулой при ${{U}_{\tau }} = 0$ (неустойчивость Рихтмайера–Мешкова [18]), и формулой при ${{U}_{n}} = 0$ (неустойчивость Кельвина–Гельмгольца [19]). Здесь ${{\rho }_{1}}$ – плотность DT-смеси, ${{\rho }_{2}}$ ‒ плотность золота после прохождения ударной волны (в соответствии с рекомендацией работы [41]).

Для произвольных значений ${{U}_{n}}$ и ${{U}_{\tau }}$ возьмем формулу

для которой $d(0) = {{d}_{0}}$, а производная $d_{t}^{'}$ равна сумме производных от формул (5) и (6).

Фигура 3 демонстрирует возможность моделирования нелинейной стадии неустойчивости Кельвина–Гельмгольца ($d{\text{/}}\lambda > 1$) c характерными наклонными гребнями (см., например, [19]).

На фиг. 4 показаны некоторые результаты расчета тестовой задачи в сравнении с формулой (7). Так как в расчете имеется небольшой начальный интервал времени, в течение которого формируется ударная волна, формула (7) показана с небольшим сдвигом по времени. В случае неустойчивости Рихтмайера–Мешкова (${{U}_{\tau }} = 0$, фиг. 4a) расчет дает качественно верную зависимость скорости роста амплитуды от времени – вначале линейный рост, а затем его постепенное уменьшение (см., например, [1]). Расчеты на разных сетках по $x$ имеют тенденцию к сходимости: расчет с 32 разностными интервалами в одной волне отличается от расчета с 64 интервалами заметно меньше, чем от расчета с 16 интервалами. Удвоение числа узлов сетки поперек границы раздела не влияло заметно на скорость роста амплитуды.

На фиг. 4б приведен результат расчета для скорости скольжения ${{U}_{\tau }} = 400$ км/с, примерно равной скорости скольжения в безударной волне основной задачи. Видно, что при той же начальной амплитуде точность воспроизведения формулы (7) намного выше. При увеличении начальной амплитуды точность воспроизведения формулы (7) уменьшается, что, скорее всего, связано с уменьшением точности самой формулы (7) относительно решения уравнений гидродинамики.

В заключение этого раздела приведем результаты расчета локального начального возмущения в виде равнобедренного треугольника, основание которого назовем длиной волны $\lambda $, а высоту – амплитудой возмущения ${{d}_{0}}$, ${{d}_{0}}{\text{/}}\lambda = 0.1$. На фиг. 5а для случая неустойчивости Рихтмайера–Мешкова (${{U}_{\tau }} = 0$) приведены профили границы раздела в три момента времени. Точки излома начального возмущения немного размазаны, но возмущение в целом остается почти прямолинейным.

При наличии скольжения (${{U}_{\tau }} = 400$ км/с, фиг. 5б) возмущение разрушается на более коротковолновые возмущения разной амплитуды и характерными для неустойчивости Кельвина–Гельмгольца наклонными гребнями. Скорость роста максимальной амплитуды (разница между максимальным и минимальным значениями $y$), как и в случае синусоидальных возмущений, заметно выше, чем в случае ${{U}_{\tau }} = 0$.

3. СПОНТАННЫЕ ВОЗМУЩЕНИЯ

Спонтанные возмущения границы раздела возникают сразу, когда она начинает двигаться под действием растущего давления горючего, нагреваемого пучком протонов. На фиг. 6 приведен фрагмент границы раздела в момент прекращения действия пучка $t = 0.016$ нс. Граница заметно сдвинулась относительно начального положения $r = 0.064$ мм. Профиль границы имеет вид почти прямой линии, на котором возмущений не видно. Чтобы увидеть их, на фиг. 6 показан профиль по координате $z$ производной $ - dr{\text{/}}dl$, где $l$ – длина профиля (отсчитываемая от точки B на фиг. 1) в сторону возрастания $z$, профиль границы раздела предполагается заданным в параметрическом виде $r = r(l)$, $z = z(l)$.

Чтобы оценить роль ошибки округления компьютера (которая определяется количеством двоичных разрядов для хранения мантиссы числа) в формировании спонтанных возмущений, расчет проводился как для обычного формата чисел типа double на языке C, так и для расширенного формата “long double” с увеличенной на 11 разрядов мантиссой (т.е. с ошибкой округления, уменьшенной в ${{2}^{{11}}} \approx 2000$ раз). Видно, что погрешность округления не играет заметной роли. Амплитуда спонтанных возмущений, как и расположение трех очагов возмущений, изменились незначительно.

При формировании детонационной волны ($t \approx 0.03$ нс) возникает локальное возмущение, которое анализировалось нами в работе [17] применительно к задаче с отражением детонационной волны от плоскости симметрии. В настоящей работе мы ограничимся демонстрацией этого возмущения на фиг. 7 через достаточно большой временной интервал после его возникновения. Профили локального возмущения похожи на профили треугольного локального возмущения для тестовой задачи на фиг. 5.

На фиг. 8 показана граница раздела при $t = 0.08$ нс. Превращение нестационарной детонационной волны в квазистационарную безударную волну еще не завершилось (безударная волна присутствует в качестве предвестника детонационной волны, подробнее см. [14]). Профиль границы раздела $r(z)$ имеет волнообразный характер вплоть до $z \approx 0.13$ мм (точка с максимальным давлением, которую естественно рассматривать в качестве положения волны горения, находится заметно правее ${{z}_{w}} \approx 0.165$ мм). На фиг. 8а показан фрагмент границы раздела с волнообразным профилем при $0.1 \leqslant z \leqslant 0.128$ мм, что достаточно далеко от локального возмущения на фиг. 7.

Отсутствие характерных для неустойчивости Кельвина–Гельмгольца наклонных гребней объясняется достаточно быстрым падением скорости горючего за волной горения: в данном случае скорость горючего обращается в нуль при $z \approx 0.1$ мм (см., также, [17]).

Профиль $r(z)$ на фиг. 8а для сетки с ${{h}_{z}} = 0.125$ мкм (сплошная линия) имеет любопытную структуру. Стрелками показаны четыре точки, где профиль имеет локальные максимумы, в промежутках между которыми имеются небольшие области, где локальная амплитуда возмущения близка к нулю, что можно трактовать как пропуск одной волны. Для наглядности этот профиль был обработан так, чтобы получить локальную амплитуду $A(z)$. Точки локального экстремума профиля $r(z)$ соединялись прямыми, что давало две ломаные линии: одну, построенную по точкам максимума, и другую, построенную по точкам минимума. Амплитуда $A(z)$ определяется как расстояние между ломаными по оси $r$ в точке $z$.

На фиг. 8б приведены функции $A(z)$, полученные обработкой профилей на фиг. 8а для двух сеток с ${{h}_{z}} = 0.125$ (сплошная линия) и 0.25 мкм (штриховая линия). Для сетки с ${{h}_{z}} = 0.125$ мкм видны четыре максимума и три минимума между ними, которые упоминались выше в связи с фиг. 8а. Укажем на близость максимальных значений $A(z)$ и расстояний между ними для разных сеток в левой части приведенного фрагмента.

Выполнить расчет для больших значений $t$ на сетке с шагом ${{h}_{z}} = 0.125$ мкм не удается из-за уменьшения шага по времени, вызванного искривлением сетки в окрестности локального возмущения. Приводимые ниже результаты получены на сетке с ${{h}_{z}} = 0.25$ мкм.

Рассмотрим процесс превращения нестационарной волны горения, имеющей вид детонационной волны с безударным предвестником, в квазистационарную безударную волну (см. [14]). На фиг. 9 показаны зависимости от времени максимального давления вдоль оси симметрии ${{p}_{{ax}}}(t)$ и вдоль границы раздела ${{p}_{b}}(t)$. Так как оба максимума расположены внутри волны горения, эти функции можно рассматривать как давления в волне горения.

Функции ${{p}_{{ax}}}(t)$ и ${{p}_{b}}(t)$ вначале растут, что указывает на нестационарность волны горения, а при $t \approx 0.1$ нс начинается переход обеих функций к примерно постоянным значениям, соответствующим квазистационарной волне горения.

Укажем на колебательный характер функции ${{p}_{b}}(t)$ при $t \gtrsim 0.1$ нс, что скорее всего связано с неустойчивостью стационарной волны Чепмена–Жуге (см., например, [20]). Как показано в [14], в квазистационарной безударной волне с хорошей точностью выполняется условие Чепмена–Жуге.

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

На фиг. 10 показаны соответствующие профили вблизи волны горения, координату которой ${{z}_{w}}$ определим как $z$-координату точки с максимальным давлением на границе раздела. Один профиль при $t = 0.08$ нс, когда неустойчивость волны горения еще не проявляется в виде колебательного режима функции ${{p}_{b}}(t)$, и другой профиль при $t = 0.15$ нс, когда колебательный режим функции ${{p}_{b}}(t)$ имеет место уже на достаточно большом интервале времени. Видно, что при $t = 0.08$ нс спонтанных возмущений вблизи волны горения нет, а при $t = 0.15$ нс спонтанные возмущения присутствуют.

На фиг. 11 показан результат описанной выше обработки профиля $r(z)$ при $t = 0.15$ нс. Интервал $0 < z < 0.1$, где располагается локальное возмущение большой амплитуды (см. фиг. 7), исключен из обработки. При $z > 0.41$ колебания настолько малы, что локальные экстремумы функции $r(z)$ не определяются (тем не менее эти колебания можно увидеть на профиле $dr{\text{/}}dl(z)$, см. фиг. 10).

На фиг. 11 показаны как амплитуда $A(z)$ (нижняя линия), так и длина волны $\lambda (z)$ (верхняя линия), под которой понимается расстояние между соседними локальными максимумами профиля $r(z)$. При замене локальных максимумов на локальные минимумы функция $\lambda (z)$ меняется незначительно.

На фоне чередующихся максимумов и минимумов функции $A(z)$ выделяются два достаточно больших интервала, не имеющих минимумов $A(z)$. На фиг. 11 эти интервалы отмечены стрелками и цифрами. Если отбросить пики профиля $\lambda (z)$, то длина волны в области между двумя интервалами заметно меньше, чем в области левее интервала 1. Укажем также на наличие областей, где локальные экстремумы амплитуды, отмеченные на фиг. 11 точками, располагаются на примерно одинаковом расстоянии друг от друга.

4. ЗАДАВАЕМЫЕ ВОЗМУЩЕНИЯ

Ограничимся рассмотрением взаимодействия задаваемых возмущений границы раздела с квазистационарной волной горения, которая состоит из головной части с линейным профилем давления и скорости (см. [14]), и основной части, включающей в себя точку с максимальным давлением вдоль границы раздела. В головную часть спонтанные возмущения не проникают из-за условия Чепмена–Жуге, а в основной части они слишком малы (их можно увидеть только на профиле $dr{\text{/}}dl$, см. фиг. 10), чтобы повлиять на развитие задаваемых возмущений.

Начальное возмущение задавалось перед уже сформировавшейся квазистационарной волной при $t = 0.1$ нс, начиная с точки ${{z}_{0}} \approx 0.33$ мм, расположенной на небольшом расстоянии от головной части волны в области радиационного предвестника, где все функции, включая $r$-координату границы раздела, незначительно отличаются от своих значений при $t = 0$. Для всех узлов $i$ вдоль границы раздела, начиная с ближайшего к ${{z}_{0}}$, $r$-координата мгновенно изменяется на величину $\Delta {{r}_{i}} = 0.5{{d}_{0}}\sin (2\pi {{z}_{i}}{\text{/}}\lambda )$, где ${{d}_{0}}$ – амплитуда, $\lambda $ – длина волны начального возмущения.

Длина волны выбирается равной 8 мкм, что для сетки с начальным шагом ${{h}_{z}} = 0.25$ мкм дает $\Delta n = 32$ разностных интервала на длине волны, как и в расчетах для тестовой задачи из разд. 2. Если учесть, что шаг сетки вдоль длины границы раздела при $t = 0.1$ нс несколько больше ${{h}_{z}}$, получается немного меньшее значение $\Delta n \approx 30.5$.

Результаты расчета при $t = 0.14$ нс представлены на фиг. 12 для двух значений начальной относительной амплитуды ${{d}_{0}}{\text{/}}\lambda $. На фиг. 12а для ${{d}_{0}}{\text{/}}\lambda = 0.01$ показаны профиль $r(z)$ границы раздела и поле давления в горючем. Упоминавшиеся выше области течения показаны стрелками и цифрами: головная часть волны горения (область 1) при $0.4 \lesssim z \lesssim 0.48$ мм (правая граница области располагается правее показанного фрагмента течения), основная часть волны горения при $0.37 \lesssim z \lesssim 0.4$ мм (область 2) и волна разрежения при $z \lesssim 0.37$ мм (область 3).

Точка начала колебаний сместилась вправо относительно исходного значения ${{z}_{0}}$, что связано с движением горючего по направлению движения волны горения в самой волне и сразу за ней в волне разрежения. В головной части волны горения амплитуда возмущения практически не растет. Ее рост начинается в основной части волны горения, где скорость скольжения максимальна. Наличие характерных наклонных гребней указывает на неустойчивость Кельвина–Гельмгольца как на основной механизм развития неустойчивости на нелинейной стадии.

На фиг. 12б показан результат описанной выше обработки профилей $r(z)$ для ${{d}_{0}}{\text{/}}\lambda = 0.01$ и 0.1. Для ${{d}_{0}}{\text{/}}\lambda = 0.1$ амплитуда возмущения заметно больше, чем для 0.01, и начинает расти уже в головной части волны, где появляются характерные наклонные гребни неустойчивости Кельвина–Гельмгольца.

Введем в рассмотрение две $z$-координаты характерных точек на границе раздела: рассмотренную ранее координату ${{z}_{w}}$ точки с максимальным давлением, которую можно выбрать в качестве координаты волны горения, и координату $z_{*}^{{}}$, в которой заканчивается головная часть волны. Так как волна горения близка к стационарной, разность $z_{*}^{{}} - {{z}_{w}}$ можно полагать не зависящей от времени.

Естественно предположить, что спонтанные возмущения не проникают в область ${{z}_{w}} \leqslant z \leqslant z_{*}^{{}}$ из-за условия Чепмена–Жуге. Если предположить также, что спонтанные возмущения вдоль волны горения (см. [20]) не меняют ее структуру и квазистационарный режим движения, то приведенные на фиг. 12б функции $d(z)$ на интервале ${{z}_{w}} \leqslant z \leqslant z_{*}^{{}}$ можно считать универсальными, т.е. функция $d\left( {{{z}_{w}} + z{\kern 1pt} '} \right)$ не зависит от $t$ на интервале $0 \leqslant z{\kern 1pt} ' \leqslant z_{*}^{{}} - {{z}_{w}}$. При необходимости можно приближенно вычислить скорость волны и получить универсальную функцию амплитуды возмущения от относительного времени $t - t_{*}^{{}}$, где $t_{*}^{{}}$ – момент прохождения правой границы головной части волны через $z_{*}^{{}}$.

ЗАКЛЮЧЕНИЕ

Исследование ограничено двумерными возмущениями границы раздела сред. Использовалась вычислительная технология, основанная на явном выделении границы раздела в виде одной из линий регулярной сетки. Для построения регулярных сеток в областях с криволинейными границами используется барьерный метод, надежно генерирующий сетки из выпуклых четырехугольников. Демонстрируется возможность применения технологии для расчета нелинейной стадии развития неустойчивости с отношением амплитуды возмущения к длине волны больше единицы. Возможность моделирования неустойчивости Рихтмайера–Мешкова и Кельвина–Гельмгольца тестируется на упрощенной задаче с однородным вдоль границы раздела возмущением.

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

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

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

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

Изучено взаимодействие задаваемого синусоидального возмущения с квазистационарной безударной волной горения, которая состоит из головной части с линейным профилем давления и основной части, являющейся окрестностью точки с максимальным давлением на границе раздела. Рассмотрены возмущения с отношением начальной амплитуды к длине волны 0.01 и 0.1. В обоих случаях с ростом амплитуды возмущения возникают характерные наклонные гребни, что указывает на неустойчивость Кельвина–Гельмгольца как на основной механизм развития неустойчивости на нелинейной стадии. Волна горения не разрушается.

Если предположить, что спонтанные возмущения не проникают в область волны горения из-за условия Чепмена–Жуге, а возможные спонтанные возмущения вдоль волны горения не меняют структуры волны и ее квазистационарность, то из полученных профилей амплитуды колебаний можно выделить универсальную часть, не зависящую от времени.