Вопросы атомной науки и техники. Серия: Математическое моделирование физических процессов , 2024, № 1

Г л а в н ы й
р е д а к т о р
Шагалиев Р. М.
З а м е с т и т е л и
г л а в н о г о
р е д а к т о р а :
Алексеев А. В., Тишкин В. Ф.
О т в е т с т в е н н ы й
с е к р е т а р ь :
Соколовская Е. В.
Ч л е н ы
р е д к о л л е г и и :
Бартенев Ю. Г., Бетелин В. Б., Бочков А. И., Вронский М. А.,
Дрёмов В. В., Залялов Н. Н., Кибзун А. И., Козелков А. С.,
Козманов М. Ю., Куркин А. А., Мартынов А. П., Петров И. Б.,
Прилуцкий М. Х., Смирнов Н. Н., Соколов С. С., Старостин Н. В.,
Степаненко С. А., Храмченков М. Г., Четверушкин Б. Н.,
Шестаков А. А., Янилкин Ю. В.
Адрес редакции и издателя: 607188,
г. Саров Нижегородской обл., пр. Мира, 37
тел.: (83130)28406, е-mail: sokol@vniief.ru.
Адрес сайта журнала: http://vant.vniief.ru/
c⃝ФГУП "РФЯЦ-ВНИИЭФ", 2023

ФГУП "РОССИЙСКИЙ ФЕДЕРАЛЬНЫЙ ЯДЕРНЫЙ ЦЕНТР — ВНИИЭФ"
ВОПРОСЫ
АТОМНОЙ НАУКИ И ТЕХНИКИ
СЕРИЯ
Математическое моделирование
физических процессов
НАУЧНО-ТЕХНИЧЕСКИЙ СБОРНИК
ВЫПУСК 1
Саров — 2024
Издается с 1978 г.
С О Д Е Р Ж А Н И Е
Мустафин А. М., Веселова Н. Н., Лебедев С. Н. Численное решение уравнения теплопроводности
в
двумерных
задачах
с
фазовыми
переходами
на
адаптивно-встраиваемых
сетках . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
Залялов А. Н., Иванов Н. В., Широков А. Е. Стохастическая методика численного моделирования радиационных поясов Земли
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Титова В. Б., Володина Н. А., Ширшова М. О., Халдеев Е. В., Сустаева Ю. М. Особенности
распространения детонации в каналах малого сечения для состава на основе тэна . . . . . . . . 30
Кошутин Д. А., Шестаков А. А. Учет кинетических эффектов в приближении лучистой
теплопроводности для расчета задач переноса излучения . . . . . . . . . . . . . . . . . . . . . . . 37
Титова В. Б., Володина Н. А., Ширшова М. О., Кирюхина М. Н., Богданов Е. Н., Становов А. А. Верификация уравнений состояния продуктов взрыва пластифицированного
октогена на результатах опытов по разгону лайнеров . . . . . . . . . . . . . . . . . . . . . . . . . 42
Шестаков А. А. 50 лет задачам Флека . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
Бугаенко А. А., Крутова И. Ю. Численное моделирование на основе линеаризованной системы
уравнений
газовой
динамики
восходящих
закрученных
потоков
с
учетом
силы
Кориолиса . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
Данилов А. С., Гордеев Д. Г., Шумилина О. Н., Арапов И. Н. Алгоритмы и особенности
программной реализации в пакете УРС-ОФ полуэмпирической модели широкодиапазонного
уравнения состояния РОСА-МФИ
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
Сведения об авторах
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

