Сравнение методов численного моделирования распада метастабильного состояния

Размещено на http://www.allbest.ru/

Сравнение методов численного моделирования распада метастабильного состояния

М.В. Чушнякова, И.И. Гончар, Р.А. Кузякин

Аннотация

В работе представлены результаты компьютерного моделирования распада метастабильного состояния, выполненного с помощью двух методов: 1) численного решения уравнения Фоккера-Планка и 2) численного решения стохастических дифференциальных уравнений. Проведено сравнение временной эволюции скоростей распада и их квазистационарных значений. Обсуждаются преимущества и недостатки этих двух методов моделирования.

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

Что может быть общего между делением атомных ядер в реакторе под действием нейтронов, диссоциацией молекул, ростом коллоидной частицы в насыщенном растворе и проводимостью джозефсоновского перехода? А между тем, для описания всех этих явлений используется одна и та же модель: тепловой распад метастабильного состояния. Это один из случаев броуновского движения (в обобщённом смысле) в поле внешней силы [1-5]. Соответствующая модель представляет собой невзаимодействующие броуновские частицы, находящиеся в потенциальной яме. Эта яма отделена барьером от области, где частиц нет совсем или очень мало. Частицы преодолевают барьер благодаря тепловым флуктуациям. Основной количественной характеристикой этого явления служит зависящая от времени скорость распада , определение которой выражается формулой

Здесь - число частиц, покинувших потенциальную яму за промежуток времени , - начальное число частиц в яме, - число частиц, покинувших яму к моменту времени .

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

Определить значение можно, моделируя процесс динамически с помощью двух подходов. В первом подходе распад описывается диффузионным уравнением Фоккера-Планка (УФП) [6, 7]

Здесь - плотность вероятности, которая зависит от координаты , импульса и времени, а потоки и определяются выражениями

Здесь и - инерционный и фрикционный параметры, которые могут зависеть от координаты [7], - обобщённая сила (часто говорят «движущая» сила, driving force). - коэффициент диффузии в пространстве импульсов. Когда найдено решение УФП, скорость распада (СР) вычисляется по методу “поток, делённый на населённость” (flux over population):

Здесь индекс «» соответствует барьеру потенциальной энергии.

Во втором подходе динамика процесса описывается стохастическими уравнениями ланжевеновского типа (УЛ) [4, 8-10]

Здесь и имеют смысл координаты и импульса броуновской частицы; случайная сила представляет собой белый шум:

Скорость распада полученная в данном подходе (будем называть её ланжевеновской СР), вычисляется по формуле (1), технические детали этих вычислений обсуждаются в [11].

Наряду с динамическими подходами, для вычисления квазистационарной скорости распада уже много лет используется приближённая формула Крамерса [1]

Здесь коэффициент затухания , ; индекс «» соответствует дну потенциальной ямы. Точность формулы Крамерса исследовалась неоднократно [2, 4, 10, 12], в некоторых случаях она даёт результат, отличающийся от точного (динамического) на 20-40%.

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

Для решения УФП (2) применяется следующая численная схема

Здесь , и - шаги сеток по времени, координате и импульсу соответственно;, , , , и т.д. Таким образом, УФП сведено к системе линейных алгебраических уравнений с неизвестными , , . Матрица коэффициентов при неизвестных является трёхдиагональной, и система решается методом прогонки [13].

Стохастические УЛ (6) моделируются с помощью схемы Эйлера-Маруямы [14]

Здесь верхние индексы маркируют последовательные моменты времени , все величины в правой части формулы (12) берутся в момент времени ; - случайные числа, имеющие нормальное распределение с нулевым средним и дисперсией 2. Они получены по методу Марсаглиа [14]. Для этого применялся генератор однородно распределённых случайных чисел, построенный нами с использованием работы [15]. Поскольку в основе этого генератора лежат целые числа удвоенной длины (long int на языке «С»), период его оценивается как , что перекрывает наши потребности. Заметим, что при использовании стандартного генератора, основанного на целых числах обычной длины, мы иногда сталкивались с проявлением периодичности.

Результаты динамических расчётов скорости распада показаны на рис. 1.

Рис. 1. Динамические скорости распада и как функции времени для трёх значений параметра , которые указаны на рисунке. Показана также крамерсова скорость

Скорости деления здесь соответствуют трём разным значениям коэффициента затухания. Его удобно характеризовать безразмерной величиной . Расчёты, представленные на рис. 1, выполнены при . Отношение , которое определяет применимость формулы Крамерса (8), составляет 4.04.

На стадии релаксации кривые, соответствующие УФП и УЛ, не совпадают, они сдвинуты друг относительно друга. Это обусловлено тем, что при решении УФП фактически определяется скорость преодоления барьера (см. уравнение ). Тогда как при решении УЛ моделирование отдельной траектории прекращается, когда броуновская частица достигает поглощающей границы, находящейся заметно за барьером. Понятно, что броуновские частицы в среднем позже достигают точки, расположенной дальше от . Именно поэтому кривые со сплошными символами на рис. 1 достигают квазистационарного значения позже. К тому же в этих двух подходах используются неодинаковые начальные условия: при решении УЛ все броуновские частицы покоятся в основном состоянии, тогда как при решении УФП предполагается некоторое распределение по координатам и импульсам.

