Эффективные численные методы и суперкомпьютерные алгоритмы для математического моделирования динамики многокомпонентных сплошных сред тема диссертации и автореферата по ВАК РФ 00.00.00, доктор наук Храпов Сергей Сергеевич

• Храпов Сергей Сергеевич
• доктор наук
• 2026

Оглавление диссертации доктор наук Храпов Сергей Сергеевич

Введение

Глава 1. Численный метод CSPH-TVD для моделирования нестационарных течений мелкой воды

1.1. Основные уравнения модели мелкой воды

1.1.1. Рельеф местности

1.1.2. Силы действующие на жидкость в рамках модели мелкой воды

1.1.3. Источники и стоки воды

1.2. Основные этапы СБРЫ-ТУЭ метода

1.2.1. Лагранжев этап — модифицированный БРЫ-подход

1.2.2. Эйлеров этап — модифицированный ТУЭ-подход

1.2.3. Условие устойчивости СБРЫ-ТУЭ метода

1.3. Результаты тестирования численной схемы СБРЫ-ТУЭ

1.3.1. Ш-тесты на плоском дне

1.3.2. Ш-тесты на неоднородном дне

1.3.3. Моделирование сдвиговых течений

1.3.4. Двумерная модель течений на сложном нерегулярном дне

1.4. Выводы первой главы

Глава 2. Суперкомпьютерное моделирование самосогласованной

динамики наносов, поверхностных и грунтовых вод

2.1. Параллельный СиЭЛ-алгоритм СБРЫ-ТУЭ метода

2.1.1. Общая структура расчетного модуля

2.1.2. Особенности распараллеливания алгоритма для СРШ

2.1.3. Эффективность параллельного СиЭЛ-алгоритма для ОРШ

2.1.4. Структура программного комплекса «<ЕсоС18-81ти1а1юп»

2.2. Моделирование динамики поверхностных вод при решении прикладных задач

2.2.1. Численное моделирование ливневого паводка

2.2.2. Численное моделирование формирования цунами

2.2.3. Численные модели гидрологических режимов ВАП

2.2.4. Построение зон затопления в период сезонных паводков

2.3. Совместная динамика наносов и поверхностных вод

2.3.1. Математическая модель динамики поверхностных вод и наносов

2.3.2. Численная модель и программная реализация

2.3.3. Численное моделирования гидродинамической аварии

2.3.4. Численное моделирования русловых деформаций

2.4. Самосогласованная динамика грунтовых и поверхностных вод

2.4.1. Математическая модель динамики грунтовых вод

2.4.2. Численная схема и программная реализация

2.4.3. Одномерные расчеты плоскопараллельных течений грунтовых вод

2.4.4. Двумерные расчеты совместной динамики поверхностных

и грунтовых вод

2.5. Основные выводы второй главы

Глава 3. Эффективные вычислительные методы суперкомпьютерного моделирования газодинамических течений

3.1. Разработка и реализация эффективных численных методов моделирования газодинамических течений

3.1.1. Основные этапы газодинамической версии численного метода СБРЫ-ТУЭ

3.1.2. Обобщение численного метода СБРЫ-ТУЭ для моделирования двумерных газодинамических течений

3.1.3. Реализация численного метода МиБСЬ

3.1.4. Условия устойчивости численных методов СБРЫ-ТУЭ / МиБСЬ

3.2. Обоснование и тестирование численных методов СБРЫ-ТУЭ /

МШСЬ

3.2.1. Анализ сходимости и вычислительной эффективности

3.2.2. Распад произвольных газодинамических разрывов

3.2.3. Трансзвуковые потоки газа с тангенциальными разрывами скорости

3.3. Суперкопьютерные алгоритмы численных методов СБРЫ-ТУЭ / МШСЬ

3.3.1. Разработка и реализация эффективных параллельных алгоритмов ОрепМР-СиЭЛ для суперкомпьютеров с СРШ

3.3.2. Анализ эффективности распараллеливания суперкомпьютерных ОрепМР-СиЭЛ-алгоритмов

3.4. Математическое моделирование динамики газовых компонентов

астрофизических объектов

3.4.1. Нелинейная динамика акустической неустойчивости в аккреционных дисках

3.4.2. Динамика газовых струй и дисков в молодых звездных объектах

3.4.3. Ударно-волновые структуры в струях из активных ядер галактик

3.5. Основные выводы третьей главы

Глава 4. Математическое моделирование неравновесных газодинамических течений

4.1. Математическая модель неравновесного колебательно-возбужденного и химически активного газа

4.1.1. Уравнения динамики неравновесного многокомпонентного газа

4.1.2. Колебательно-поступательная релаксация и химические реакции

4.1.3. Коэффициенты вязкости, теплопроводности и диффузии

4.2. Газодинамические неустойчивости в неравновесных средах

4.2.1. Акустическая и тепловая неустойчивости

4.2.2. Влияние химического энерговыделения на акустическую

и тепловую неустойчивости

4.2.3. Анализ решений дисперсионного уравнения, границы устойчивости и области запрещенных частот

4.2.4. Неустойчивости тангенциального разрыва скорости

4.2.5. Неустойчивость границы раздела между средами с различными значениями степени неравновесности

4.3. Численное моделирование одномерных газодинамических течений в неравновесных средах

4.3.1. Динамика акустической неустойчивости и образование ударно-волновых импульсов

4.3.2. Влияние времени колебательной релаксации на ударно-волновые импульсы и инициирование тепловых взрывов

4.3.3. Ударно-волновые структуры в дозвуковых и сверхзвуковых потоках неравновесного газа

4.3.4. Ударно-волновые импульсы в химически активном газе и образование детонационных волн

4.4. Численное моделирование двумерных газодинамических течений

в неравновесных средах

4.4.1. Двумерные ударно-волновые структуры в дозвуковых и сверхзвуковых потоках неравновесного газа

4.4.2. Влияние неравновесности среды на динамику неустойчивости Кельвина -Гельмгольца

4.4.3. Ударно-волновые и вихревые структуры в сверхзвуковых струях неравновесного газа

4.5. Основные выводы четвертой главы

Глава 5. Реализация на GPU лагранжевых методов для моделирования динамики бесстолкновительных и газовых компонентов галактик

5.1. Особенности параллельной реализации численных моделей N-тел HaGPU

5.1.1. Уравнения движения в модели N-тел и численный алгоритм316

5.1.2. Параллельный алгоритм динамики N-тел на GPU

5.1.3. Результаты моделирования динамики N-тел на GPU

5.1.4. Анализ вычислительной эффективности параллельных алгоритмов на multi-GPU платформах

5.2. Параллельный CUDA-алгоритм метода сглаженных частиц

5.2.1. Математическая модель столкновения двух галактических гало

5.2.2. Численный метод сглаженных частиц

5.2.3. Параллельный алгоритм динамики SPH-частиц на GPU

5.2.4. Основные результаты моделирования динамики SPH-частиц на GPU