C O N T E N T S
Mustafin A. M., Veselova N. N., Lebedev S. N. Numerical solution of the radiative heat transfer
equation in 2D problems with phase transitions on adaptive meshes . . . . . . . . . . . . . . . . . . .
4
Zalialov A. N., Ivanov N. V., Shirokov A. E. Stochastic method of numerical simulation of radiation belts of the Earth
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Titova V. B., Volodina N. A., Shirshova M. O., Khaldeev E. V., Sustaeva Y. M. Specific features
of detonation propagation in small cross-section channels for the PETN-based composition . . . . . . 30
Koshutin D. A., Shestakov A. A. Account for kinetic effects in approximation of radiating heat
transfer to compute radiation transfer problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
Titova V. B., Volodina N. A., Shirshova M. O., Kityukhina M. N., Bogdanov E. N., Stanovov A. A.
Verification of the equations of state fir the explosion products of the plasticized hmx using the
experimental results on liners acceleration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
Shestakov A. A. 50 years of fleck problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
Bugaenko A. A., Krutova I. Yu. Numerical Simulation on the basis of linearized system of gasdynamic equations of rising swirl flows with the account for Сoriolis forces . . . . . . . . . . . . . . . 57
Danilov A. S., Gordeev D. G., Shumilina O. N., Arapov I. N. Algorithms and specific features
of program realization of semi-empirical model of wide-range ROSA-MFI equation of state in
URS-OF package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
Information about authors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
Ответственный за выпуск Е. В. Соколовская
Редактор
Е. Н. Старченко
Корректоры
Т. А. Меркушева
Е. А. Окатьева
А. В. Федоренко
Н. Ю. Костюничева
Дата выхода в свет 25.03.2024
Формат 60×84/8
Усл. печ. л. ∼9, 1
Уч.-изд. л. ∼8, 4
Тираж 1000 экз.
Зак. тип. 108-2024
8 статей
Учредитель: ФГУП "РФЯЦ-ВНИИЭФ"
Свидетельство о регистрации ПИ №ФС77-29789 от 04 октября 2007 г.
выдано Роскомнадзором
Оригинал-макет подготовлен
в Математическом отделении ФГУП "РФЯЦ-ВНИИЭФ"
Отпечатано в ИПЦ ФГУП "РФЯЦ-ВНИИЭФ"
607188, г. Саров Нижегородской обл., ул. Силкина, 23
c⃝ФГУП "РФЯЦ-ВНИИЭФ", 2024

В настоящем выпуске публикуются тексты отдельных докладов, представленных на XVI Международной конференции "Забабахинские научные чтения",
которая проходила в г. Снежинске Челябинской области с 30 мая по 1 июня
2023 года.

ВАНТ, сер. Математическое моделирование физических процессов. 2024. Вып. 1
УДК
519.63
ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ
В ДВУМЕРНЫХ ЗАДАЧАХ С ФАЗОВЫМИ ПЕРЕХОДАМИ
НА АДАПТИВНО-ВСТРАИВАЕМЫХ СЕТКАХ
А. М. Мустафин, Н. Н. Веселова, С. Н. Лебедев
(ФГУП "РФЯЦ-ВНИИТФ им. академика Е. И. Забабахина",
г. Снежинск Челябинской области)
Предлагается методика численного решения двумерного уравнения теплопроводности с учетом фазовых переходов на адаптивной сетке с использованием схемы "Ромб",
включающая построение адаптивной разностной сетки, и исследуется эффективность
применения разработанной методики. Для решения задач о распространении тепла при
фазовых переходах, называемых задачами типа Стефана, используется метод сквозного счета на основе применения разностного аналога δ-функции. Дано описание схемы
"Ромб" на адаптивно-встраиваемой сетке. Приведены численные решения задач о распространении тепловой волны в неоднородной области и о воздействии теплового пучка на лед. На задаче с тепловой волной наглядно показаны преимущества применения
алгоритмов адаптации в области фронта тепловой волны в двумерном пространстве.
Численные результаты как первой, так и второй задачи наглядно демонстрируют преимущества применения адаптивно-встраиваемых сеток совместно с методом сквозного
счета по сравнению с использованием традиционных подходов в случае решения задач
типа Стефана.
Ключевые слова: задача Стефана, разностная схема "Ромб", адаптивная разностная
сетка, двумерное уравнение теплопроводности, фазовый переход, тепловая волна.
Введение
В работе исследуется обобщение одного из подходов к численному решению уравнения теплопроводности в трехфазных средах с двумя нестационарными межфазовыми границами. Задачи такого
типа принято называть задачами Стефана [1]. Рассматривается общий случай динамики температурных полей твердого вещества в результате нагрева. Сложность построения вычислительного
алгоритма для задачи Стефана заключается в том, что граница раздела фазовых состояний подвижна. Существует множество методов решения данной проблемы [2]: метод энтальпий, метод ловли
фронта фазового перехода в узле сетки, квазистационарный подход, метод выпрямления фронта,
использование δ-функции и сглаживание коэффициентов [3]. Однако далеко не все они применимы
в случае многомерных задач.
Данная работа рассматривает метод сквозного счета со сведением нескольких уравнений для разных фаз к одному общему уравнению теплопроводности посредством введения одного из аналогов
δ-функции. Точность выбранного метода повышается за счет использования вместо классической
разностной сетки адаптивной [4, 5], предполагающей динамическое сгущение сетки в окрестности
фронта фазового перехода.
Ранее был предложен подход к решению задачи Стефана в одномерном случае, основанный на
методе сквозного счета и применении схемы "Ромб" [6, 7] к решению уравнений. Исследования [8]
показали высокую эффективность такого подхода к решению одномерных задач с фазовыми переходами.
– 4 –

