Математическое моделирование нестационарных течений с ударными волнами в пространственно–неоднородных средах тема диссертации и автореферата по ВАК РФ 00.00.00, кандидат наук Горкунов Сергей Владимирович

• Горкунов Сергей Владимирович
• кандидат наук
• 2023

Оглавление диссертации кандидат наук Горкунов Сергей Владимирович

Введение

Глава 1. Описание математических моделей

1.1 Общие положения

1.2 Физико-математическая модель движения газа с учетом

влияния фрагментов на турбулентные характеристики потока

1.2.1 Модель движения газовой среды

1.2.2 Влияние фрагментов оболочки на течение газа

1.2.3 Модель диффузионного горения

1.3 Физико-математическая модель движения твердой оболочки как термоупругой среды

1.3.1 Модель движения оболочки

1.3.2 Пластические и упругие деформации металлов

1.4 Уравнения состояния

1.4.1 Уравнение состояния идеального газа

1.4.2 Уравнение состояния JWL

1.4.3 Уравнение состояния твердой оболочки

1.5 Континуальная модель фильтрационного течения с поверхностями фазового перехода в пористой среде

Глава 2. Пакет прикладных программ

2.1 Решатель SPH_solveг

2.1.1 Исходные данные для SPH_solveг

2.1.2 Описание численного метода SPH

2.1.3 Граничное условие типа свободная поверхность

2.1.4 Алгоритм параллельного вычисления SPH_solveг

2.2 Описание решателя GK_solveг

2.2.1 Структура исходных данных для GK_solveг

2.2.2 Метод Годунова-Колгана

2.2.3 Эффективное использование расчетной области при расчете параметров ударной волны

2.2.4 Алгоритм параллельного вычисления GK_solveг

2.3 Сопряжение решателей SPH_solveг и GK_solveг

2.4 Решатель для расчета течения с фазовыми переходами в

пористой среде

Глава 3. Валидация и верификация расчетных модулей

3.1 Верификация 8РЫ_зо1уег

3.1.1 Распад разрыва в упругой среде

3.1.2 Истечение газа в вакуум

3.2 Верификация решателя СК_эо1уег

3.2.1 Галилеевская инвариантность

3.2.2 Регулярное отражение ударной волны от жесткой стенки

3.2.3 Разлет газа в вакуум от косой стенки

3.3 Валидация расчетных модулей

3.3.1 Распространение ударной волны в канале с препятствиями

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

3.4 Верификация решателя для расчета течения с фазовыми переходами в пористой среде

Глава 4. Примеры применения программного комплекса для

решения практических задач

4.1 Моделирование эволюции ударных волн в условии промышленной застройки

4.1.1 Постановка задачи

4.1.2 Результаты расчетов

4.2 Распространение ударных волн в протяженных каналах с препятствиями

4.2.1 Постановка задачи

4.2.2 Эффективная ударная адиабата канала с препятствиями

4.2.3 Задача Римана для распада начального разрыва

давления и плотности в канале с препятствиями

4.2.4 Течение газа в сильно загромождённых каналах

4.3 Исследование коротковолновой неустойчивости контактной границы жидкость-газ в пористой среде

Заключение

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

Введение диссертации (часть автореферата) на тему «Математическое моделирование нестационарных течений с ударными волнами в пространственно–неоднородных средах»

Введение

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

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

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

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

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

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

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

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

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

3. Реализовать разработанный метод в виде алгоритма численного расчета для высокопроизводительных вычислительных систем, создать программный комплекс.

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

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

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

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

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

3. Результаты исследования нелинейной стадии развития неустойчивости границы раздела газ/жидкость в пористой среде.

4. Результаты математического моделирования эксперимента о прохождении ударной волны через канал с периодически расположенными препятствиями. Автомодельные решения задачи о распаде произвольного разрыва на границе канала с препятствиями и гладкого

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

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

Программный комплекс для высокопроизводительных вычислительных систем использовался при выполнении научных проектов, поддержанных грантом из федерального бюджета в форме субсидии (Соглашение № 075-15-2020-785), Российским научным фондом (проекты № 16-19-00188, № 16-11-10195, № 21-11-00126), а также в рамках Государственного задания (Рег. N НИОКТР АААА-А21-121012190116-5).

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

Основные положения, выносимые на защиту:

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

2. Проведен анализ устойчивости границы раздела газ/жидкость в пористой среде при вытеснении жидкости газом. Установлено, что скорость роста амплитуды возмущения на нелинейной стадии увеличивается с уменьшением длины волны при использовании закона Дарси. Призна-

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

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

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

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

Апробация работы. Основные результаты работы докладывались на следующих конференциях и семинарах: 8-я, 11-я, 12-я Международная научная школа молодых ученых «Вихри и волны в сложных средах» (г. Москва, 2017, 2020, 2021); 11-я, 14-я «Всероссийская школа-семинар «Астрофизика и физическая механика классических и квантовых систем (АФМ)» (Москва, 2017, 2020); Всероссийская конференция «Физика взрыва: теория, эксперимент, приложения» (Новосибирск, 2018); 10-th International conference for promoting the application of mathematics in technical and natural sciences «AMITANS 2018» (Albena, 2018); XVII международная научная конференция и школа молодых учёных «Физико-химические процессы в атомных системах» (Москва, 2019); 9-я всероссийская научная конференция с международным участием «Механика композиционных материалов и конструкций, сложных и гетерогенных сред» (Москва, 2019); Ежегодная научная конференция отдела горения и взрыва (Москва, 2021); V, VI, VII Международная конференция «Лазерные, плазменные исследования и технологии - ЛАПЛАЗ» (Москва, 2019, 2020, 2021).

Личный вклад. Все представленные в диссертационной работе результаты численных расчетов были получены автором лично либо при его непосредственным участии. Автор лично участвовал в разработке численных методов, алгоритмов и пакета прикладных программ, созданного при выполне-

нии данной работы. Постановка задач и обсуждение результатов проводились совместно с научным руководителем.

Публикации. Основные результаты по теме диссертации изложены в 19 печатных изданиях [1—19]. Из них семь [1—7] опубликованы в журналах, рекомендуемых ВАК, в том числе шесть [1—5; 7] в индексируемых международными базами данных Scopus и/или Web of Science. Две [5; 7] из этих 6 публикаций относятся к иностранным журналам из первого квартиля Q1. Зарегистрирована 1 программа для ЭВМ [9]. Результаты работы содержатся также в 11 тезисах докладов [8; 10—19].

Объем и структура работы. Диссертация состоит из введения, 4 глав, заключения. Полный объём диссертации составляет 141 страницу, включая 63 рисунка и 7 таблиц. Список литературы содержит 156 наименований.

Глава 1. Описание математических моделей

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

1.1 Общие положения

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

