Астрономический журнал, 2023, T. 100, № 1, стр. 6-18

2.1. Постановка задачи

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

Лабораторные эксперименты по исследованию пространственного распределения пучков ускоренных протонов проводились в Ц-НИИМАШ группой В.С. Беляева и А.П. Матафонова на созданной ими 10 ТВт пикосекундной лазерной установке “Неодим” [16]. Описание установки и результатов эксперимента, а также численного моделирования эксперимента в аксиально-симметричной МГД модели с полоидальным магнитным полем изложено в работе [15].

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

Масштабы физических процессов в астрофизических джетах и масштабы в джетах, возникающих при воздействии лазерного импульса на мишень, сильно отличаются. Однако применение критериев подобия позволяет получить соответствие ряда характеристик астрофизического и лабораторного джета. Обоснование этого метода можно найти в работах [21, 22]. Ранее мы использовали методы подобия для исследования замагниченных джетов в полоидальном магнитном поле [15]. С помощью масштабирования и применения критериев подобия было получено удовлетворительное сходство ряда параметров джетов, истекающих из молодых звезд, и джета, наблюдаемого на установке “Неодим”. В настоящей работе мы моделируем образование и коллимацию джета с тороидальным магнитном полем.

2.2. Область моделирования

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

Для моделирования задачи выбрана инерциальная цилиндрическая система координат $\left( {r,\phi ,z} \right)$ с центром в середине мишени, ось $z$ перпендикулярна мишени. В задаче предполагается осевая симметрия в распределении всех макроскопических величин $\rho $, $T$, $v$, $B$, $\partial {\text{/}}\partial \phi = 0$, однако вычисляются все три компонента скорости $v$ и магнитного поля $B$. Моделирование проводится в половине координатной плоскости $r - z$. В задаче используется равномерная разностная сетка, координаты узлов сетки лежат в пределах $0 \leqslant R \leqslant {{R}_{{\max }}}$, $0 \leqslant z \leqslant {{Z}_{{\max }}}$. Моделирование проводится на равномерной сетке $(r,z)$ с числом ячеек ${{N}_{R}} \times {{N}_{Z}} = 257 \times 513$, или ${{N}_{R}} \times {{N}_{Z}} = 513 \times 1026$.

Центральная область мишени, к которой однородно и мгновенно прикладывается энергия лазерного излучения, представляет собой однородный цилиндр с радиусом ${{R}_{d}}$ и толщиной ${{Z}_{d}}$. При этом ${{R}_{d}} \ll {{R}_{{\max }}}$, ${{Z}_{{\max }}}$ и ${{Z}_{d}} \ll {{R}_{{\max }}}$, ${{Z}_{{\max }}}$. Схема области моделирования приведена на рис. 1.

2.3. Система уравнений

Для численного моделирования процессов образования и коллимации джетов в лабораторном эксперименте используется стандартная система уравнений магнитной гидродинамики с учетом конечной проводимости [20] в отсутствие гравитации.

Переменные в системе уравнений имеют следующие значения: $\rho $ – плотность вещества, $v$ – скорость, $B$ – индукция магнитного поля, $\varepsilon $ – внутренняя энергия единицы массы вещества, ${\mathbf{J}}$ – плотность электрического тока, $\sigma $ – эффективная электропроводность среды, $c$ – скорость света. Для описания термодинамических свойств вещества используется адиабатическое уравнение состояния идеального газа: Здесь $\gamma $ – показатель адиабаты. Для текущей задачи выбрано значение $\gamma = 5{\text{/}}3$, которое соответствует показателю адиабаты одноатомного нерелятивистского газа. В системе уравнений (1)–(5) учтен закон Ома в виде: Магнитная вязкость ${{\eta }_{m}} \equiv {{c}^{2}}{\text{/}}(4\pi \sigma )$, проводи-мость $\sigma $ и показатель адиабаты $\gamma $ считаются постоянными во всей области расчета и не изменяются во времени.