5.3. Суперкомпьютерный программный комплекс «SPH + N-body» для моделирования самосогласованной динамики многокомпонентных галактик

5.3.1. Структура и оптимизация суперкомпьютерного кода «SPH

+ N-body»

5.3.2. Анализ вычислительной эффективности параллельного кода «SPH + N-body» на GPUs

5.4. Основные выводы пятой главы

Глава 6. Суперкомпьютерное моделирование самосогласованной динамики многокомпонентных галактических систем

6.1. Численная модель спиральной галактики Млечный Путь

6.1.1. Исходные параметры динамической численной модели Галактики по данным наблюдений

6.1.2. Математическая модель динамики спиральных галактик

6.1.3. Критерий гравитационной устойчивости

6.1.4. Результаты численного моделирования

6.1.5. Спиральная структура Млечного Пути: наблюдения Gaia DR 2 и численные эксперименты

6.1.6. Динамика газа в центральных областях Галактики и ограничения на параметры балджа

6.2. Спиральная структура карликовых галактик

6.2.1. Численные модели карликовых галактик

6.2.2. Эволюция звездно-газовых дисков карликовых галактик

6.2.3. Критерий устойчивости Тоомре для карликовых галактик

6.2.4. Спиральная структура карликовых галактик

6.3. Численное моделирование столкновения галактических гало

6.3.1. Математическая модель сталкивающихся галактических гало

6.3.2. Результаты численного моделирования

6.4. Динамика галактик с крупномасштабным противовращением

6.4.1. Описание модели ретроградной аккреции

6.4.2. Замещение газового диска галактики и активность галактических ядер при ретроградной аккреции

6.4.3. Образование полярных колец при ретроградной аккреции межгалактического газа

6.5. Имитационное моделирование эволюции взаимодействующих многокомпонентных галактик

6.5.1. Лобовое столкновение двух богатых газом спиральных галактик

6.5.2. Взаимодействие массивной спиральной галактики с карликовой дисковой галактикой

6.6. Основные выводы шестой главы

Заключение

Список публикаций по теме диссертации

Список литературы

Введение диссертации (часть автореферата) на тему «Эффективные численные методы и суперкомпьютерные алгоритмы для математического моделирования динамики многокомпонентных сплошных сред»

Введение

Основные направления и актуальность исследований

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

Выделим две основные проблемы, указывающие на актуальность диссертационного исследования:

1. При построении численного решения, описывающего динамику жидкости или газообразной среды в природных, технических и астрофизических системах, возникает проблема устойчивости, сбалансированности и точности численных методов, связанная с наличием в расчетной области крупномасштабных и мелкомасштабных движущихся границ типа «газ - вакуум» или «жидкость - сухое дно». Аналогичные проблемы с устойчивостью и точностью численных методов возникают при моделировании сильных разрывов, таких как ударные и детонационные волны, гидравлические скачки, тангенциальные и контактные разрывы. Учет внешних сил, характерный пространственный масштаб неоднородности которых сравним с размером вычислительной сетки и/или пространственным масштабом неустойчивых волновых структур, также может приводить к снижению точности, сбалансированности и устойчивости численного решения. Релаксационные процессы в неравновесных колебательно-возбужденных и химически активных сплошных средах, приводящие к развитию новых газодинамических неустойчивостей, образованию ударно-волновых и вихревых структур с малыми пространственно-временными масштабами, требуют модификации условий устойчивости численных методов. При моделировании динамики поверхностных вод на реалистичном рельефе местности в областях затопления/осушения, содержащих разрывы и излома профиля дна, известные численные методы интегрирования уравнений мелкой воды могут приводить к некорректным (нефизичным) решениям с отрицательными глубинами и аномально высокой скоростью течения. В неравновесных течениях колебательно-возбужденного и/или химически активного газа, а также в газовых подсистемах астрофизических объектов (аккреционных дисках и струях) важную роль при формировании наблюдаемых ударно-волновых и вихревых структур играют газодинамические неустойчивости. Детальное исследование линейной и нелинейной стадий эволюции этих неустойчивостей при неоднородном распределении

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

2. Для исследования нелинейной динамики нестационарных течений многокомпонентных сплошных сред требуется проведение вычислительных экспериментов с высоким пространственным разрешением и минимальным временем расчетов, что требует использования высокопроизводительных вычислительных систем. Например, в прикладных задачах моделирования аварий на гидротехнических сооружениях и гидрологических режимов крупных водных объектов (морей, озер, водохранилищ, рек и междуречий с обширными поймами) необходимо рассчитывать зоны затопления/осушения на территориях с площадью 102-105 км2 с разрешением 1-100 м в течение длительного времени от суток до нескольких месяцев. Для оперативного (часы - дни) получения информации о динамике жидкости на таких территориях требуемая производительность вычислений на операциях с числами двойной точности (РР64) составляет порядка 10 Тфлопс на один вычислительный узел (суперкомпьютер). Такая плотность вычислений достижима только при использовании параллельно выполняющихся программ на гибридных вычислительных системах с графическими процессорами (СРШ). Аналогичная плотность вычислений необходима при исследовании динамики газодинамических неустойчивостей в неравновесных течениях колебательно-возбужденного и/или химически активного газа, детальной структуры ударных и детонационных волн, образующихся на нелинейной стадии развития этих неустойчивостей, когда требуется высокое пространственное разрешение численных моделей ~ 1-100 мкм при размерах энергетических установок порядка нескольких метров. Такая же проблема возникает в астрофизических приложениях, например, при численном моделировании галактик, когда количество гравитационно взаимодействующих между собой частиц превышает 106. Следовательно, разработка эффективных параллельных алгоритмов для исследования динамики сплошных сред в природных, технических и астрофизических системах на гибридных суперкомпьютерах с СРШ также является важной и актуальной задачей.

Цели и задачи диссертационной работы

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

Достижение поставленной цели предусматривает решение следующих задач.

1. Создание нового численного метода интегрирования эволюционных уравнений модели мелкой воды, основанного на совместном использовании лагранже-ва (SPH — Smoothed Particle Hydrodynamics) и эйлерова (TVD — Total Variation Diminishing) подходов. Проведение анализа вычислительных свойств разработанного комбинированного лагранжево-эйлерова метода CSPH-TVD.

2. Построение математических моделей совместной динамики наносов, поверхностных и грунтовых вод. Разработка параллельных алгоритмов CSPH-TVD метода с использованием технологий CUDA (Compute Unified Device Architecture) и OpenMP (Open Multi-Processing) для имитационного моделирования самосогласованной динамики наносов, поверхностных и грунтовых вод на суперкомпьютерах c GPUs. Изучение особенностей динамики нестационарных процессов затопления/подтопления территорий с учетом неоднородностей их характеристик и самосогласованного описания динамики поверхностных и грунтовых вод, транспорта наносов и размыва дамб.

