Содержание к диссертации
Введение
Глава 1 7
Обобщенная модель ионизационного равновесия 7
1.1 Модель ионизационного равновесия 7
1.2 Модели неидеальности 14
1.3. Микрополевая модель неидеальности 32
1.4 Интерполяционные формулы для функций ферми-дирака 37
Глава 2 41
Построение микрополевой модели 4i
2.1 Микрополевая поправка на неидеальность 41
2.2 Статистические суммы. Системы энергетических уровней лития и водорода 44
2.3 Численный алгоритм расчета модели 48
Глава 3 52
Верификация модели si
3.1 Програмный комплекс 52
3.2 Параметры неидеальности 54
3.3 Результаты расчетов 59
Основные результаты 65
Литература 66
Приложение! 70
- Модели неидеальности
- Интерполяционные формулы для функций ферми-дирака
- Статистические суммы. Системы энергетических уровней лития и водорода
- Параметры неидеальности
Введение к работе
Многократная ионизация разреженных газов описывается системой классических уравнений Саха ([5]). Для этой системы в [5] доказано существование и единственность решения, а в [6] предложен эффективный метод решения, пригодный для реализации на ЭВМ.
Классические уравнения Саха не учитывают квантового вырождения электронной компоненты, а также взаимодействия заряженных (и незаряженных) частиц. Поэтому при повышении плотности плазмы их физическая точность ухудшается. В модель требуется включать соответствующие эффекты.
Для учета вырождения электронов в многокомпонентной плазме в работах [7] и [8] была предложена модель ионизационного равновесия. Ионы при этом можно считать классическими частицами, так как эффект вырождения ионов проявляется лишь при плотностях в 1000 раз превышающих твердотельные, при этом вряд ли правомерно считать вещество плазмой.
Взаимодействие заряженных частиц существенно уже при плотностях порядка атмосферной (например, в капиллярных газовых разрядах). Оно приводит к уменьшению (сдвигу) всех потенциалов ионизации, расщеплению и уширению линий спектров, обрезанию статистических сумм атомов и ионов. Предложено много теоретических моделей для описания этих эффектов. Однако выводы этих моделей лишь физически правдоподобны и отнюдь не строги, а результаты их применения сильно различаются. Поэтому остановимся на этом вопросе подробнее.
Сначала строились классические дебаевские поправки ([9]). Однако большинство физиков считало, что эти поправки сильно завышают эффект неидеальности при высоких плотностях. (Экспериментов, с которыми можно было бы сравниться, не существовало.) К тому же в классической модели Дебая-Хюккеля нарушается электронейтральность среды в целом, а также наблюдается не единственность решений уравнений ионизационного равновесия [30]. Столь серьезные недостатки вообще ставят под сомнение правомерность использования данной модели. Позднее стала популярной модель БДХ (Дебай-Хюккель в большом каноническом ансамбле,[10]). Она приводит к гораздо меньшим эффектам неидеальности; на наш взгляд она их сильно занижает. Эти и многие другие работы обосновывались различными приближениями теории возмущений с помощью фейнмановской диаграммной техники. Обзор это направления есть в [11].
Другим подходом, который мы считаем наиболее перспективным в настоящий момент, является микрополевая модель неидеальности ([12]). В ней использована идея плазменного микрополя, давно применявшаяся в спектроскопии ([13]). Было предложено рассчитывать термодинамические поправки на неидеальность с помощью этого микрополя. Эта модель существенно лучше согласовывалась с результатами экспериментов для очень плотной плазмы, проведенных в 1980-е годы, чем [9]-[11].
Однако оказалось, что микрополевая модель [12] требует доработки. Во-первых, использовавшееся там выражение для сдвига потенциала ионизации не позволяло написать простые выражения для поправок на неидеальность к термодинамическим функциям; их можно было восстановить только очень громоздкими численными расчетами (многократным интегрированием сдвига потенциала ионизации). Во-вторых, ряд деталей модели был сделан недостаточно тщательно. На первый взгляд это казалось малосущественным. Однако, при попытке рассчитать состав высокоплотного слабо нагретого вещества, алгоритмы «срывались».
Таким образом, ставилась задача, во-первых, построить такой сдвиг потенциалов ионизации с учетом мгновенного флуктуирующего микрополя, который имел бы несложную поправку на неидеальность, пригодную для расчетов термодинамических характеристик плазмы. Во-вторых, необходим алгоритм, позволяющий производить расчеты модели ионизационного равновесия с указанной поправкой на взаимодействие на ЭВМ. В данной диссертационной работе приводится решение указанных проблем.
Первая глава работы содержит некоторые основные результаты, которые необходимы для дальнейших изысканий. Здесь представлена модель ионизационного равновесия, вводится понятие мгновенного флуктуирующего микрополя, а также приводится микрополевая модель неидеальности из [12], взятая в данной работе за основу. Помимо этого в первой главе содержится обзор наиболее известных (классических) моделей неидеальности, а также необходимая в дальнейшем интерполяционная формула для функций Ферми - Дирака, позволяющая свести интегральную зависимость к несложному дробно-рациональному выражению, удобному для численных расчетов.
Вторая глава посвящена построению микрополевого сдвига потенциалов ионизации А<р для модели ионизационного равновесия. В ней указывается соответствующая этому сдвигу поправка на неидеальность AF, необходимая для расчета термодинамических характеристик плазмы. Здесь же приводится и алгоритм, пригодный для численных расчетов ионизации по описанной модели. Помимо этого во второй главе описывается построение системы энергетических уровней для Li, позволяющей достаточно аккуратно воспроизвести статистические суммы с учетом их обрезания.
В третьей главе дается описание программного комплекса, построенного на основе предложенной микрополевой модели неидеальности, и позволяющего производить расчет ионизации и термодинамических свойств плазмы одного элемента. С его помощью проводится анализ количественных характеристик различных параметров неидеальности. Рассчитываются таблицы ионизации для Li, Al, Си, демонстрирующие работоспособность программы, достаточную устойчивость численного метода и адекватность приведенной микрополевой поправки на взаимодействие.
Основные результаты были изложены в работах [1-4]. Они докладывались на VII-x Харитоновских чтениях, проходивших в марте 2005 года в городе Саров.
Для написания программы использовался язык C++. Текст приводится в приложении 1.
Модели неидеальности
Параметр неидеальности вводится в разных работах неодинаково. Наиболее естественным определением является энергетическое: УЕ=-ЕПОТ1ЕШН, ЕШН =Ее+Ея, ЕП0Т=АЕ + ЗЕ (1.20)
В нем не учитывается та доля потенциальной энергии, которая идет на ионизацию и возбуждение, ибо она не имеет отношения к неидеальности. Для моделей, удовлетворяющих теореме вириала (1.12) - (1.13), выполняется 7Е=-2Ы 1РКИН, Ркш =Ре+Ря. (1.21) Упомянем также параметр, более соответствующий конденсированной (т.е. очень плотной) Z - зарядной полностью ионизованной плазме (без нейтральных атомов): Г = Z2 IRT (или Г = Z Xj/RT для смеси); (1.22) при газовых плотностях он неудобен. Энергетический же параметр остается разумным при любых плотностях.
Область применимости моделей неидеальности часто характеризуют диапазоном значений параметра (1.20). Например, модель идеальной плазмы получается из (1.2) при AF = 0 и обрезании сумм (1.14) на некоторых определенных номерах. Очевидно, в этой модели AE = AP = A(p = 8P = S p = 0, а условием ее применимости с точностью 1% будет YE — 0 01 (разумеется, эту величину уЕ надо вычислять с помощью более точных моделей; в самой модели идеальной плазмы по определению всегда уЕ = 0). Из-за такого жесткого ограничения данную модель используют лишь для достаточно разреженной плазмы.
Напомним некоторые свойства идеальной плазмы. Строго доказано, что ее уравнения ИР имеют единственное положительное решение для сколь угодно сложного состава и любой степени ионизации [7]. Давление и энергия всегда положительны. Это удобно в прикладных расчетах, когда зачастую приходится использовать модель далеко за границами ее применимости. Видно, что модель идеальной плазмы даже при далекой экстраполяции дает качественно правдоподобные результаты, хотя количественная точность может стать плохой.
Частично ионизованный газ состоит из подсистем нейтральных и заряженные частиц. Величина АЕ обусловлена взаимодействием только заряженной подсистемы, ав Ея входит кинетическая энергия как ионов Еи, так и нейтральных атомов ЕА. При слабой ионизации последняя может стать преобладающей. Поэтому рассматривают также параметр неидеальности заряженной подсистемы ул=-АЕ/(Ее+Еи) (1.23) Если плотность плазмы не чрезмерно велика и электроны являются невырожденными, то Ее =3/2хеТ и справедливо соотношение Гп ГЕ=(\ + хе)1(\ + хе -х0), х0 = ,0. (1.24) і
Для слабо ионизованного газа 1 -х0 «хе «1 и УЕ УП\ ПРИ ЭТОМ возможна ситуация, когда плазменная подсистема будет сильно неидеальна, а газ в целом с хорошей точностью идеален.
Описанные выше способы являются классическими (традиционными) подходами к количественному вычислению параметра неидеальности в плазме. Однако в третьей главе данной работы будет предложен другой подход к решению этого вопроса.
Модель Томаса-Ферми является квазиклассическим приближением к квантово-механической модели атома. В прецизионных квантово-механических моделях учитывается оболочечная структура внутриатомных электронов. В модели Томаса-Ферми предполагается, что оболочки «смяты» сильным взаимодействием соседних атомов, т.е. высоким давлением окружающей среды. Спектр атомных электронов при этом вместо дискретного становится сплошным.
Предполагается, что в центре атомной ячейки находится ядро. При этом плотность электронов внутри ячейки непостоянна, а зависит от расстояния от ядра г. Локальная плотность кинетической энергии электронов здесь также выражается через функции Ферми-Дирака, как и в модели ИР, только с учетом зависимости от г. Она учитывает электронно-статистическое взаимодействие электронов друг с другом и с ядром. Потенциальная энергия электронов в модели Томаса-Ферми не мала. Поэтому в рамках данной модели возможна сильная неидеальность. Однако конкретный вид этой сильной неидеальности для нас не представляет интереса, поскольку в ней не учтена оболочечная структура (а для модели ИР эти эффекты очень существенны). Тем не менее, у этой модели имеется интересный предельный случай при Епот 1ЕКИН - 0 - однородный электронный газ.
Интерполяционные формулы для функций ферми-дирака
Функции Ферми-Дирака определяются следующим образом: 7 ( )= Г т -т к -х С1-60) + ехр( -х) v Они возникают в задачах квантовой механики при изучении системы фермионов (например, электронов в атоме, металле, плазме) при усреднении различных степеней импульса по распределению Ферми-Дирака. В этих задачах индекс к принимает целые и полуцелые значения. Используя соотношение I[ (х) = klk_x (х), можно доопределить функции для всех нецелых к -\. Свойства функций Ферми-Дирака и подробные таблицы для наиболее употребительных функций приведены в работах [40-44].
В задачах квантово-статистической модели атома требуются не только эти, но и другие функции, различные их комбинации и более сложные выражения. (В данной работе, в частности, необходимо вычислить значение аргумента х для известного значения функции /1/2W.) Затабулированы лишь не многие из них. Кроме того, таблицы не удобны для численного решения задач на ЭВМ.
Итак, требуется аппроксимировать функции на бесконечном интервале -оо х оо. Вдали от особых точек большинство функций можно разложить в полином. Вблизи особых точек и на бесконечном интервале нам нужно искать не полиномиальное разложение.
Рассматриваемые функции вблизи особых точек имеют разложения следующего вида fKaelxYjakekX, - (1.61а) Jt=0 f flxHftfikx- kf39x-+ao (1.61b) к=0 Если в качестве аргумента взять У = I\iiix), то асимптотики будут выглядеть так: / аук(\ + аУ+...) при х -со, у - 0 ; (1.62а) / /?Х і/3(1 + Д.у"8/3+...)прИ х- 00, д;_ оо (i.62b) Можно искать аппроксимацию как дробно-рациональную функцию следующего вида: т п =2 ,У(5 У)- . (1.63) /=о 1=0
При этом исследование асимптотик налагает дополнительные ограничения на пределы суммирования в числителе и знаменателе и на отношение главных (т.е. либо старших, либо младших) коэффициентов сумм.
Такая аппроксимирующая функция не может удовлетворить сразу обеим асимптотикам, т.е. придется разбивать интервал — о х оо минимум на два и в каждом подбирать свои наилучшие параметры. Кроме того, поскольку структуры аппроксимирующей и исходной функций не вполне совпадают, для хорошей точности требуется довольно много свободных параметров. Например, в [44] для нахождения аппроксимационной формулы для 13/2(х) необходимо было разбивать область определения на три интервала, а для получения точности 10"8 брать до 10 варьируемых коэффициентов в каждом интервале.
Поэтому, учитывая приведенный пример, будем требовать, чтобы аппроксимирующая функция в интервале аппроксимации имела разложение, близкое к исходной функции, а вне этого интервала имела бы правильное асимптотическое поведение. Такими являются, например, следующие функции: т п z = ay(l + apy2py(} + 2bpy2n-v,0 y y, (1.64а) Р=\ /7=1 т п z = а0у ( сср+УрПУ( Ьр+гУірПГ ,УЇУ « . (1.64Ь) /?=0 р=0
Особенно удобными оказались v = l/3 и // = 1/4. Главные коэффициенты и пределы суммирования определяются с помощью разложений (1.62).
Т.к. значения функции Ферми-Дирака изменяются от очень маленьких до очень больших значений, то в качестве критерия наилучшего приближения требовалась малость относительной а не абсолютной погрешности, т.е. функция заданного класса z(x,a0,av...,an) является наилучшим приближением к F(x), если параметры a0,alt...,an выбраны так, что max z / F -11 минимален.
Очевидно, параметрами формул (11) являются не только коэффициенты a,b,a,fl, но и величина у. Пределы суммирования следует рассматривать как свободные параметры, так как они определяют число коэффициентов и тем самым - класс аппроксимирующих функций.
Предварительные расчеты показали, что прибавление каждой лишней пары свободных коэффициентов увеличивает точность формул приблизительно в 10 раз (из вида формы следует, что число коэффициентов может меняться только на четную величину). Хорошая точность обеспечивается примерно шестью свободными коэффициентами. Это число и было выбрано для дальнейших расчетов.
Коэффициенты находились итерационным процессом. Теоретически его сходимость доказать не удалось, но практически процесс всегда сходился, хотя довольно медленно. Описание процесса мы приводить не будем, поскольку он не представляет самостоятельного математического интереса.
Статистические суммы. Системы энергетических уровней лития и водорода
Построение системы энергетических уровней лития. Далее нам потребуется два квантовых числа - главное п и орбитальное / (0 / п-\). Система уровней атома водорода хорошо известна и описывается следующими выражениями:
Для детального изучения влияния неидеальности на ионизацию мы выбрали Li. Он достаточно прост для расчета, т.к. имеет всего 3 электрона. При низких температурах это металл с плотностью р = 0,54г/см3 (атомарный объем V »150 а.е.). Он плавится при t = 453.7К и испаряется при t = \620K. Потенциалы его ионизации равны: Li - 5.4 эВ, Ы+ - 75.6 эВ, Ы++- 122.4 эВ.
Соответственно в статистических суммах (2.7) суммирование производится не по одному индексу т, а по обоим указанным индексам п и /. Но поскольку для водорода имеет место вырождение по /, то статистическая сумма для Li++ сводится к суммированию только по одному индексу гг. Для нейтрального Ы статистические веса уровней по-прежнему определяются формулой (2.13). Однако система энергий уровней уже много сложнее. Она приведена в табл.2.1 (там же для сравнения приведена система уровней Н), взятой из [45]. Видно, что у Ы нет вырождения по /. Более того, зависимость от / значительна. Казалось бы, при мощности современных компьютеров можно ввести данные табл.2 непосредственно в статистические суммы и вычислять их. Однако энергии состояний с / 3 и я 11 пока экспериментально не измерены.
Т.о. для практического использования табл.2.1 в расчетах необходимо построить разумную аппроксимацию энергетической системы уровней. Заметим, что с увеличением п зависимость от / ослабевает, и сами энергии сближаются с соответствующими значениями энергий для Н. (В физике часто заменяют высоковозбужденные уровни различных атомов и ионов водородоподобными). Исходя из этого, удалось подобрать следующую неплохую аппроксимацию: (рн 22.6 Р, (2.14) п! п2 (и + 2/)3 5 где рн «13.6эВ - потенциал ионизации водорода. Ее точность вполне достаточна для вычисления термодинамических величин.
Для иона Li+ данных по энергетическим уровням в [45] нет. Однако электронная структура этого иона аналогична нейтральному Не. Поэтому для него использовался аналогичный прием, но вместо реперного столбца энергий Я брались энергии водородоподобного Не+ (они вчетверо выше, чем у Н). Т.о. для системы уровней Ы+ использовалась аппроксимация (2.14) с заменой потенциала ионизации водорода # я«13.6эВ на потенциал 2-й ионизации Не РНе+ «54.4эВ: Рш 22.6
Это хорошо воспроизводило потенциал второй ионизации Li. Разумеется, такая аппроксимация гораздо менее точна, чем для нейтрального атома Li.
Однако сами ионы Li+ появляются лишь при достаточно высоких температурах. Анализ показал, что в этих условиях влияние системы уровней (а, соответственно, и величина погрешности аппроксимации) на состав плазмы существенно меньше, чем при низких температурах.
Подстановка полученных формул в расчеты статистических сумм позволяет учесть любое число высоко возбужденных состояний. Такие численные расчеты были проведены. Они показали, что реально статистические суммы для не слишком разреженных газов очень быстро сходятся благодаря наличию формфактора, учитывающего неидеальность. Величины статистических сумм обычно не привышают 3-4 при газовых плотностях и были еще меньше при твердотельных плотностях.
Параметры неидеальности
Для системы классических уравнений Саха в [6] предложен эффективный метод решения, пригодный для реализации на ЭВМ и доказана его сходимость. Теперь удалось обобщить этот подход и для модели ИР с микрополевой поправкой на неидеальность (см. [1]). Теоретический вопрос о сходимости метода пока остается открытым в виду сложности и громоздкости модели, однако практически алгоритм сходился в очень широком диапазоне температур и плотностей.
Систему уравнений ионизационного равновесия (1.7) можно свести к одному уравнению. Обозначим относительную концентрацию электронов за х. Теперь перепишем уравнения (1.7) в несколько ином виде:
Теперь для нахождения относительных электронных концентраций достаточно воспользоваться балансным соотношением 2 bcrt = хе; i,k В виду крайней сложности системы уравнений ионизационного равновесия (1.7), а, следовательно, и уравнения (2.20), пользоваться для его решения более сложными методами типа метода Ньютона весьма затруднительно. Поэтому мы решаем его комбинацией метода деления пополам и метода простых итераций. Предлагается следующая процедура расчета. Имея на (л + 1)-ом шаге х"и хк (на первом шаге берем начальное приближение х и хк) вычисляем (простая итерация) ju(x") (для этого пользуемся интерполяционной формулой для функций Ферми - Дирака (1.65)), А рк" = А рк"(х",хк") и 8 рк =S(pk"(x",xk"). Далее по (2.22) делением пополам получаем хп+\ по (2.21) вычисляем х%+ и т.д. Несмотря на хорошую устойчивость данного метода, иногда возникают ситуации, когда в процессе вычисления искомый корень уравнения оказывается за пределами отрезка. По-видимому, это связано с тем, что итерации фактически ведутся по двум переменных х и ju. Поэтому по окончанию итерационного процесса необходимо поставить проверку электронного баланса хе = 2-( и, в i,k случае его невыполнения, положить границы отрезка равными х±0.5х. После этого заново повторить итерационный процесс.
Численные исследования показали, что областях низких температур и высоких плотностей данный алгоритм может давать сбой. Это связано с появлением здесь у функции (2.22) трех неотрицательных корней. На рис.2.2 демонстрируется поведение функции (2.22) для лития при р»0.8г/сл/ и Т&60000К, Т&3000К Т&600К, при этом значение переменной х изменяется от 0 до 1.1. Рис. 2.2. Поведение функции (2.16) для Li при р 0.8г/ам3и Т 60000/:, Т 3000/: Г 600/:. Аргумент х меняется от 0 до 1.1. 0.400 0.200 0.000 -0.200 -0.400 0.600 -0.800 -1.000 0.400 0.200 0.000 -0.200 -0.400 -0.600 -0.800 -1.000 0.400 0.200 0.000 -0.200 -0.400 -0.600 -0.800 -1.000 ГЛАВА 3. ВЕРИФИКАЦИЯ МОДЕЛИ. 3.1 ПРОГРАМНЫЙ КОМПЛЕКС
Программный комплекс реализован на языке C++ [3,4]. Он состоит из непосредственно программного ядра и прилагающейся к нему таблице потенциалов ионизации элементов, взятой из [46]. Она приведена в приложении 1. В этой таблице содержатся потенциалы ионизации всех элементов периодической системы от z = 1 (водород) до 2 = 103(лоуренсий), причем кратностей от к = 1 до z. В этой таблице потенциалы низших ионизации элементов взяты из экспериментальных данных, собранных в [45]. Для более высоких кратностей ионизации потенциалы построены в [46] на основе расчетов по релятивисткой модели Хартри-Фока-Слэтера [47] с учетом специальных интерполяционных поправок.
В работе [2] мы проводили наиболее тщательные расчеты с Ы и Н, для которых нам удалось подобрать хорошую аппроксимацию системы электронных уровней, описанную в пункте 2.2. Благодаря этому для указанных элементов удалось с достаточной точностью вычислить статистические суммы с учетом обрезания всех уровней. Однако для других элементов система уровней имеет слишком сложное поведение, и для их статистических сумм не удалось найти достаточно простого способа расчета. Тем не менее, расчет ионизации Ы и Н показал, что аккуратный учет возбужденных уровней слабо влияет на ионизацию и термодинамику. Доминирующий вклад вносит основное состояние и его обрезание