В предыдущей серии расчетов [15] магнитное поле в области моделирования имело полоидальную конфигурацию (${{B}_{z}}$, ${{B}_{r}}$). Для точного выполнения условия равенства нулю дивергенции магнитного поля $\nabla \cdot {\mathbf{B}} = 0$ в расчетах вместо полоидального компонента магнитного поля ${{B}_{z}}$ и ${{B}_{r}}$ использовался тороидальный компонент ${{A}_{\phi }}$, вектор потенциала магнитного поля ${\mathbf{A}}$, ${\mathbf{B}} = \nabla \cdot {\mathbf{A}}$. В текущих расчетах в области задано тороидальное магнитное поле ${{B}_{\phi }}$, для которого в качестве наиболее подходящей модели выбрано поле проводника с током, направленного вдоль оси $Z$, как схематично показано на рис. 1.

После приведения системы уравнений (1)–(5) к безразмерному виду в ней возникают два характерных безразмерных параметра. Первый из них $\beta $,

определяет соотношение между характерными значениями давления вещества ${{P}_{{0jet}}} = {{\rho }_{{0jet}}}c_{{s0}}^{2}{\text{/}}\gamma $ и энергии единицы объема магнитного поля. Здесь ${{\rho }_{{0jet}}}$ – плотность втекающего вещества на границе расчетной области, ${{P}_{{0jet}}}$, ${{c}_{{s0}}}$ – давление и скорость звука во втекающем веществе, ${{B}_{0}}$ – индукция магнитного поля в области, ${{V}_{{A0}}} = \frac{{{{B}_{0}}}}{{\sqrt {4\pi {{\rho }_{{0jet}}}} }}$ – альфвеновская скорость втекающего вещества на границе расчетной области.

Второй параметр – характерная безразмерная магнитная вязкость (или обратное ей характерное магнитное число Рейнольдса Re_M):

здесь $\sigma $ – постоянная эффективная проводимость плазмы; ${{L}_{0}}$ – характерный размер и ${{V}_{{A0}}}$ – характерная альфвеновская скорость, используемая при определении магнитного числа Рейнольдса Re_M. В качестве характерного размера ${{L}_{0}}$ был взят масштаб задачи ${{R}_{{\max }}}$. В качестве характерной скорости ${{V}_{0}}$ взята альфвеновская скорость втекающего вещества на границе расчетной ${{V}_{{A0}}}$.

Для решения системы уравнений использовалась осесимметричная разностная МГД схема с конечной проводимостью гибридного типа, основанная на методе локальных итераций [23], и методе коррекции потоков [24, 25]. Для решения астрофизических задач эта схема впервые была использована в работах [26, 27], и получила дальнейшее развитие в работах [28, 29].

2.4. Граничные и начальные условия

Система уравнений (1)–(5) численно решалась в цилиндрической области $(0 \leqslant z \leqslant {{Z}_{{\max }}}$, $0 \leqslant r \leqslant {{R}_{{\max }}})$. Предполагалось, что на внешних границах есть магнитное поле, но отсутствуют токи, и были заданы свободные граничные условия, т.е. для любой величины $f$ $\partial f{\text{/}}\partial n = 0$. Здесь $\partial f{\text{/}}\partial n$ – производная по нормали к поверхности, т.е. на внешней границе $f(r,z) = {\text{const}}$. В соответствии с предположением об осевой симметрии, граничные условия на оси $z$ ($r = 0$, $0 \leqslant z \leqslant {{Z}_{{\max }}}$) записываются как

