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

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

ФГУП "РОССИЙСКИЙ ФЕДЕРАЛЬНЫЙ ЯДЕРНЫЙ ЦЕНТР — ВНИИЭФ"
ВОПРОСЫ
АТОМНОЙ НАУКИ И ТЕХНИКИ
СЕРИЯ
Математическое моделирование
физических процессов
НАУЧНО-ТЕХНИЧЕСКИЙ СБОРНИК
ВЫПУСК 2
Саров — 2024
Издается с 1978 г.
С О Д Е Р Ж А Н И Е
Веселова Н. Н., Лебедев С. Н. Схема "Ромб" для решения уравнения теплопроводности на
квадратных адаптивно-встраиваемых сетках . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
Мжачих С. В., Лапшина Ю. Н. Трехэтапный алгоритм расчета локально комонотонного
кубического сплайна класса C1 с приближенным проецированием . . . . . . . . . . . . . . . . . . 18
Шестаков А. А., Старцев А. Н. Тестовая задача с аналитическим решением уравнений
переноса излучения и энергии в двумерном осесимметричном случае . . . . . . . . . . . . . . . . 28
Колесов Г. Н., Дубинов А. Е. Моделирование динамики сильноточного релятивистского
электронного пучка в камере дрейфа с неоднородным магнитным полем . . . . . . . . . . . . . . 35
Глазырин И. В., Ершова А. В., Михайлов Н. А. Учет спонтанных магнитных полей в трехмерной программе "Фокус"
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
Новиков К. А. Математическое моделирование переноса примеси в поверхностных водах в
расчетном коде GeRa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
Арапова Е. Ю., Тихонов А. В. Формирование слоя пирамид для перехода от гексаэдральных
к тетраэдральным сеткам
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
Сведения об авторах
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

C O N T E N T S
Veselova N. N., Lebedev S. N. A "Romb" scheme for solving the radiative heat transfer equation
on adaptive square meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
Mzhachikh S. V., Lapshina Yu. N. Three-step algorithm for calculating a locally comonotonic
cubic C1-class spline with approximate projection
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Shestakov A. A., Startsev A. N. A test problem with analytical solution of radiation and energy
transfer equations in two-dimensional axisymmetric case
. . . . . . . . . . . . . . . . . . . . . . . . 28
Kolesov G. N., Dubinov A. E. Modeling the dynamics of a high-current relativistic electron beam
in a drift chamber with an inhomogeneous magnetic field
. . . . . . . . . . . . . . . . . . . . . . . . 35
Glazyrin I. V., Ershova A. V., Mikhailov N. A. Accounting for spontaneous magnetic fields in
"Focus" 3D program . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
Novikov K. A. Mathematical modeling of impurity transport in surface waters using the GeRa
computation code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
Arapova E. Yu., Tikhonov A. V. Formation of a layer of pyramids for the transition from hexahedral to tetrahedral meshes
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
Information about authors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
Ответственный за выпуск Е. В. Соколовская
Редактор
Е. Н. Старченко
Корректоры
Т. А. Меркушева
Е. А. Окатьева
А. В. Федоренко
Н. Ю. Костюничева
Дата выхода в свет 25.06.2024
Формат 60×84/8
Усл. печ. л. ∼9, 9
Уч.-изд. л. ∼7, 1
Тираж 1000 экз.
Зак. тип. 903-2024
7 статей
Учредитель: ФГУП "РФЯЦ-ВНИИЭФ"
Свидетельство о регистрации ПИ №ФС77-29789 от 04 октября 2007 г.
выдано Роскомнадзором
Оригинал-макет подготовлен
в Математическом отделении ФГУП "РФЯЦ-ВНИИЭФ"
Отпечатано в ИПЦ ФГУП "РФЯЦ-ВНИИЭФ"
607188, г. Саров Нижегородской обл., ул. Силкина, 23
c⃝ФГУП "РФЯЦ-ВНИИЭФ", 2024

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

Н. Н. Веселова, С. Н. Лебедев
Постановка задачи
В двумерной области D рассматривается краевая задача для двумерного уравнения переноса
энергии в приближении лучистой теплопроводности [8]:
∂ε
∂t = −1
ρdiv ⃗W + f.
(1)
Здесь div ⃗W = 1
rd
∂rdWr
∂r
+ ∂rdWz
∂z

, где r, z — декартовы или цилиндрические координаты, при
этом для декартовых координат d = 0, для цилиндрических d = 1 (в этом случае ось Oz является осью симметрии); ⃗W = −κ∇T — вектор теплового потока, где κ = κ (ρ, T) — коэффициент
теплопроводности; ε (ρ, T) — удельная внутренняя энергия, где T (t, r, z) — температура, ρ (t, r, z) —
плотность; f (t, r, z) — источник.
Уравнение (1) рассматривается с начальными условиями
T
t0, r, z