3. Обобщение численного метода CSPH-TVD для интегрирования полных уравнений газодинамики и разработка параллельных алгоритмов (CUDA и OpenMP) для моделирования нестационарных течений газа на суперкомпьютерах c GPUs. Исследование вычислительных свойств обобщенного метода CSPH-TVD и его параллельных реализаций в одномерном и двумерном приближении, сравнение с другими хорошо известными и апробированными методами типа MUSCL (Monotone Upwind Scheme for Conservation Laws).

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

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

6. Построение математических и численных моделей динамики бесстолкнови-тельных и газовых компонентов галактик. Разработка параллельных алгоритмов лагранжевых методов (SPH, N-body) и реализация их в виде комплекса программ для суперкомпьютеров с GPUs. Исследование точности и вычисли-

тельной эффективности разработанного программного обеспечения.

7. Разработка алгоритмов и методов суперкомпьютерного моделирования (SPH + N-body) для описания самосогласованной динамики многокомпонентных гра-витирующих систем галактик. Численное имитационное моделирование столкновения галактических систем, динамики спиральных галактик и эволюции их спиральной структуры.

8. Построение численных моделей динамики многокомпонентных (газ + звезды + темная материя) галактик, позволяющих объяснить наблюдаемую морфологию/кинематику этих гравитирующих систем и возникающих в них ударно-волновых структур.

Научная новизна.

1. Созданы новые численные методы семейства CSPH-TVD для численного интегрирования уравнений мелкой воды и полной системы уравнений газодинамики, основанные на комбинированном лагранжево-эйлеровом подходе (SPH + TVD) и обладающие высокой вычислительной эффективностью.

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

3. Разработаны и реализованы новые параллельные алгоритмы для методов CSPH-TVD, MUSCL (Monotonia Upstream-centered Scheme for Conservation Laws), SPH и N-body с прямым методом вычисления гравитационных сил (PPM — Particle-Particle Method), основанные на технологиях параллельных вычислений OpenMP-CUDA и GPUDirect/HostCopy для суперкомпьютеров с GPUs, обладающие высокой производительностью вычислений и эффективностью multi-GPU распараллеливания. Разработан новый иерархический сеточный метод сортировки на основе параллельного каскадного алгоритма нахождения частичных сумм в CUDA блоке, который позволяет повысить производительность вычислений при расчете газодинамических сил в SPH-методе.

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

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

газодинамических неустойчивостей (акустической, тепловой и тангенциального разрыва скорости) и определены границы устойчивости звуковых, энтропийных и вихревых мод в зависимости от степени неравновесности среды, моделей времени колебательно-поступательной (УТ — У1Ьга1юпа1-ТгапБ1а1юпа1) релаксации, нагрева и охлаждения.

6. Впервые на основе газодинамического моделирования с использованием численных методов СЭРЫ-ТУБ/М^СЬ исследована динамика звуковых волн, генерируемых источником возмущений, в неравновесной колебательно-возбужденной среде на линейной и нелинейной стадиях развития акустической неустойчивости. Проведена валидация разработанных численных моделей динамики неравновесного газа на основе сравнения результатов вычислительных экспериментов с аналитическими решениями линейной модели. Показано, что на нелинейной стадии развития акустической неустойчивости происходит формирование пилообразной системы слабых ударных волн (УВ), которые в неравновесной колебательно-возбужденной среде оказываются неустойчивыми и распадаются с образованием системы квазистационарных ударно-волновых импульсов (УВИ) большой интенсивности.

7. Впервые на основе разработанных численных моделей нелинейной динамики неравновесных сред показано, что в колебательно-возбужденном газе с высокой степенью неравновесности из-за взрывного развития тепловой неустойчивости (теплового взрыва) могут квазипериодически образовываться сильные ударные волны с очень малым скачком плотности на фронте, так называемые тепловые ударные волны (ТУВ), которые оказываются нестабильными и быстро эволюционируют к устойчивой квазистационарной системе ударно-волновых импульсов (УВИ). Ударно-волновые импульсы также образуются при дозвуковом натекании неравновесного колебательно-возбужденного газа на твердую стенку. В этом случае сначала образуются нестабильные слабые ударные волны, которые также быстро распадаются с образованием УВИ. Таким образом, предложены (обнаружены) три новых механизма формирования ударно-волновых импульсов (УВИ). На основе анализа результатов численных экспериментов доказано, что структура ударно-волновых импульсов определяется только параметрами неравновесной среды (в основном степенью неравновесности, моделями времени УТ-релаксации и нагрева/охлаждения) и не зависят от параметров и типа начальных возмущений, следовательно, УВИ являются ударно-автоволными импульсами (УАВИ). Показано, что в неравновесной среде даже с незначительным химическим энерговыделением УАВИ могут инициировать взрывное развитие тепловой неустойчивости (тепловой взрыв) за фронтом УВ, приводящее к образованию детонационных волн высокой интенсивности. Анализ результатов вычислительных экспериментов и структуры УВ, описываемых неравновесной ударной адиабатой, позволил определить новые критерии

существования квазистационарных ударно-волновых структур в колебательно-возбужденном газе. На основе разработанных аналитических и вычислительных методов анализа математических моделей неравновесных газодинамических потоков с тангенциальными разрывами скорости впервые показано, что колебательная неравновесность среды значительно усиливает неустойчивости Кельвина-Гельмгольца и волноводно-резонансных мод сверхзвуковой струи. А на нелинейной стадии эволюции этих неустойчивостей образуются более интенсивные вихревые и ударно-волновые структуры.

8. Соискателем проведено несколько больших серий вычислительных экспериментов самосогласованной эволюции многокомпонентных (газ + звезды + темная материя) галактических систем в различных условиях. Обработка и анализ результатов численного моделирования позволили сделать ряд новых выводов о динамике и морфологии различных галактических систем, согласующихся с данными астрономических наблюдений. При столкновении галактических сфероидальных систем продемонстрированы эффекты существенного обмена газом между двумя сталкивающимися гало галактик и формирования конусов горячего газа из-за сверхзвукового натекания на потенциальную яму с образованием ударных волн. Для нашей галактики Млечный Путь показано образование долгоживущей 3 млрд. лет) спиральной структуры с преобладанием неосе-симметричных трех- и четырех рукавных мод и центрального звездного бара размером около 3 кпк, приводящего к формированию максимума на кривой вращения в окрестности 1 кпк, что хорошо согласуется с данными наблюдений о кинематике и морфологии звездно-газового диска Галактики. Показано, что в карликовых галактиках (dS) с большим содержанием газа в галактическом диске могут формироваться бар и спиральная структура только за счет внутренних процессов, связанных с развитием гравитационной неустойчивости. Впервые получены ограничения на толщину дисков, при которых возможно формирование глобального спирального узора в dS-галактиках. Показано, что при аккреции на галактический звездно-газовый диск ретроградного газа из межгалактического пространства, обладающего противоположным угловым моментом, всего за 1 миллиард лет может сформироваться галактика с противовращающимися звездным и газовым дисками, а часть межгалактического газа образует вокруг этой галактики полярное газовое кольцо. Обнаружено, что при лобовом столкновении галактических дисков с большим содержанием газа могут образовываться газовые мосты с характерной мелкомасштабной волокнистой структурой, наблюдаемые в астрофизических объектах типа Taffy. Впервые продемонстрирована сильная зависимость пространственной структура и термодинамических свойств газа в этих объектах от геометрии столкновения и параметров исходных галактик. Предложены новые возможные механизмы формирования компактных и ультракомпактных галактик (cE/UCD) в окрестности массивной