На границе $(z = {{Z}_{{{\text{min}}}}},$ ${{R}_{d}} \leqslant r \leqslant {{R}_{{\max }}})$ были заданы фиксированные граничные условия, $\rho = {{\rho }_{0}}$ и ${{v}_{0}} = 0$, чтобы запретить втекание вещества в область моделирования через эту границу. На противоположной границе $(z = {{Z}_{{\max }}},$ $0 \leqslant r \leqslant {{R}_{{\max }}})$, были заданы “свободные” граничные условия $\partial {\text{/}}\partial n = 0$, чтобы позволить веществу вытекать из области. Втекание вещества через эту границу в область моделирования также запрещено. На цилиндрической границе $({{Z}_{{\min }}} \leqslant z \leqslant {{Z}_{{\max }}},$ $r = {{R}_{{\max }}})$ также были заданы свободные граничные условия $\partial {\text{/}}\partial n = 0$. Для изучения влияния граничных условий была проведена серия расчетов с различными размерами области моделирования.

Для расчетов с полоидальным магнитным полем использовался векторный потенциал магнитного поля $A = {{A}_{\phi }}$, который был определен и задан во всей области расчета. Для однородного магнитного поля $B$, сонаправленного с осью $z$, векторный потенциал записывается следующим образом: ${{A}_{\phi }} = rB$, ${{A}_{r}} = 0$, ${{A}_{z}} = 0$, тороидальный компонент поля ${{B}_{\phi }} = 0$. Для расчетов с тороидальным магнитным полем ${{B}_{\phi }}$ было задано по формуле поля проводника с током, направленного вдоль оси $Z$, как схематично показано на рис. 1 и 2. Таким образом, ${{A}_{\phi }} = 0$ во всей области, для $0 \leqslant r \leqslant {{R}_{d}}$ ${{B}_{\phi }} = {{B}_{0}}r{\text{/}}{{R}_{d}}$, для ${{R}_{d}} \leqslant r \leqslant {{R}_{{\max }}}$ ${{B}_{\phi }} = {{B}_{0}}{\kern 1pt} {{R}_{d}}{\text{/}}r$, где ${{B}_{0}}$ – постоянная величина. При этом автоматически гарантировано, что поле на оси $Z$ ${{B}_{\phi }} = 0$.

В начальный момент времени $t = 0$ магнитное поле в области представляло собой “чистое” тороидальное поле. Начальная плотность вещества во всей области моделирования, за исключением мишени, была однородной и имела малое значение: $\rho = {{\rho }_{0}}$. Начальная скорость вещества во всей области моделирования, за исключением мишени, ${{v}_{z}} = {{v}_{r}} = 0$. Вещество не имело углового момента, ${{v}_{\phi }} = 0$.

В области мишени $0 \leqslant R \leqslant {{R}_{d}}$, $0 \leqslant z \leqslant {{Z}_{d}}$ было задано значение плотности $\rho \gg {{\rho }_{0}}$ и сверхзвуковая скорость вылетающего из мишени вещества. Большинство расчетов было проведено для $M = 3$, $M = 6$. Были исследованы случаи, когда вылетающее вещество не несет магнитное поле, и когда вылетаюшее вещество несет тороидальное магнитное поле, по величине и структуре аналогичное магнитному полю в области моделирования.

2.5. Джет в гидродинамическом случае

В качестве тестирования МГД модели традиционно использовалось гидродинамическое приближение, когда было выключено магнитное поле, и положено $B = 0$ во всей области. Это позволило оценить корректность выбора размеров области моделирования, гидродинамических параметров и граничных условий. В тестовом случае начальная плотность вещества во всей области моделирования была однородна и имела малое значение: $\rho = {{\rho }_{0}}$. Начальная скорость вещества во всей области моделирования, за исключением мишени, была нулевой ${{v}_{z}} = {{v}_{r}} = 0$. Вещество не имело углового момента, ${{v}_{\phi }} = 0$.