= T 0 (r, z) .
(2)
На внешней границе L области D задаются граничные условия вида
αT + β ⃗W · ⃗n = ϕ (t,⃗r) ,
t > 0,
⃗r ∈L,
(3)
где ⃗n — внешняя нормаль к L. С помощью коэффициентов α, β и функции ϕ (t,⃗r) задаются различные граничные условия.
Для численного решения задачи (1)—(3) строится неявная конечно-разностная схема на основе
подхода "Ромб" [7]. Для аппроксимации распределения физических величин и их производных по
пространству используется квадратная адаптивная сетка.
Построение адаптивной сетки
Для пространственной аппроксимации задачи используется адаптивно-встраиваемая дробная счетная сетка. Это сетка, которая строится путем многоуровневого дробления основной сетки и может
быть адаптирована в процессе счета под особенности течения, конструкции и т. п.
Вся рассчитываемая система рассматривается как одна счетная область. Граница области разбивается на четыре участка Γ1, Γ2, Γ3, Γ4 (рис. 1, а) в порядке обхода контура области по часовой
стрелке.
Рис. 1. Начальная регулярная разностная сетка (а) и адаптивно-встраиваемая квадратная разностная сетка (б)
– 4 –

Схема "Ромб" для решения уравнения теплопроводности. . .
За основную сетку в области берется квадратная разностная сетка, имеющая регулярную структуру. Данная разностная сетка является начальной сеткой на каждом временном шаге. Процесс
построения адаптивно-встраиваемой дробной сетки производится в цикле итераций путем дробления имеющихся ячеек сетки на четыре части согласно установленному набору критериев. Каждая
вновь образованная ячейка имеет свой уровень дробления и информацию о номерах ячеек, с которыми она граничит.
Например, уровень ячейки, принадлежащей начальной сетке, полагается
нулевым; ячейки, образованные после одного дробления ячеек начальной сетки, имеют уровень 1
и т. д.
Построенная таким образом сетка имеет нерегулярную структуру (рис. 1, б).
Схема "Ромб" на адаптивной сетке
Аппроксимация уравнений. Счетная область покрывается основной квадратной регулярной начальной разностной сеткой. Границы
Γ1 и Γ3 разбиваются на I интервалов, а границы
Γ2 и Γ4 –– на J интервалов. В результате процесса адаптации получаем нерегулярную квадратную адаптивную разностную сетку, соответствующую некоторым критериям. Каждая ячейка этой сетки имеет свой номер, уровень дробления, список всех соседей по границам, координаты своих узлов (рис. 2).
Проведем построение разностной схемы для
уравнения (1).
Записывая уравнение теплопроводности в потоковой форме, интегрируя по ячейке адаптивной сетки и применяя неявную аппроксимацию
по времени, получаем конечно-разностную схему
Рис. 2. Ячейка сетки для аппроксимации уравнений
εn+1 −εn
τ
= −1
ρDIV(h) ⃗W n+1 + fn+1;
⃗W n+1 = −κn+1GRAD(h)T n+1,
(4)
где DIV(h) и GRAD(h) — разностные аналоги дифференциальных операторов.
Для решения системы (4) используется подход на основе схемы "Ромб".
Вводятся нормальные составляющие вектора теплового потока на ребрах ячейки (см. рис. 2):
Wξ1 = ⃗W · (−⃗n1), Wη2 = ⃗W · ⃗n2, Wξ3 = ⃗W · ⃗n3, Wη4 = ⃗W · (−⃗n4), а также нормальные компоненты
вектора ⃗W в центре ячейки: Wξ, Wη. Здесь |⃗n1| = |⃗n3| = ∆r и |⃗n2| = |⃗n4| = ∆z, где ∆r, ∆z ––
размеры ячейки в направлении r и z соответственно.
Разностный оператор дивергенции с учетом квадратной ячейки имеет следующий вид:
DIV(h) ⃗W =
1
∆V
Z
∆V
div ⃗WdV =
1
∆V
4
X
k=1

rd ⃗W · ⃗n

k ,
где ∆V — объем ячейки: ∆V = π (r3 + r1) ∆S при d = 1; ∆V = ∆S при d = 0; ∆S = (z3 −z1) (r2 −r1) —
площадь ячейки; rd
k = (rk + rk+1) /2 — радиус-вектор середины ребра.
Разностный оператор градиента GRAD(h)T получается при интегрировании grad T =
∂T
∂r , ∂T
∂z