спиральной галактикой типа Млечного Пути при многократном прохождении спутника (карликовой дисковой галактики) через диск этой массивной галактики.

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

Результаты диссертационного исследования, в частности, «Программный комплекс для моделирования самосогласованной динамики поверхностных вод и наносов на основе СБРЫ-ТУЭ метода» [В7, В9, В12, В17, В19, В20, В21] успешно внедрен в коммерческих компаниях России (ООО «РИЦ «ТЕЛЕНОВО», АО «Волговодпроект», ООО «Сергиевское карьерное управление») и Казахстана (ТОО «СЕОБАТ») при проведении научно-технической экспертизы гидротехнических проектов, связанных с моделированием зон затопления/подтопления территорий в период сезонных паводков, размыва дамб хвостохранилищ и русловых деформаций при добыче речного песка.

Полученные результаты использовались при выполнении НИР и НИОКР в рамках грантов, в которых соискатель являлся руководителем: РНФ 23-21-00401 «Математическое моделирование газодинамических неустой-чивостей в акустически активных средах», РФФИ 18-47-340007 «Суперкомпьютерное гидродинамическое моделирование русловых процессов, размыва и последствий прорыва заградительных дамб на водных объектах Волгоградской области», РФФИ 16-07-01037 «Суперкомпьютерное моделирование динамики жидкости и газа в природных и технических системах на основе лагранжево-эйлерова СБРЫ-ТУЭ метода высокого порядка точности», РФФИ 15-45-02655 «Прогнозирование гидрологического режима территории на основе нестационарных моделей динамики поверхностных вод», РФФИ 13-07-97056 «Геоинформационный портал для поддержки научных исследований в области экологии и рационального природопользования», РФФИ 10-07-97017 «Геоинформационная система для математического моделирования нелинейной динамики поверхностных вод суши», ФЦП «Старт» Проект № 11-3-Н1.4-0566 «Разработка программного комплекса для моделирования нелинейной динамики жидкости и газа».

Соискатель был официальным исполнителем в грантах РНФ 23-71-00016 и РФФИ 18-47-340003, 16-02-00649, 15-47-02642, 15-02-06204, 14-08-97044, 13-0197062, 11-07-97025, 09-02-97021.

Диссертационная работа выполнялась в рамках госзаданий Министерства науки и высшего образования РФ: «Создание программного обеспечения для моделирования физических сред и природных явлений» (проект N 2.852.2017/4.6,

2017-2019 гг.), «Программное обеспечение, суперкомпьютерные вычисления, математическое и компьютерное моделирование» (N 0633-2020-0003, 2020-2022 гг.), «Разработка методов моделирования сложных систем на основе интеграции прямых вычислительных экспериментов и интеллектуального анализа данных» (FZUU2026-0005, 2026-2028 гг.).

Основные результаты получены на оборудовании ЦКП «Суперкомпьютерный центр коллективного пользования ВолГУ», часть расчетов проводилась с использованием ресурсов ЦКП Сверхвысокопроизводительные вычислительные ресурсы МГУ им. М.В.Ломоносова (РНФ 23-71-00016, RFMEFI62117X0011).

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

Методология исследования включает теорию гидро- и газодинамики, общие принципы математического моделирования, методы решения уравнений гиперболического типа, эйлеровы и лагранжевы методы численного интегрирования дифференциальных уравнений в частных производных, методы разработки алгоритмов, структурное и объектно-ориентированное программирование, стандарты и технологии параллельного программирования (OpenMP, CUDA C/C++, GPUDirect и HostCopy), геоинформационные технологии для обработки и визуализации пространственных данных.

Достоверность результатов и выводов диссертационной работы определяется применением строгих математических моделей и хорошо апробированных численных методов/алгоритмов, сопоставлением результатов с полученными ранее, валидацией численных методов на основе данных натурного эксперимента или анализа математических моделей, сравнением численных и аналитических решений в некоторых предельных случаях, а также четким физическим смыслом полученных результатов и согласованностью их с современными представлениями о предмете исследования. Основные результаты, выносимые на защиту, опубликованы в рецензируемых высокорейтенговых журналах, входящих в Перечень ВАК, Белый список и/или индексируемых в международных базах Web of Science (WOS), Scopus, MathSciNet в том числе в изданиях, имеющих квартили Q1/Q2, уровни У1/У2 и категории К1/К2, докладывались на

научных семинарах и конференциях российского и международного уровней.

Апробация работы. Основные результаты диссертации докладывались на: Всероссийской научной конференции с международным участием «Параллельные вычислительные технологии (ПаВТ'2025)» (Москва, 2025); XI Международной конференции «Суперкомпьютерные дни в России 2025» (Москва, 2025); 8-м семинаре по численному моделированию в МГД и физике плазмы: методы, инструменты и результаты (Москва, 2025); Международной научной конференции «Кибер-физические системы: проектирование и моделирование» (CyberPhy-2024) (Казань, 2024); Всероссийской астрономической конференции ВАК-2024 (САО РАН, Нижний Архыз, 2024); 6th International Conference on Control Systems, Mathematical Modeling, Automation and Energy Efficiency (Липецк, 2024); International conference «Active galaxies at different scales and wavelengths» (САО РАН, Нижний Архыз, 2024); 11 Международном симпозиуме «Неравновесные процессы, плазма, горение и атмосферные явления», NEPCAP-2024 (Сочи, Адлер, 2024); Всероссийских научных конференциях «Современная звездная астрономия 2021, 2023» (Москва, 2021, Волгоград, 2023); 5th International Conference «Creativity in Intelligent Technologies and Data Science (CIT&DS) 2023» (Волгоград, 2023) Международных научных конференциях ФизикА.СПб 2021, 2023, 2025 (Санкт-Петербург); IV International conference «Supercomputer Technologies of Mathematical Modelling» (SCTeMM'19)» (Moscow, 2019); Конференции, посвященной 100-летию со дня рождения А.Г. Масевич (ИНАСАН, Москва, 2018); International Conference «Applied Mathematics, Computational Science and Mechanics: Current Problems» (Воронеж, 2017-2018 гг.); Научной конференции в рамках Летней Суперкомпьютерной Академии (Москва, 2017); Конференции «Современная звездная астрономия 2017» (Екатеринбург, 2017); XII Всероссийской школе-конференции молодых ученых и специалистов «Управление большими системами» (Волгоград, 2015); Научной конференции «Russian Supercomputing Days» (Moscow, 2016-2019 гг.); XII Межрегиональной научно-практической конференции «Проблемы устойчивого развития и эколого-эконо-мической безопасности регионов» (Волжский, 2016); Национальном Суперкомпьютерном Форуме (Переславль-Залесский, 2013-2015 гг.); ИнтерКарто / Ин-терГИС 20 «Устойчивое развитие территорий: картографо-геоинформационное обеспечение» (Белгород, 2014); Международных научных конференциях Интер-Карто/ИнтерГИС 15-18 «Устойчивое развитие территорий: теория ГИС и практический опыт» (Пермь, 2009, Ростов-на-Дону, 2010, Барнаул, 2011, Сен-Дье-де-Вож, Франция, 2012); Международной научной конференции «Математические методы в технике и технологиях» (Волгоград, 2012); Conference «Dynamics and evolution of disc galaxies» (Pushchino, Russia, 2010); Научно-практической конференции «Актуальные проблемы современной физики и математики» (Элиста, 2010); Международная конференция «Progress in the Study of Astrophysical