В области мишени $0 \leqslant R \leqslant {{R}_{d}}$, $0 \leqslant z \leqslant {{Z}_{d}}$ было задано значение плотности $\rho \approx 300{\kern 1pt} {{\rho }_{0}}$ и сверхзвуковая скорость вылетающих из мишени частиц $v = M{\kern 1pt} {{c}_{s}}$, ${{c}_{s}} = \gamma \frac{{{{P}_{0}}}}{{{{\rho }_{0}}}}$ – скорость звука в невозмущенной области, $M$ – число Маха. Для тестового гидродинамического случая $M = 3$. Для простоты тестовых расчетов скорость вылетающего вещества в начальный момент времени имела только $z$-компоненту ${{v}_{z}}$; ${{v}_{r}} = 0$.

На рис. 3 показано начальное распределение плотности и скорости вещества, а также результаты моделирования выброса вещества в гидродинамическом случае. Цветной заливкой обозначен логарифм плотности вещества, стрелками показаны векторы скоростей, длина стрелки пропорциональна модулю скорости. Под воздействием тепловой энергии лазера вещество мишени начинает лететь вдоль оси $Z$ в направлении трекового детектора, попутно расширяясь в направлении оси $R$. Образуется ударная волна, которая движется к границе камеры, оставляя за собой небольшое устойчивое истечение конусовидной формы, параметры которого сохраняются с течением времени. Вблизи мишени наблюдается небольшая коллимация вещества, обусловленная наличием фонового газа. Остальное вещество рассеивается в области моделирования.

Форму истечения также можно увидеть на рис. 4, на котором показано распределение плотности вдоль оси $R$ на разных расстояниях $z$ от поверхности мишени в момент времени $t = 20$. Плотность потока вблизи мишени, $z = 0.1$ максимальна у оси $Z$ и распределена в узком конусе, ей соответствует тонкая сплошная линия. На расстоянии $z = 1.2$ распределение становится более широким, но еще сохраняет форму конуса, ей соответствует штрихпунктирная линия. Пунктирная линия соответствует плотности в середине камеры, $z = 2.5$. Вблизи детектора, $z = 5.0$ распределение плотности становится практически равномерным, что показывает штриховая линия.

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

2.6. Джет в тороидальном магнитном поле

Проведенные тестовые расчеты для гидродинамического случая показали корректность выбора размеров области моделирования, гидродинамических параметров и граничных условий. На следующем этапе моделирования было включено тороидальное магнитное поле ${{B}_{\phi }}$, которое задавалось во всей области следующим образом: при $0 \leqslant r \leqslant {{R}_{d}}$ ${{B}_{\phi }} = {{B}_{0}}{\kern 1pt} r{\text{/}}{{R}_{d}}$, при ${{R}_{d}} \leqslant r \leqslant {{R}_{{\max }}}$ ${{B}_{\phi }} = {{B}_{0}}{\kern 1pt} {{R}_{d}}{\text{/}}r$, где ${{B}_{0}}$ – постоянная величина. Были проведены расчеты для разных значений параметра $\beta = \frac{{8\pi {{P}_{\infty }}}}{{B_{0}^{2}}}$ и числа Маха $M$. Параметр магнитной вязкости был одинаковым во всех расчетах ${{\tilde {\eta }}_{M}} = \frac{1}{{{\text{R}}{{e}_{M}}}} = {{10}^{{ - 5}}}$. Остальные параметры соответствовали гидродинамическому случаю, т.е. начальная плотность вещества во всей области моделирования однородна и имеет малое значение: $\rho = {{\rho }_{0}}$. Начальная скорость вещества во всей области моделирования, за исключением мишени, нулевая ${{v}_{z}} = {{v}_{r}} = 0$. Вещество не имеет углового момента, ${{v}_{\phi }} = 0$.