Численное решение уравнения теплопроводности в двумерных задачах. . .
В настоящей работе данная методика обобщается на случай двумерных задач. Для этого строится
двумерная разностная схема "Ромб" с учетом особенностей адаптивной разностной сетки, связанных, прежде всего, с ее нерегулярностью в двумерном пространстве. Рассматриваются две задачи:
о распространении тепловой волны в неоднородной области и о воздействии теплового пучка на лед.
Для каждой из них приводятся полученные численные результаты и их анализ.
Постановка двумерной задачи
Вещество, заполняющее область G в двумерном пространстве XOY (рис. 1), в начальный момент времени находится в твердом состоянии при некоторой постоянной температуре T0.
На одной из границ области устанавливается мощный источник тепла, например, с температурой Th, причем Th > Ts > Tf,
где Ts, Tf — соответственно температура кипения (испарения) и плавления (кристаллизации) данного вещества. В общем случае температура Th является функцией от времени и
координат границы.
В поставленных условиях вещество начнет плавиться и в последующем испаряться.
Обозначим границу плавления как Гf и при
Рис. 1. Схематичное изображение задачи
этом заметим, что температура в точках, принадлежащих границе, всегда будет равна температуре
плавления Tf. Границу испарения обозначим Гs с условием, что температура в точках, принадлежащих границе, всегда будет равна температуре кипения Ts.
Таким образом, рассматриваемая часть вещества имеет три зоны с соответствующими фазами
состояния — твердым, жидким и газообразным, а также две границы фазовых переходов — плавления и испарения. Переход вещества из твердого в жидкое состояние на границе осуществляется
при затрачивании удельной теплоты Lf. При испарении затраты соответствуют удельной теплоте парообразования Ls. В текущей постановке предполагается, что перенос тепла осуществляется
только за счет теплопроводности, а движением вещества пренебрегается.
Процесс распространения тепла в веществе описывается дифференциальным уравнением вида
ρi
∂εi
∂t = div (χi grad (T)) + f,
(1)
где i = {s, l, g} — индекс состояния вещества (твердое, жидкое, газообразное); T = T (x, y, t) —
температура; εi (ρi, T) — удельная внутренняя энергия; ρi — плотность; χi (ρi, T) — коэффициент
теплопроводности; f (x, y, t) — плотность тепловых источников.
На границах фазовых переходов (плавления и испарения) вводятся соответственно условия Стефана [1] для потоков:
χl grad (T)|Гf−0 −χs grad (T)|Гf+0 = Lfρ∂ξ1
∂t ;
χg grad (T)|Гs−0 −χl grad (T)|Гs+0 = Lsρ∂ξ2
∂t ,
(2)
где ξ1 (x, y, t), ξ2 (x, y, t) — функции, описывающие изменение положения межфазовых границ с
течением времени.
Таким образом, для каждой зоны со своим фазовым состоянием вещества решается уравнение со
своими коэффициентами. Согласно методике, ранее разработанной для одномерных трехфазных задач теплопроводности [8], суммарную энергию, включающую в себя энергию фазовых превращений,
можно записать с помощью функции Хевисайда ℏ(T):
– 5 –

А. М. Мустафин, Н. Н. Веселова, С. Н. Лебедев
˜ε (T) = ε + Lfℏ(T −Tf) + Lsℏ(T −Ts) ,
ℏ(T) =