Disks: Collective and Stochastic Phenomena and Computational Tools» (Волгоград, 2003); IAU Colloquium 184: AGN Surveys (Бюракан, Армения, 2001).

Результаты обсуждались на научных семинарах (ВолГУ, ВолГТУ, Калм-ГУ, САО РАН, ИПТМУ РАН, МГУ, ЮФУ).

Публикации. Список основных публикаций по теме диссертации содержит 95 наименований, включая 74 научные работы и 21 РИД (результаты интеллектуальной деятельности). Основные результаты диссертации опубликованы в 74 научных работах (из них 9 без соавторов), входящих в Белый список, перечень ВАК и индексируемых в международных базах данных WoS/Scopus/Math-SciNet, в том числе 6 статей в изданиях с квартилем Q1, 7 статей — Q2, 2 статьи в журналах уровня У1, 5 статей — У2 и 22 статьи из перечня ВАК, включая 4 статьи в журналах категории К1. Имеется 21 свидетельство о регистрации программ для ЭВМ и баз данных, из них 10 без соавторов.

Список литературы диссертационного исследования доктор наук Храпов Сергей Сергеевич, 2026 год

дф - -

дд хх г дд хх г

д(пр) д(пр) - -

дд хх, дд хх,

дШ гк

дхг '

У^Пк рк

дШ г

гк

дхг

фк

(3.14)

где Ш,к = Ш (|хг — Хк|, к) = ЬШ (|хг — Хк|, к). Заметим, что функция сглаживающего ядра Ш,к обладает следующими свойствами:

Шк = Шк,,

дШ г

гк

дШкг дШа

дхг

дхк

дхг

= 0.

В численном методе СБРЫ-ТУЭ учитывается взаимодействие ¿-той частицы только с частицами, находящимися в соседних ячейках, в отличии от традиционного БРЫ-подхода [430], в котором ¿-ая частица взаимодействует со всеми частицами в радиусе 2к. В соответствии с идеологией явных сеточных методов для устойчивости алгоритма необходимо, чтобы в ¿-ую ячейку за время интегрирования поступала информация только из соседних ячеек ^ ± 1).

Поскольку рассматриваемый алгоритм является комбинированным лагран-жево-эйлеровым методом, то требование взаимодействия ¿-ой частицы только с соседними ^ ± 1) распространяется и на лагранжев этап. Поэтому в соотношениях (3.14) суммирование по индексу к ведется в пределах ¿ ± 1.

Заменяя входящие в уравнение (3.8) пространственные производные конечными суммами (3.14), содержащими аналитически вычисляемые производные от сглаживающего ядра Ш, получим

Фг