В области мишени $0 \leqslant R \leqslant {{R}_{d}}$, $0 \leqslant z \leqslant {{Z}_{d}}$ задано значение плотности $\rho \approx 300{\kern 1pt} {{\rho }_{0}}$ и сверхзвуковая скорость вылетающих из мишени частиц $v = M{{c}_{s}}$, где ${{c}_{s}} = \gamma \frac{{{{P}_{0}}}}{{{{\rho }_{0}}}}$ – скорость звука в невозмущенной области, $M$ – число Маха. Скорость втекающего в расчетную область через границу мишени вещества в начальный момент времени имеет только $z$-компонету ${{v}_{z}}$; ${{v}_{r}} = 0$. Магнитное поле во втекающем веществе задано по такой же формуле, что и начальное поле в области: для $0 \leqslant r \leqslant {{R}_{d}}$ ${{B}_{\phi }} = {{B}_{0}}{\kern 1pt} r{\text{/}}{{R}_{d}}$, для ${{R}_{d}} \leqslant r \leqslant {{R}_{{\max }}}$ ${{B}_{\phi }} = {{B}_{0}}{\kern 1pt} {{R}_{d}}{\text{/}}r$, где ${{B}_{0}}$ – постоянная величина.

2.6.1. Тороидальное поле $\beta = {{10}^{{ - 4}}}$, $M = 3$. Результаты моделирования выброса вещества для параметра $\beta = {{10}^{{ - 4}}}$ представлены на рис. 5. На левой панели показана эволюция течения вещества. Цветная заливка соответствует логарифму плотности вещества, стрелками обозначены векторы скоростей. Можно видеть, что вещество мишени под воздействием энергии лазерного излучения начинает лететь вдоль оси $Z$. Сначала течение немного расширяется во всех направлениях, но затем коллимируется магнитным полем. Образуется ударная волна вытянутой формы, которая движется от мишени к границе камеры, оставляя за собой устойчивое истечение. Влияние внешнего газа на коллимацию потока в данном случае незначительно.

Форму течения также можно увидеть на рис. 6, где показано распределение плотности вдоль оси $R$ на разных расстояниях от поверхности мишени в момент времени $t = 12$. Тонкая сплошная линия соответствует плотности потока у границы мишени $z = 0.1$. Штрихпунктирная линия соответствует расстоянию $z = 1.2$, пунктирная соответствует плотности по середине камеры, $z = 2.5$, штриховая линия соответствует плотности вблизи границы камеры, $z = 5.0$. Можно видеть, что распределение плотности в потоке плазмы зависит от расстоянии от мишени. Плотность потока вблизи мишени максимальна у оси $z$ и распределена в узком конусе. По мере движения вещества к границе камеры поток сохраняет форму конуса, умеренно расширяясь к границе. С помощью рис. 5 и 6 можно оценить угол отклонения. Он приблизительно равен $ \approx {\kern 1pt} 11^\circ {\kern 1pt} - {\kern 1pt} 13^\circ $. И, соответственно, угол раствора конуса $\theta \approx 22^\circ {\kern 1pt} - {\kern 1pt} 26^\circ $. Таким образом, на трековом детекторе можно будет наблюдать четкое сплошное пятно джета.

По правой панели рис. 5 можно проследить эволюцию магнитного поля в области моделирования. Цветная заливка соответствует логарифму ${{B}_{\phi }}$. В начальный момент времени поле в области увеличивается при значениях $0 \leqslant r \leqslant {{R}_{d}}$, а затем, при ${{R}_{d}} \leqslant r \leqslant {{R}_{{\max }}}$ спадает к границе области. После выброса вещества первоначальное магнитное поле сначала деформируется, но остается достаточно сильным для коллимации вещества. По мере формирования джета поле начинает восстанавливаться. Начальное распределение индукции магнитного поля ${{B}_{\phi }}$ вдоль оси $R$ показано сплошной линией на рис. 7. Его структура в начальный момент времени ${{t}_{0}} = 0.0$ полностью соответствует полю проводника с током, направленного вдоль оси $Z$, схематично представленного на рис. 2. В момент времени ${{t}_{1}} = 1.2$ поле гораздо меньше по абсолютному значению и слабо спадает вдоль $R$ к границе области (пунктирная линия). В момент времени ${{t}_{2}} = 12.0$ структура ${{B}_{\phi }}$ снова возвращается к полю проводника с током, но с меньшими абсолютными значениями (штриховая линия). Таким образом, тороидальное магнитное поле сохраняется в области на протяжение всего расчета и играет главную роль в коллимации потока.