0,
T < 0;
1,
T > 0,
где ε — некая обобщенная функция, определяющая внутреннюю энергию в зависимости от фазового
состояния вещества.
При подстановке суммарной энергии в исходное уравнение (1) функция Хевисайда обращается в
δ-функцию Дирака. Таким образом, уравнение, описывающее трехфазную задачу Стефана (1), (2),
приобретает вид
ρ eC (T) ∂T
∂t = div (χ grad (T)) + f,
eC (T) = ∂ε
∂T + Lfδ (T −Tf) + Lsδ (T −Ts) ,
(3)
где ρ = ρ (T) — функция, задающая плотность вещества; χ = χ (T) — функция, задающая коэффициент теплопроводности. В случае линейной зависимости энергии от температуры внутренняя
энергия определена выражением ε (T) = CV T, в котором CV = CV (T) — функция, определяющая
коэффициент удельной теплоемкости.
Уравнение (3) является общим для всей области определения, объединяющим в себе уравнения
для подобластей с разными фазами состояния вещества [3, 9].
Для рассматриваемого вещества
вводятся разрывные функции для плотности, коэффициентов теплоемкости и коэффициентов теплопроводности, так как в зависимости от фазы они меняют свое значение. Поскольку смена фазы
определяется температурой, то и эти величины будут зависеть от интервала температур для каждой
фазы. То есть
CV (T) =









Cs,
T < Tf;
Cl,
Tf ≤T ≤Ts;
Cg,
T > Ts;
χ (T) =









χs,
T < Tf;
χl,
Tf ≤T ≤Ts;
χg,
T > Ts;
ρ (T) =









ρs,
T < Tf;
ρl,
Tf ≤T ≤Ts;
ρg,
T > Ts,
где Cs, Cl, Cg — удельные коэффициенты теплоемкости; χs, χl, χg — удельные коэффициенты
теплопроводности; ρs, ρl, ρg — плотности соответственно для твердого, жидкого и газообразного
состояний исследуемого вещества.
В численных расчетах δ-функция Дирака заменяется аналогом (аппроксимацией) вида
δ (T) =
1
a√πe−(T/a)2,
где a = ∆/ (2√π), ∆— параметр, определяющий полуширину "колокола" функции [10] на оси
температур.
В выбранном подходе предлагается от разрывных функций характеристик вещества, зависящих
от температурных интервалов, перейти к непрерывным гладким функциям. В области действия
δ-функции (T ≤|∆|) производится сглаживание перехода от одного значения к другому в разрывных
функциях теплоемкости, теплопроводности и плотности. Вводится функция сглаживания
γ (T) = 1
2 cos
π (T + ∆)
2∆

+ 1
2.
Для каждого фазового перехода определяется свой температурный диапазон воздействия функции сглаживания и δ-функции. В случае трехфазной задачи задаются диапазоны ∆f и ∆s — плавления и испарения соответственно.
Введя понятие теплового потока согласно закону Фурье [11], запишем уравнение (3) в потоковой
форме:
∂T
∂t = −
1
ρ eC (T)
div

⃗W

+
1
ρ eC (T)
f,
⃗W = −χ grad (T) .
(4)
– 6 –

Численное решение уравнения теплопроводности в двумерных задачах. . .
В области G = [ax ≤x ≤bx] × [ay ≤y ≤by] при t ∈[0, tmax] для уравнения (4) решается начальнокраевая задача
T (x, y, 0) = ϕ (x, y) ,
(x, y) ∈G;
αT + β

⃗W · ⃗n

= ψ (x, y, t) ,
(x, y) ∈Γ,
где Γ — внешняя граница прямоугольной области G; ϕ (x, y) — заданное начальное распределение
температуры; ψ (x, y, t) — заданная функция на границах области; α и β — параметры, с помощью
которых можно задавать граничные условия различного типа; ⃗n — внешняя нормаль к границе Γ.
Адаптивная сетка
Пусть на плоскости XOY определена равномерная разностная сетка с шагами hx (по оси OX) и
hy (по оси OY )
ωhxhy = ωhx × ωhy,
ωhx = {xj = jhx, j = 0, 1, 2, . . . , J} ,
ωhy = {yi = ihy, i = 0, 1, 2, . . . , I} ,
расширенная во временную область с некоторым шагом по времени τ:
ωhxhyτ = ωhx × ωhy × ωτ,
ωτ = {tn = nτ, n = 0, 1, 2, . . .} .
Взамен такой исходной сетки вводится адаптивная сетка [12], схематично изображенная на рис. 2.
Данная сетка является нерегулярной и не может быть записана в привычном матричном виде.
Пусть количество ячеек исходной регулярной сетки составляет IJ.
Процесс адаптации исходной сетки имеет несколько (L) уровней. При переходе на каждый следующий уровень происходит
дробление ячеек сетки, удовлетворяющих условию адаптации, на четыре ячейки меньшего размера. Пусть Pl−1 — число ячеек, подвергшихся дроблению при переходе с (l −1)-го на l-й уровень
(l = 1, 2, . . . , L). Тогда число ячеек, участвующих в счете после L уровней адаптации, K = IJ +
+
L
X
l=1
3Pl−1. Геометрические размеры k-й ячейки (k = 0, 1, 2, . . . , K −1) описываются формулами

