CC BY
28
УДК 533.9.01
И.М. Данилов, В.В. Демченко
Московский физико-технический институт (государственный университет)
Математическое моделирование волны термоядерного горения, инициированной сферическим игнитором в мишенях инерциального термоядерного синтеза
Представлены результаты работы программного пакета, который моделирует распространение волны термоядерного горения по DT-мишени. Проведено сравнение результатов моделирования модельной задачи с аналитическим решением.
Ключевые слова: контролируемый термоядерный синтез, численное моделирование, сферический игнитор.
I. Введение
Управляемый термоядерный синтез (УТС) привлекателен многими аспектами, но один из самых важных для практического применения — это неисчерпаемость термоядерного топлива на нашей планете. В качестве топлива можно использовать тяжелые изотопы водорода, например, дейтерий.
В [1] сделаны оценки минимальной энергии лазерного импульса с учётом тепловыделения. Из этих оценок следует, что для термоядерного синтеза (ТС) выгоднее всего использовать DT-плазму, сжатую до высоких плотностей, многократно превышающих плотность твёрдого водорода.
В статье [2] был сформулирован подход быстрого зажигания (fast ignition) мишени. Вещество мишени сначала сжимается длинным лазерным импульсом до большой плотности при относительно малой температуре. Затем небольшая часть сжатого вещества быстро нагревается до необходимой для зажигания температуры. В рамках концепции fast ignition зажигание предполагалось осуществить на конечной стадии сжатия с помощью релятивистских электронов, ускоренных пе-таваттным лазерным импульсом (третьим по счету) в плотном веществе внутри мишени, куда лазерное излучение должно было проникнуть, распространяясь внутри канала, созданного вторым лазерным импульсом.
В данной работе расоматривается двухстадийный механизм инициирования волны термоядерного (ТЯ) горения в мишенях лазерного термоядерного синтеза (ЛТС). Этот процесс протекает следующим образом: на первом этапе формируется область высокоплотной, но относительно холодной плазмы, на втором этапе внутри этой области при помощи дополнительного энергетического импульса образуется высокотемпературная подобласть, называемая игнитором. Для реализации второго этапа — создание игнитора — существует несколько вариантов: во-первых, при помощи потока быстрых электронов, образующих-
ся при воздействии на вещество лазерного излучения; во-вторых, используя пучки многозарядных ионов, или с помощью протонов, ускоренных лазерным импульсом релятивистской интенсивности в специальной мишени, расположенной в непосредственной близости от ТЯ мишени. При реализации любого из перечисленных методов быстрого поджига возникает одна существенная проблема: форма игнитора и основной плазмы мишени перестает быть сферически симметричной, что может сказаться на эффективности выгорания ТЯ топлива, в частности, увеличить критическую энергию игнитора, необходимую для поджига мишени. Приведённая выше модель инициирования позволяет существенно снизить требования к профилированию и стабильности основного лазерного импульса, при этом коэффициент усиления мишени может не претерпевать особых изменений.
При исследовании процесса ТЯ горения мишени определённый интерес представляют физические явления, происходящие на начальной стадии в высокотемпературной плотной плазме.
Представленная работа связана с разработкой и тестированием пакета программ, моделирующего ТЯ процесс горения в сферически симметричном случае. В частности, изучаются физические процессы, протекающие в разные моменты времени как на самой ранней стадии, так и на более поздней. Осуществлено сравнение аналитического решения с результатами численного моделирования, исследована сходимость численных решений при измельчении сетки.
II. Математическая модель
Рассматривается осесимметричная задача ТЯ горения плотной высокотемпературной дейтери-ево-тритиевой плазмы. Используется гидродинамическое одножидкостное двухтемпературное приближение, включающее процессы электронной и ионной теплопроводности, электрон-ионной столкновительной релаксации и выделение энергии за счёт ТЯ реакций.
Для описания ТЯ горения плазмы в переменных Эйлера (г; в) сферической системы координат выбрана следующая математическая модель, которая имеет следующий вид:
др 1 d(r2pU)
о, ~п о
dt
dr
O.
dpU | 1 d[r2{pU2+P)\ |
—77 ~ h
1 9 (sin 9pV)
r sin в дв
1 d(sin вpUV)
де
dr
dpV | 1 d(r2p\JV) |
—77 ~
r sin в дв 2 P r
1 d[sin в(pV2 + P)]
(І)
(2)
dt
dr
r sin в cos 0P r sin в
дв
(З)
дрЕ | 1 9[г2£/(/эР + Р)] ,
^ ~ h
dt
dr
sin в
дв
+■
dr
+
1 %2(We}r] , 1 d\sm.6{We)e\
r sin в дв
1 d[sin в(Wi )e ]
+
1%Щ|
r2 dr
dpsj 1 d{r~U p£j)
dt r2 dr
+«i4^ +
r sin в дв
1 d(sin вV pєi)
(4)
r sin в дв 1 d(sin вV)
+
r2 dr
r sin в
дв
-n (T 7M I 1 dp(^}r] 1 c)[sin^(Wi)e]
Че*1 e r2 dr rsintf 90 ’
(5)
где p — плотность; U — радиальная г-компонента скорости; V — тангенциальная 0-компонента скорости; £i — ионная удельная внутренняя энергия; ее — электронная удельная внутренняя энергия; E = ee + ei + (U2 + V2 )/2 — полная удельная энергия; P = (y — 1)p(ee + £i) — давление; 7 = 5/3 — показатель адиабаты; Pi = (7 — 1)pei — ионное давление; We,i — тепловые потоки за счёт электронной и ионной теплопроводности:
We = ке grad Te,
Wi = ki grad Ti,
где Te,Ti — соответственно электронная и ионная температуры;
k1
кє = \J k\ + Щ;
9.S3 ■ 1O19 ■ T5j\ эрг
k2
In Aej ^ см • с • кэВ 2,54716 ■ 1O17 ■ p ■ Te эрг
1пЛ^
ln Леі = 0,5ln[1 + 3,59 ■ 1O
см с кэВ
4 T2 + S.S62 ■ 1O~5 ■ p4/3
1пЛґ = ln[1 + 5.27 ■ (O.39 + 51 ■ Ti/(p)1/3]; ki = 2.51967 ■ 1O18 ■ (Ti)5/2/lnЛі[
эрг
см ■ с ■ кэВ lnЛі = O.5ln[1 + 4.3993 ■ 1O4 ■ Tf/p];
Qei — коэффициент электронно-ионного обмена: Qei = 2.S15 ■ 1O24x
p2 ■ ln[l + 3,59 • 104 • (T2 + 8,862 • lO^5 • p4/3)/p] Te12 + 1,262 • IQ-3/?
[-
эрг
‘см3 • с • кэВ ^т.я — энерговыделение за счёт ТЯ реакций:
дТЯ = 8,483 • 1029х
2 (1 + 0,232 •^/4)ехр(-20/^1/3) эрг Р / I------------------------т~ I- я -I’
туз _ у/1 + 9141.10-5.^13/4 СМ • с
Тг] - кэВ; [Те] - кэВ; [р] - г/см3.
Исходная система уравнений включает уравнение неразрывности вещества (1), законы количества движения (2)-(3) (для радиальной и тангенциальной компоненты импульса соответственно), законы сохранения для полной энергии (4) и изменение ионной составляющей плазмы (5).
В физической модели предполагается, что энергия а-частиц локально передаётся электронной составляющей плазмы, а нейтроны не взаимодействуют с веществом мишени и свободно её покидают.
Используемая на первом полушаге схема расчёта предполагает, что гидродинамическая часть уравнения, имеющая гиперболический тип, аппроксимируется явно. Рассматривается двухмерный случай и используется явная схема первого порядка точности, для которой выполняется условие Куранта-Фридрихса-Леви. Вместе с гидродинамикой учитывается ТЯ выход. В результате мы получаем окончательное значение плотности, компоненты вектора скорости и предварительные значения электронной и ионной температур. На следующем полушаге производится учёт членов с электронной и ионной теплопроводностью по неявной схеме (схема пространственного покомпонентного расщепления с итерациями по нелинейным коэффицентам электронной и ионной теплопроводности). В итоге, учитывая граничные условия для электронной и ионной температур, получаем краевую разностную задачу относительно них, которая решается итерациями с использованием одномерных прогонок по лучам. На третьем полушаге учитывается электронно-ионная релаксация.
III. Постановка задачи
В данной работе рассматривается ТЯ горение специальной сферически семметричной мишени. В ее центральной части, состоящей из плотной
r
r
r
r
].
плазмы, создают область повышенной температуры (см. введение). На рис. 1 представлена структура расчётной области. Игнитор — это сфера наименьшего радиуса 35 мкм, в которой находится горячая и плотная плазма. Следующая область, имеющая форму шарового слоя с внутренним радиусом 35 мкм и внешним 300 мкм, — плотная и холодная плазма, занимающая большую часть мишени. Пространство между границами мишени и границей расчётной области заполнено относительно холодной и разреженной плазмой.
Рис. 1
Воспользуемся аналитическим решением вспомогательной задачи. Аналитическое решение полной задачи, состоящей из системы нелинейных уравнений, в настоящее время получить затруднительно, поэтому отладка программы проводится на модельной задаче, для которой существует аналитическое решение.
IV. Сравнение аналитического и численных решений на модельной задаче
В качестве модельной задачи для отладки первого этапа разностной схемы, на котором учитывается гидродинамика, выбран случай гидродинамического распада разрыва с образованием волны сжатия и волны разрежения.
В начальный момент времени наблюдаем разрыв на расстоянии 35 мкм от центра мишени — граница игнитора. В игниторе температура 12 кэВ, за его границей и до конца мишени — 1 кэВ. Мы наблюдаем распад разрыва с образованием ударной волны и волны сжатия (рис. 2).
На ударной волне должны выполняться три основных закона сохранения: массы, количества движения и энергии. Запишем соотношения Рэн-кина-Гюгонио, отражающие эти законы:
Р\(и\ — До) = ро(ио — До), Р\ + Р\(и\ — Б,)2 = Ро + ро(и, — До)2,
(и і — Во){рі[є\ +
(Ц. — До)2
— (и — До){ро[єо +
(и — До)2
]+рі}
]+Ро|.
На контактной границе скорость и давление связаны следующими соотношениями:
V = П2, Р1 = Р2.
В волне разрежения сохраняются значения инварианта на характеристике V + С:
и3 + ^К = и,+ ^
73 — 1
73 — 1’
где Сз и С2 — скорость звука в игниторе и в области постоянного течения между волной разрежения и контактной границей соответственно.
При помощи этих соотношений строим характеристики и + С в волне разрежения.
Закон сохранения массы выполняется для сеток с разным количеством узлов с различной степенью точности. Радиус расчётной области в модельной задаче составляет 60 мкм. Масса, распределенная по сфере этого радиуса:
V = 4 • п • 0,63/3 = 0,90478 • 10~6 см3,
т = V • 100 г/см3 = 90,478 • 10~6 г.
Численное моделирование для равномерных сеток с различными шагами точности дало следующее распределение массы: 500 узлов на каждый луч, всего 4 луча — 78,78651 • 10~6 г, 1000 узлов на каждый луч, всего 8 лучей — 87,25661 • 10~6 г, 2000 узлов на каждый луч, всего 16 лучей — 89,31799 • 10~6 г, 4000 узлов на каждый луч, всего 32 луча — 89,99397 • 10~6 г. В результате математических преобразований [6] получаем уравнение для г = Р1 /Р3:
X2^2п - а0у2ХЯг‘+2 + 2а^(р + v)XZn+1-
-[2 + (р + V)2а0]ХЯп - v2Z2 + 2^(^ + у)2-
-(р + V )2 + 1 = 0.
После подстановки числовых значений:
5
73=70=3, А‘ = °, а0 = 4, X2 = 144,
п = 5, V =60 получаем алгебраическое уравнеие 10-й степени:
144z10—2880z7+5760z6—2904z5—60z2+120z—59 = 0.
(6)
При решении уравнения (6) получаем 10 корней, из которых 8 корней являются комплексными, а из оставшихся один корень положительный. Найденное решение:
Р1/Р3 = г = 0,8707787311,
где Р3 — давление в игниторе, а Р1 = 4,639-1017 дин/см2 давление в области постоянного течения между контактным разрывом и ударной волной. Затем находим скорость распространения ударной волны:
Д0 = 8,027 • 107 см/с,
скорости перемещения контактного разрыва и плотности в областях постоянного течения:
и1 = 4,818 • 107 см/с,
V! = П2,
Р1 = 250,181 г/см3, р2 = 66,027 г/см3. результатов представлено на
Сравнение рис. 3, 4, 5.
У. *10кэВ У. ‘10|4Д1Н,см
— Давление(501) — Давление(1001)
— Давление(2001) — Давление(4001)
— Давление(4001 с учетом) Аналитическое решение
— давление (8000 с учетом на области 1.2млм)
Рис. 4
Т=10пс
X. *100мкм
— Э.т.(500) — Э.т.(ЮОО)
— Э.т.(2000) — Э.т.(4000)
— Э.т. (4000 с учетом) — Э.т. (8000 с учетом на области 1.2млм)
-—аналитическое решение
На рис. 3 показано распределение плотности по радиальному направлению. «Плотность (501)», «плотность (1001)», «плотность (2001)», «плотность (4001)» — результаты численного моделирования модельной гидродинамической задачи на области 60 мкм для четырёх сеток с разным шагом. Самый нижний график — результат моделирования при использовании сетки с 500 узлами, для следующего графика использовалась сетка с 1000 узлами, затем представлен график для сетки с 2000 узлов, график, больше всех приближающийся к аналитическому решению, получен на сетке с количеством узлов 4000. Также на рис. 3 показано распределение плотности для ре-шеня полной задачи с учётом ТЯ выхода с применением сетки в 4000 узлов на области 60 мкм — «плотность (4001 с учётом)». «Плотность (8000 с учётом на области 1,2 мм)» — график плотности для задачи с учётом ТЯ выхода, решаемой на области 1,2 мм, число узлов применяемой сетки — 8000. «Аналитическое решение» — результаты аналитического решения, полученные в п. 4, где от 0 • 100 мкм до 0,225 • 100 мкм — невозмущённая область, в точке абсциссы 0, 225 • 100 мкм голова волны разрежения, в 0, 29 • 100 мкм — её хвост, [0,29; 0,395) • 100 мкм — область постоянного течения, в точке 0,395 • 100 мкм — контактный разрыв, далее следует область постоянного течения (0,395; 0,43) • 100 мкм, 0, 43 • 100 мкм — ударная волна, распространяющаяся вправо.
На рис. 4 нарисовано распределение давления по радиальному направлению. Здесь представлены четыре результата численного моделирования модельной гидродинамической задачи на области 60 мкм — «Давление (501)», «Давление (1001)», «Давление (2001)», «Давление (4001)». Самый нижний график — результат моделирования при использовании сетки с 500 узлами, для следующего графика использовалась сетка с 1000 узлами, далее — сетка с 2000 узлов, график, больше всех приближающийся к аналитическому решению, получен для сетки с числом узлов 4000. При сравнении данных, полученных при использовании различных сеток, выяснилось, что точность решения (обусловленная его приближением к аналитическому) возрастает в два раза при увеличении количества узлов в два раза. Также на рис. 4 показано распределение давления для ре-шеня полной задачи с учётом ТЯ выхода с применением сетки в 4000 узлов на области 60 мкм — «Давление (4001 с учётом)». «давление (8000 с учётом на области 1,2 млм)» — график давления для задачи с учётом ТЯ выхода, решаемой на области 1,2 млм, число узлов применяемой сетки — 8000. «Аналитическо решение» — результаты аналитического решения, полученые в п. 4, где от 0 • 100 мкм до 0,225 • 100 мкм — невозмущённая область, в точке абсциссы 0,225 • 100 мкм голова волны разрежения, в 0,29 • 100 мкм — её хвост, [0, 29; 0,395) • 100 мкм — область постоянного течения, в точке 0, 395 • 100 мкм — контактный раз-
рыв, далее следует область постоянного течения (0,395; 0,43) • 100 мкм, 0,43 • 100 мкм — ударная волна, распространяющаяся вправо.
На рис. 5 изображено распределение электронной температуры по радиальному направлению. Здесь представлены четыре результата численного моделирования модельной гидродинамической задачи на области 60 мкм — «500», «1000», «2000», «4000». Соответствие графиков и используемых сеток такое же, как и на предыдущих рисунках. При сравнении данных, полученных при использовании различных сеток, выяснилось, что точность решения (обусловленная его приближением к аналитическому решению) возрастает в два раза при увеличении количества узлов в два раза. «4001 (с учётом)» — распределение давления для решеня полной задачи с учётом ТЯ выхода с применением сетки в 4000 узлов на области 60 мкм. «8000 (с учётом на области 1,2 мм)» — график температуры для задачи с учётом ТЯ выхода, решаемой на области 1,2 млм, число узлов сетки — 8000. «Аналитическое решение» — результаты аналитического решения, полученные в п. 4.
Из графиков видно, что схема действительно первого порядка точности, так как при двукратном увеличении количества узлов на той же расчётной области точность возрастает в два раза.
V. Результаты численного моделирования полной задачи
При развитии ТЯ горения в мишени создаются необходимые условия для возникновения детонационной волны. Во-первых, это происходит за счёт повышения плотности за фронтом ударной волны. Во-вторых, при торможении а-частиц происходит разогрев электронной компоненты плазмы (кулоновское взаимодействие), а за счёт эле-тронно-ионной релаксации происходит нагрев ионной компоненты.
На 17,5 пс тепловая волна входит в волну сжатия (в область постоянного течения за фронтом ударной волны) и начинает разогревать плазму за фронтом ударной волны, усилилась ТЯ реакция, в результате чего повышается давление, сформировалась ударная волна — это подтверждает максимальное энерговыделение (рис. 6). Волна сжатия сильно ослабляется тепловой волной, опережающей её.
На 22,5 пс наблюдаем два максимума у графика плотности: первый максимум — это результат действия ослабленной волны сжатия (максимум едва сместился в направлении от центра иг-нитора), его значение уменьшилось; второй максимум — результат действия детонационной волны. Наблюдаем начало процесса переноса вещества по направлению к центру игнитора — появился локальный максимум у графика скорости, этот максимум находится под первым максимумом плотности (рис. 7).
Рис. 6. Представлены графики на момент времени 17,5 пс. «Плотность» — по оси ординат плотность в единицах г/см . «Скорость» — радиальная скорость распространения вещества в 10 см/с. «EQ*1 000 000 (энерговыделение ТЯ)» — энерговыделение ТЯ горения, умноженное на 106 для соблюдения масштаба, в 1014 кэВ/см3. «Ионная температура» — температура ионов плазмы в кэВ. «Электронная температура» — температура электронов плазмы, в кэВ. «Давление /100» — давление в плазме, деленное на 100 в 1016 дин/см2
Рис. 7. Представлены графики на момент времени 22,5 пс. «Плотность» — по оси ординат отложена плотность в единицах г/см3. «Скорость» — радиальная скорость распространения вещества в 107 см/с. «EQ*1 000 000 (энерговыделение ТЯ)» — энерговыделение ТЯ горения, умноженное на 106 для соблюдения масштаба, в 10 кэВ/см . «Ионная температура» — температура ионов плазмы в кэВ. «Электронная температура» — температура электронов плазмы, в кэВ. «Давление /100» — давление в плазме, деленное на 100 в 1016 дин/см2
Рис. 8. Нейтронный выход по прошествии 300 пс
На 37,5 пс детонационная волна сошлась в центре игнитора, её энергия перешла в кинетическую энергию ударной волны и началось сильное удельное энерговыделение — от центра начала распространяться новая ударная волна. Этот процесс из-за малых объёмов в центре мишени мало влияет на ход ТЯ реакции. Реакция поддерживается за счёт повышенного давления за фронтом детонационной волны и ионной теплопроводимости. На рис. 8 показан выход ТЯ реакции.
VI. Заключение
1. В результате сравнения аналитического и численного решений модельной задачи, проведённого в п. 4, было установлено, что:
1) схема имеет заявленный первый порядок аппроксимации;
2) утечка массы незначительна, особенно для сеток с количеством узлов на луч больше 2000;
3) в начале реакции происходит распад разрыва с образованием ударной волны и волны разрежения.
2. Анализ результатов численного моделирования полной задачи показал, что после распада раз-
рыва происходит зарождение и протекание само-поддерживающейся ТЯ реакции.
Литература
1. Прохоров А.М., Анисимов С.И., Паши-нин П.П. Лазерный термоядерный синтез // УФН. — 1976. — Т. 119. — С. 401-422.
2. Tabak M., Hammer J., Yanovsky V.P. [et al.] Fast Ignition // Phys. Plasmas. — 1994. — V. 1. — P. 1626-1673.
3. Бессараб А.В., Гаранин С.Г., Кирилов Г.А., Кочемасов Г.Г., Муругов И.М., Мхитарьян Л.С., Пунин В.Т, Сухарев С.А. Эксперименты на установках «Искра 4» и «Искра 5», развитие диагностик и результатов последних двух лет // Физика плазмы. — 1998. — Т. 24. — С. 197-225.
4. Зельдович Я.Б., Райзер Ю.П. Физика ударных волн и высокотемпературных гидродинамических явлений. — М.: Наука, 1966.
5. Демченко В.В. Распад произвольного гидродинамического разрыва: учебное пособие. — М.: МФТИ, 1998.
6. Демченко В.В. Нелинейные уравнения и системы уравнений: учебное пособие. — М.: МФТИ, 2003.
7. Белоцерковский О.М., Демченко В.В., Косарев В.И., Холодов А.С. Численное моделирование некоторых задач лазерного сжатия оболочек // ЖВМиМФ. — 1978. — Т. 18, № 2. — С. 420-444.
8. Самарский А.А., Губинская Н.Б., Гуськов С.Ю., Демченко В.В., Змитренко Н.В., Ильин Д.В., Левковский А. А, Розанов В.Б., Шерман В.Е. Математическое моделирование волны термоядерного горения при ИТС, инициированной несферическим смещённым относительно центра мишени игнитором // Исследовано в России. — 2006. — С. 1137-1147.
9. Андреев А.А., Гуськов С.Ю., Ильин Д.В. [и др.] Двустадийный процесс инициирования волны термоядерного горения в DT-топливе // ЖЭТФ. — 2003. — Т. 123. — С. 783-787.
10. Basko M.M. // Plasma Phys. — 2003. — V. 45. — P. 125-132.
11. Басов Н.Г., Воловски Е., Денус С., Заха-ренков Ю.А., Карнаухов А. А, Мурз В., Склизков Г.В., Фарны Ю., Шиканов А.С. Диагностика нейтральной компоненты плазмы микросфер, сжимаемых лазером // Физика плазмы. — 1988. — Т. 14. — С. 77-111.
12. Анисисмов С.И., Прохоров А.М., Фортов В.Е. Применение мощных лазеров для иссле-
дования вещества при сверхвысоких давлениях // УФН. — 1984. — Т. 142, в. 3. — С. 395-431.
13. Басов Н.Г, Розанов В.Б. [и др.\. О диагностике плазмы в условиях ЛТС по продуктам вторичных термоядерных реакций // Физика плазмы. — 1986. — Т. 12. — С. 916-1117.
14. Выговский О.Б., Гуськов С.Ю., Ильин Д.
B, Левковский А. А, Розанов В.Б., Шерман В.Е. Влияние ядер отдачи на скорость термоядерной реакции дейтерий-тритиевой плазмы // Физика плазмы. — 1990. — Т. 16. — С. 1509-1512.
15. Потапов А.В., Лобода Е.А., Пронин В.А., Политов В.Ю., Антонова Л.В., Липин А.В., Как-шин А.Г., Чефонов О.В. Спектры рентгеновского излучения лазерных мишеней в экспериментах на установке Сокол П // ФП. — 2007. — Т. 33. —
C. 3--15.
16. Торризи Л., Гаммини С., Селона Л., Краса И., Ласка Л., Воловски И. Лазерная генерация плазмы на установке в Национальном институте ядерной физики // Физика плазмы. — 2006. — Т. 32. — С. 558--564.
Поступила в редакцию 05.03.2008.