2.6.2. Тороидальное поле $\beta = 1{{0}^{{ - {\mathbf{3}}}}}$, $M{\mathbf{ = 3}}$. В следующем расчете мы уменьшили магнитное поле. На рис. 8 представлены результаты моделирования для случая с параметрами $\beta = {{10}^{{ - 3}}}$, $M = 3$. Левая панель отображает эволюцию течения вещества, цветной заливкой показан логарифм плотности. Вещество мишени под воздействием тепловой энергии лазера летит вдоль оси $Z$ в направлении детектора. Сначала течение расширяется во всех направлениях, но затем коллимируется магнитным полем. Образуется конусообразная вытянутая ударная волна, которая движется от мишени к границе камеры, оставляя за собой устойчивое истечение. Влияние внешнего газа на коллимацию потока незначительно.

Рисунок 9 показывает распределение плотности вдоль оси $R$ на разных расстояниях от поверхности мишени в момент времени $t = 6$, по которому можно определить форму течения. Тонкая сплошная линия соответствует плотности потока у границы мишени $z = 0.1$. Штрих-пунктир соответствует расстоянию $z = 1.2$, пунктир соответствует плотности по середине камеры, $z = 2.5$, штриховая линия соответствует плотности вблизи границы камеры, $z = 5.0$. Видно, что распределение плотности в потоке плазмы зависит от расстоянии до мишени. Плотность потока вблизи мишени максимальна у оси $z$ и распределена в узком конусе. Если в середине камеры поток сохраняет форму конуса, умеренно расширяясь с увеличением $z$, то у границы $z = 5.0$ распределение плотности образует четкую кольцевую структуру. С помощью рис. 8 и 9 можно оценить угол отклонения. Он приблизительно равен $ \approx {\kern 1pt} 15^\circ {\kern 1pt} - {\kern 1pt} 18^\circ $. Соответственно, угол раствора конуса $\theta \approx 30^\circ {\kern 1pt} - {\kern 1pt} 36^\circ $. Таким образом, на трековом детекторе можно будет наблюдать джет в виде четкой кольцевой структуры.

Правая панель рис. 8 отражает эволюцию магнитного поля в области моделирования. Цветная заливка соответствует логарифму ${{B}_{\phi }}$. В начальный момент времени поле в области увеличивается при значениях $0 \leqslant r \leqslant {{R}_{d}}$, а затем, при ${{R}_{d}} \leqslant r \leqslant {{R}_{{\max }}}$ спадает к границе области. После выброса вещества первоначальное магнитное поле деформируется, но остается достаточно сильным для коллимации вещества. По мере формирования джета поле так же начинает восстанавливаться, но остается неоднородным вдоль $Z$.

На рис. 10 показаны начальное распределение индукции магнитного поля ${{B}_{\phi }}$ вдоль оси $R$ и его эволюция во времени. Структура поля в начальный момент времени ${{t}_{0}} = 0.0$ полностью соответствует полю проводника с током, направленного вдоль оси $Z$, как схематично изображено на рис. 2. В момент времени ${{t}_{1}} = 1.5$ уменьшилось по абсолютному значению и слабо спадает вдоль $R$ к границе области (штриховая линия). В момент времени ${{t}_{2}} = 6.0$ структура ${{B}_{\phi }}$ снова возвращается к полю проводника с током, но с меньшими абсолютными значениями (штрихпунктирная линия). Таким образом, и в этом случае тороидальное магнитное поле сохраняется в области на протяжении всего расчета и играет главную роль в коллимации потока.