Обычно рассматриваются четыре закона сохранения:

— закон сохранения массы;

— закон сохранения импульса;

— закон сохранения энергии;

— закон сохранения отдельных компонентов, составляющих систему.

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

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

В задачах, связанных с метанием твердых оболочек продуктами детонации, разрушение оболочки приводит к изменению параметров потока газовой фазы [20—22]. Нарушение сплошности оболочки существенно меняет поле течения не только количественно, но и качественно [23—25]. По мере развития процесса происходит активное смешение продуктов детонации с воздухом, что может сопровождаться значительным энерговыделением за счет экзотермических реакций.

1.2 Физико-математическая модель движения газа с учетом влияния фрагментов на турбулентные характеристики потока

1.2.1 Модель движения газовой среды

Согласно [26] отношение масштабов турбулентности (максимального к минимальному) составляет Поэтому при численном решении уравне-

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

Для описания течения газа постоянного состава рассмотрим систему уравнений нестационарного движения для гомогенной сжимаемой среды в трехмерной декартовой системе координат. Уравнение неразрывности, три уравнения для компонент импульса и уравнение для полной энергии используются в виде:

д_

~дъ

Рд

Рд

Рд Ед

+

д

дх^

рд

Рдид,гид,з + ^Рд - тд,ьз = ° (1Л)

(рд Ед + Рд )ид,] — ид,гтд,г,з + Яд,з где Рд - плотность газа, Рд - давление газа, - символ Кронекера,и-]-ая компонента скорости в декартовых координатах, Ед- полная энергия газа, равная сумме внутренней (ед) и кинетической энергий. Тензор напряжений определяется как ), д^д^)-]-компонента теплового потока. Здесь и далее индексом д обозначается газовая фракция.

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

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

Турбулентные потоки характеризуются рядом масштабов длины и времени, при этом малые масштабы не могут быть смоделированы с помощью доступных в настоящее время вычислительных ресурсов. Таким образом, для реализации стратегии моделирования, учитывающей все масштабы физической значимости, широко используется модель больших вихрей (Large Eddy Simulation - LES). В ней полностью разрешены большие масштабы течения, в то время как малые масштабы (которые не могут быть смоделированы) являются неразрешенными (их также называют подсеточными масштабами или subgrid-scale - SGS). Разделение масштабов в турбулентных течениях впервые было проанализировано Колмогоровым [28; 29]. Для получения системы уравнений модели больших вихрей необходимо провести осреднение уравнений (1.1). В рамках существующих подходов к статистическому описанию турбулентности предполагается, что осреднение по временным и по пространственным масштабам эквивалентны. В силу более простой записи далее представлено осреднение по времени, тогда как в численных расчетах удобно использовать общий масштаб осреднения А = (AxAyAz)3, где Ах, Ay, Az определяют локальный размер расчетной сетки по х,у и z- направлениям.

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

/ = 7 + /'• (1.2)

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

как:

7 =

АТ

'АТ

I' = I - }■

(1.3)

(1.4)

Причем временной масштаб осреднения АТ достаточно большой, чтобы = 0.

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

7 + 9 = 7 + 9, 1'9 = 7д,

к/ = к/, (к= еопэ^,

7 = 7,

т ¡л а л 97 д7

к = 0, (к = еопэ^, — = —,

оз оз

7' = 0,

19 = 79-

Кроме того, применяется взвешенное осреднение (осреднение по Фавру

[30—32]):

7 = 7 + 7'',

где / - средняя величина, - пульсации определяются как:

(1.5)

' = ч-

7" = 7 - 7-

Осреднение по Фавру имеет следующие свойства:

(1.6) (1.7)

91" = 0,

р/ = р/ = р/ + Р Г ,

7 = 7 + РХ

* Р '

7" = - — = 0. р

Для получения системы уравнений модели больших вихрей нужно осред-нить рд, Рд, согласно (1.2), а щ,дсогласно (1.5) и, следуя подходу [27; 33], подставить в (1.1). В результате получается следующая система:

д_

Ж

р9

рд ид.г Р~9 Е9

+

дхп

рдид.г

+ Ь • Р — т • ■ + т

I иг. п1 а ал. п >

р дид.г ид.З 1 Ь г.З1 д тд .ЬЗ • ~д.г.] р дид ¿Ед + ид.]Рд + Чд.] — ид.г тд .¡.] + ® П] + ^а.

д.]

= 0,

(1.8)

где переменная р определяет осредненную плотность газа, - скорость газа по 1-му направлению, Рд - осредненное давление газа, Ед - осредненная полная энергия на единицу массы, определяющаяся как Ед = ед + + када

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

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

Осредненный вязкий тензор напряжений и осредненный вектор теп-

лового потока определяются как:

^=* (^ + Ш)- 2 ^£, (1.9)

—ВТ

<Ш = -кя ^. (ио)

Здесь, Тд- средняя температура газа, рд- осредненная вязкость газа при заданной температуре газа, кд- тепловая проводимость при заданной температуре газа.

Тензор вязких напряжений с учетом турбулентности:

т71з = Р - ^ ЪА (1.11)

Поток энтальпии с учетом турбулентности:

Я7з = Ря(Еаид,г - Едщ,) + (ид,3Рд - Щ^эРд). (1.12)

Диссипация энергии в турбулентном потоке:

®д-з = (из,з тд,ьз — ид,з тд,ьз )• (ОЗ)

Тензор вязких напряжений определяется через скорость "деформации":

1 —- 2

т^ = -2р^(Бд^ - -+ 6,,3када. (1.14)

Скорость "деформации" Sgjj определяется как:

^=2( Ц+Э • (1.15)

Подсеточная турбулентная вязкость ^ определяется из соотношения:

гн = С у/к~А, (1.16)

где А - локальный масштаб осреднения, С - константа модели подсеточной турбулентности.

Уравнение, учитывающее перенос, диссипацию и генерирование кинетической энергии турбулентности, записывается в виде:

| (т,*) + ^ (рЛ.^.) = ^ + рг - от. (1.17)

Здесь, Р^98 и 08кд8 определяют генерирование и диссипацию подсеточной кинетической энергии соответственно

г) / Bksgs\ ksgs 3/2

^=ъ ж), =ср—, (1Л8)

где С£ - константа модели подсеточной турбулентности.

Значения Cv и С£ можно задать равными 0.067 и 0.916 в соответствии с [34]. Эти значение получены в рамках модели Localized Dynamic k-Equation Model (LDKM) [35—37].

Плотность потока подсеточной полной энтальпии:

ъ = Pgpt, S" (1Л9)

Здесь, Рrt -турбулентное число Прандтля, Hg - осредненная полная энтальпия.

1.2.2 Влияние фрагментов оболочки на течение газа

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

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

