Содержание к диссертации
Введение
1 Численные методы решения многомерных задач нестационарной газовой динамики в многокомпонентных системах 9
1.1 Лагранжевы методы 12
1.2 Лагранжевы методы с перестройкой сетки 14
1.2.1 Метод Годунова на подвижных криволинейных сетках 15
1.3 Методы с перестройкой связей между лагранжевыми узлами 17
1.3.1 Метод «сглаженных частиц» SPH 18
1.4 Методы «частиц в ячейках» 21
1.5 Эйлеровы методы 25
1.5.1 Алгоритмы отслеживания контактных и свободных границ тел на эйлеровой сетке 27
1.5.2 Алгоритм адаптивного уточнения сетки 30
2 Метод индивидуальных частиц для расчета газодинамических многокомпонентных нестационарных течений с большими деформациями 32
2.1 Общая схема расчетной процедуры 34
2.2 Предварительный шаг 36
2.3 Основной шаг 41
2.4 Дробление/объединение частиц 42
2.5 Алгоритм определения ориентации контактных/свободных границ тел 44
2.6 Граничные условия 45
2.6.1 Граничные условия на границах эйлеровой сетки . 45
2.6.2 Граничные условия на внутренних поверхностях раздела 49
2.7 Интегрирование по времени 51
2.8 Параллельная реализация для многопроцессорных ЭВМ с распределенной памятью 52
2.8.1 Многопроцессорные ЭВМ и параллельные вычисления 52
2.8.2 Параллелизм расчетного алгоритма 59
2.9 Примеры тестовых расчетов 62
3 Численное моделирование процесса высокоскоростного удара 69
3.1 Постановка задачи 70
3.2 Результаты численного моделирования 71
3.3 Обсуждение результатов 78
4 Численное моделирование воздействия интенсивных пучков тяжелых быстрых ионов на вещество 82
4.1 Постановка задачи 83
4.2 Расчет энерговклада в ячейках сетки 85
4.3 Энерговклад в смешанных ячейках 86
4.4 Параллельная реализация 87
4.5 Численное моделирование воздействия пучка ионов урана на тонкую свинцовую фольгу 94
Заключение 104
Литература 106
- Алгоритмы отслеживания контактных и свободных границ тел на эйлеровой сетке
- Алгоритм определения ориентации контактных/свободных границ тел
- Многопроцессорные ЭВМ и параллельные вычисления
- Численное моделирование воздействия пучка ионов урана на тонкую свинцовую фольгу
Введение к работе
Настоящая диссертация посвящена численному моделированию процессов при высоких плотностях энергии на современных высокопроизводительных параллельных ЭВМ.
Актуальность. Исследования физических процессов и свойств вещества при высоких плотностях энергии имеют ряд особенностей. Автомодельные решения существуют для ограниченного числа теплофизических процессов и газодинамических течений, как, например, для случаев из-энтропического сжатия и расширения, одномерного сжатия стационарной ударной волной и т.п. [1]. В задачах мощного импульсного энерговыделения (менее Ю-5 с) газодинамическое течение, как правило, неодномерно и нестационарно, и анализ получаемой в таких условиях информации требует привлечения аппарата численного моделирования [2, 3]. Использование при компьютерном моделировании современных эффективных разностных схем совместно с моделями реальных свойств веществ позволяет описать результаты лабораторного эксперимента с высокой степенью точности и достоверности. Отметим в этой связи, что большое число практически важных задач, как, например, разработка противометеоритной защиты космических аппаратов [4], недоступно для изучения в лабораторных условиях. Появившиеся в последнее время новые экспериментальные методы физики экстремальных состояний, в частности, использование интенсивных пучков тяжелых ионов для нагрева вещества [5, 6], и полученные с их помощью
результаты требуют тщательного теоретического исследования. В такой ситуации численное моделирование с помощью кодов, ранее проверенных при решении тестовых задач лабораторного эксперимента, является единственно доступным инструментом исследования. Вышеперечисленные обстоятельства определяют актуальность настоящей работы.
Целью работы является развитие численных алгоритмов решения системы уравнений газовой динамики в трехмерной постановке с учетом специфики параллельных вычислений, реализация разработанных алгоритмов для современных многопроцессорных вычислительных систем, численное моделирование процессов при высоких плотностях энергии, таких как высокоскоростное пробивание и воздействие интенсивных пучков тяжелых ионов на вещество.
Диссертация состоит из введения, четырех глав и заключения с основными результатами работы.
Во введении описана структура диссертации.
Первая глава содержит обзор современных численных методов решения многомерных задач нестационарной газовой динамики в многокомпонентных системах, применимых для моделирования задач физики высоких плотностей энергии. Обзор ограничивается лишь теми методами, на основе которых разработаны не просто программы для решения той или иной задачи, а были созданы пакеты прикладных программ, позволяющие проводить расчет широкого круга задач, с возможностью быстрого перехода от задачи к задаче.
Вторая глава посвящена подробному описанию модифицированного метода индивидуальных частиц для расчета газодинамических многоком-
понентных нестационарных течений с большими деформациями в трехмерной постановке. Подробно рассматривается расчетная схема, алгоритмы определения положения и ориентации внутренних границ раздела сред, с дальнейшей постановкой на них особых граничных условий. Также описывается специфика параллельной реализации метода для многопроцессорных ЭВМ с распределенной памятью, приводятся результаты тестовых расчетов.
Третья глава содержит результаты многомерного численного моделирования процесса высокоскоростного удара. Описывается постановка задачи, особенности, возникающие при ее решении. Обсуждаются результаты численного моделирования в сравнении с экспериментальными данными. Показано влияние используемой в расчетах модели уравнения состояния вещества на результаты моделирования.
Четвертая глава посвящена проблеме численного моделирования трехмерных газодинамических процессов, происходящих при воздействии на вещество интенсивных пучков тяжелых быстрых ионов. Предлагается численный алгоритм учета динамики энергетических потерь ионов в веществе, реалистичного трехмерного пространственного распределения параметров и временной зависимости интенсивности импульса пучка. Приводятся результаты численного моделирования процесса нагрева тонкой свинцовой фольги ионами урана.
В заключении сформулированы основные результаты работы.
Практическая ценность. Разработанные в работе модели, методы, алгоритмы и газодинамические программные коды являются современным и эффективным инструментом для исследования практически важных задач
физики высоких плотностей энергии. Практическая ценность работы определяется использованием полученных результатов для решения прикладных задач в ИПХФ РАН, ИТЭС ОИВТ РАН, Институте тяжелых ионов (GSI, Германия).
Научная новизна состоит в следующем: на основе двумерного метода индивидуальных частиц разработан численный алгоритм для решения трехмерных задач газовой динамики в многокомпонентных системах который реализован для режима параллельных вычислений. С его помощью изучено влияние моделей уравнений состояния вещества на результаты математического моделирования высокоскоростного пробивания и выполнено численное моделирование экспериментов по исследованию свойств металлов, нагреваемых интенсивными пучками тяжелых ионов.
На защиту выносятся следующие положения:
Модифицированный метод индивидуальных частиц для решения трехмерных задач нестационарной газовой динамики в многокомпонентных системах для моделирования процессов физики высоких плотностей энергии;
Численный алгоритм расчета энерговклада интенсивных пучков тяжелых быстрых ионов в трехмерных газодинамических процессах, с учетом динамики энергетических потерь ионов в веществе, реалистичного трехмерного пространственного распределения параметров и временной зависимости интенсивности импульса пучка;
Результаты численного моделирования процесса высокоскоростного удара, в которых показано влияние используемого в расчете уравнения состояния на результаты расчетов;
Результаты трехмерного численного моделирования нагрева интенсивным пучком ионов урана тонкой свинцовой фольги в трехмерной постановке.
Апробация работы. Результаты исследований докладывались и обсуждались на Всероссийской научной конференции «Высокопроизводительные вычисления и их приложения» (Черноголовка, 2000), Международной конференции «Математическое моделирование» (Самара, 2001), третьей Всероссийской молодежной школе «Суперкомпьютерные вычислительно-информационные технологии в физических и химических исследованиях» (Черноголовка, 2001), Международной конференции «Уравнения состояния вещества» (Приэльбрусье, 2002 и 2004), Международной конференции «Воздействие интенсивных потоков энергии на вещество» (Приэльбрусье, 2003 и 2005), Международной конференции «Parallel Computational Fluid Dynamics 2003» (Москва, 2003), Научно-координационном совещании-симпозиуме «Проблемы физики ультракоротких процессов в сильнонеравновесных средах» (Новый Афон, 2003 и 2004), Международном семинаре «Супервычисления и математическое моделирование» (Саров, 2004), а также на научных семинарах Межведомственного суперкомпьютерного центра, ИПХФ РАН, ИММ РАН и ИТЭС ОИВТ РАН.
Публикации. По теме диссертации опубликовано 18 научных работ [7-24].
Благодарности
Автор выражает глубокую признательность своим научным руковожу
дителям Берзигиярову Парвазу Куламовичу и Ломоносову Игорю Влади
мировичу за терпеливое и своевременное руководство, предоставленную
творческую свободу и всяческую поддержку на всем протяжении работы.
Сотрудникам лаборатории Уравнений состояния вещества: Острику А.В.,
» Шутову А.В., Султанову В.Г. и коллективу отдела Экстремальных состо-
яний вещества ИПХФ РАН в целом — за многочисленные полезные обсуждения результатов диссертации, помощь при её подготовке, а также за приятное и плодотворное совместно проведенное время и хорошее отношение. Родителям — за веру, терпение и поддержку.
Алгоритмы отслеживания контактных и свободных границ тел на эйлеровой сетке
Перед рассмотрением особенностей и возможностей группы методов использующих идеологию частиц следует указать, что основополагающий метод этой группы — метод «частиц в ячейках» Харлоу [78] (РІС — «Particle in Cell») по имеющимся в литературе сведениям был исторически первым численным методом, позволившим разработать пакет достаточно универсальных программ для решения широкого круга нестационарных многомерных задач механики сплошных сред [79].
Происхождение этого метода отличается от происхождения других методов тем, что при его разработке основное внимание первоначально уделялось не столько решению дифференциальных уравнений в частных производных, сколько моделированию основных физических процессов с помощью рассмотрения сплошной среды как дискретных движущихся частиц [80]. При этом указанные частицы фактически являются лагранже-выми и используются для определения как параметров текущей жидкости так и для описания движения контактных границ. В то же время для определения значений гидродинамических полей эта группа методов использует эйлерову сетку. Во всех методах этой группы продвижение решения на один слой по времени производится в несколько этапов (в исходном варианте Р1С-метода — в два этапа). Как было впервые отмечено в работе [81] — фактически имеет место расщепление по физическим процессам. При этом для всех этих методов характерным является то обстоятельство, что конвективный перенос полевых переменных осуществляется путем использования дискретных частиц, а расчет нужных пространственных производных осуществляется на фиксированной эйлеровой сетке. Таким образом можно сказать, что здесь недостаток лагранжевых методов связанный с искажением лагранжевой сетки преодолевается не за счет локального изменения топологии как в методах предыдущей группы, а за счет «усреднения» параметров по ячейкам эйлеровой сетки, которое автоматически учитывает изменения связей между лагранжевыми частицами при их перетекании из ячейки в ячейку. В то же время дискретное представление элементов среды позволяет естественным образом учитывать контактные и свободные границы внутри расчетной области: точное знание типа вещества каждой частицы делает возможным проводить расчет даже при сильном перемешивании, когда само понятие контактной границы теряет всякий смысл.
Однако, несмотря на все усовершенствования, нужно указать, что по экономичности (как в отношении требуемой расчетной памяти так и по объему вычислений) эта группа методов заметно уступает чисто лагранжевым методам, на задачах, решение которых под силу последним. По существу, метод частиц увеличивает на единицу размерность задачи [80]. Еще один недостаток с точки зрения «безавостности» состоит в том, что при расширении вещества в волне разгрузки среднее число частиц на одну ячейку может стать столь малым, что в объеме образуются искусственные пустоты — эйлеровы ячейки с нулевой плотностью не содержащие частиц. Следует также отметить недостаток, связанный с дискретным представлением конвективных потоков: при использовании малого числа частиц на ячейку или высоких градиентов плотности в решении возникают сильные нефизические осцилляции. Во избежание этого приходится заранее учитывать (если это возможно) положение областей сильного разрежения, первоначально создавая в них большое число частиц, что приводит к завышенным требованиям к количеству используемых аппаратных ресурсов. Также в известной степени спасает положение введение в рассмотрение частиц имеющих конечные размеры и форму [82]. При этом переход частицы из ячейки в ячейку происходит более плавно, перенос всей ее массы осуществляется за несколько шагов по времени. Однако, при этом усложняется задача расчета ячеек, содержащих несколько веществ и локализации контактных границ.
Вообще, проблема локализации границ контакта в этой группе методов достаточно актуальна. Во первых, в том случае когда для расчета давления используются переменные с эйлеровой сетки, возникает необходимость определения давления в ячейках, содержащих несколько веществ с различающимися уравнениями состояния. Вторая группа затруднений связана с необходимостью обеспечения непрерывности нормальной к контактной границе компоненты скорости частиц. В третьих, существуют сложности в алгоритмах вычисления геометрических характеристик контактных границ (вектора нормалей, касательных плоскостей и т.д.).
За последние 30 лет на базе основополагающей идеи, которая была заложена в исходный метод РІС, как за рубежом так и в нашей стране была проделана большая работа как по совершенствованию самого этого метода, так и по разработке существенно отличающихся от него методов. Так, например, модификации метода активно и успешно применялись для расчетов процессов горения и детонации Мейдером [83].
Алгоритм определения ориентации контактных/свободных границ тел
После вычисления нового положения и размеров частиц на предыдущем этапе может возникнуть такая ситуация, когда в одной эйлеровой ячейке окажутся центры нескольких частиц с одинаковыми идентификаторами (т.е. частицы одного и того же сорта), и информация о течении в данной ячейке будет избыточной. В областях же расширения вещества из-за сильного увеличения объема и линейных размеров частиц некоторые ячейки могут оказаться пустыми (т.е. в них не будет содержаться ни один центр ни одной частицы).
Во избежание этих ситуаций, а также с целью увеличения степени однородности данных хранящихся в памяти ЭВМ, вводится этап дробления и объединения частиц. После того как вычислены значения координат центров частиц: xn+1, уп+1 и zn+l и их геометрические размеры: dxn+1, dyn+1
Дробление частиц (приводится для двумерного случая). Исходная частица, граница которой - сплошная линия, а центр обозначен черным квадратом, дробится вдоль границ ячеек (координатные линии проходящие через узел (f, j)). Здесь же черными кружками обозначены центры четырех вновь образованных частиц, а пунктиром — их границы. и dzn+1 происходит дробление частиц по границам-плоскостям эйлеровых ячеек (рисунок 2.4). Их параметры (удельные) равны параметрам исходной частицы за исключением координат центров и геометрических размеров, которые пересчитываются из очевидных соотношений.
Объединение частиц с одинаковыми идентификаторами, центры которых расположены в одной ячейке, производится попарно с удовлетворением условий сохранения «центра объемов», объема, массы, импульса и полной энергии (рисунок 2.5). В результате параметры объединенной частицы определяются как: J+1"Объединение частиц (приводится для двумерного случая). Исходные частицы с номерами 1 и 2, центры которых обозначены черными кружками. Граница частицы полученной в результате их объединения показана пунктиром; ее центр — черным квадратом.
Отношения размеров новой частицы полагаются равными соответствующим отношениям размеров эйлеровой ячейки (для кубической ячейки dx = dy = dz = \/V).
Внутри расчетной области для постановки граничных условий проскальзывания материалов друг относительно друга вдоль контактных границ, или на свободной поверхности при расчете течений сред с нешаровым тензором напряжений, необходимо знать ориентацию границ в пространстве. В описываемом методе для этого используется следующая процедура. Единичный вектор нормали N, задающий ориентацию границы, определяется для каждого узла эйлеровой сетки, имеющего маркер контактной или свободной границы. В первом случае вектор нормали равен нормированной сумме единичных векторов, ориентированных из рассматриваемого узла в направлении тех из 26-ти (8-ми в двумерном случае) ближайших узлов, которые содержат маркер какого-либо одного тела. Во втором случае, аналогично, производится сложение векторов, ориентированных на ближайшие «вакуумные» узлы.
Точность определения вектора нормали к поверхностям можно повысить за счет усреднения направлений нормалей в соседних узлах на контактной/свободной границе.
После предварительного шага в узлах эйлеровой сетки определены значения давления, компонент скорости, тип узла и если узел является узлом на свободной или контактной поверхности — единичный вектор нормали к соответствующей поверхности. Граничные условия на эйлеровой сетке можно разделить на две группы: условия на внешних границах расчетной области и условия на внутренних свободных и контактных поверхностях.
На границе области интегрирования используется общепринятый способ отражения: вдоль границы вводится дополнительный слой фиктивных ячеек, в узлах которых задаются значения типа узлов, компонент скорости и давления [3].
Так, в случае неотражающей жесткой границы или плоскости симметрии в фиктивных узлах задаются условия антисимметричного отражения нормальной компоненты скорости и симметричного отражения остальных величин (например для нижней границы в направлении х, фиктивные ячейки имеют индексы г = —1;г = 0 — ячейки на границе области и т.д.)
Многопроцессорные ЭВМ и параллельные вычисления
В настоящее время рейтинг самых мощных суперкомпьютеров ТОР500 [115] возглавляет система BlueGene, состоящая из 32768 процессоров PowerPC 440 0.7 GHz, с общей производительностью 70720 GFlops. В России, по данным отечественной версии рейтинга супер-компьютеров ТОР50 [116], самым мощным является кластер МВС 1000М в Межведомственном суперкомпьютерном центре РАН, состоящий из 276 двухпроцессорных узлов на базе процессоров PowerPC 970 2.2 GHz с памятью 4 GB RAM на узел, объединенных сетями Myrinet и 2xGigabit Ethernet [117].
Следует сказать, что наибольшее распространение получили именно кластерная архитектура с распределенной памятью. К концу 1990 года кластерная схема проектирования суперкомпьютеров была слабо распространена, системы такого типа чаще всего собирались в академических организациях, сами при этом являясь целью и объектом исследований. Об использовании их для тяжелых расчетов, запуске на них практических приложений в то время не было и речи. Большинство таких кластеров было с небольшим числом процессоров, и как результат, в ноябрьскую версии списка ТОР500 1999 года вошло всего 7 кластерных систем. Но уже в ноябрьском списке ТОР500 2004 года кластеры занимают доминирующее положение: 294 позиции из 500 (рисунок 2.10). Основные достоинства кластерных систем — сравнительно низкая стоимость, хорошая масштабируемость. Сфера их применения постоянно растет и уже сейчас они находят самое широкое применение не только в науке и ВПК, но в финансовой и производственной областях.
На протяжении уже более чем 20-ти лет параллельное программирование является одним из главных направлений исследований в области вычислительной техники. За это время было разработано множество подходов и языковых средств для программирования параллельных вычислительных систем [118]. Тем не менее, до сих пор прикладные программисты тяготеют к использованию низкоуровневых расширений традиционных языков, таких как HPF («High Performance FORTRAN») [119], и библиотекам с обменом сообщениями: PVM («Parallel Virtual Machine») [120], BSP («Bulk Synchronous Parallelize») [121] и MPI («Message Passing Interface») [122].
MPI или интерфейс передачи сообщений — библиотека реализующая идеологию обмена сообщениями [122], являющаяся стандартом «de facto» среди средств обеспечения межпроцессорных коммуникаций в области параллельных вычислений. Библиотека разработана и постоянно развивается участниками международного сообщества MPIF [123] («Message Passing Interface Forum») в которое входит более 40 организаций. Первый стандарт был опубликован в июне 1994 года. К настоящему моменту времени в свет выпущено расширение стандарта MPI-2. Доступны реализации стандарта MPI (как коммерческие, так и свободно-распространяемые), практически для всех известных в настоящее время платформ, что и определяет его высокую популярность среди разработчиков параллельного программного обеспечения.
Библиотека MPI используется совместно с языками последовательного программирования, такими как FORTRAN-77, С, C++ и поддерживает все базовые типы данных принятые в них, также как и позволяет пользователю определять свои собственные типы данных. Интерфейс содержит порядка 120 процедур для организации и обеспечения: Индивидуальных взаимодействий; Коллективных взаимодействий; Задания пользовательских топологий процессов и др. На этапе выполнения MPI-программа представляет собой совокупность автономных процессов, функционирующих под управлением своих собственных программ и взаимодействующих посредством стандартного набора библиотечных процедур для передачи и приема сообщений. Таким образом, в самом общем случае MPI-программа реализует MPMD-модель программирования («Multiple Program - Multiple Data» — множество потоков управления при множестве потоков данных). На практике, однако, чаще всего ограничиваются SPMD-моделью («Single Program - Multiple Data» — один поток управления при множестве потоков данных). В данной модели все процессы исполняют в общем случае различные ветви одной и той же программы. Такой подход обусловлен тем обстоятельством, что задача может быть достаточно естественным образом разбита на подзадания, каждое из которых может решаться отдельно и независимо от других, но одинаковым со всеми образом. Отметим, что с теоретической точки зрения любое множество MPMD-программ может быть объединено в одну SPMD-программу. С практической же точки зрения SPMD-подход оказывается более предпочтительным при параллельной реализации уже имеющихся последовательных программ, так как не требует значительной переделки кода.
Численное моделирование воздействия пучка ионов урана на тонкую свинцовую фольгу
Процесс высокоскоростного удара является одной из наиболее сложных задач для многомерного газодинамического моделирования. В ходе ее решения возникает ряд характерных особенностей, которые выдвигают высокие требования к используемым численным схемам. Сложная динамично развивающаяся структура течения, наличие ударных волн высокой интенсивности, внутренних контактных и свободных границ, положение которых заранее не известно и также подлежит определению, высокие градиенты давления, плотности и других параметров и большие деформации среды делают, зачастую, использование таких расчетных схем как, например, метод конечных элементов, проблематичным. Практическая важность теоретических исследований в этой области, в частности, обусловлена необходимостью разработки и оценки эффективности конструкционных средств противометеоритной защиты космических летательных аппаратов [4]. Так как интересный для исследования диапазон скоростей соударения практически не реализуем в натурных экспериментах, математическое моделирование остается единственно возможным методом изучения явления.
В данной главе приводятся результаты трехмерных расчетов процесса высокоскоростного пробивания свинцовым шариком свинцовой преграды (со скоростью 6.6 км/с). Данная задача представляет интерес по многим причинам. На момент времени 30 мкс методом теневой рентгеновской съемки зафиксировано пространственное распределение плотности вещества ударника и преграды [126], что дает возможность провести непосредственное сравнение результатов расчета и эксперимента. Высокая скорость удара приводит к плавлению ударника и преграды с последующим испарением при расширении. Следует также отметить, что при двумерном численном моделировании данной задачи [127] был получен парадоксальный результат: наилучшее описание распределения экспериментальной плотности достигалось с использованием наиболее простого по физической модели уравнения состояния свинца.
В двух- и трехмерной постановках моделируется взаимодействие свинцового сферического ударника (диаметр 1.5 см, масса 20 г) со свинцовой пластиной (толщина 0.63 см) при нормальной скорости соударения 6.6 км/с. Расчет проводился с пространственным разрешением 2 ячейки на 1 мм вплоть до момента времени 30 мкс. Схема эксперимента показана на рисунке 3.1.
В расчетах использовались широкодиапазонное УРС типа Ми-Грюнайзена в калорической форме [128], не содержащее информации о положении фазовых границ плавления и кипения и описывающее в обобщенной форме свойства метастабильной жидкости в области растяжения, и полуэмпирическое многофазное УРС в табличной форме [129]. Качество термодинамического описания с помощью калорического и многофазного УРС в области параметров, соответствующих процессам высокоскоростного удара и последующего расширения, иллюстрируется рисунком 3.2. Как следует из рисунка 3.2, данные уравнения состояния с хорошей точностью описывают имеющуюся совокупность экспериментальных данных по ударному сжатию и изэнтропическому расширению свинца. Учет эффекта испарения в многофазном уравнении состояния для области параметров, соответствующих эксперименту [126] приводит к более высокой конечной скорости расширения (до 15%) по сравнению с калорическим У PC.
На рисунках 3.3 показаны распределения плотности в осевом сечении на разные моменты времени для вариантов с калорическим широкодиапазонным и табличным многофазным уравнение состояния свинца.
Расчетное распределение плотности в сечении вдоль оси симметрии для варианта с многофазным УРС на момент времени 30 мкс подробно показано на рисунке 3.4 (различными тонами показаны плотности в диапазоне 0—1 г/см3).
