Доклады Российской академии наук. Математика, информатика, процессы управления, 2023, T. 512, № 1, стр. 89-95

1. ВВЕДЕНИЕ

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

Диффузионное приближение широко применяется в задачах высокотемпературной газовой динамики для описания вклада излучения [1, 2]. Его использование позволяет понизить размерность решаемой задачи о лучистом переносе энергии за счет того, что во многих практически важных случаях позволяет эффективно определить поток W и плотность энергии излучения U,

где c – скорость света, не прибегая к дорогостоящему расчету интенсивности энергии излучения I в зависимости от направления полета фотонов Ω. Например, в задачах, связанных с генерацией короткоживущей плазмы под действием импульсного источника энергии (лазер, электроразрядный генератор) на конденсированное вещество, анизотропия излучения в основном заметна в начальной стадии процесса, в сравнительно узкой зоне вблизи фронта излучения, и мало влияет на следующее за ним течение плазмы. В подобных ситуациях диффузионное приближение дает вполне приемлемое качественное и количественное описание лучистого теплообмена.

Диффузионная модель переноса энергии излучением включает следующие уравнения [1–3]. Все уравнения записаны в системе единиц из работы [1].

где U – плотность энергии излучения; T – температура среды; ${{U}_{p}} = \frac{{4\sigma {{T}^{4}}}}{c}$ – энергия равновесного излучения тела при локальной температуре материала T (σ – постоянная излучения Cтефана-Больцмана); æ = 1/l – осредненный по спектру коэффициент поглощения; l – длина свободного пробега фотона, t – время, с – скорость света; $W$ – поток энергии излучения.

При этом вклад лучистой энергии в удельную внутреннюю энергию среды Е описывается уравнением

Начальные и граничные условия системы (1.1)– (1.3) задаются для конкретной задачи. Несмотря на широкое распространение в практике решений задач динамики излучающего газа [3] использование диффузионного приближения является не простой алгоритмической проблемой.

Аппроксимации системы уравнений (1.1)–(1.2) посредством неявных схем приводят к системам сеточных уравнений, которые необходимо решать с помощью тех или иных итерационных схем. Соответствующие алгоритмы в параллельной реализации показывают резкое снижение эффективности параллельных вычислений, особенно при использовании систем с экстрамассивным параллелизмом [4]. Локально-итерационные методы решения сеточных систем уравнений имеют хорошую перспективу, однако, несмотря на достигнутый в последние годы прогресс в их разработке (см., напр., [5]), требуются дальнейшие исследования по их адаптации (в частности, повышение эффективности при обменах данными в подходах типа геометрического параллелизма) для расчетов на перспективных гетерогенных вычислительных системах, в которых в качестве ускорителей используются графические платы. Большой проблемой при моделировании нелинейных процессов, зачастую снижающей эффективность использования итерационных методик в системах распределенных вычислений, является разбалансировка нагрузки компьютера.

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

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

Модифицированное уравнение диффузии и алгоритм его решения рассмотрим в следующем параграфе.

2. МОДИФИЦИРОВАННОЕ УРАВНЕНИЕ ДИФФУЗИИ

Уравнения диффузионной модели (1.1)–(1.2) в приближении $\frac{{\partial \;U}}{{\partial \;t}} = 0$ описывают квазистационарное поле излучения, которое соответствует “практически мгновенному” установлению лучеиспускательной способности и коэффициента поглощения излучения, зависящих от температуры и плотности вещества, которые, в свою очередь, могут зависеть от времени. Это приближение справедливо, когда характерное время изменения состояния вещества в расчетной области (время гидродинамических процессов) много больше времени прохождения этой области светом. В дальнейшем изложении не рассматривается зависимость излучательной способности и коэффициента поглощения от частоты излучения. В случае необходимости ее учета все полученные результаты легко обобщаются для многогрупповой по спектру модели [3].

Остановимся на физических и математических аспектах построения диффузионной модели (1.1)–(1.2).

Уравнение для одномерного, зависящего от времени переноса излучения имеет вид [1, 2]:

где I – удельная интенсивность излучения, z – пространственная координата, μ – косинус угла относительно оси z.

В прикладных расчетах [1–3] широко используются приближения, основанные на угловых моментах интенсивности излучения (2.1). Первые два угловых момента I для плоского слоя определяются как:

где U – плотность энергии излучения, W – поток излучения. Нулевой момент уравнения (2.1) удовлетворяет уравнению [1, 2]: которое является выражением энергетического баланса поля излучения.

Первый момент уравнения (2.1) находится путем умножения его на μ, и интегрирования по полному телесному углу:

Здесь появляется третий угловой момент, P = = $\frac{{2\pi }}{c}\int\limits_{ - 1}^{ + 1} {d\mu \mathop \mu \nolimits^2 I\left( \mu \right)} $, который представляет собой давление лучистого потока. Соответственно, уравнение (2.3) выражает баланс импульса фотонов.

Учитывая уравнение баланса материальной энергии (1.3) для определения T, уравнения (2.2) и (2.3) имеют три неизвестных: U, W и P. Известный и самый простой в реализации способ замыкания получается в предположении зависимости давления излучения от плотности энергии:

Уравнение (2.3) при этом преобразуется к виду:

Уравнение энергетического баланса (2.2) остается неизменным. Уравнения (2.2), (2.4) составляют основу так называемого приближения P₁. В работе [6] было показано, что в пределе оптически толстой среды решение в модели P₁ имеет те же разложения нулевого и первого порядка, что и решение исходного уравнения переноса.

Уравнения диффузионной модели (1.1)–(1.2) соответствуют приближению P₁, если в (2.4) пренебречь производной по времени от потока излучения W. Нестационарность поля излучения в диффузионном приближении учитывается членом $\frac{{\partial U}}{{\partial t}}$ в (2.3). Исключая из (1.1)–(1.2) поток W, получим нестационарное уравнение диффузии вида [2]:

Здесь второе слагаемое записано в дивергентной форме, т.к. коэффициент æ в общем случае может зависеть от температуры и плотности среды и, тем самым, от пространственной координаты. Начальные и граничные условия задаются для конкретной задачи. Уравнение (2.5) является уравнением параболического типа. С целью построения эффективного вычислительного алгоритма заменим (2.5) на уравнение гиперболического типа с малым параметром $\tau \ll 1$, имеющим размерность времени, при второй производной $\frac{{{{\partial }^{2}}U}}{{\partial {{t}^{2}}}}$:

По форме уравнение (2.6) совпадает с предложенным и исследованным в [8, 9] и других работах гиперболизированным квазилинейным параболическим уравнением. Поэтому для его численного решения можно применять явную трехслойную разностную схему, удовлетворяющую общим условиям устойчивости [7]. Такая схема предложена и изучена в работах [8, 9].

Малый параметр при второй производной τ в данном случае можно получить из системы (2.2), (2.4), отбросив производные по времени от медленно меняющихся параметров U_p и æ, которые зависят от термодинамического состояния среды. Параметр τ имеет ясный физический смысл – это время прохождения светом расстояния порядка длины свободного пробега фотона τ = 1/(сæ) = l/с. Однако в сеточной модели минимальное расстояние равно шагу сетки и τ = h/c (h = δl, δ ~ 10^–2). Поэтому с целью стабилизации схемы и увеличения допустимого шага интегрирования по времени предлагается увеличить значение этого параметра, или, что эквивалентно, уменьшить эффективную скорость света с → αс, 0 < α < 1. Значение α следует выбирать таким образом, чтобы, прежде всего, выполнялось условие устойчивости. Ограничение на параметр α по устойчивости имеет и другой аспект. При моделировании радиационно-гидродинамических процессов корректирующий фактор α должен подчиняться физическому ограничению: модифицированная скорость света с' = αс должна быть много больше характерных скоростей гидродинамических процессов, учитываемых в конкретной радиационно-гидродинамической модели. Поэтому обоснованный выбор значения параметра α является основной проблемой рассматриваемой нами модификации параболического или эллиптического уравнения: он должен обеспечивать значительный вычислительный эффект, позволяя выполнять расчет по трехслойной явной схеме с разумным с вычислительной точки зрения шагом по времени.

Выпишем аппроксимирующую (2.6) трехслойную явную схему:

Здесь Δt – шаг по времени при решении эволюционной радиационно-газодинамической задачи; ${c \mathord{\left/ {\vphantom {c 3}} \right. \kern-0em} 3}\mathop {({1 \mathord{\left/ {\vphantom {1 \unicode{230} }} \right. \kern-0em} \unicode{230} }\,U{{x}_{k}})}\nolimits^j {{{\mathbf{x}}}_{k}}$ – разностная аппроксимация дивергенции потока cdiv l/3 grad U, k – индекс пространственной координаты; $U_{i}^{j}$ – значение плотности энергии излучения U на j-слое по времени в центральном элементе разностного шаблона.

Можно предложить следующую оценку, подтвержденную рядом вычислительных экспериментов

где h – пространственный размер ячейки сетки, V – характерная скорость процесса.

Из общей теории устойчивости разностных схем [7] можно получить критерий устойчивости трехслойной явной схемы типа (2.7), который имеет вид [8]

что значительно более приемлемо, чем условие устойчивости явной двухслойной схемы для параболического уравнения [7]