по ячейке разностной сетки
1
∆S
Z
∆S
Z
grad TdS =
1
∆S
Z
⃝T⃗ndl, где ⃗n — единичная нормаль к границе
– 5 –

Н. Н. Веселова, С. Н. Лебедев
ячейки. В случае квадратной ячейки (рис. 2)
GRAD(h)T =
1
∆S

(T2 −T4) √α2, (T3 −T1) √α1

,
где ∆S — площадь ячейки; α1 = (r3 −r4)2 = ∆r2; α2 = (z3 −z4)2 = ∆z2.
Вводя простые итерации по нелинейности коэффициента теплопроводности κ
κµ = κ
ρn+1, T µ

,
где µ — итерационный индекс, получаем систему разностных уравнений
εµ+1 −εn
τ
+
1
∆m
h
(RWη)2 −(RWη)4 + (RWξ)3 −(RWξ)1
iµ+1
= fn+1,
W µ+1
ξ
+ κµ
∆S α1 (T3 −T1)µ+1 = 0,
(5)
W µ+1
η
+ κµ
∆S α2 (T2 −T4)µ+1 = 0,
где ∆m — масса ячейки; Rk — радиус-вектор середины ребра k
k = 1, 4


, Rk =
rk + rk+1
2
d
,
r5 = r1, d = 0, 1.
Все величины, записанные без нижнего индекса, относятся к центрам ячеек.
С помощью метода стабилизирующей поправки [9], итерации которого совмещены с итерациями
по нелинейности коэффициента теплопроводности [10], систему (5) сводим к двум одномерным системам. Счет на одной µ-итерации разбивается на два этапа. На первом этапе решается система
при j + 1/2 = const, где j = 0, J −1 — индекс основной регулярной разностной сетки вдоль границ
Γ2 и Γ4:
εµ+1/2 −εn
τ
+
1
∆m
h
(RWη)2 −(RWη)4
iµ+1/2
+
1
∆m
h
(RWξ)3 −(RWξ)1
iµ
= fn+1;
W µ+1/2
η
+ κµ
∆S α2 (T2 −T4)µ+1/2 = 0.
(6)
На втором этапе решается система при i + 1/2 = const, где i = 0, I −1 — индекс основной регулярной разностной сетки вдоль границ Γ1 и Γ3:
εµ+1 −εn
τ
+
1
∆m
h
(RWη)2 −(RWη)4
iµ+1/2
+
1
∆m
h
(RWξ)3 −(RWξ)1
iµ+1
= fn+1;
W µ+1
ξ
+ κµ
∆S α1 (T3 −T1)µ+1/2 = 0.
(7)
Решение систем разностных уравнений. Если бы разностная сетка имела регулярную структуру, то решение одномерных систем (6) и (7) свелось бы к стандартному для метода "Ромб" алгоритму [11], использующему потоковую прогонку [5].
Однако алгоритм в стандартном виде не
применим к нерегулярной сетке, а в рассматриваемом случае эти системы необходимо решать с
учетом адаптации основной сетки. Пример сеточной структуры, для которой необходимо решать
системы разностных уравнений, приведен на рис. 3.
Потоковая прогонка предполагает рекуррентный расчет прогоночных коэффициентов для прямого направления от ячейки к ячейке, начиная с левой границы, для обратного направления — от
правой границы к левой.
Напрямую использовать потоковую прогонку на адаптивной разностной сетке невозможно.
В
этом случае вдоль координатных линий основной сетки необходимо организовать передачу прогоночных коэффициентов от ячейки к ячейке.
– 6 –

Схема "Ромб" для решения уравнения теплопроводности. . .
Рис. 3. Пример структуры основной сетки на линии i + 1/2 = const
В соответствии с направлением решения проводится упорядочение ячеек. На рис. 4 приведен
пример упорядочения ячеек на двух этапах, где указаны номера ячеек для первого (см. рис. 4, а) и
второго (см. рис. 4, б) этапов.
Для упорядочения ячеек используется рекурсивная функция, в которой перебор ячеек осуществляется в порядке возрастания индекса основной сетки i для первого этапа и j — для второго этапа
от ячейки с нулевым уровнем по всем ее вложенным ячейкам. В массив упорядочения записывается
номер ячейки, принадлежащей данной адаптивной сетке. В результате получается массив упорядоченных ячеек с номерами k = 0, . . . , Kj −1, где Kj — число ячеек нерегулярной адаптивной сетки,
содержащихся в ячейках основной сетки вдоль выбранной координатной линии.
Рассмотрим решение на адаптивной сетке системы разностных уравнений (6) первого этапа.
Ввиду нелинейности удельной внутренней энергии систему (6) решаем итерационным методом
Ньютона. Линеаризуя внутреннюю энергию, получаем
 ∂ε