Уравнения движения воздуха и фрагментов на стадии, когда расстояние между фрагментами много больше их характерных размеров, записываются в следующем виде [22; 38—41]:

Ф 9 Р 9и9,г

фд Рдид,гид,3 + фдЬг¿Рд - фдТд^ фд р дид^Ед + Фдид, зРд + фдЯд,з - фд ид,г Тд ,ьз фд Рдуд,к (ид,з + Уд,з,к) 0

(6*Р*- Т*)Ц

(6,3Р*и* - Т*и*)Ц Шк

где Фд - объемная доля газовой фазы, рд - плотность газовой фазы, Рд - давление газовой фазы, 6^ - символ Кронекера, ид^ - г-ая компонента скорости газовой фазы, Ед - полная энергия газовой фазы, т^- - тензор турбулентных напряжений, - тепловой поток в ^'-ом направлении, Удк - массовая доля вещества к-ого вида, Уд^,к - скорость диффузии к-ого вещества в j-ом направлении,Шк - прирост вещества к-ого вида за счет химических реакции, Ff,i - импульс передаваемый в газовую фракцию от фрагментов, рf - межфазный массоперенос, Qf - поток тепла от фрагментов , Wf - работа, совершаемая фрагментами над газовой фракцией, Sf,к - межфазный перенос массовой доли компонентов для к-го вещества, верхний индекс * обозначает значения на границе раздела двух фаз.

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

ед = ХХ^к, (1.21)

к=1

где Мс обозначает количество компонентов в газовой фазе, ед^к - внутренняя энергия вещества к-ого вида, вычисляемое согласно (1.72).

Тензор турбулентных напряжений:

* (щ-+%)+^дх., (1.22)

д_ дъ

Ф

Ф 9 Р9ид,г Ф д РдЕ д дх^

Ф д Р дУд,к

рf

Ff,i

+

Qf + Wf

Sf,к

(1.20)

где Рд - коэффициент вязкости газа, Л - объемная вязкость, в соответствии с [42] можно считать Л = -2/-рд.

Для газов вязкость является возрастающей функцией температуры. Существуют разные модели для описания этой зависимости [42]. Наиболее распространенной является:

Рд = Р° (Ю , (1.23)

где ро вязкость газа при температуре Т0, п обычно принимает значение близкое к 0.7 [42].

Тепловой поток :

дТ ^

Яд,3 = к9+ Р^ 19,кЬд,кУд^к, (1.24)

3 к=1

где кд,к - энтальпия вещества к-го вида газовой фазы, кд- теплопроводность газовой фазы.

Скорость диффузии вещества к-го вида с учетом закона Фика определяется как:

Уа1к = - ^, (1.25)

9V , Яг ■ V

где Одк - коэффициент диффузии вещества к-го вида.

В процессах, рассматриваемых в работе, фрагменты состоят из твердой фракции и у них отсутствуют фазовые переходы, поэтому:

Р/ = (1.26)

к = (1.27)

В силу того, что объемная доля фрагментов мала и рассматриваемые процессы достаточно быстрые, можно пренебречь теплообменом между фрагментами и газовой фазой. Поэтому:

^ = 0. (1.28)

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

1 я

Рк = Д^ ^ |_2 ^^9д 1и1>г>п - - Щ^пЦ . (1.29)

п=1

Обмен энергией вследствие работы сил сопротивления:

1 ГП

^ = \_2 Т)С1 - и9,%п | (и К,п - ид,г,п)и/,г,п\ . (1.30)

п=1

Изменение скорости фрагмента:

(!и у I 1 п о . .. , . ,

= ДУ2Г°С199 1щ>г - ид>г1(и^г - (1.31)

в уравнениях (1.29) - (1.31) ту - масса фрагмента, гу - радиус фрагмента, АУ - объем, по которому происходит осреднение (обычно это объем ячейки разностной сетки), Су - коэффициент сопротивления фрагмента, и^ - г компонента скорости фрагмента, N - количество фрагментов в объеме А У.

Турбулентность имеет ключевое значение в модели горения. Именно турбулентная диффузия обеспечивает доставку кислорода в область, занятую реакционноспособными продуктами. В рамках модели больших вихрей уравнение (1.20) с учетом (1.26) - (1.28) записываются в виде:

фд Рд фд РЩ,г

д фд РЩ,г д дхj фд Р дид,гид,з + фдЬг¿Рд - фдТд+

дъ фд Р^ _ ф дРуд,к _ фд Рдид,гЕд + Фдид^Рд + Фд ^ - Ф дУ, д ^ + _ + у) + Щг;дк+

+Ф Тдз-

+ф^т

0

дФ

дхл

дф дхл

Ы к

0

+

0

(1.32)

Осредненный вектор теплового потока определяется как:

^ ___^__^^

= -кд Ид,кУд^к + ^ ^, (1.33)

3 к=1 к=1

где к подсеточный поток тепла.

Скорость диффузии вещества к-го вида:

ЪГ* = - ^. (1.34)

Подсеточный конвективный поток вещества к-го вида:

¥дТк = Р(ид,г¥д,к - Щ,гУд,к)• (1.35) Подсеточный диффузионный поток вещества к-го вида:

Список литературы диссертационного исследования кандидат наук Горкунов Сергей Владимирович, 2023 год

- СТ,

щ

-1 гз

3-1

Р' 3Г г, 3

(К - <)(г? - гв) + КР - Лг* - г*))^?-

2

(2.21)

2.1.3 Граничное условие типа свободная поверхность

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

1. поиск граничных БРИ частиц с последующим взаимодействием с искусственными БРИ частицами вакуума;

2. пересчет плотности частиц итерационно или через увеличенный радиус взаимодействия.

Поиск граничных частиц с последующей разгрузкой представлен в [86]. Для определения границы для каждой БРИ частицы вычисляются взвешенная суммарная длина векторов до ближайших соседей Л г, взвешенный вектор п г указывающий направление отсутствия частиц, и взвешенное количество соседей иг.

С> г = ^ , Л г = ^ Гг] Шг] , Шу = ш( qг з), Ягз

Пг = Г

где т,, р,, -- масса и плотность БРН-частицы, г^ — расстояние между г и ] частицами, ш(д) — весовая функция, обеспечивающая больший вклад ближайшим соседям, определяемая как:

ш(д) = I +2-'1.7)12 0 ^ 2 .