Данное различие особенно заметно проявляется при использовании разностных сеток большого размера. Именно для такого случая подробной пространственной сетки, требующей для своей обработки больших вычислительных мощностей, используются высокопроизводительные вычислительные системы.

Как показывает теоретический анализ линеаризованной задачи [10], при малом значении τ различие между гиперболизированной и исходной системами незначительно. Однако из опыта решения больших задач, имеющих реальный физический смысл, известно, что оценки, полученные в линейном приближении, не всегда выполняются.

В следующем параграфе рассмотрим результаты вычислительного эксперимента применительно к решению задачи радиационной газовой динамики, связанной с выбором параметра τ и сравнением решения гиперболизированного (2.8) и исходного (2.7) уравнений диффузии.

3. ВЫЧИСЛИТЕЛЬНЫЙ ЭКСПЕРИМЕНТ

Вычислительные эксперименты проводились на одномерных задачах, поскольку явная схема для многомерного случая будет отличаться только способом аппроксимации пространственных производных, и ее реализация не вызывает принципиальных трудностей. Для выбора параметра τ следует использовать линейный размер ячейки пространственной сетки h.

Для тестирования разработанной методики использовались точные решения из работ [11], [12]. В этих работах приведены решения системы уравнений энергии и диффузии излучения в плоском одномерном случае без учета рассеяния

при специально выбранных начальных и граничных условиях, а также с использованием уравнения состояния и коэффициента поглощения вида:

В каждом случае задача решалась с помощью как исходной (2.5), так и гиперболизированной системы (2.6). Для решения гиперболизированного уравнения использовалась трехслойная схема (2.7) .

На рис. 1 и 2 приведено решение для двухобластной задачи с разрывным коэффициентом поглощения, рассмотренной в ([11], задача 4).

Решение системы (1.1)–(1.3) при данных начальных и граничных условиях:

Расчет выполнялся на отрезке [0, 2], шаг пространственной сетки h = 0.01. На графиках представлены рассчитанные профили плотности энергии излучения U и внутренней энергии Е (сплошные линии, левая шкала) в системе единиц [11], где с = 3, при двух значениях малого параметра при второй производной: τ = h/c и τ = = h/(0.01c) на момент времени t_max = 1. Соответствующие этим значениям графики 1 и 2 на обоих рисунках визуально не отличимы друг от друга и от точного решения, полученного в [11]. Пунктирными линиями показана погрешность разностного решения в норме С, им соответствует правая шкала. Относительная погрешность для U и для Е при τ = h/c в равномерной норме имеет порядок 10^–4, см. табл. 1. Погрешность в интегральной норме на порядок меньше. При увеличении τ в 100 раз погрешность увеличилась в ~10 раз. Тем самым экспериментально подтверждены оценки (2.8) и (2.9) для явных схем решения соответственно систем (2.5) и (2.6). При этом шаг по времени, выбранный из условия (2.8), составлял 1.4 ×10^–4 и 1.4 × 10^–3 соответственно, тогда как допустимый шаг по времени для явной схемы (2.9) для исходного параболического уравнения Δt_max = h²/(2æ_max) = 2.5 × 10^–5. Таким образом, благодаря использованию модифицированного уравнения диффузии и трехслойной разностной схемы было достигнуто увеличение шага интегрирования по времени более чем в 50 раз при сохранении высокой точности разностного решения.

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

4. ЗАКЛЮЧЕНИЕ

Метод гиперболизации уравнения диффузионного переноса излучения имеет хорошую перспективу в отношении решения задач радиационной газовой динамики на современных и перспективных многопроцессорных вычислительных системах. Тестовые решения показали эффективность гиперболизации для задач в приближении “серого” вещества, однако при соответствующих доработках, обусловленных необходимостью детализации выбора стабилизирующего фактора в зависимости от спектрального состава излучения, метод вполне применим в многогрупповом спектральном приближении. Также стабилизированные трехслойные схемы для модели диффузионного переноса будут полезны в вычислениях на многомерных пространственных сетках сложной структуры, таких как, например, неструктурированные или локально-адаптивные иерархические сетки.

Важным направлением дальнейшего развития методики является исследование возможности выполнения расчетов с переменным шагом по времени и/или переменным малым параметром при второй производной по времени. Это позволит существенно повысить эффективность вычислений при сильных изменениях параметров среды и при использовании динамически перестраиваемых сеток.

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

Как уже отмечалось, данный подход позволяет существенно сократить время расчета. Это является актуальным для решения многомерных задач радиационной газовой динамики с учетом спектральной зависимости поля излучения. Именно моделирование таких процессов является наиболее ресурсоемким и требует привлечения для своего решения вычислительных систем высокой и сверхвысокой производительности.