(

0

\

рг

Е

9г

дШ

дШк , , рк —--1--фк 0

рг дхг

0 гк

дхг

\

рг

рк

пг + Пк дШ гк (оп)г дШ г

2

дхг

+

рг

фк

гк

дхг

/

(3.15)

где величина Ш= Ш (|х0 — х01, к), вычисляемая в центрах ячеек, используется для аппроксимации внешних потенциальных сил (стационарных или нестационарных). Аппроксимация (3.15) обеспечивает консервативность метода СБРЫ-ТУЭ на лагранжевом этапе и выполнение условия хорошей сбалансированности в присутствии неоднородных внешних сил.

Отметим, что сглаживающее ядро Wik в (3.15) можно интерпретировать как удельный потенциал газодинамического взаимодействия частиц, а максимальный радиус их взаимодействия равный 2h не следует путать с размерами частиц £,(t).

Заметим также, что аппроксимация градиента давления (для одномерного случая dp/dx) в классическом SPH-методе [430] существенно отличается от представленной в формуле (3.15) и имеет вид

— « д(x-1) V мк( 0(xi;t) + 0(xk;t) + П Л dWik

dx, n ^ k\д(xi,t)2 0(xk,t)2 %k) dxi '

где nik — искусственная вязкость, введение которой необходимо для устранения многопотоковости (пролета частиц сквозь друг друга), 0 и 0 — определяются в соответствии с формулами (3.11) и (3.13). Данная аппроксимация градиента давления обеспечивает консервативность традиционного метода SPH. Однако в случае p = const и д = const представленное соотношение приводит к по-

дР п

явлению отличной от нуля силы давления, что противоречит условию —— = U

dx

вытекающему из постоянства давления. Таким образом, в классическом SPH-подходе при неоднородном распределении плотности и постоянном давлении появляются нефизичные силы, что приводит к неадекватному описанию контактных разрывов и нарушению условия хорошей сбалансированности.

С учетом (3.15) система уравнений (3.8)-(3.9) сводится к системе обыкновенных дифференциальных уравнений относительно t и для ее численного интегрирования необходимо использовать схемы Рунге-Кутты [283, 536], удовлетворяющие TVD-условию [316]. Ограничимся схемой Рунге-Кутты 2-го порядка точности. На шаге «предиктор» методом ломаных находим промежуточные значения U* и x* в момент времени tn+1 = tn + Tn

U* = un + Tn Фi (un, xn),

(3.16)

x* = xn + T un + xi = xi + 'n 2 '

где к = (i — 1, i, i + 1), знак «~» над вектором консервативных переменных Ui означает, что центр масс частицы xi смещен относительно исходного положения x0. На шаге «корректор» по наклону интегральной кривой в точке U* и x* вычисляем приращение значений Ui и xi на временном слое tn+1 = tn + Tn:

__ Un + U * T

u "+1 = + Tf Ф, (U k ,xk

xn + x* T un + un+1

xn+1 = xi + xi + Tn ui + ui

x Л --\

i 2 2 2

(3.17)

Таким образом, рекуррентные соотношения (3.16) и (3.17) позволяют для момента времени ¿?+1 вычислить интегральные характеристики частиц (НИП+1), движущихся под действием газодинамических и внешних сил. В процессе этого движения частицы деформируются (рис. 3.1), однако в рамках рассматриваемого СБРЫ-ТУЭ алгоритма нет необходимости явным образом вычислять £г(£).

Отметим, что для уменьшения ошибок округления при численной (компьютерной) реализации алгоритма в рекуррентных соотношениях (3.16) и (3.17) вместо жг(£) лучше использовать переменную £г(£) = жг(£) — ж0 (£? = 0), характеризующую смещение центра масс частиц относительно исходного положения ж0. В этом случае величину жг — жк, используемую при расчете Ш¿к и ее производной, можно представить в виде жг — = Н(г — к) + — .

Эйлеров этап численного метода CSPH-TVD

На данном этапе в момент времени £п+1 необходимо определить связь между величинами иП+1 и ИП+1, т.е. в соответствии с основным положением метода СБРЫ-ТУЭ надо найти изменение интегральных характеристик частиц, вызванное их перемещением в исходное состояние. Для этого воспользуемся хорошо известным в механике сплошных сред соотношением [65]

( Г ч , Г Г д„, , д Л

И(ж, ¿) (ж = у И(ж,£) + —Е(ж,^)| (ж, (3.18)

где и = Е)т, Е = Ии = (^и, £и2,Еи)т — вектор плотности потока

массы, импульса и энергии соответственно, непосредственно связанный с движением газа. Проинтегрируем уравнение (3.18) по времени в пределах от до £п+1 с учетом введенных ранее в § 3.1 обозначений. В результате левая часть уравнения примет вид

¿У (И (£) = Н (Ип+1 — И?) . (3.19)

При вычислении правой части (3.18) необходимо зафиксировать положение частицы в пространстве, например, в исходном состоянии (ж0, £0), тогда

¿„+1 Х0+1/2 Х0+1/2 ¿„+1

У — У И(ж,¿) (ж(£ + У У Е(ж,£) (¿(ж =

Х0-1/2 Х0-1/2

(3.20)

п ¿+1/2 — Е г—1/2,

Н (ИП+1 — И?) + Т? (Е¿+1/2 — Ег—1/2) ,

¿„+1

- 1 [

где Е(ж) = — Е(ж, ¿) — среднее значение потока на временном промежут-Тп 7

ке [tn,tn+1], F¿±1/2 = F(ж0±1/2) — усредненное по времени значение потока на границах эйлеровых ячеек (ж0±1/2 = x0 ± h/2). Приравнивая выражения (3.19) и (3.20), приводя подобные и сокращая на h, получим

un+1 = иn+1 - | (F¿+1/2 - Fi-1/2) . (3.21)

Для нахождения потоков F¿±1/2 можно применять различные численные методы, основанные на эйлеровом подходе и удовлетворяющие принципу не возрастания полной вариации численного решения (TVD-условие) [316]. В данной работе используются реализуемый в схемах годуновского типа подход и приближенные методы решения задачи Римана (Riemann Problem — RP) о распаде произвольного разрыва: метод Лакса-Фридрихса (LF) [235], метод Хар-тена-Лакса-ван Лира (HLL) [315] и модифицированный метод HLLC [573].

Задача Римана решается отдельно для каждой из границ эйлеровых ячеек. При этом в качестве начальных условий для RP необходимо задать слева (L) и справа (R) от рассматриваемой границы значения параметров потока, которые определяют величину скачка и могут быть получены на основе кусочно-полиноминальной реконструкции функции U(x,t). От порядка реконструкции зависит точность численного алгоритма. Здесь ограничимся рассмотрением случая кусочно-линейной реконструкции, обеспечивающей численной схеме второй порядок точности по пространству [589].

Используемая в методе CSPH-TVD модификация TVD-подхода заключается в том, что пространственная кусочно-линейная реконструкция величины U(x,t) внутри эйлеровых ячеек производится относительно положения «жидких» частиц, находящихся в момент времени tn+i на расстоянии Ах'П+1 = £'n+l от центов ячеек, где 0 < l < 1 (см. п. 1.2.2). Значение l = 0 соответствует стандартному TVD-подходу, а в газодинамической версии метода CSPH-TVD будем использовать l = 1. Пример кусочно-линейной реконструкции величины U(x,t) показан на рис. 1.6 (см. п. 1.2.2).

В газодинамических численных методах лучше использовать кусочно-линейную реконструкцию физических (примитивных) переменных V = (g,u,p)T [228], а уже на основе этой аппроксимации определять значения консервативных переменных U слева и справа от границы ячейки при решении задачи Римана. Если же сначала проводить линейную аппроксимацию величин U внутри ячеек и на их основе рассчитывать значения физических переменных V, то в областях сильного разрежения (p ^ 0, д ^ 0) с большими градиентами газодинамических величин или на границе с вакуумом может нарушаться позитивность (устойчивость) метода из-за появления отрицательных значений давления при его вычислении по формуле p = (7 — 1) [E — 0.5(gu)2/g]. Поэтому для обеспечения позитивности метода CSPH-TVD будем использовать кусочно-линейную

реконструкцию функции У(х,^) внутри ¿-й ячейки в момент времени £п+1

У(х, 1п+1) = Уп+1 + (х — хП+1) ХП+1 , х е [х0—1/2, х0+1/^ , (3.22)

где ХП+1 — вектор наклонов (угловых коэффициентов) линейного распределения У(х, €) внутри ¿-й ячейки, хП+1 = х0 + <^п+1 — положение центра масс частицы в момент времени tn+1 при ее смещении относительно исходного состояния на расстояние <^п+1 = £г(1п+1).

Наклоны кусочно-линейного распределения (3.22) должны удовлетворять ТУЭ-условию [316]. Для этого воспользуемся функцией-ограничителем [316,

561]

( Vп+1 _ V"п+1 Vп+1 _ Vп+1 п+1 _ г I г+1 г г г—1

ХП+1 = С

К h ' 1 h

где кг = 1 + ЙД1 — <^гп+1, ^ = £г/Ь. Величину (кг + кг—1)/2 можно интерпретировать как коэффициент деформации ¿-й частицы в момент времени tn+1. В качестве функций ограничителей, удовлетворяющих ТУЭ-условию и подавляющих нефизические осцилляции вблизи разрывов, могут применяться различные функции (1.45)—(1.47) (см. п. 1.2.2).

Выражения для У(х,^) слева (Ь) и справа (Я) от границы х0±1/2 с учетом (3.22) можно представить в виде:

Vl = V+1 + ( h _ д+Л xn+1, Vr = vn+1 _(2 + tn+1) xy+1, (3.23)

где i = i для границы x0+1/2 и i1 = i_1 для границы x0_1/2. С учетом (3.23), зная значения физических переменных Vl,r = (ql,r,ul,r,Pl,r), можно легко определить значения консервативных переменных Ul,r = (qlir, qlir ulr, Elr), где 2

QL,R ul,r , PL,R ™ .

Elr =--1--. Тогда используя методы приближенного решения за-

2 y _ 1

дачи Римана (LF, HLL, HLLC) можно вычислить значения потоков Fi±1/2, входящих в уравнение (3.21). В методе LF [235] для вектора потоков на границах ячеек имеем

F¿±1/2 = Fl+Fr + S. , (3.24)

где Fl,r = F (UL,R), S. = max(|SL|, |SR|) — скорость распространения единственного разрыва, разделяющего две области с постоянными значениями UL, Ur (см. рис. 1.7 в п. 1.2.2). Для минимальной SL и максимальной Sr скоростей распространения волн внутри ячейки справедливы оценки [243, 272]

Sl = min (ul _ Cl, ur _ Cr) , Sr = max (ul + Cl, ur + Cr) , где clrr = \JY Pl,r/Ql,r — адиабатическая скорость звука в газе.

В методе ЫЬЬ [315] выражение для потоков на границах ячеек имеет вид

«±1/2

Е

=

5яЕь — 5ьЕя + 5ь5я (ия — иь)

— 5

ь

Е

при 0 < 5ь,

при 5ь < 0 < (3.25)

при 0 >

В данном случае решение задачи Римана содержит два, распространяющихся со скоростями 5ь и 5я разрыва, которые разделяют три области с постоянными значениями иь, и*, и я (см. рис. 1.7 в п. 1.2.2). Значение величины и* между разрывами 5ь и определяется выражением

и

Еь — Ея + 5я ия — 5ьиь

5Я — 5

ь

Существуют и другие модификации метода ЫЬЬ, например, метод ЫЬЬС, в котором помимо ударных волн и волн разряжения учитываются контактные и тангенциальные разрывы [573]. В методе ЫЬЬС рассматривается три разрыва, которые распространяются со скоростями 5ь, 5*, 5я и соответственно разделяют четыре области газодинамического течения с постоянными значениями иь, и*ь, и*я, и я (см. рис. 1.7 в п. 1.2.2). Величина 5* в одномерном приближении определяет скорость контактного разрыва. В этом случае выражение для потоков на границах ячеек может быть представлена в виде

г±1/2

Еь,

при 0 < 5ь,

=

Е*ь = Еь + 5ь(и*ь — иь), при 5ь < 0 < 5*, Е*я = Ея + 5я(и*я — ия), при 5*< 0 < 5Я,

(3.26)

Ея, при 0 > 5я.

Для величин и*к с индексом К = (Ь, Я) имеем

и

и к — Ек + Р*к О

, Б* = [0, 1, 5*]т,

— 5*

Р*к = (1 — 1)рк + (5к — ик)(5* — ик),

(3.27)

где 1 = 0 для стандартной реализации ЫЬЬС в соответствии с [573], 1 = 1 для модифицированной версии ЫЬЬС в методе СБРЫ-ТУЭ. Выражение для скорости контактного разрыва 5* получается из уравнения р*ь = р*я

(1 — 1)(ря — Рь) + (5ь — иь) — £я«я (5я — ия)

—

£ь(£ь — иь) — ^я(5я — ия)

(3.28)

С учетом (3.27) ЫЬЬС-потоки Е*к в (3.26) можно представить в виде

5*(5к — Ек) + О* =-^-~-. (3.29)

Ок — О*

3.1.2. Обобщение численного метода CSPH-TVD для

моделирования двумерных газодинамических течений

Обобщим алгоритм СБРЫ-ТУЭ для двумерных (2Э) течений газа. Будем исходить из интегральных законов сохранения массы, импульса и энергии для идеального газа в декартовой системе координат г = (ж, у):

-^УУ^хЫу = 0, (3.30)

JJ £и-ж-у = — JJ — /Х^ -ж-у, (3.31)

*Л —у=—// ^

— 11 £г> -ж-у = — II ( —--£/у ) -ж-у, (3.32)

Е-ж-у = — 11 [ —---1-----£и/х — £у/у ) -ж-у, (3.33)

где Е(£) — «2О-объем» «жидкой частицы», деформирующийся произвольным

^М2

образом в процессе ее движения, V = (и,г>) — скорость газа, Е =--Ъ £,

дф дф ^ ^

/х = —и = ——--компоненты удельной потенциальной внешней силы.

дж ду

Система уравнений (3.30)-(3.33) замыкается уравнением состояния идеального газа (3.4).

Введем характерное значение плотности £о, а также пространственный гд и временной ¿о масштабы задачи. Далее будем использовать только безразмерные величины (£, г, £, V, р, Е, ф) переопределив их следующим образом:

£ г £ ^о р£2 „ Е£д . ф£о

£ = —, г = —, £ = —, V = -, р = -2 , Е = -2 , Ф = ^Г •

¿0 Го £о го £ог2 £ог2 Г

о

Рассмотрим двумерную реализацию численного алгоритма СБРЫ-ТУЭ. Как было отмечено в первой главе (п. 1.2.1), на лагранжевом этапе можно использовать два пространственных шаблона аппроксимации — «звезда» и «крест»

(см. рис. 1.4). В первом варианте для расчета сил в (i, j)-ячейке используются все соседние ячейки, включая угловые, а во втором — угловые ячейки не учитываются.

Воспользуемся стандартными процедурами дискретизации сплошной среды, применяемыми в численных схемах, основанных на эйлеровом и лагранже-вом подходах. Покроем расчетную область равномерной эйлеровой (неподвижной) сеткой с пространственным шагом h, где hi,j = h = const (i и j пространственные индексы) и x0+1 = x0 + h, y0+1 = y0 + h. Совместим, в начальный момент времени, частицы с ячейками эйлеровой сетки и введем временные слои tn+1 = tn + тп с неравномерным шагом тп.

В двумерном случае численный алгоритм CSPH-TVD также содержит лаг-ранжев и эйлеров этапы, подробное описание которых было представлено в п. 3.1.1. Здесь рассмотрим отличия, возникающие при построении численного алгоритма для двумерных газодинамических течений. Частицы будем характеризовать массой Mi,j (t) = gdxdy, импульсом Pi,j (t) = / / gv dxdy и

(t)

Si

энергией (£) = JJ Edxdy.

^ (*)

Лагранжев этап ОБРИ-ТУВ метода

Как отмечалось выше, необходимость перехода от лагранжева этапа к эй-

леровому и обратно требует введения в численный алгоритм величин:

/ \

1

AiJ (t) = h2

A(x, t)dxdy

(3.34)

\Si,j (t)

/

являющихся аналогом средних значений функции А = (д, дп, ду, Е,р, ф) в ячейках при конечно-объемной аппроксимации на неподвижной сетке.

Подставляя в систему уравнений (3.30)-(3.33) соотношение (3.34) и используя вспомогательную переменную р = л/2р получим:

dUij dt

Q

i,3

(

0

\

Pi

Pi

dp

dxi dp

3 dyi

gi

Зф

dxi

дф

Pi,

u.

i,3

dx.

+

dp | д(up)

dxi дф

+ v.

dy3

dp dyi

\

-(gu)ij^--(gv)

dxi

'ij

+

дф dyj

d (vp)

dyi

/

(3.35)

2

T

где U = (g, gu, gv,E) — вектор консервативных переменных, u«j =

v«j =

(£v)«j g«J

(£u)«j g«J

и

скорости центра масс частиц. Значение величины ^«j можно

выразить через компоненты вектора консервативных переменных и в виде:

^«j = л/2(7 - 1Ь/E«j

(gu)2? (gv)2

2g«j 2g«,j

Для аппроксимации пространственных производных в (3.35) воспользуемся аналогичной процедурой на основе SPH-подхода со сглаживающим ядром W(q) (см. п. 3.1.1). Выражения для нормировочных постоянных сглаживающего ядра и переменной q в (1.26)—(1.28) (см. § 1.2) зависят от пространственного шаблона численного алгоритма — «звезда» или «крест» (см. рис. 1.4). В шаблоне «звезда» имеем q = \J(Дж)2 + (Ay)2/h, а в шаблоне «крест» — q = |Ax|/h или q = |Ay|/h в зависимости от направления, вдоль которого вычисляются силы. Здесь Дж = (ж« — ж«/) и Ay = (yj — yj/) — разность соответствующих координат между (i,j)-частицей и соседними (¿', j')-частицами.

Введем следующие обозначения: W = W(q) h2, W = W(qo) h2,

dW = dW(q) dq = dW(q) (ж« — ж«/) дж^ dq дж« dq q h2 '

dW = dW(q) dq = dW(q) (y- — yj) dyj dq dyj dq q h2 '

тогда правую часть уравнения (3.35) Q«j в методе CSPH-TVD с шаблоном «звезда» можно записать в виде

Q

«j

0

ÖXj

dW g.

+--Wj'" о

•+1 jti ( dW + dW0

«'=«—1 j'=j —1 ^ «+1 j+1 Г

^«j E E

+

«j

dy? ^

Ф«

dW0

«j

dyj

«'j ' л ' «'j '

«—1 J/=J—1 ^

«+1 j+1 Г . / ч o

«'=«—1 j '=j —1 \ «+1 j+1 ,

E E W

u«j + u«/,j/ dW (gu)«,

2

d^«

+

dW ^-

d^«

+

«'j '■

v«j + v«/,j/ dW i (gv)«jA dWl

«'=«—1 j'=j—1

dy«

+

-ф.

«'j

dy« J

(3.36)

2

При использовании в методе СБРН-ТУЭ шаблона «крест» БРН-аппрокси-мация (3.35) примет вид

д

¿3

¿+1

¿'=¿-1 3+1

Рь3 ^ { Р3

О

дШ д.

+

¿,3

3 г ,3 дх.

дхг рг,з дШ дз

Фъ

дШ0

+

ф.

дШ

о

дУз ' Рг/113 ду

г+1

Е

3'=3 -1

3

¿'=¿-1 I

3+1 г

'

Щ,з + Щ'з дШ (ди)^

2

дх:.

+

+ Уз дШ (ду)Ьз

з'=з -1

2

дy:

+

Зфг',з

ф3

дШ

о

3 дх., дШ

+

Рьз

дуг

(3.37)

Для численного интегрирования по времени системы ОДУ (3.35) с учетом (3.36)-(3.37) воспользуемся схемой Рунге-Кутты 2-порядка точности (см. п. 3.1.1). В двумерном приближении схему «предиктор-корректор» (3.16)-(3.17) можно представим в виде:

и_ из + тп (ип,рп,Сп) ,

р* _ т ипз + и*3 _ 1 п

2

уп. + V*.

/-* _ т ¿3 ¿3 Ц3 _ 1 п'

(3.38)

2

ип. + и * ■

ттп+1 _ "^¿3 + и 3

¿3 2

+ 2 д

Р * гт ип + иП+1 рп+1 _ + Тп ¿л ¿,3

ри 2+2 2

М,

т Уп + Уп+1

л п+1 _ + тп ¿^ ¿,3 ^ 2+2 2

(3.39)

2 2 2 ^3

где величины р^ (г) _ х^ (г) - х0 и (¿,3 (¿) _ У^ (¿) - У0 (Рп _ О, (п _ О) определяют смещения центра масс частиц в эйлеровых ячейках относительно исходного положения (х°,у0).

Эйлеров этап численного метода CSPH-TVD

На эйлеровом этапе вычисляются потоки массы, импульса и энергии, обусловленные перемещением газа через границы ячеек, в момент времени гп+у2 _ гп + тп/2. На этом этапе находится приближенное решение задачи Римана. Разность втекающих и вытекающих в ячейки потоков позволяет определить изменения характеристических «жидких» частиц, рассчитанных на предыдущем лагранжевом этапе в момент времени гп+1 _ гп + тп:

ТТп+1 _ тт п+1 ' п (ъ

Т ¿,3 _ Т ¿,3 - X ^ ¿+1/2,3

гп+1

п

Т,

¿-1/2,3

3 - ¿,3+1/2 - С¿,3-1/^ , (3.40)

где Е = Ии = (ди, ди2, диг>,Еи)Т, С = Иг> = (дг>, диг>, дг>2, Еи)Т — векторы адвективных потоков физических величин (массы, импульса и энергии),

¿п+1 ¿п+1

- 1 [ - 1 [

Е(х,у) = — Е(х,у,£) и С-(х,у) = — С(х,у,£) — средние значе-

тп «/ тп «у

ния потоков на временном промежутке [£п,£п+1], а Е¿±1/2,^' = Е(х0±1/2, у0±1/2) и СС¿,¿±1/2 = СС(х0±1/2, у0±1/2) — значения этих потоков на границах эйлеровых ячеек (х0±1/2 = х0 ± Ь/2 и у0±1/2_= У0 ± _/2).

Для вычисления Е¿±1/2,7 и С-¿,¿±1/2 можно использовать приближенные методы решения задачи Римана из п. 3.1.1 (3.24)-(3.26), которые легко адаптируются на случай двумерных течений. Кусочно-линейную реконструкцию сеточных функций будем проводить, также используя вектор физических переменных V = (д,и,г>,р)Т. Обобщая выражение (3.22) для двумерного случая в пределах (г, 3)-й ячейки (ж0_1/2 < х < ж0+1/2, у0_1/2 < у < у0+1/2) получим

У(х, у, £п+1) = VП+1 + (х _ х0 _ ХП+1 + (у _ у0 _ С^) Y:;+1 , (3.41)

где 1 и Y;г+1 — векторы наклонов линейного распределения У(х,у,£) внут-

ри (г,3)-й ячейки вдоль х- и у-направлений, соответственно. Наклоны кусочно-линейного распределения (3.41) должны удовлетворять ТУЭ-условию [316, 561]:

=

=

-Vп+1 _ -Vп+1

у ¿+1,.7 у ^

(х) г,

п+1 _ Vп+1

У М+1 У ^

«у _

■Vп+1 _ Vп+1

м_г-1,7

(х) г,

-у п+1 _ ^^^п+1