[0 2 <д

После вычисления соответствующих величин проверяется условие:

а, < 38(\цг— 0.07)+ 5. (2.22)

В случае если частица признается граничной, для нее вводится виртуальная частица на расстоянии к, в направлении п. Далее решается задача Римана о разгрузке в вакуум с виртуальной БРН-частицей для нахождения скорости истечения в вакуум рассматриваемой частицы. Найденную скорость необходимо подставить в (2.18) в качестве и*д. При этом вклад от зеркальной частицы учитывается с весом (тах^ев1 — &г). В, обозначает множество, содержащее г частицу и её соседей.

В [86] не представлено обоснование критерия (2.22). Для оценки работоспособности описанного подхода были рассмотрены случаи типичного расположения граничных частиц, представленные на рис. 2.1. При этом считается, что все частицы имеют одинаковую массу т плотность р и расположены в узлах кубической сетки с длиной грани ячейки к.

а) б) в) г)

Рисунок 2.1 — Варианты граничных частиц. Красным цветом показана целевая

частица.

Виртуальная БРН-частица добавляется так, чтобы компенсировать недостающие частицы в области вакуума. В случае если алгоритм построения

виртуальной частицы работает исправно, то для ¿-той частицы должно выполняться:

^ W' (njji) + Wl (Fiji) = S = const, (2.23)

где первое слагаемое соответствует сумме вкладов реальных частиц, а второе слагаемое соответствует вкладу от виртуальной частицы. Согласно [86] вклад виртуальной частицы определяется как:

Wl(h~h) = (maxjeBi(&j) - &i)W'(h~h). (2.24)

В таблице 1 представлены значения S согласно (2.23) для случаев, когда частица не является граничной (столбец 1), изображенных на рис. 2.1. Результаты, представленные в таблице 1, показывают, что значения S, полученные согласно (2.24), отличаются для различных вариантов. Было замечено, что значительно лучшие результаты можно получить, определив вклад виртуальной частицы как:

Wl(hF) = 0.b(maxjeBi(fy) - at)W'(h,h). (2.25)

Таблица 1 — Сумма вкладов от зеркальной и реальных частиц.

Частица не граничная рис. 2.1 (а) рис. 2.1 (б) рис. 2.1 (в) рис. 2.1 (г)

Е7 W'(ггз,Ь) -2.553 -1.918 -1.399 -1.522 -1.631

Wl(h,h) (2.24) 0 -1.211 -2.137 -1.486 -1.566

S (2.24) -2.553 -3.129 -3,536 -3.008 -3.197

Wl(hh) (2.25) 0 -0.605 -1.068 -0.743 -0.783

S (2.25) -2.553 -2.523 -2.467 -2.265 -2.414

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

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

Другим способом учесть разгрузку граничных БРЫ-частиц является пересчет плотности частиц. В этом случае вместо интегрирования уравнения сохранения массы (2.18) плотность на новом временном слое рассчитывается через положения и массы окружающих частиц:

Пересчет плотности может использоваться для частиц из внутренней области, плотность которых меняется мало. Для частиц, граничащих со свободной поверхностью, вычисление плотности по (2.26) приведет к неправильным результатам из-за отсутствия частиц. Согласно [92—94] использование нормализации:

приводит к более точным результатам и позволяет расширить применение на граничную область. Тем не менее, (2.27) корректно использовать только в случае мало сжимаемых веществ, для которых можно считать hi = const.

При моделировании сжимаемых материалов значения hi на новом временном слое также являются неизвестными величинами. В [95] предложено на каждом шаге итерационно методом Ньютона или делением отрезка пополам совместно решать (2.8) и (2.26).

Для исключения необходимости поиска итерационного решения в [91; 96] предложено использовать приближенный подход и для определения плотности в (2.26) подставлять hi.

N

3=1

(2.26)

(2.27)

(2.28)

где

N

р* = ^m3W(гijhi), h* = Csmoothh,

3=1

(2.29)

Csmooth - положительная константа. В [91] Csmooth G [2,3], в [96] Csmooth = 2. В данной работе используется Csmooth = 2. Сравнение и эффективность изложенных методов приводится в разделе 3.1.2.

2.1.4 Алгоритм параллельного вычисления SPH solver

Выделением расчетных областей и распределением нагрузки занимается главный процесс (master). Подчиненные процессы (slave) получают полный объем данных для проведения моделирования от master-процесса. Декомпозиция расчетной области master-процессом схематично представлена на рис. 2.2. В первую очередь происходит разделение расчетной области на слои вдоль одной оси. В случае необходимости выделенные слои могут быть разделены на центральную часть и некоторое количество секторов.

Расчётная область

Рисунок 2.2 — Декомпозиция расчетной области таэ1ег-процессом.

Представим время выполнения одной итерации рассматриваемого алгоритма в последовательном и параллельном варианте:

Fcompiexit y{N

Тр Fcompiexity(N/P) + ^decomposition + T'community,

где N - количество SPH-частиц; Т1 - время итерационного шага при последовательном решении; Тр - время итерационного шага при параллельном решении; р - количество slave-процессов; Fcompiexitу(п) - сложность алгоритма; TdeœmPosition - время, затраченное master-процессом на разбиение расчетной области; Tcommunity - время, затраченное на обмен данными средствами MPI.

Эффективность алгоритма распараллеливания будем оценивать, основываясь на характеристиках ускорения и коэффициента масштабируемости [97; 98].

Ускорение (speedup):

s (р N ) — Т■ Тр

Коэффициент масштабируемости (efficiency):

«PN> - » -

Распределение нагрузки между slave-процессами - линейная операция относительно количества SPH-частиц. Поэтому можно считать, что Tdecomp0sition — a0N, где а0 - некоторая константа.

Предполагая, что толщина слоев, представленных на рис. 2.2, одного порядка с их радиусом, можно считать, что время обмена данными средствами MPI составляет:

2

Т — (N\2

Т community а1р I I 5

\Р/

где а1 - некоторая константа.

Согласно [99] сложность метода сглаженных частиц с учетом использования cell-linked list для поиска соседей составляет Nln(N). Тогда можно записать ускорение и коэффициент масштабируемости в виде:

N ln(N )р

S (pN ) —---—, (2.30)

aip2( f) + a Np + N f) N ln(N )

p2( f) 2/3 + ао Np + N f)

E(P,N) —-r^-—-; г■ (2.31)

ai p\ p

На рис. 2.3 представлено полученное ускорение $(р^) при моделировании задачи о расширении цилиндра со стальной оболочкой с использованием разного количества процессоров. В расчетах использовалось постоянное количество БРИ частиц N * = 10015231.

6 5

К*4

2 1

-1- у'* +

-f

1 2 3 4 5 6