∂T
ν T ν+1
τ
+
1
∆m
h
(RWη)2 −(RWη)4
iν+1
+
1
∆m
h
(RWξ)3 −(RWξ)1
iµ
=
= εn −εν
τ
+
 ∂ε
∂T
ν T ν
τ + fn+1;
(8)
W ν+1
η
+ κµ
∆S α2 (T2 −T4)ν+1 = 0,
где T ν=0 ≡T µ.
Для замыкания системы (8) используются следующие соотношения [7], связывающие температуру
и компоненты вектора теплового потока в центре ячейки и на ее ребрах:
Рис. 4. Упорядочение ячеек на первом (а) и втором (б) этапах расчета уравнений
– 7 –

Н. Н. Веселова, С. Н. Лебедев
Tk =
Tk,2 + Tk,4
2
ν+1
+ δk

W ν+1
η,2
−W ν+1
η,4

k ;
W ν+1
η,k
= 1
2 (Wη,2 + Wη,4)ν+1
k
+ θk

T ν+1
k,2
−T ν+1
k,4

,
k = 0, . . . , Kj −1.
(9)
Здесь δ, θ — параметры схемы:
δν
k = 1
4
∆S
α2˜κν
k + εδ
;
θν
k = 0,
где εδ = 10−15 — малая константа, ограничивающая сверху значение параметра δ; ˜κ — аппроксимационный коэффициент теплопроводности, ˜κν
k = κµ
k +
Λκµ
k
8
1 + 2Kν
k

, Λ
κµ
k


= κн −2κµ
k + κв, Kν
k =
= (Rk,2 + Rk,4) κµ
kτα2
2∆S∆mεν
T
. В определении κн, κв участвуют все ячейки, которые граничат с ячейкой k
соответственно снизу и сверху.
К системе уравнений (8), (9) добавляются разностные граничные условия
αiT µ+1/2 + βiW µ+1/2
η
= ϕi,
i = 0, I.
(10)
Исключая неизвестные функции, определенные в центрах ячеек адаптивной сетки, имеем замкнутую систему разностных уравнений
(A0)ν
k T ν+1
k,4
+ (B0)ν
k T ν+1
k,2
+ (C0)ν
k (Wη)ν+1
k,4 + (D0)ν
k (Wη)ν+1
k,2 = (F0)ν
k ;
(A1)ν
k T ν+1
k,4
+ (B1)ν
k T ν+1
k,2
+ (C1)ν
i (Wη)ν+1
k,4 + (D1)ν
k (Wη)ν+1
k,2 = 0,
(11)
где
(C0)ν
k = −2δν
k −
2τ
∆mεν
T
Rk,4;
(D0)ν
k = 2δν
k +
2τ
∆mεν
T
Rk,2;
(A1)ν
k = −(B1)ν
k = −2˜κµα2
∆S ;
(A0)ν
k = (B0)ν
k = (C1)ν
k = (D1)ν
k = 1;
(F0)ν
k = 2
εν
T
εn −εν + εν
T T ν + τfn+1

−
2τ
∆mεν
T

R3W µ
ξ,3 −R1W µ
ξ,1

k ,
k = 0, . . . , Kj −1.
Система (11) решается с помощью одного из вариантов потоковой прогонки. Для этого предполагается связь неизвестных функций, справедливая для всех ребер ячеек сетки, в виде
Xk,4T n+1
k,4
+ Yk,4W n+1
k,4
= Zk,4.
(12)
Исключая из (11) и (12) неизвестные функции, приходим к формулам прямой прогонки
Xk,2 = (C1B0 −C0B1)ν
k ˜Xн + (B1A0 −B0A1)ν
k ˜Yн;
Yk,2 = (C1D0 −C0D1)ν
k ˜Xн + (D1A0 −D0A1)ν
k ˜Yн;
(13)
Zk,2 = (C1F0)ν
k ˜Xн −(F0A1)ν
k ˜Yн + (A1C0 −A0C1)ν
k ˜Zн,
k = 0, 1, . . . , Kj −1,
где ˜Xн =
Xн
|Xн| + |Yн|, ˜Yн =
Yн
|Xн| + |Yн|, ˜Zн =
Zн
|Xн| + |Yн| — прогоночные коэффициенты, которые
упорядоченная ячейка с номером k принимает от соседних ячеек через нижнее ребро. Если ребро
является граничным, то из граничного условия (10) следует: ˜Xн = αν
0; ˜Yн = βν
0; ˜Zн = ϕν
0.
Встречная прогонка ведется в обратном направлении.
Предполагая связь неизвестных функций в виде
Xk,2T n+1
k,2
+ Yk,2W n+1
k,2
= Zk,2
(14)
– 8 –