Диссертация посвящена математическому моделированию динамики процессов высокотемпературной релаксации в веществе при облучении мощными ионными пучками. В работе рассматривается ряд актуальных применений системы нелинейных дифференциальных уравнений в частных производных вида.
DE р-^ = -PeV • V + V • (KeVTe) ± д, {Те — ТО + qe (r, 0, (v.l) DEP.
Целью работы является развитие методов математического моделирования взаимодействий ускоренных ионов с веществом, разработка эффективных алгоритмов и программ для численного решения систем нелинейных дифференциальных уравнений в частных производных для описания процессов высокотемпературной релаксации, сопровождающей эти взаимодействия, построение математических моделей ряда важных с практической точки зрения процессов таких, как взаимодействие интенсивных ионных пучков с тонкими образцами, формирование треков в высокотемпературных сверхпроводниках, процессов схлопывания кавитационного пузырька в дейтерированном ацетоне.
За последние двадцать лет диапазон технологических применений источников концентрированных потоков энергии и, в частности, ионных ускорителей, существенно расширился [1]. В этой связи можно упомянуть различные типы термической обработки материалов, получение тонких пленок и покрытий, а также сверхчистых материалов, направленную модификацию физико-механических и химических свойств металлов и сплавов. Особое значение эти исследования приобретают в связи с дальнейшей перспективой миниатюризации компьютерных технологий, возникшей благодаря открывшейся в последнее время возможности получения капельных наноструктур (кластеров) проводящих материалов, а также создания в поверхностном слое полупроводников так называемых квантовых точек. В качестве отдаленной цели уже сформулирована и развивается концепция квантового компьютера, создание которого будет связано с развитием принципиально новых технологий, включающих, по-видимому, и различные методы обработки материалов ускоренными ионными.
Требования к ускорителям, предназначенным для генерации необходимых ионных пучков, определяются как поставленными технологическими задачами, так и особенностями взаимодействия ионов с материалами. В частности, для исключения активации мишеней необходимо ограничить энергию ионов величиной, находящейся ниже порога ядерных реакций. Для протонов, которые обычно всегда присутствуют в ионном пучке из-за различных загрязнений, это соответствует Е < 400 кэВ. Далее, важно также обеспечить существенную плотность выделяющейся энергии. Например, эксперименты показывают, что для того чтобы воздействие пучка на механические и антикоррозионные свойства поверхности металла стало заметным, необходимо, чтобы энерговклад пучка превышал энергию сублимации [2, 3], а для уменьшения влияния процессов тепловой релаксации важно добиться еще и достаточно высокой мощности поступающей энергии. Оценки показывают, что для практических целей обработки металлов часто бывает достаточно обеспечить длительность ионного импульса At ~ 10~7 с при энергии ионов Е ~ 300 кэВ и потоке энергии 1−10 Дж/см2. Кроме этого, могут предъявляться также определенные требования к однородности плотности потока по сечению пучка, высокой долговечности установки, хорошей повторяемости параметров пучка (обычно нестабильность не должна превышать 10—20%).
С точки зрения вышеперечисленных критериев весьма перспективным представляется использование источников ионов со взрывной ионной эмиссией [4]. В этом случае при энергии ионов Е ~ 100 — 300 кэВ и радиусе пучка г ~ Ю-2 см удается получить плотность выделившейся энергии в интервале Qo ~ 10~2 — 102 Дж/см2, что отвечает плотности выделяемой мощности Ро ~ Ю5 — 108 Вт/см2 при длительности импульса Ai ~ 10~8 -К)" 7 с.
Физические модели взаимодействия пучков ионов с веществом.
Картина взаимодействия мощных ионных пучков с мишенью весьма сложна. Она содержит большое разнообразие самых различных физических явлений: от тепловых и механических, до газодинамических и плазменных. Хотя многие особенности этих процессов уже ясны, в настоящее время общая физическая модель, описывающая в широком диапазоне параметров такие взаимодействия, пока отсутствует. Характерная черта воздействия мощного пучка ионов на металл состоит в том, что зона структурных изменений в этом случае более чем на два порядка превосходит длину свободного пробега ионов (эффект дальнодействия). Например, дислокационные структуры, образующиеся при ионном облучении, могут располагаться на расстояниях до 100 мкм от поверхности при средней длине свободного пробега ионов менее 1 мкм. Удовлетворительного объяснения этого явления в настоящее время не существует. В частности, гипотеза, согласно которой дефекты образуются вблизи поверхности, а затем мигрируют на большие глубины [5, 6], противоречит общепринятым представлениям о том, что дефекты не могут преодолевать границы зерен вещества. Учет механизмов образования в облучаемых материалах волн сжатия, переходящих в процессе нелинейной эволюции в ударные, представляется весьма привлекательным. При этом дислокации могут образовываться за счет энергии упругих волн, возбуждаемых в образце [7, 8, 9]. Однако и в этом случае имеются экспериментальные данные, которые, по-видимому, не укладываются в рамки такого описания. Мы имеем в виду образование дефектов на расстоянии 100 мкм от источника размером 10 нм, поскольку в этом случае расходящаяся сферическая акустическая волна должна быстро терять плотность энергии, необходимую для образования дефектов [10].
Развитие новых радиационных технологий обработки материалов, наблюдающееся в последние годы, сделало актуальной задачу разработки моделей расчета прохождения потоков ионов через вещество. Такие модели должны учитывать явления генерирования и эволюции тепловых и термомеханических полей, а также процессы фазовых превращений, приводящих как к полезной модификации приповерхностного слоя, так и к нежелательному разрушению облучаемой поверхности. В общем случае эти процессы описываются сложной системой уравнений сплошной среды, решаемых на современных ЭВМ.
При низких интенсивностях пучка, когда средняя тепловая энергия е ~ кТ молекул среды существенно меньше теплоты испарения АН, вылет вещества из мишени может осуществляться только за счет высокоэнергетических хвостов распределений — реализуется тепловой активационный механизм испарения. Скорость распространения фронта испарения V в таком режиме описывается формулой Аррениуса [11]:
V ~ - ехр{-АН/кТ). т.
Здесь параметры, а — постоянная решетки и г — период собственных колебаний молекул, связанные со скоростью звука с соотношением а/т ~ с. При более высоких мощностях пучка, когда испарение еще не носит характера взрыва, все еще можно говорить о существовании у образца поверхности, разделяющей паровую и конденсированную фазы. Очевидно, что скорость перемещения этой поверхности детерминирована законом сохранения энергии, а учет испарения в уравнении теплопроводности может быть произведен путем наложения соответствующих условий на границе.
раздела конденсированной и паровой фаз. В простейшей форме граничное условие на фронте испарения имеет вид [12] дТ.
K~ + F~LRVf = О, где К — коэффициент теплопроводности конденсированной среды (теплопроводностью паровой фазы здесь пренебрегаем), F — поглощенная в волне испарения плотность энергии падающего пучка, L — удельная теплота плавления, R — плотность конденсата, Т — его температура, V/ - скорость перемещения фронта испарения. Если Т превышает температуру плавления, то испарение происходит не с твердой, а с жидкой (кипящей) поверхности при давлении насыщенного пара, причем температура пара и вещества с разных сторон поверхности совпадают и равны температуре кипения. Расчет температурных полей в этом случае может быть выполнен с использованием уравнения теплопроводности в рамках задачи Стефана, учитывающей скорость движения границы раздела твердой и жидкой фаз вещества на основе уравнения теплового баланса [13]. При этом температура на границе раздела, как и в случае границы жидкость-пар, меняется непрерывно, а ее производная претерпевает скачок, описываемый соотношением теплового баланса: дт (хь-о, г) дт{хь + о, t).
——~-= QKVb, ох ох где Q — удельная теплота фазового перехода, Kl, Ks ~ теплопроводности в жидкой и твердой фазах соответственно, Vj, — скорость движения границы: Vb = dxb/dt. Динамика фазового перехода типа плавления рассматривается подробно в разделе 1.2.
При дальнейшем повышении плотности энергии падающего пучка механизм равновесного испарения нарушается. В настоящее время описание процессов испарения такого типа производится в рамках модели с кнудсе-новским слоем [И]. Критерием перехода от модели с кипящей поверхностью к модели с кнудсеновским слоем является, очевидно, условие.
Р > Р.
1 i сг> где Р — давление вблизи поверхности, Pcr ~ критическое давление перехода жидкость-пар. Расчеты, с использованием уравнения Больцмана показывают, что в кнудсеновском слое пара, непосредственно прилегающем к поверхности конденсата и имеющем толщину порядка нескольких длин свободного пробега, устанавливается температура Тк несколько меньшая, чем температура Т жидкой фазы:
Тк = 0,65 Т и плотность.
РК = Pv, где ру — плотность насыщенного пара. Кроме этого, на границе раздела следует задать условия, отвечающие законам сохранения энергии, импульса и массы.
Даже этот метод описания испарения может стать совсем неадекватным задаче в том случае, когда средняя энергия молекул е в приповерхностной области намного превосходит теплоту испарения. В этих условиях резкой границы между конденсированным и испарившимся веществом не существует, так что уравнения среды в переходной области следует решать совместно с уравнениями состояния. Подобные граничные условия изучались в задачах инерционного ядерного синтеза (см., например, [14]), когда генерируются интенсивные ударные волны, т. е. фактически происходит взрыв поверхностного слоя вещества. Критерии перехода от модели двухфазной границы к этому режиму испарения сформулированы в работе [15].
Фазовые переходы в рамках задачи Стефана.
Как известно, моделирование фазовых переходов типа плавления — затвердевания! испарения — конденсации приводит к задаче Стефана. Соответствующие математические модели характеризуются наличием подвижной, заранее неизвестной границы S фазового перехода (см., например, [23]). Формулировка задачи Стефана впервые была дана Г. Ламе и Б. П. Клапейроном в 1831 году, для частного случая сосуществования жидкой и кри-сталической фаз при одинаковой температуре [24]. В 1889 Стефан опубликовал четыре работы, в которых была сформулирована эта проблема в общей формеи которая, фактически, принята и в настоящее время [25]. В соответствии с этим подходом для межфазовой границы принимается следующее условие.
Г{х3 + О, t) дТ{х3 — 0- t) soV-дх—q-дх-= LPsoiVs• (v.3).
Здесь Vs = d^s/dt — скорость движения границы: S, Ksoi и Кцд — коэффициенты теплопроводности материала для твердой и жидкой фаз, L и psoi — удельная теплота плавления (или энтальпия фазового, перехода) и плотность соответственно. Условие (v.3) имеет ясную физическую интерпретацию. Действительно, согласно закону Фурье, поток тепла j пропорционален градиенту температуры, j =- —К grad Т.
Поэтому, левая сторона (v.3), представляющая собой разность потоков тепла поступающего на S и выходящего из нее, описывает тепло, поглощенное единицей площади границы за единицу времени. Выражение в правой, стороне (v.3), очевидно, представляет собой тепло, связанное с плавлением или затвердеванием материала, который за единицу времени пересекла единица площади перемещающейся границы раздела S.
Полная математическая формулировка задачи Стефана включает, кроме (v.3), условие, учитывающее тот факт, что фазовый переход происходит при постоянной температуре,.
T|s = TV • (v.4) где Т* - температура плавления, а также — закон сохранения энергии: дТ.
PC = -divj + q (x, t).
Здесь g (x, t) описыват мощность внешнего источника тепла, а С — коэффициент теплоемкости. В основополагающей работе Стефана рассматривался частный случай.^(х, t). = 0, так что вся динамика теплопереноса полностью определялась лишь температурным градиентом в среде.
Если рассматриваемая задача включает достаточно реалистичные начальные и граничные условия, очень часто довольно сложные, то задачу Стефана практически никогда не удается решить аналитически и приходится прибегать к численным расчетам на ЭВМ. Некоторые конкретные примеры граничных условий рассматриваются ниже.
Обычно соотношения (v.3) — (v.4) используются в численных алгоритмах в явном виде [26]. В 1953 году А. Н. Тихоновом и А. А. Самарским был предложен другой подход [27], согласно которому условия (v.3) — (v.4) включаются непосредственно в уравнение сохранения энергии. В результате получаем следующую обобщенную формулировку задачи Стефана: рС + L 5{Т — Г*)) ^ + v grad т) = div (K grad Т) + д (х,?), (v.5) где L 5(T — T*) dT/dt дополнительный вклад в теплоемкость тепла, израсходованного на фазовый переход, v grad Т учитывает возможное температурное изменение, обусловленное конвекцией (в задачах, рассматриваемых в данной главе, не возникает необходимости учитывать этот член). Таким образом, главная идея этого подхода сводится к предложению учета удельной теплоты плавления L в качестве дополнительной компоненты теплоемкости рС, которая, однако, дает вклад только в точке фазового перехода, когда Т = Т.
В дальнейшем Самарский и его ученики превратили эту идею в эффективные численные алгоритмы (см., например [28, 29]). Но даже в этом подходе уравнение (v.5) обычно считается только следствием условия (v.3).
В разделе 1.2 будет показано, что на самом деле условие (v.5) дает правильное описание фазовых переходов даже в том случае, когда соотношения (v.3) и (v.4) не применимы.
Процессы трекообразования в раках модели температурного пика (МТП).
Латентные треки ускоренных ионов в высокотемпературных сверхпроводниках способны играть роль центров пиннинга магнитных потоков, существенно увеличивая при этом плотность критического тока [33, 34, 35]. Не смотря на большую практическую важность этих процессов, удовлетворительной теории формирования треков в высокотемпературных сверхпроводниках пока не существует. Теоретически возможны несколько различных механизмов, среди которых модель температурного пика (термопика) [36], модель ионного пика [37], а также другие более изощренные модели (см., например, [38]). Модель ионного пика объясняет формирование трека созданием положительного ионного облака вокруг траектории иона, которое «взрывается» из-за электростатического отталкивания — происходит так называемый кулоновский взрыв. Согласно модели термопика, материал плавится, в пределах цилиндра вокруг траектории ускоренного иона, если температура решетки превышает температуру плавления. Последующий процесс быстрого остывания приводит к формированию аморфной области, которая может наблюдаться методами электронной микроскопии высокого разрешения1.
Первая попытка разработки модели термопика (МТП) для описания процесса формирования треков в высокотемпературных сверхпроводниках была предпринята в [39]. Физическими причинами, которые привели авторов к этому, послужили следующие факты:
• нанодифракция в УВагСизС^-я.показывает, что материал в области треков находится в аморфном состоянииэтом смысле укоренившиеся термины «латентный» или «скрытый» по отношению к рассматриваемым в диссертации трекам не вполне корректны.
• электронная микроскопия показала наличие искажения решетки вокруг треков, соответствующее расширению материала внутри треков [39, 40].
Хотя оба эти факта действительно могут быть интерпретированы естественным образом как последствия процессов плавления материала (сопровождающегося его расширением) и последующего затвердевания, вполне возможны также и другие объяснения. Поэтому исследования с использованием почти той же самой экспериментальной-техники привели авторов работы [34] к другому заключению о том, что механизм формирования трека в УВазСизО^ основан на процессах ионизации [37].
Первое описание формирования треков в высокотемпературных сверхпроводниках с помощью МТП, предпринятое в [39], не учитывало теплоту плавления и поэтому предсказывало существенно большие размеры треков, чем экспериментально измеренные. Чтобы каким-то образом интерпретировать это различие, в [39] была предложена интересная гипотеза «эпитаксиального восстановления», согласно которой не вся расплавленная область переходит в аморфную фазу, а только ее внутренняя часть, в то время как внешняя часть «восстанавливается», т. е. возвращается в прежнее кристаллическое состояние.
В работе [41] для объяснения эволюции размеров треков с изменением энергии налетающего иона для облучаемых сверхпроводников УВагСизОу-^ и Bi2Sr2CaCu208 был предложен феноменологический подход, также основанный на концепции термопика. Хотя эта модель успешно объясняла ряд наблюдаемых особенностей трекообразования, она содержала параметры, не зависящие от физических свойств материалов и, поэтому, ее можно рассматривать только как полезное предварительное исследование проблемы.
Более детальная модель процесса формирования треков в УВагСизОу-^, основанная на системе связанных уравнений для температур электронов и атомов, была предложена в работе [42] по аналогии с моделью термопика [31], развитой в Казне (Франция) для описания треков в аморфных металлах и полупроводниках. В качестве свободного параметра в этой версии МТП принята средняя длина свободного пробега электрона относительно атомов решетки, Л =• /Der. Здесь De — температуропроводность возбужденных электронов вблизи траектории иона, которая предполагается постоянной (для данного материала) и приблизительно равной 1 — 2 см2/с [43]. Параметр г — время электронно-атомной релаксации — может быть оценен с использованием результатов фемтосекундых лазерных экспериментах [44, 45]. Другие количественные параметры, используемые в модели, это известные теплофизические характеристики облучаемого материала, такие как теплопроводность электронов и атомов Ке и К^ теплоемкость Се и Сг, плотность р твердой и жидкой фаз, температура плавления Тт, удельная теплота плавления Qf.
Значение параметра, А ~ 18 нм, найденное в [42] для УВагСизОу-^, оказалось близким к соответствующей величине, полученной для аморфных металлов и полупроводников, а время электронно-атомной релаксации г находилось в хорошем согласии с лазерными фемтосекундными экспериментами. Однако, простые аналитические оценки выполненные в [42], показали, что экспериментально наблюдаемую зависимость радиусов треков от выделившейся при прохождении иона энергии можно объяснить только в том случае, если предположить приблизительно линейную зависимость т от Те (такая зависимость следует, в частности, из теории Аллена [46]). В этом отношении описание формирования трека в УВагСизС^-я оказалось существенно отличным от предложенного в Казне, где величина параметра г считалась не зависящей от температуры. В диссертации т (Те) зависимость учтена в явном виде, т. е. функция т (Те) включена в систему уравнений, описывающих формирование трека.
Охлопывание кавитационного пузырька в D-ацетоне.
С тех пор как в 2002 году в работе [78] было, сообщено о регистрации продуктов термоядерных реакций D (d, p) T и D (d, n)3He на, уровне 105 с-1 при акустической кавитации в СзБеО, эта область исследований стала объектом пристального внимания не только акустиков, но и специалистов по ядерной физике, а также физике плазмы. Последние экспериментальные данные по наблюдению продуктов ядерных реакций при схлопывании кавитационного пузырька в D-ацетоне опубликованы в [77]. Они по-прежнему свидетельствуют в пользу этой новой, чрезвычайно простой возможности осуществления реакций термоядерного синтеза2.
Описание процесса схлопывания пузырька, предложенное Р.И. Нигма-тулиным с соавторами [80], основано на численном решении следующей системы уравнений в форме законов сохранения массы, энергии и импульса для газообразной и жидкой фаз в приближении сферической симметрии коллапсирующей кавитационной полости: I + Я (Риг2) = 0, + ?| (ри^).+% = о, ч где р, и, р, Т, е и, А — плотность массы,-скорость, давление, температура, плотность полной энергии и коэффициент теплопроводности среды.
В решаемую систему входят также уравнения состояния и условие Герца-Кнудсена-Ленгмюра [79] на границе фазового раздела, учитывающее коэффициент конденсации а. Для ацетона, а — 1. Это означает, что доля молекул пара, которые при ударе о поверхность раздела фаз отражаются обратно, пренебрежимо мала. По этой причине общая масса пара при сжатии пузырька быстро уменьшается. Кроме того, в случае поглощения молекулы границе раздела передается импульс, который приблизительно в два раза меньше, чем импульс, передаваемый при отражении. Все это сильно снижает сопротивление паров схлопыванию пузырька и, по-видимому,.
7 2Результаты, приведенные в [?], имели широкий научный резонанс и стимулировали весьма оживленную дискуссию (см. обзор [84].) является одной из основных причин создания высоких температур и давлений при акустической кавитации в D-ацетоне.
При сверхвысоких давлениях в [80] используются уравнение состояния Ми-Грюнайзена [81], потенциальные функции Борна-Майера [81], которые описывают сверхсжатые жидкости, и современные экспериментальные данные об ударной адиабате в ацетоне [82]. Учитываются диссоциация и ионизация молекул в момент схлопывания пузырька, и соответствующие им потери энергии, а также влияние молекулярной, электронной и ионной теплопроводности.
Численное решение полученной системы уравнений было найдено в [80] с помощью метода Годунова [83]. В результе было установлено, что значения температуры ионов Т и плотности р в пузырьке на стадии схлопывания достигают очень больших значений: Т ~ 107 -г Ю8 К, р ~ 100 г/см3 на расстояниях г ~ 5 -г 25 нм от центра пузырька в течение интервалов времени At ~ Ю-13 -fЮ-12 с [78].
Выход термоядерных нейтронов в расчете на один пузырек может быть рассчитан по формуле где П£> - объемная концентрация ядер дейтерия- < va (y) > - эффективная скорость реакций D — D-слияния. Он оказывается в хорошем согласии с экспериментом [84].
В данной модели, в соответствии с гипотезой Р. И. Нигматулина, для диссоциированной фазы берутся те же уравнения состояния, что и для недиссоциированного ацетона («замороженная» ударная адиабата). Хотя подобное предположение и противоречит экспериментальной ударной адиабате [82], оно может быть оправдано быстротечностью завершающей стадии процесса схлопывания. Поскольку характерные времена протекания процесса в данном случае существенно короче, чем в экспериментах, описанных в [82], можно допустить, что существенная перестройка пространственных положений фрагментов, по сравнению с их расположением в молекулахне успевает произойти. Таким образом, физический смысл этой гипотезы сводится к приближенному учету термодинамической неравновесности процесса.
Кроме этой гипотезы, в теоретическом описании [80] используются два правдоподобных, хотя и не доказанных строго, предположения о том, что в многопузырьковой системе в акустическом поле во время схлопывания формируется дополнительный положительный импульс давления в несколько десятков атмосфер и о том, что электроны плазмы не успевают нагреться за время схлопывания 10~12 с). Последнее предположение можно назвать гипотезой перегрева ионной компоненты неравновесной плазмы.
Диссертационная работа устроена следующим образом:
В первой главе диссертации предпринята попытка построения реалистической модели термоупругих процессов в материалах, испытывающих интенсивную ионную бомбардировку. С этой целью развита модель пространственно-временной динамики энерговыделения в тонких пленках, отвечающая имеющимся источникам со взрывной ионной эмиссией, а также произведены расчеты последующих процессов термоупругой релаксации. В результате анализа используемой математической модели получен важный вывод о невозможности постановки соответствующей задачи Стефана для описания движения границы раздела твердой и жидкой фаз (не смотря на ее повсеместное использование в подобных расчетах). Предложен другой, адекватный рассматриваемой задаче метод решения, который затем используется, существенно во второй главе диссертации при моделировании процессов трекообразования.
Во второй главе исследуется процесс формирования треков в высотем-пературных сверхпроводниках на основе модели температурного пика. Модель, используемая в этой работе, построена по аналогии с МТП развитой в Казне (Франция), однако все основные предположения МТП обсуждаются здесь гораздо более основательно, чем это было сделано ранее, причем как с физической, так и математической точек зрения. Дополнительно учтена зависимость скорости электрон-атомной релаксации от температуры электронов, а также предлагается более точное описание процесса энерговыделения. Невозможность описания процесса трекообразования в рамках традиционной постановки задачи Стефана потребовала разработки и тщательного тестирования новых методов численного решения системы связанных уравнений модели.
Третья глава посвящена проверке гипотезы перегрева, а также связанному с этим развитию модели в направлении явного учета процессов охлаждения ядер электронами в сверхплотной сильно неравновесной плазме, образующейся при схлопывании кавитационного пузырька в D-ацетоне. Поскольку в данной модели предполагается, что электроны характеризуются определенной температурой Те, расчеты энергии электронов Ее> фактически сводятся к расчетам их теплоемкости Се (Те). Информация о параметре электрон-ионной связи д, как и в процессе трекообразования, эквивалента знанию времени т электрон-ионной релаксации. Основные проблемы, связанные с предлагаемым в диссертации обобщением модели [82] связаны именно с необходимостью расчета параметров g и Се{Те), так как информация об остальных параметрах уже содержится в той или иной форме в работах [82] и [51].
Для исследований процессов термоупругой релаксации описаных в первой главе диссертации и для изучения процессов трекообразования описаных во второй главе, разработан комплекс программ на языке FORTRAN.
Для расчетовов теплоемкости электронов в дейтерированном ацетоне в процессах схлопывания кавитационного пузырька разработан комплекс программ в среде SCILAB 3.0-RC1 (May 17, 2004), которая. является бесплатно распространяем аналогом MATLAB. Тексты программ находятся в приложении.
3.1.5 Выводы.
Оценки с использованием аналитической формулы (3.4), полученной Спитцером [92], показывают, что электронная компонента плазмы действительно может оставаться относительно холодной на протяжении заметной части времени существования сверхплотного состояния вещества в схлопы-вающемся кавитационном пузырьке в CsDgO. Расчеты, выполненные на основе такого модельного допущения, предсказывают в этом случае возможность осуществления термоядерных реакций D-D слияния со скоростью порядка одного события на одно схлопывание пузырька при температуре D-ацетона около 273 К [78, 84]. Оптимизма относительного этого новой возможности реализации термоядерного синтеза прибавляют также два дополнительных обстоятельства. Во-первых, согласно [95], более тщательное численное моделирование показывает, что оценки, выполненные с помощью аналитических формул типа (3.4), использовашихся в данной работе, дают некоторое занижение величины времени релаксации. Во-вторых, расчеты в рамках гидродинамической модели [78, 84], говорят о том, что дальнейшее понижение температуры D-ацетона приводит к уменьшению времени существования термоядерной плазмы почти до значения Ю-13 с, в то время как выход нейтронов при этом увеличивается.
3.2 Теплоемкость электронов в плазме, образующейся при охлопывании кавитационного пузырька в D-ацетоне.
3.2.1 Фомулировка модели и результаты расчета.
Предполагая локальную компенсацию электрического заряда электронов и ионов, среднее число электронов с энергией Е будем описывать распределением Ферми для невзаимодействующих частиц:
W = е (Е-«(т))/квт + i ' (3−9) где Т — температура электронов, К в — постоянная Больцмана, /л (Т) — химический потенциал электронов, находящихся в электрическом поле ионов. Поскольку для ацетона температура диссоциации Т^ существенно меньше температуры ионизации Ti (Td ~ 3 • 103 К, Т- ~ 1,2 • 105 К), будем считать, что ионизация происходит позже диссоциации и, следовательно, энергии связи электронов совпадают с потенциалами ионизации отдельных атомов [96].
Из того, что при Т = 0 все электроны находятся в связанном состоянии следует, что значение химического потенциала ц при нулевой температуре равно наибольшей энергии связи электрона в атомах С, D и О, т. е. 871,1 эВ. Выделяя в выражении для химического потенциала эту величину отдельным слагаемым,.
КТ)=М0) + Дм (Т), выражению (1) нетрудно придать вид е-(а+А")/квт + i ' е (ек-л цуквТ + х ' соответственно для связанных и свободных состояний. Здесь.
Ei = fi (0) — Ei энергия связи, отвечающая г—му электронному уровню, р1.
Sk —.
2 т энергия свободного электрона. Полученные формулы позволяют трактовать химический потенциал как своего рода потенциальную яму, при выходе из которой электрон становится свободным (затратив на выход энергию, равную энергии связи). Зависящая от температуры поправка к постоянной части химического потенциала А/л (Т) находится из условия сохранения полного числа электронов г к где суммирование выполняется с учетом кратности вхождения атома данного сорта в молекулу СзБбО. Заменяя обычным образом [93] сумму по состояниям свободных электронов интегралом (здесь v — объем, приходящийся на одну молекулу вещества) л/2г>т3/2 v^dey/2v (mKBTf'2 J.zdz.
TT2/i3 J о e (e-An)/KBT+1 JQ ez~Av/KBT + 1 ' получим уравнение для определения величины поправки Afi (T).
Ущ (-, Л * Т-1 = const Т3/2 Г /ldl-, (3.11).
Z-J г у е-{?г+Ац)/КвТ + IJ JQ ez-An/KBT + I' ^ > где щ — кратность вхождения атомов в формулу вещества.
На рис. 3.2 показана. расчетная зависимость А/л (Т) для СзБбО. Тот факт, что поправка Ац (Т) имеет отрицательное значение объясняется большим статистическим весом свободных состояний, приводящим к понижению температуры начала эффективной ионизации £i/К в, которую можно было бы ожидать без учета этого обстоятельтва. На рис. 3.3 показано среднее число электронов,.
NU) = 1 — (Nitb), 82.
Т [106 К].
Рис. 3.2. Расчетная зависимость Ац{Т) для СзБ60 при сжатиях р/р0 = 1, 10, 50 и 100 (кривые 1, 2, 3 и 4 соответственно). 1.
0.8 ^ 0.6 0.4 0.2.
0 0.5 1 1.5 2 2.5.
Т [106 К].
Рис. 3.3. Среднее число электронов, выбиваемых с уровней? i — 11.3, 138.1, 490, 871.1 (кривые 1, 2, 3 и 4 соответственно).
Рис. 3.4. Температура 90% ионизации различных электронных уровней для степени сжатия р/ро — 1 и 100.
Рис. 3.5. Температура 90% ионизации в зависимости от сжатия р/ро для самого глубокого электронного уровня? i — 871.1 эВ. освобождающихся с i—го уровня в зависимости от температуры. Видно, что представленная здесь картина последовательного возбуждения все более глубоких уровней вполне аналогична той, которая имеет место при ударной ионизации или в равновесной плазме. Рисунки 3.4 и 3.5 иллюстрируют некоторые важные с физической точки зрения характеристики плазмы, образующейся при схлопывании кавитационного пузырька в D-ацетоне. Так, на рис. 3.4 показана зависимость температуры почти полной ионизации электронных уровней от их энергии связи. Рис. 3.5 дает информацию о зависимости температуры почти полной ионизации всех атомов плазмы от степени сжатия.
Теплоемкость электронов, отнесенная к одной молекуле вещества, может быть найдена прямым дифференцированием полной энергии свободных электронов: где суммирование выполняется по всем связанным состояния электронов. На рис. 3.6 представлена рассчитанная таким образом теплоемкость электронов в области сильной ионизации при сжатиях р/ро =1,10 и 100, переведенная в боле’е удобные для практического использования единицы. Видно, что она не описывается линейной зависимостью С = 7 Т, справедливой для вырожденного электронного газа в широкой потенциальной яме с плоским дном при температурах меньше или порядка температуры Ферми, когда постоянная (Зоммерфельда) 7 может быть рассчитана по формуле где N — полное число электронов в молекуле [93].
3.2.2 Интерполяционные формулы.
Для моделирования процессов схлопывания кавитационого пузырька важно иметь простые, но в тоже время, достаточно точные интерполяци.
6 4 2 О О.
0.1.
0.2.
0 0.2 0.4 0.6 0.8.
1.2 1.4 1.6 1.8 2.
Т [10е К].
Рис. 3.6. Теплоемкость электронов для молекулы D-ацетона при разных сжатиях р/ро, рассчитанная в соответствии с рассматриваемой моделью. онные формулы расчета зависимости Су (Те, р/ро). Их нахождение представляет собой отдельную вычислительную задачу, которой посвящен настоящий раздел. С этой целью в табл. 3.1 приведены соответствующие рис. 3.6 значения теплоемкости Су для некоторых характерных температур электронов. Точки в табл. 3.1 выбирались следующим образом. С помощью программы 1, помещенной в приложении 2, вычислялись значения Су для Те (:Tmin, Tmax) К с точностью Ю-4, где Tmin =9000, 9500 и 10 000 К, Ттах = 2 • 106 К для р/ро =1, 10 и 100. Таким образом, для принятого разбиения.
Cvi = Cy (Ti), Ti = Tmin т • г, % — 0,., п, т = (Ттах — Tmin)/n, где п — число делений рассматриваемого интервала температур. Потом выбирались предварительные значения элементов таблицы (Т^Сук), к? (0, п) в точках, для которых первая или вторая производная Су равны нулю, а также некоторые дополнительные промежуточные точки, в которых значения теплоемкости составляли определенный процент относи.
1 1 1.
1 ШАлААа^ Р/Р0=1.
0 0.5 1 1.5 2.
IML Л/VA/W^ /V p/pQ=10 0 0.5 1 1.5 i i i 2 p/pn=100.
0.3 о.
UJUVVY.А /.
О 0.5 1 1.5 2.
T [10е К].
Рис. 3.7. Относительная ошибка интерполяции кубическим сплайном теплоемкости при промежуточных температурах, получающаяся при использовании табличных значений теплоемкости. тельно скачка между экстремальными значениями на соответствующих температурных интервалах. Таким образом, число табулируемых точек для каждого табличного значения степени сжатия p/pQ было фиксированным. На границах интервалов температур в точках Tmin и Ттах вычислялись приближенные значения производной теплоемкости по формулам Cy{Tmin) = (CVi — Cvo)/t и С’у (Ттах) = (CVn — Cvn-i)/T. Эти значения приведены в табл. 3.2.
Далее по выбранным точкам строился кубический сплайн для теплоемкости с вышеуказанными значениями производной на границах, с целью сравнения полученных таким образом величин Су к с точными. Это позволило определить характерную для табл. 3.1 относительную ошибку интерполяции теплоемкости при промежуточных, не охватываемых таблицей, значениях температур электронов. Далее предварительные элементы.
Заключение
.
В дисертационной работе проведены исследования процессов термоупругой релаксации и процессов трекообразования в материаллах облучаемых ускоренными ионными пучками. Разработаны модели, алгоритмы и программы для математического моделирования этих процессов и для вычисления неизвестных параметров моделей. Основные результаты работы можно сформулировать в следующих пунктах:
1. Разработаны, тестированы и применены новые эффективные алгоритмы расчета, использующие разностные схемы второго порядка точности для решения уравнений теплопроводности с учетом фазовых переходов и источников энерговыделения.
2. Разработан комплекс математических моделей и программ для расчета процессов тепловой релаксации в образцах, подвергаемых облучению ускоренными ионами. С его помощью были выполнены исследования взаимодействий ионов с веществом и получены следующие результаты: а) Изучен процесс формирования и релаксации термоупругих напряжений в материалах, испытывающих интенсивную ионную бомбардировку. Сформулирована модель пространственно-временной динамики энерговыделения в тонких пленках, отвечающая имеющимся источникам со взрывной ионной эмиссией. Установлено, что максимальные механические напряжения, возникающие вследствие теплового расширения вещества, не превышают ни предел текучести железа, ни его временное сопротивление разрыву и таким образом не являются ответственными за модификацию образца. б) В результате численных экспериментов и аналитических оценок получен вывод о невозможности постановки традиционной задачи Стефана для описания движения границы раздела твердой и жидкой фаз в том случае, когда в уравнении модели присутствует источник энерговыделения. в) В рамках концепции температурного пика рассмотрены процессы трекообразования в высокотемпературных сверхпроводниках. С учетом имеющихся теоретических и экспериментальных данных разработана модель процесса энеговыделения, описывающая его пространственно-временную динамику. Уточнена модель электронной подсистемы с использованием теории Аллена и результатов фемтосекундных лазерных экспериментов. г) Продемонстрирована применимость концепции температурного пика для описания процессов трекообразования в иттриевых сверхпроводниках. Установлена бифуркационная зависимость решений от параметра температуропроводности электронов, связанная с процессом электронной закалки.
3. Разработана математическая модель электрон-ионной релаксации в плазме, образующейся при схлопывании кавитационного пузырька в D-ацетоне и создан комплекс программ для расчета неизвестных параметров модели. Вычислено время охлаждения ядер электронами, а также теплоемкость электронов во всем интервале температур и плотностей образующейся плазмы. Решена задача интерполяции теплоемкости электронов с точностью, достаточной для математического моделирования процесса релаксации.
Благодарности.
Выражаю огромную благодарность своим научным руководителям Б. Ф. Костенко и Э. А. Айряну за постановку задач, поддержку и постоянное внимание. Выражаю глубокую признательность профессорам Е. П. Жидкову и И. В. Пузынину за проявленный интерес к работе и полезные дискуссии. Особо хочу поблагодарить О. И. Стрельцову, Я. Бушу, И. Покорны за ценные консультации. Так же хочу поблагодарить соавторов и коллег по работе: И. В. Амирханова, И. Н. Гончарова, Т. А. Стриж, Д. В. Подгайного и Т. JI. Бояджиева. Искренне благодарен своей супруге Юдите за неоценимую помощь и поддержку.
Отдельно хочу поблагодарить дирекцию Лаборатории информационных технологий Объединенного института ядерных исследований за предоставленные хорошие условия для работы.