5(р, Л/ * )/Б( 1,Л/*)

Рисунок 2.3 — Зависимость времени вычисления от ускорения для разных р. Точки 1 - время вычисления модуля 8РИ_зо1уег полученное на разных р, линия

2 - аналитическая зависимость (2.30).

2

2.2 Описание решателя GK solver

Решатель СК_эо1уег предназначен для моделирования газодинамических процессов, таких как эволюция ударных волн в окружающей среде с учётом влияния движущихся фрагментов оболочки и препятствий на формирование поля течения газовой фазы. В решателе учитывается возможность диффузионного горения в окружающей среде [6].

2.2.1 Структура исходных данных для GK solver

Исходными данными для решателя могут быть результаты 3И моделирования решателя 8РИ_зо1уег (облако продуктов разложения, фрагменты

метаемой оболочки). В этом случае происходит сопряжение по временному слою БРИ частиц и структурированной расчетной сетки. После сопряжения происходит продолжение численного счета решателем СК_эо1уег.

2.2.2 Метод Годунова-Колгана

Решатель СК_эо1уег предназначен для трехмерного моделирования эволюции ударных волн в открытых пространствах, содержащих препятствия. Численное интегрирование уравнений, представленных в параграфе 1.2, производится методом Годунова-Колгана [100; 101] на прямоугольной сетке в переменных Эйлера.

Для описания метода введем Qn = (ии™, и™, рп, Рп) - вектор параметров на временном слое п. Для методов семейства Годунова разностная схема ячейки с индексами г,],к для отыскания значений на п + 1 слое представляется в виде:

ттп+1 = ттп _

г,],к ^ 1,з,к

Аъ

_ (рт+0.5 _ _ _

л г I Г г+0.5,^',^ * I-0.5,],к + * +0.5,к * 1,3 -0.5,к + * 1,],к+0.5 * 1,],к-0.5 I , 4 7

(2.32)

5 _ гш+0.5 + гш+0.5 _ -рп+0.5 + -рп+0.5 _ *т+0.5

где

ип = и =

/р N

?их риу риг

\Е /

рп+0.5 рп+0.5 рп+0.5 рп+0.5 рп+0.5

Г г+0.5,],к, Г г-0.5,],к, Г г,]+0.5,к, Г 1,3 -0.5,к, Г г,],к+0.51

Р п+0.5 Гг,з,к-0.5

(2.33)

потоки через грани рас-

четной ячейки с индексом г,],к.

Метод Годунова предполагает постоянные параметры в пределах одной ячейки. Для вычисления потоков через границы ячеек используется решения задачи Римана. Способы решения задачи Римана представлены в [102; 103].

Вслед за [101] обозначим Ргт как функцию потоков, определяемую из решателя задачи Римана. Тогда в случае использования схемы Годунова [104] можно записать потоки через границы в виде:

рп+0.5 _ р (^п ^п )

г %+0.5,з,к _ г гт(Ч1^,к, 4{+1^,к),

рп+0.5 _ г (^п ^п )

г %л+0.5,к _ Г гт(Ч1^,к, ^{¿+1,к) ,

гп+0.5 _ г (^п ^п )

г %,з,к+0.5 _ г гт(Ч1^,к,

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

ОЧ±0.5^,к _ ± 0.5АхОЧ^,к, ]±0.5,к _ з,к ± 0.5Ау ОЧ^^,

],к±0.5 _ Цп,з,к ± 0.5А,гОЧ^^•

Тогда потоки через границы ячеек будут рассчитываться по соотношениям:

трп+0.5 _ тр (глп глп \

Г г+0.5, к _ Г гт (Чг+0.5^,к, 4(1+1)-0.5,],к) ,

-лт+0.5 _ ^ (^п ^п )

г ¡¿+0.5,к _ г гт\Ч1,^+0.5,к, ^¿,(^+1)-0.5,к),

трп+0.5 _ тр (ГЛп глп \

Г г,о,к+0.5 _ Г гт (Ч'1,^,к+0.5, Ч^,(к+1)-0.5) •

Для вычисления Qп±0.5jк, Яп,э±{).ь,к, Яп^,к±0.5 Колган предложил использовать принцип минимальных производных, который можно представить в виде:

д, . д, ч ,А /¿-0.5, если 1/¿-0.51 < 1 /¿+0.51

А;г _ тгптоа(А;г-0.5,А/¿+0.5) _ \ •

А/,+0.5, если | /¿-0.51 > | /¿+0.51

где А/,-0.5 _ /г - /г-1, А/г+0.5 _ /г+1 - /г.

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

ттто(1(а,1)) = <

a, если |а| ^ |&| и аЬ > 0

b, если |а| > |&| и аЬ > 0 . 0. если аЬ ^ 0

Схема Годунова-Колгана имеет второй порядок аппроксимации по пространству и первый порядок по времени.

2.2.3 Эффективное использование расчетной области при расчете

параметров ударной волны

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

Известно, что уменьшение размера ячейки по пространству приводит к существенному повышению точности. Однако такой путь повышения точности связан со значительным ростом затрат вычислительных ресурсов. Поскольку для расчета распространения ударных волн на высокопроизводительных вычислительных системах с использованием параллельных вычислений обычно применяются явные схемы, то уменьшение шага по пространству в 2 раза приводит к шестнадцатикратному (24) росту количества действий и, следовательно, к такому же росту потребности во времени вычисления. Использование неравномерной по размеру ячеек расчетной сетки позволяет повышать точность расчета локально, в интересующих местах или в областях с сильно изменяющимися параметрами по пространству. Разумное использование такого подхода позволяет получить расчет с хорошей точностью за приемлемое время. На рис. 2.4 представлены варианты расчетных сеток.

В трехмерной постановке для ячеек в форме параллелепипеда существуют три основных варианта построения неравномерной сетки:

- сгущение/разрежение пространственных шагов вдоль одного направления [105]. Расчетная сетка с использованием такого подхода представлена на рис. 2.4 (б). В этом случае разностная сетка представляет собой совокупность параллелепипедов, примыкающих друг к другу, с изменяющимися длинами

№ К

а) б) в) г)

Рисунок 2.4 — Варианты расчетных сеток. а) равномерная сетка; б) сетка со сгущением/разрежением пространственных шагов; в) AMR сетка; г) перестраивающаяся сетка.

граней. Изменение шага происходит сразу во всей плоскости, перпендикулярной направлению изменения шага. Поэтому, если сгущение/разрежение ячеек требуется только в одной локальной части пространства, то для его создания ячейки будут сгущены/разрежены во всем соответствующем поперечном сечении, что приводит к большему увеличению числа ячеек, чем этого требуют пространственные особенности. Координаты узлов сетки со временем не меняются. Именно такой подход используется в коммерческом коде FLACS [106]. Следует отметить, что такое использование разностной сетки с переменным размером ячейки обычно снижает порядок разностной схемы до первого, что может приводить к снижению точности разностного метода в практических расчетах; - AMR (adaptive mesh refinement) метод [107—109]. Вложение сеток производится путем разделения крупных ячеек на несколько более мелких, как показано на рис. 2.4 (в). Такое деление, как правило, охватывает не одну ячейку, а некоторую область. В этом случае разностная сетка представляет собой совокупность областей в форме параллелепипеда, примыкающих друг к другу, при этом в одних областях шаг меньше, чем в других. Обычно области с меньшим шагом оказываются вложенными в области с большим шагом, как матрешки. Такая организация разностной сетки (в отличие от сгущающихся/разрежающихся структурированных сеток) позволяет организовать мелкий шаг лишь в ограниченной трехмерной области, а не во всем поперечном сечении расчетной области. Следует отметить, что такое построение разностной сетки требует использования довольно сложных алгоритмов иерархической организации данных и вычисления потоков через границы ячеек; также для такой сетки могут потребоваться и сложные алгоритмы перемещения области с вложенной сеткой; т.е. работа с вложенными сетками требует существенной переработки кода.

Наиболее эффективны такие сетки в задачах, где требуется точно рассчитать контактную поверхность без выделения ее явным образом. Типовым примером такой задачи может служить задача о диффузионном горении. - перестройка сетки (пример такой сетки приведен на рис. 2.4 (г)); в этом случае разностная сетка представляет собой совокупность параллелепипедов, примыкающих друг к другу, причем по направлениям может существовать как сгущение/разрежение ячеек, так и равномерное их разбиение. Последнее даже предпочтительнее, поскольку позволяет избежать возможного снижения точности (порядка точности) разностной схемы. Перестраивающаяся сетка заведомо эффективна в тех случаях, когда со временем процесс развивается по пространству, охватывая все новые и новые области; в этом случае нет необходимости проводить расчет сразу на всей области, поскольку параметры среды где-то и не меняются вовсе. В начале моделирования выбирается сравнительно небольшая область и расчет ведется только в ней; когда течение развивается настолько, что вот-вот выйдет за границы расчетной области, размеры сетки увеличиваются на заранее заданную величину, чтобы позволить течению развиваться дальше, по-прежнему не взаимодействуя с границами расчетной области; и расчет ведется до тех пор, пока течение в своем развитии во времени опять не приблизится к границе расчетной области; и в этом случае опять происходит перестройка сетки. Такой порядок действий «расчет - контроль за подходом течения к границе расчетной области на заданное расстояние - перестройка сетки - возобновление расчета на новой сетке» повторяется до тех пор, пока не будет достигнут заранее заданный критерий остановки расчета. Наиболее характерным примером задачи, где может быть использована перестраивающаяся сетка, является расчет распространения ударной волны. В этом случае критерием перестройки сетки является приближение ударной волны к границе расчетной области. При перестройке сетки обеспечивается сохранение массы, импульса и энергии. Для этого сначала решается чисто геометрическая задача о том, какие ячейки из старой (более мелкой) сетки войдут в каждую ячейку новой (более крупной) сетки. Для ячейки новой сетки будем использовать индекс %п. Для множества ячеек старой сетки , которые входят (полностью или частично) в ячейку новой сетки с индексом гп, будем использовать индекс г80, где в = 1..М{п. Можно определить объемную долю ячейки старой сетки , которая войдет в ячейку гп новой сетки. Будем обозначать такую объемную долю как а8. У%п, обозначим объемы новых и старых ячеек сетки.

При перестройке сетки в ячейке гп новой сетки для продолжения расчета следует задать следующие значения плотности:

1 "

Р = ^Е (Р Л a (2.34)

составляющих скоростей u,v,w:

Уг ,

%n s=1

1 "

^n = aS(2.35)

K^n e=i

1 "

Уг. = -~Е(Уго Р г30Уг3о ^ ), (2-36)

Р г П^гП £=1

1 w

Win = ^Ргг^ а). (2-37)

п п =1

полной энергии Е:

1 м

Егп = —^Е а). (2.38)

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

2.2.4 Алгоритм параллельного вычисления GK solver

Для организации параллельных вычислений используется подход, когда все процессы параллельно выполняют копии одной программы. Такая модель программирования называется SPMD (Single Program, Multiple Data) [110—112]. Ввод/вывод данных задачи осуществляется процессами независимо. Такой подход позволяет добиться равномерного распределения затрат оперативной памяти по всей вычислительной системе. В этом случае эффективность расчета зависит от времени, затраченного на обмен информацией между процессами.

Время затрачиваемое на обмен информацией определяется количеством пересылок, которые должен выполнять каждый процесс, и объемом информации, отправляемой за одну пересылку.

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

абв Рисунок 2.5 — Варианты декомпозиции расчетной области по (а) одной (б) двум

(в) трем осям.

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

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

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

Рисунок 2.6 — «Перекрестные пересылки» вычислительных областей, которые рассчитываются на разных вычислительных процессах.

После этого обмен данными происходит по следующему алгоритму. Во-первых, процессы обмениваются информацией на поверхностях, перпендикулярных оси X, затем оси У и после этого оси 2. На обрабатываемой оси процесс отправляет информацию в процесс, который обслуживает от него "левую"область и получает данные от этого процесса. После этого аналогичные операции выполняются с процессом, который находится «справа».

Время выполнения одной итерации рассматриваемого алгоритма в последовательном и параллельном варианте можно записать в виде:

Fcompiexity (N

Fp Fcompiexity (^/P) + FCommunity,

где N - количество расчетных ячеек, F1 - время итерационного шага при последовательном решении, Тр - время итерационного шага при параллельном решении, р - количество slave процессов, FCOmpiexity(п) - сложность алгоритма, F community - время, затраченное на обмен данными средствами MPI.

Сложность метода Годунова-Колгана зависит линейно от количества ячеек. В [113] показано, что в рассматриваемом алгоритме при декомпозиции расчетной области по трем осям FCOmmunity = ai VN2 (tfp - 1), тогда:

Ti = aoN,

Fp = ^ + ai^ (tfp - 1)

где <Хо, ai - некоторые положительные константы.

В этом случае ускорение (speedup) запишется в виде:

S (р N) = ^ = р

Тр 1 + а (ур - 1)'

ao^N

Коэффициент масштабируемости (efficiency):

в (:pN) = s w) = 1

р 1 + о (Ур - 1)'

аоУЫ

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

2.3 Сопряжение решателей SPH solver и GK solver

Целью сопряжения решателей является продолжение численного счета существенно другим численным методом, используя другой решатель в определенный момент времени. Для смены численного метода необходимо заполнить расчетную сетку решателя СК^о1уег параметрами среды, определенными распределением БРН-частиц в решателе SPH_so1veг. Для этого предлагается найти объемные доли БРН-частиц, попадающих в каждую ячейку. Будем считать,

что БРН-частица занимает объем куба со стороной (*)3, где т и р - масса и

плотность БРН-частицы соответственно. Обозначим оу - объемную долю ^'-ой БРН-частицы в к-ой расчетной ячейке. Тогда можно представить объемную долю газовой фазы в расчетной ячейке с индексом к в виде:

Ф дк = 1 - Е а-к,о, (2.39)

где Ы8 - множество БРН-частиц, описывающее твердую фракцию.

Плотность в к-ой ячейке расчетной области определяется как:

Рк =

Ф

^(.Р 3&к, 3) + Р Агг ( Фдк - ^ ^к

3=1

3=1

где - множество БРЫ-частиц, описывающее газовую фракцию, р^ ность ]-ой БРЫ-частицы, р а^ - плотность окружающей среды.

(2.40)

плот-

Вычисление составляющих скорости в к-ой ячейке расчетной об-

ласти:

Пк =

Ф д к Р к

^ / Ъ

№ ак + РА гПАгг ( Фд к - ^ ^к,з =1 =1

Ук =

Фдк Рк

^(РЯ,ак,з) + РАггЪАгг ( Ф дк - ^

ак,з

=1 =1

(2.41)

(2.42)

Ык =

Ф д к Р к

Щ I ^

^(Р№ ак,з) + РА , г™ А г г ( Ф дк - ^

ак,з

=1 =1

(2.43)

где и^^^,Wj - х,у,г компоненты скорости ]-ой БРЫ-частицы, имг^Ак-^мг -компоненты скорости окружающей среды.

Аналогично можно записать выражение для полной энергии в к-ой ячейке расчетной области:

Ек =

Ф д к Р к

^(РзЕ3 ак + Р А г г-ЕА г г ( Ф дк - Е ак,

3=1

=1

(2.44)

где Е^ - полная энергия ]-ой БРЫ-частицы, Еа¡г - удельная полная энергия окружающей среды.

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

1

1

1

1

1

2.4 Решатель для расчета течения с фазовыми переходами в

пористой среде

Расчетный модуль предназначен для решения задачи (1.90) описанной в параграфе 1.5. Рассматривается пространство, занятое пористой средой. Внутри пористой среды находится область О в Я2, насыщенная жидкостью. Граница этой области состоит из двух незамкнутых линий Г1 и Г2, не имеющих общих точек. Эти линии заданы в виде однозначных функций ^ = й ¡(х), где г = либо в параметрическом виде ^ = й (£),х = (£,) где £ -- параметр, отсчитываемый от некоторой точки на линии.

Будем считать, что функция й ¡(х), является периодической с периодом Тх и выполняется равенство

5 ¡(х) = зг(х + Тх), (2.45)

а соответствующее решение уравнения Лапласа

АР = 0 (2.46)

так же является периодической функцией с тем же периодом Тх.

Выделим в области О, область О1, с границей Га, которая определяется функцией <§1(х) при XI < х < х1 +Тх, границей Гс, которая задается функцией й 2(х) при х1 < х < х1 +Тх и границами Гъ и Г^. Граница Гъ имеет координату х = х2 и координату ^ € (я 1(х2),з2(х2)), а для точек границы Г^ выполняется условие х = х1 иг € (в 1(х1), 52(х1)).

Пусть на границе Га давление постоянно и равно Р1, а на границе Гс постоянно и равно Р2.

На границах Гъ и Г^ давление и скорость фильтрации неизвестны, но равны в силу предположения о периодичности решения, т.е.

Р (г ,х1 ) = Р (г ,х1 +Тх) (2.47)

и

УР (х х)х=х1 = УР (г х)х=хл+тх. (2.48)

В рассматриваемом случае возможно непосредственное применение численного метода граничных элементов [114], поскольку область О1 имеет

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

Заменим границу Г области О1 ломаной линией Гп, состоящей из Мп отрезков (панелей) Г длиной Ц, соединяющих точки (х^,у^) и (х^,у^+1). Панель Г можно описать параметрически следующим образом: х = х(£),у = у(£), - 0.5Ц < £ < 0.5Ц.

Каждой из панелей Г поставим в соответствие источник потенциала двойного слоя, плотность интенсивности которого ф^ изменяется от края панели до ее середины по линейному закону:

ф (£) = 14 - (Щ - Щ'+О0 < £ ^ 0-5Ь

Ф \щ + (щ - и-)-°.Ы3 ^ £ ^ 0 '

Здесь £ - координата, которая отсчитывается от середины панели с номером ] и принимает значения от -^/2 до Ц/2, а Uj - значение ф в центре панели Г.

Зависимость (2.4) обеспечивает непрерывное распределение интенсивности источников потенциала ф вдоль ломаной линии Г п.

Решение представим в виде суперпозиции решений Фj(х^^^щ,Щ+1) для потенциала, создаваемого в точке (х,у) каждой панелью по отдельности. Тогда функцию распределения давления можно записать как

м

Р (х1 У) = ЕФ (х,У,из-ъЩ ,из+1)- (2.49)

j=1

Решение

Фj(х,у,Щ-1,Щщ+1)= Ф(х,ух(£)ф(£))ф0(£) (2.50)

является потенциалом двойного слоя и решением уравнения Лапласа вне Г;

(£)) = 2п [^ФЩфФШШФ(£))) • (2.51)

где г(х,у,х^(£)ф(£)) - вектор, проведенный из точки (хj(£)ф(£)) в точку (х,у), а п(хíj(£),ф(£)) - орт вектор нормали к Г в точке (хj(£)ф(£)).

Функция Фj(х,y,Uj-1,Uj,Uj+1) линейно зависит от переменных ир -1),щ,ир + 1). Эти величины определяются с помощью граничных условий.

Для этого потребуем выполнения граничных условий в середине каждой панели, то есть в контрольных точках с координатами (х]+1/2,Уз+1/2), где х]+1/2 = 0.5(xj+1 + х^),у+1/2 = 0.5(у^+1 + у^). В результате получим линейные уравнения относительно неизвестных Uj.

м

и. ^

р (хЗ+1/2-У]+1/2) = 2 + (хЗ+1/2-УЗ+1/2-Ц—1-Ц ЦЗ+1) (2.52)

j=l

Запишем полученные уравнения применительно к сформулированным выше граничным условиям.

Пусть на границе Га находятся панели с номерами от 1 до га. Тогда

м

и. ^

р1 = ~2 + (Ъ+1/2Ы+1/2Щ-иЩ Ц]+1)- (2.53)

j=1

при 1 ^ % ^ га.

Если на границе Гс находятся панели с номерами от ц + 1 до ъс, то

м

и. ^

р2 = 2 + ^З+^УЗ+ФЦЗ — Щ Ц]+1)- (2.54)

j=1

при %Ь + 1 ^ % ^ %с.

Построим ломаную линию, соответствующую участкам границы Гь и Г'й таким образом, чтобы панели на этих границах имели одинаковый размер и панель с номером %а + т такую же у-координату, что и панель с номером М — т + 1, где т ^ ц — %а. Тогда условие равенства давлений в точках с одной и той же координатой у на линиях Гь и Г' и можно записать в виде

Щ

м

2 + (хк+1/2- Ук+1/2 ,Щ—1 Цк Цк+1) =

3=1

м

и + (х з+1/2-1 УЗ+1/2,цЗ—1 -Ц] ЦЗ+1)- (2.55)

=1

при ъа + 1 ^ % ^ ц и ц = М — (г — ъа) + 1.

Кроме того, на Гь и Г' должно иметь место равенство скоростей (градиентов давления) в точках с равной у-координатой. Условие (2.55) обеспечивает

равенство у-компонент градиента давления, поэтому достаточно записать условие для равенства х-компонент:

м д

^ дхФ (х-Ук+1/2-ик—1-ик Цк+1)х=хк+1/2 =

=1

м д

= Е дхФ (х-У3+1/2-Ц — Щ -Щ+1)х=х0+1/2 . (2.56)

=1

Система уравнений (2.53), (2.54), (2.55) и (2.56) обеспечивает выполнение граничных условий в контрольных точках, состоит из М линейных уравнений и решается относительно неизвестных ит.

Метод граничных элементов был выбран именно ради возможности вычислять скорость движения границы аналитически. Если бы использовался метод конечных разностей с разностной сеткой, заполняющей всю область то производные пришлось бы находить приближенно с помощью сеточных функций давления. При использовании представленного подхода производные найденного приближенного решения вычисляются аналитически и в расчет не вносится дополнительная ошибка. При расчете движения границы Ь предполагается, что каждая панель двигается в течении малого времени т (шага интегрирования по времени) с постоянной скоростью V*, направленной по нормали. Метод вычисления нового положения границ панели изменяется незначительно. Вычисления проводятся только для одного периода изменения координаты х. Метод граничных элементов применялся при исследование на устойчивость фазового перехода в пористых средах [13; 18].

Глава 3. Валидация и верификация расчетных модулей

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

При работе над данным разделом диссертации использованы следующие публикации автора, в которых согласно Положению о присуждении ученых степеней МИФИ отражены основные результаты, положения и выводы исследования:

— Трехмерное моделирование образования ударных волн и разгона металлических оболочек высокоэнергетическими соединениями / В. А. Шар-гатов, С. В. Горкунов // Горение и взрыв. - 2021. - Т. 14. - № 2. - С. 92-99. - DOI 10.30826/CE21140210.

- Stability of finite perturbations of the phase transition interface for one problem of water evaporation in a porous medium / V. A. Shargatov, S. V. Gorkunov, A. T. Il'ichev // Applied Mathematics and Computation. - 2020. V - 378. № 125208.

3.1 Верификация SPH solver

3.1.1 Распад разрыва в упругой среде

Рассматривается задача о распаде разрыва продольной компоненты скорости в упругопластической среде. В качестве среды используется сталь со следующими характеристиками: начальная плотность ро = 7850 кг/м3, модуль объемного сжатия К = 156 ГПа, модуль сдвига С = 83 ГПа, предел текучести УО = 0.3 ГПа. Среда описывается уравнением состояния вида: Р(р) = К(1 - ро/р).

В начальный момент времени в области х < 0 заданы следующие параметры: Р = 0, р = ро, их = 300 м/с„ иу = 0 иг = 0, = 0, а в области х ^ 0: Р = 0, р = р0, их = 0, иу = 0, и = 0, Бу = 0. Здесь Р - давление, р - плотность, их, иу, иг - х,у,х компоненты скорости, - компоненты девиатора напряжений.

При моделировании использовалось по 12 частиц вдоль осей У иZ, и 1000 частиц вдоль оси X. На рис. 3.1 представлены результаты расчета в сравнении с аналитическим решением на момент времени £ =0.63 мс. Полученные в ходе расчета профили плотности и скорости с приемлемой точностью описывают аналитическое решение.

300-

200-

8100

8050

СП

8000

ад

7950

с!

7900

7850

3 100-

- 2

0

х, т

- 4

-2

0

х, т

аб Рисунок 3.1 — Распределение (а) плотности, (б) скорости вдоль оси X в резуль-

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

0

4

2

4

2

4

3.1.2 Истечение газа в вакуум

Рассмотрим следующую задачу. В начальный момент времени в области в виде шара радиуса Я0 находится газ под высоким давлением. Снаружи шара находится вакуум. Обозначим индексом 1 параметры газа в начальный момент времени. Газ описывается уравнением состояния идеального газа.

В расчете использовались следующие параметры: плотность р! = 1800 кг/м3, давление Р! = 30 ГПа, радиус шара Я0 = 0.1 м, показатель адиабаты у = 3. Моделирование методом БРИ выполнялось в трехмерной постановке с равномерным укладыванием частиц в начальный момент времени в узлах

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

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

— Расчет с поиском и разгрузкой частиц на границе согласно (2.22)- (2.25).

— Моделирование выполнено с пересчетом плотности, используя больший радиус сглаживания согласно (2.26).

— Для оценки качества расчетов, выполненных SPH методом, получено решение одномерным кодом [115—117].

На рис. 3.2, 3.3 представлены полученные распределения плотности и модуля скорости соответственно.

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

3.2 Верификация решателя GK _ solver

Верификация решателя GK_solver выполнена на основе признанной научным сообществом методики [118; 119]. Принципиальным моментом методики является использование для тестов только тех задач, которые имеют аналитическое решение. Тем самым устраняется традиционное сомнение в сходимости именно к решению дифференциальной задачи, а не к какому-либо другому. Такое сомнение возникает, если сходимость проверяется путем измельчения сетки. Для верификации отобраны тесты, которые могут быть применены для задач, решаемых в переменных Эйлера в трехмерной постановке.

аб Рисунок 3.2 — Распределение плотности от радиуса в момент времени (а) 20

мкс, (б) 50 мкс. Линия 1 - метод SPH без использования граничных условий,

2 - вычисление плотности согласно (2.26), 3 - расчет с поиском и разгрузкой

частиц на границе (2.22)- (2.25), 4 - решение одномерным SIN кодом.

аб Рисунок 3.3 — Распределение модуля скорости от радиуса в момент времени (а)

20 мкс, (б) 50 мкс. Линия 1 - метод SPH без использования граничных условий,

2 - вычисление плотности согласно (2.26), 3 - расчет с поиском и разгрузкой

частиц на границе (2.22)- (2.25), 4 - решение одномерным SIN кодом.

3.2.1 Галилеевская инвариантность.

В тесте рассматривается галилеевская инвариантность двумерных и трехмерных распределений. Цель теста: проверка сохранения распределения плотности при переходе в подвижную систему координат. Известно, что численный метод Годунова-Колгана приближенно инвариантен на двумерном