˜hx

k = hx
2Lk ;

˜hy

k = hy
2Lk ,
(5)
где hx и hy — пространственные шаги исходной регулярной сетки; Lk — уровень адаптации k-й
ячейки.
Решение разностных уравнений осуществляется шагами по времени. Перед расчетом
следующего шага имеем поля температуры и
плотности, определенные на некоторой сетке
n-го временного слоя. Для следующего шага
(n + 1) по времени строится новая адаптивная
сетка, на которую отображаются температура и плотность. Исходные или уже раздробленные ячейки, удовлетворяющие условиям
адаптации, делятся.
Раздробленные ячейки
с предыдущего временного слоя, не удовлетворяющие условию адаптации, объединяются. В таком случае температура определяется
по суммарной энергии ячеек, образующих исходную ячейку.
Рис. 2. Адаптивная разностная сетка
– 7 –

А. М. Мустафин, Н. Н. Веселова, С. Н. Лебедев
Условием адаптации при переходе с нулевого уровня на первый выступает попадание центра ячейки в ε-окрестность заданной точки: Oε (D) =

(x, y) :
q
(x −xD)2 + (y −yD)2 < ε

. Для перехода
на второй уровень дробления центр ячейки должен попадать Oε/2 (D). В данной работе точка D
принадлежит множеству точек границ фазовых переходов.
Также иногда удобно использовать в качестве условия адаптации минимальное пороговое значение
градиента температур Tgrad, превышение которого означает, что ячейка подлежит дроблению. При
переходе на каждый следующий уровень адаптации пороговое значение увеличивается в два раза.
Схема "Ромб" на адаптивной сетке
Рис. 3.
Аппроксимация уравнения в ячейке
сетки
Аппроксимация системы уравнений в пределах одной ячейки. Для решения системы уравнений (4) используем схему "Ромб" [7].
Заменив производные по времени и пространству на первые разностные производные, запишем систему в разностном
виде (индексы центров ячеек (j + 1/2, i + 1/2) здесь и
далее опущены):
T n+1 −T n
τ
= −1
ρ eC
DIV(h)

⃗W n+1
+ 1
ρ eC
fn+1;
⃗W n+1 = −χn+1GRAD(h)
T n+1

,
(6)
где DIV(h) и GRAD(h) — разностные аналоги дифференциальных операторов div и grad соответственно.
Оператор дивергенции получается при интегрировании div

⃗W

по ячейке (рис. 3) разностной сетки
(согласно теореме Гаусса—Остроградского):
DIV(h)

⃗W

=
1
∆S
4
X
k=1

⃗W · ⃗nk

.
Оператор градиента GRAD(h) (T) определяется путем интегрирования grad (T) =
∂T
∂x , ∂T
∂y

по
площади ячейки сетки:
1
∆S
Z
∆S
Z ∂T
∂x dS =
1
∆S
I
γ
T⃗nxdl =
1
∆S

T2˜h2
y −T4˜h2
y

=
1
∆S (T2 −T4) ˜h2
y;
1
∆S
Z
∆S
Z ∂T
∂y dS =
1
∆S
I
γ
T⃗nydl =
1
∆S

−T1˜h2
x + T3˜h2
x

=
1
∆S (T3 −T1) ˜h2
x.
В случаях, когда имеет место нелинейность коэффициента теплопроводности, вводятся простые
итерации по индексу µ. В результате с учетом замены m = ρ∆S (масса ячейки) система (6) примет
вид
T µ+1 −T n
τ
= −1
eCm
(Wy3 −Wy1 + Wx2 −Wx4)µ+1 + 1
eCρ
fn+1;
W µ+1
x
+ χµ
x
˜h2
y
∆S (T2 −T4)µ+1 = 0;
W µ+1
y
+ χµ
y
eh2
x
∆S (T3 −T1)µ+1 = 0,
(7)
где введены нормальные составляющие вектора теплового потока на ребрах ячейки Wy1 = ⃗W ·(−⃗n1),
Wy3 = ⃗W ·⃗n3, Wx2 = ⃗W ·⃗n2, Wx4 = ⃗W ·(−⃗n4) и нормальные компоненты Wx, Wy вектора ⃗W в центре
– 8 –