Из рис. 1 видно, что квазистационарные значения скоростей распада, полученные с помощью УФП () и УЛ (), хорошо согласуются друг с другом и с : типичное различие при данных значениях параметров составляет меньше 2%. Это соответствует области большой и средней диссипации (), так как именно в этом приближении выведена формула Крамерса (8). При меньшем затухании согласие с аналитической формулой разрушается, однако, как показано на рис. 2а, согласие между и сохраняется. Таким образом, вопрос о том, какой из подходов предпочтительнее, решается исключительно на основе анализа затрат компьютерного времени.

Рис. 2. В зависимости от безразмерной величины , характеризующей затухание, показаны: а) отношение квазистационарных динамических скоростей и к аналитической ; б) затраты компьютерного времени (в часах) при двух методах моделирования; в) шаг сетки по времени при решении УФП.

На рис. 2b представлены результаты такого анализа: показано, сколько времени требуется для получения скоростей распада с приемлемой точностью обсуждаемыми двумя методами. Из рисунка видно, что это время при решении УФП сильно зависит от коэффициента затухания (или от ). Дело в том, что при уменьшении коэффициента затухания приходится уменьшать шаг сетки по времени. В противном случае нерегулярно осциллирует, что не имеет никакого физического смысла. Какой шаг использовался в расчётах, показано на рис. 2с.

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

При моделировании уравнений Ланжевена точность скорости распада определяется количеством вовлечённых траекторий (т.е. количеством броуновских частиц). При большей диссипации их общее число приходится увеличивать, так как всё меньше частиц достигает поглощающей границы к заданному моменту времени. В этом методе каждый раз приходится оценивать величину погрешности, с которой определена любая физическая величина - в нашем случае . Однако эта особенность одновременно является и преимуществом метода: можно сделать «пристрелочный» расчёт с небольшим числом траекторий и быстро получить результат, пусть и с большой погрешностью. Это особенно ценно для отладки, а также при подготовке длительных расчётов.

Библиографический список

1. Kramers H. A. Brownian motion in a field of force and the diffusion model of chemical reactions // Physica. 1940. Vol. 7. P. 284-304. Doi: 10.1016/S0031-8914(40)90098-2.

2. Bьttiker M., Harris E. P., Landauer R. Thermal activation in extremely underdamped Josephson-junction circuits // Physical Review B. 1983. Vol. 28. P. 1268. Doi: 10.1103/PhysRevB.28.1268

3. Hanggi P., Talkner P., Borkovec M. Reaction-rate theory: fifty years after Kramers // Review of Modern Physics.1990. Vol. 62. P. 251-341. Doi: 10.1103/RevModPhys.62.251.

4. Gonchar I. I., Kosenko G. I. Is the Kramers formula applicable for describing the decay of highly excited nuclear systems? // Soviet Journal of Nuclear Physics. 1991. Vol. 53. P. 133-142.

5. Jiang Z., Smelyanskiy V. N., Isakov S. V. et al. Scaling analysis and instantons for thermally assisted tunneling and quantum Monte Carlo simulations // Physical Review A. 2017. Vol. 95. P. 012322. Doi: 10.1103/PhysRevA.95.012322.

6. Risken, H. The Fokker-Planck equation. Methods of Solution and Applications. Berlin: Springer. 1996. 472 p. Doi: 10.1007/978-3-642-61544-3.

7. Adeev G. D., Gontchar I. I., Pashkevich V. V. et al. Diffusion Model of Formation of Fission-Fragment Distributions // Soviet Journal of Particles and Nuclei. 1988. Vol. 19. P. 529.

8. Abe Y., Gregoire C., Delagrange H. Langevin approach to nuclear dissipative dynamics // Journal de Physique Colloques (Paris). 1986. Vol. 47. P. 329. Doi: 10.1051/jphyscol:1986436.

9. Frobrich P., Xu S.Y. The treatment of heavy-ion collisions by Langevin equations // Nuclear Physics A. 1988. Vol. 477. P. 143-161. Doi: 10.1016/0375-9474(88)90366-1.

10. Gontchar I. I., Chushnyakova M. V., Aktaev N. E. et al. Disentangling effects of potential shape in the fission rate of heated nuclei // Physical Review C. 2010. Vol. 82. P. 064606.

Doi: 10.1103/PhysRevC.82.064606.

11. Гончар И. И., Крохин С. Н. Прецизионные вычисления скорости деления возбуждённых атомных ядер // Вестник Омского Университета. 2012. Т. 4. С. 84-87.

12. Edholm O., Leimar O. The accuracy of Kramers' theory of chemical kinetics // Physica A. 1979. Vol. 98. P. 313-324. Doi: 10.1016/0378-4371(79)90182-1.

13. Турчак Л. И., Плотников П. В. Основы численных методов. М.: Физматлит. 2003. 304 с.

14. Kloeden P. E., Platen E., Schurz H. Numerical solution of SDE through computer experiments. Berlin: Springer-Verlag Berlin Heidelberg. 1994. 270 p. Doi: 10.1007/978-3-642-57913-4.

15. Marsaglia G. Xorshift Random Number Generators //Journal of Statistical Software. 2003.

Doi: 10.18637/jss.v008.i14.

Размещено на Allbest.ru