Содержание к диссертации
Введение
2 Теоретическая часть 13
2.1. Математическая формулировка класса задач 13
2.2. Схема Годунова в подвижных сетках 14
2.2.1. Интегральные уравнения сохранения для подвижных сеток 14
2.2.2. Вывод численной схемы для двумерного случая 16
2.2.3. Решение задачи произвольного разрыва 22
2.3. Алгоритм 30
2.3.1. Общий алгоритм решения задач методом Годунова в подвижных сетках. 30
2.3.2. Анализ точности по пространству. Дробление или огрубление сетки. (Пункт 2 алгоритма.) 32
2.3.3. Анализ формы области. Разбиение-области на подобласти. (Пункт 3 алгоритма.) 36
2.3.4. Поиск пересечений границ подобластей, изменение типов границ.( пункт 5 алгоритма) 38
2.4. Вычисление энерговклада от различных источников энергии 41
2.4.1. Решение уравнений химической кинетики 41
2.4.2. Численная схема решения уравнения теплопроводности 44
2.4.2.1. Анализ устойчивости численной схемы решения уравнения теплопроводности 46
2.4.3. Расчет энерговклада пучков ионов с учетом реального распределения мощности по радиусу пучка и фокусировки в мишени 47
2.5. Тестирование алгоритма 52
2.5.1. Распространение плоской ударной волны и волны разгрузки в идеальном газе, (одномерный тест) 52
2.5.2. Расчет нерегулярного отражения ударных волн.( двумерный тест)...55
2.5.3. Тестирование численной схемы решения уравнений химической кинетики 57
2.5.4. Тестирование численной схемы решения задачи теплопроводности.61
2.5.4.1. Одномерный тест 61
2.5.4.2. Двумерный тест 64
3. Результаты численного моделирования 66
3.1. Высокоскоростное соударение 67
3.1.1. Соударение свинцового шара со свинцовой пластиной 67
3.2. Задачи с химической кинетикой 73
3.2.1. Вычисление критического диаметра ТНТ 73
3.2.2. Двумерный расчет детонации водородовоздушной смеси 77
3.2.3. Численное моделирование разгона стальной цилиндрической оболочки заполненной прессованным ТНТ, при различных способах подрыва заряда 80
3.3. Взаимодействие ионных пучков с конденсированными мишенями 97
3.3.1. Моделирование развития гидродинамической неустойчивости при разгоне металлических фолы пучком протонов 97
3.3.2. Взаимодействие пучков тяжелых ионов с мишенями 102
3.3.2.1. Режимы проникания пучка тяжелых ионов в материал 104
3.3.2.2. Взаимодействие мощного ионного пучка с мишенью из разнесенных свинцовых пластин 108
3.3.2.3. Сравнительные расчеты параметров сжатого вещества для сплошных и полых мишеней 114
3.3.2.4. Оценка оптимальных профиля мощности пучка и размеров полости, обеспечивающих максимальные плотности энергии сжатого вещества 120
Основные результаты: 126
Основные обозначения 127
Литература 129
- Интегральные уравнения сохранения для подвижных сеток
- Численная схема решения уравнения теплопроводности
- Вычисление критического диаметра ТНТ
- Взаимодействие пучков тяжелых ионов с мишенями
Введение к работе
Предметом исследования данной работы являются нестационарные явления и процессы, сопровождаемые достижением в веществе высоких плотностей энергии. В работе рассматриваются такие задачи как, инициирование и развитие детонации конденсированных взрывчатых веществ, динамика многокомпонентных химически реагирующих газовых смесей, взаимодействие пучков легких и тяжелых ионов с конденсированными мишенями. Математическая модель эти процессов и явлений представляет собой систему уравнений в частных производных. Аналитические решения этой системы уравнений найдены лишь для одномерных случаев со специфическими граничными условиями. В общем случае решение этих уравнений требует численного интегрирования.
Основными целями данной работы являлись создание вычислительной методики, способной эффективно решать вышеперечисленные задачи и получение с ее помощью конкретных решений, практически важных задач.
В настоящее время для решения уравнений Эйлера разработано множество численных методов. Эти методы можно классифицировать на методы конечных разностей и методы конечных объёмов. Методы конечных объёмов основаны на численном решении уравнений Эйлера, записанных в интегральном виде, что приводит к автоматическому выполнению законов сохранения массы импульса и полной энергии, т.е. консервативности по этим переменным. Далее здесь будут рассматриваться только методы конечных объемов. Представленная в данной работе методика ориентирована на расчет течений с большими деформациями с выделением контактных разрывов между различными материалами, поэтому обзор существующих численных методик так же ограничен методиками, позволяющими проводить расчет таких течений.
Наибольшее распространение для задач такого рода получили т.н. методы частиц, сочетающие в себе черты эйлерова и лагранжева подходов.
Первый метод частиц в ячейке- РІС, был предложен Харлоу в 1957 году [1]. Метод впоследствии породил множество модификаций. Так, например, модификации метода активно и успешно применялись для расчетов процессов горения и детонации Мейдером [2]. Далее были разработаны методы крупных частиц FLIC [3-4], и метод гладких частиц SPH, предложенный в [5] активно развиваемый в последнее время в США [6-21]. В Новосибирске создан метод индивидуальных частиц [22].
В методе частиц в ячейке область решения покрывается эйлеровой сеткой, но сплошная среда аппроксимируется совокупностью точечных частиц постоянной массы, которые движутся через эйлерову сетку. Частицами через ячейки переносятся кинетическая, внутренняя энергии, масса. На эйлеровой сетке вычисляются параметры полей давления, плотности, внутренней энергии. Метод позволяет исследовать сложные явления в динамике многокомпонентных сред, хорошо отслеживает контактные разрывы, выдерживает произвольные деформации вещества в расчетной области. Однако предложенный Харлоу метод обладает рядом недостатков. Дискретность изменения массы при переходе из одной ячейки в другую ведет к возникновению колебаний в решении, амплитуда которых определяется количеством частиц. Увеличение количества частиц естественным образом щ ведет к пропорциональному увеличению вычислительных затрат. Для устойчивости счета метод требует введения явных диссипативных членов с искусственной вязкостью. Метод плохо описывает области разрежения, из-за пропорционального разряжению уменьшения числа частиц.
Метод крупных частиц вместо совокупности частиц в ячейке использует как частицу саму эйлерову ячейку, чем и определяется название метода.
Лагранжева сущность метода сводится к вычислению потоков через границы эйлеровых ячеек. Уменьшение числа частиц до одной на ячейку ведет соответственно к увеличению эффективности метода по сравнению с PIC . Этот метод также требует введения искусственной вязкости.
Особенностью SPH метода является отсутствие эйлеровой сетки и вообще численной сетки как таковой. В SPH методе присутствуют только частицы с координатами их центров. Частицы в этом методе не являются точечными. Плотность в частице распределена симметричным образом по радиусу частицы, который меняется пропорционально при сжатии и разрежении. Значения поля параметров в заданной точке пространства вычисляются по координатам соседних частиц и функциям распределения параметров в них. Метод хорошо описывает поведение материала при сильных деформациях, от сильного сжатия до сильного разрежения. К недостаткам метода можно отнести зависимость результатов расчета от выбора функции распределения по радиусу частицы, правила определения соседства частиц. Как следствие существует необходимость «правильного» выбора вида функции для конкретного типа задач.
Метод индивидуальных частиц базируется на РІС методе и отличается от него тем, что частицы не точечные. Для трехмерного случая частицы имеют форму параллелепипеда с гранями параллельными плоскостям эйлеровой сетки. Метод включает в также себя процедуры объединения и дробления ячеек, контролирующие наличие определенного числа частиц в одной ячейке сетки. Как и РІС метод требует применения искусственной вязкости.
В отличие от методов частиц, другие методы конечных объемов, используют только ячейки численной сетки и связанную с ними информацию для определения численного решения. Поэтому для остальных методов очень важно создать численную сетку таким образом, что бы было возможно получить решение с заданной точностью, и минимизировать при этом затраты машинного времени.
В работе Дж. Ф.Томпсона [24] приводится обзор по методам расчёта сеток в вычислительной аэродинамике. Рассматриваются различные методы построения сеток, их преимущества и недостатки. Кроме того, автор пишет, какими должны быть методы расчета сеток (МРС). Отмечается, что по существу, расчёт сетки представляет собой расстановку узлов сетки в физической области, таким образом, чтобы обеспечивалась удобная связь между узлами сетки и возможность представления всех физических явлений, происходящих в этой области, минимальным числом узлов, в соответствии с требуемой точностью. МРС должны позволять вести расчёты физических явлений без ограничения на формы областей, что даст возможность разработать общие алгоритмы расчета, в которых границы областей задаются в качестве входных данных. Границы области могут быть подвижными, причем либо их движение может задаваться извне, либо определяться в процессе решения задачи как следствие его эволюции. Узлы сетки также могут изменять свои положения, чтобы следить за развитием явлений по положению градиентов функций (искомых решений физической задачи). Таким образом, МРС должны обеспечить возможность: Упорядоченной расстановки узлов (с тем, чтобы соседние узлы сетки можно было легко находить, а данные в этих узлах эффективно обрабатывать и запоминать при помощи ЭВМ; + Построение сеток, смещение узлов которых не нарушает гладкости сеточных линий; 4 Аппроксимации непрерывных функций их значениями в узлах сетки с требуемой степенью точности и оценки погрешности этой аппроксимации; Ф Определение расстановки узлов сетки по оценке погрешности аппроксимации и изменения способа расстановки узлов.
В работе [24] также обсуждаются принципы построения сеток, конфигурации расчетных областей (где рассматриваются также методы разделения области на части); проводится анализ погрешности аппроксимации разностных выражений, анализ источников погрешностей (в том числе погрешностей привносимых сеткой); рассматриваются МРС с помощью конформных отображений, эллиптических, параболических и гиперболических систем уравнений; алгоритмы построения пространственных, ортогональных, адаптивных сеток, использование подвижных конечных элементов.
В заключение Дж. Ф. Томпсон пишет, что наилучшим МРС для произвольных областей является метод разделения области на части, согласно которому для каждой подобласти рассчитывается своя сетка. Одним из наиболее перспективных направлений он считает разработку методов расчета адаптивных сеток, динамически связанных с решением. В связи с этим он указал также, что применение практически любого вычислительного алгоритма не приводит к затруднениям, если расчетная сетка выбрана правильно, в соответствии с градиентами функции, описывающей искомое решение. Основной вывод, который автор делает в статье, заключается в том, что наибольшего успеха при численном решении уравнений с частными производными можно достичь не тщательной разработкой способов разностной аппроксимации или алгоритмов решения сеточных уравнений, а разработкой адаптивных сеток, динамически связанных с решением.
Этот вывод нашел подтверждение в [25], где описывается адаптация сеток для задач газовой динамики. В [25] проведено сравнение численных расчетов распространения пламени от раскаленного шара в обтекающую его горючую смесь газов проведенных одним и тем же методом на равномерной и адаптивной сетках. Было показано, что в случае адаптивной сетки устранялись численные колебания температуры в областях с большими градиентами, характерными для данного метода при его использовании на равномерной сетке.
Основные выводы авторов о применении адаптивных сеток для данного рода задач, следующие: + адаптация сетки с учетом первой производной от зависимой переменной требует задавать априори максимально допустимое изменение между узлами сетки; + резкие изменения шага, возникающие при адаптации сетки, не приводят, по-видимому, к росту погрешности аппроксимации; + адаптация сетки с учетом первой производной зависимой переменной позволяет решить известную проблему сеточного числа Рейнольдса, возникающую при расчете конвективно-диффузионных зон больших градиентов. В [26] показано, что при использовании адаптивных сеток для решения нестационарных уравнений в частных производных могут возникать неустойчивые решения для движения ячеек сетки, когда система уравнений диссипативна. Используя метод линейного возмущения, авторы нашли простой критерий определения устойчивости системы и показали, как создавать устойчивые дифференциальные системы определения скоростей ячеек адаптивной сетки.
Разновидностью адаптивных методов, позволяющих значительно повысить точность расчетов течений с бесконечно большими градиентами в решении (решения с разрывами) являются методы, позволяющие явно выделять поверхность разрыва в решении. Один из таких методов приведен в [27]. Здесь решались двумерные задачи, для которых поверхность разрыва представлялась в виде ломаной линии. Все поле течения вычислялось на равномерной эйлеровой сетке по методу Мак-Кормака, за исключением ячеек, которые пересекала ломаная. Для этих ячеек пересчет параметров производился с учетом величин по обе стороны разрыва. Движение самого фронта и параметров на нем вычислялось из решения задачи Римана со вторым порядком точности. Этим методом были решены тестовые задачи по распространению сферически симметричной ударной волны и проведены сравнения с аналитическими решениями. Было получено очень хорошее соответствие решений. Далее приведены результаты расчетов положения выделяемого ударного фронта, образующегося при обтекании клина сверхзвуковым потоком и параметры течения, а также решена задача о возникновении неустойчивости Кельвина-Гельмгольца на границе раздела движущихся с различными скоростями жидкостей разной плотности. Расчеты проводились как на грубой, так и на мелкой сетке. Было показано, что даже грубая сетка позволяет получать приемлемые результаты, благодаря выделению поверхностей разрыва.
Выделение поверхностей разрыва удобно проводить на подвижных сетках, кода поверхность разрыва является одновременно поверхностью уровня разностной сетки.
Примером такого метода может служить разработанный в России метод Годунова в подвижных сетках с выделением поверхностей разрыва. В [28] представлена реализация метода Годунова со вторым порядком точности по пространству. Проведено сравнение метода Годунова и обобщенного метода Годунова второго порядка при расчете двумерных стационарных течений. Показано, что метод второго порядка обеспечивает существенное повышение точности численных решений (вследствие снижения схемной вязкости) и в то же время столь же надежен, как и оригинальный метод Годунова, не использующий искусственной вязкости.
В то же время использование подвижных регулярных сеток в некоторых случаях может приводить к сильным искажениям сетки, которые в свою очередь будут приводить к сильным искажениям искомого решения. Альтернативой в этом случае может быть применение нерегулярных сеток, которые к тому же обладают высокой степенью адаптивности. В [29] приведен пример использования сетки такого рода. Сетка построена на треугольных ячейках. На сетке реализован метод конечных объемов. Параметры в ячейке сетки вычисляются с использованием потоков через окаймляющие ее ребра, поэтому, чтобы легко пересчитывать параметры внутри ячейки и в тоже время иметь возможность добавлять или исключать вершины и ребра треугольников, не изменяя записей, относящихся к другим ячейкам, информация о параметрах в ячейке хранится в списках, привязанных к ребрам. Разработанный алгоритм был применен к расчету равновесия плазмы со свободными границами, где была продемонстрирована его способность выделять свободные границы и быстро приспосабливаться к решению. Подобные программы использовались в [30, 31] с большим успехом для решения уравнений Эйлера.
Предлагаемая в данной работе методика использует в своей основе метод Годунова [23] для четырехугольных криволинейных подвижных сеток.
Методика включает в себя возможность разреза в процессе расчета каждой численной области, покрытой одной сеткой, вдоль одной из линий сетки на две сетки, в каждой из которых сетка строится независимо. Это делает методику более адаптируемой к изменяющейся форме области счета, приближая, по степени адаптации к форме области, данный метод к методам на ^ нерегулярных сетках. По сути метод использует преимущества регулярных сеток для расчета течений внутри каждой из сеток и в тоже время позволяет вводить нерегулярность в тех местах, где она необходима, через разрезы сеток.
Для получения более ортогональных сеток имеется алгоритм автоматического изменения расстановки узлов по контуру сетки. v Внутри каждой из сеток имеется возможность адаптации сетки к решению путем сгущения узлов сетки в области больших градиентов.
Пространственное разрешение в каждой из сеток поддерживается автоматически на заданном уровне, путем добавления или удаления узлов сетки независимо по каждому из направлений сетки. ^ Разработанная методика реализована в виде модулей, написанных на языке программирования Фортран-77.
Разработанная методика позволяет проводить численное моделирование задач в двумерной постановке для плоских и осесимметричных течений. Задачи решаются в гидродинамическом приближении для невязких сжимаемых сред, ф движение которых описывается уравнениями Эйлера, модифицированных членами, учитывающими вклад дополнительных источников энергии.
Во второй главе дана общая математическая формулировка моделируемых процессов, описываемых системой дифференциальных уравнений в частных производных. Подробно описан метод Годунова в подвижных криволинейных сетках решения уравнений гидродинамики. Приведены алгоритмы адаптации сеток и изменения типов граничных условий на сетках. Описаны алгоритмы вычисления энерговклада от различных источников энергии. Приведены тестовые расчеты, демонстрирующие адекватность примененных моделей и правильность реализации алгоритмов.
В третьей главе приведены результаты численного моделирования, полученные разработанным кодом. Приведены решения задач высокоскоростного соударения тел, детонации конденсированных взрывчатых веществ, динамики химически реагирующих газовых смесей, взаимодействия пучков легких и тяжелых ионов с конденсированными мишенями.
В тексте нумерация объектов, таких как формулы, рисунки и таблицы своя в каждом из параграфов, в случае ссылки на объект из другого параграфа перед номером ссылаемого объекта добавляются через точку номера соответствующих параграфов.
Интегральные уравнения сохранения для подвижных сеток
Для вывода интегральных уравнений необходимых для построения схемы Годунова уравнения (2.1.1) удобно переписать в виде: гиперплоскостями t=const т.е. для каждой из компонент вектора F мы имеем уравнение здесь подынтегральное выражение есть дивергенция четырехмерного вектора g: g = (F,F Поэтому согласно теореме Остроградского-Гауса выражение (3) равносильно следующему: где v - орт внешней нормали к поверхности Г. Если 1- орт оси t и п - орт внешней нормали к сечению a)(t) гиперплоскостью t=const то и, следовательно, соотношение (3) принимает вид: Рассмотрим область Г ограниченную гиперповерхностями r(ti) rfe) , где sin(v,t)=0 и cos(v,t[)= -1, cos(v,t2)= 1 и поверхностью rt, где sin( vfytO и запишем (6) в виде: Тогда для поверхности rt можно ввести в рассмотрение Dn скорость перемещения поверхности Г в направлении нормали п, тогда вектор D„ п+ 1 ортогонален вектору v, следовательно выполняется соотношение: Для двумерного случая уравнения (2.2.1.10) для каждой из компонент вектора F имеют вид: рассмотрим счетную область Q(t). В фиксированный момент времени t = tn это криволинейный четырехугольник с границами ГД Г", Г/, и ГУ (см. рис.1). Ячейку сетки с индексами (ij) в момент времени t„ обозначим R"p. Она определяется координатами четырех углов с координатами определяемыми уравнениями аналогичными (10). Объём области fljj, в трёхмерном пространстве (x,y,t)cR xlRe определяется таким образом двумя последовательностями точек (2) и (3) (см. рис. 2). Уравнений Эйлера (1) для объёма Q"tj запишем в виде где F обозначает вектор переменных F=(p,pu,pv,pe) имеющих скорость и = (и,и ) и Q?m,ij, m= 1»— 4, обозначают потоки по нормали к поверхностям боковых граней четырёх сторон fig. Эти боковые грани fi\j являются поверхностями, соединяющими соответствующие углы Ry и Rjj1; Например левая боковая грань Si определяется точками (см. рис. 3). Остальные стороны справа спереди и сзади обозначены на рис. 3 Sr,Sf, и Sb, соответственно.
Поскольку поверхности граней, в общем случае, неплоские, существует некоторый произвол в их определении Необходимые для определения потоков в уравнениях (1) направление нормали к грани и значение скорости D грани в направлении нормали предложено определять следующим образом. Рассмотрим левую грань S1 ,которая отделяет 0.пи от 0."Ы]. Скорость грани в направлении нормали является функцией пространства и времени. В численной схеме будем использовать некую усредненную скорость: где П(5/) - площадь проекции грани 5/ на плоскость t=t„. Эта площадь есть мера изменения площади ячейки R"j связанная с движением ребра V"mi.m,P"mii.m. Перемещение ребра в направлении нормали получается при делении площади проекции на средний размер ребра, обозначаемый d\. В расчетах средняя длина ребра определяется как расстояние между усреднёнными координатами его концов Если поверхность грани плоская, то скорость D, определяемая (5) совпадает с локальной скоростью движения этой грани в каждой точке. Данный численный метод является непосредственной аппроксимацией щ интегральных уравнений (4) и относится к методам конечных объёмов в подвижных сетках. Он имеет вид где F"j, -усредненное по ячейке R"j приближение решения F, Qi", Q", Qf и Qbn - приближение для суммарных потоков между ячейками сетки, называемые to численными потоками. Запишем формулу для вычисления потока Q" Все параметры в этой формуле используются, как значения на грани движущейся со скоростью D„. Остальные потоки вычисляются аналогично. Для того чтобы их вычислить нужно знать положение внешней нормали к , поверхности. Естественно выбрать для грани Si, нормаль определяемую точками (6), которые уже использовались при определении скорости D поверхности Si. Эта нормаль Л/ и соответствующее касательное направление заданы формулами
Численная схема решения уравнения теплопроводности
Воспользуемся уравнением теплопроводности Фурье : где с-теплоемость, р-плотность, Т-температура, dV- элемент объёма, к-коэффициент теплопроводности, n-направление нормали к поверхности, dS-элемент площади поверхности. Запишем это уравнение для произвольной ячейки сетки с индексом і: где J-количество граней с площадью Sy отделяющих і-ю ячейку с объёмом Vi от соседних ячеек с индексами j. Сконструируем схему по принципу "предиктор-корректор". Заменим производные конечными разностями: Л А где индекс отмечает значение для верхнего слоя по времени, ха-расстояние л между центрами соседних ячеек, Т-средняя за шаг по времени температура, ys коэффициент учитывающий поправку на разницу в направлении между нормалью к поверхности грани и вектором соединяющим центры ячеек, At-шаг , вводя обозначение Эта формула используется для вычисления новых значений температур по неявной схеме. На шаге "предиктор" для правой части уравнения вместо температур на верхнем временном слое используются температуры с нижнего слоя. На стадии "корректор" используется эта же формула, где в правой части используются значения для верхнего слоя с предыдущей итерации. Итерации продолжаются до удовлетворения условия сходимости: Полученная схема является консервативной и имеет второй порядок точности по пространству и времени: где Т(х,і)-численное решение, Т (хД)-точное решение, dx2,dy2-mar по пространству и по времени. Для анализа устойчивости полученной численной схемы используется метод Фурье. Для упрощения анализ проводится на равномерной прямоугольной сетке с постоянным шагом по каждому направлению и постоянными параметрами к,с,р. В этом случае Ау одинаково для всех ячеек вдоль каждого из направлений, а именно: Пусть для ячейки с индексом к: где і - мнимая единица, j обозначает направление (j=1,2,3 для 3-мерного случая) . Используя уравнение (5) запишем:
То есть решения, получаемые по данной схеме всегда устойчивы. Но при X О в решении могут возникать затухающие колебания. Для получения решений без колебаний шаг по времени должен удовлетворять условию: Численное определение энерговклада пучка в сечении мишени вдоль оси симметрии построено несколькими способами и включает в себя вычисление энерговклада пучка с детальной трассировкой траекторий ионов в мишени, многогрупповой подход с построением оптимальных групп и упрощенный подход, заключающийся фактически в сведении многогруппового подхода к одной группе. Вычисление энерговклада пучка с детальной трассировкой траекторий ионов основано на предположении, что ионы движутся по прямолинейным траекториям. Пучок представляется набором частиц с координатами и скоростями, заданными в одном из сечений пучка. Запишем сначала уравнения движения отдельного иона в трехмерном пространстве в декартовой системе координат {x,y,z}. Пусть ионы движутся по прямолинейным траекториям в направлении оси х и в сечении с координатой х0, заданы все координаты ионов {yQ,z0}u производные от траекторий ионов вдоль оси х {y 0,z0}, тогда траектории ионов описываются функциями от х: Те же уравнения в цилиндрической системе координат {х,г,ф}: Поскольку для интересующего нас осесимметричного случая поворот на любой угол р не должен иметь значения, то используя для начального положения иона систему декартовых координат с осью ъ\, направленной вдоль радиус вектора г, перепишем первое из уравнений (2) относительно повернутой вокруг оси х на угол р0 системы координат {x,yi,zi}, сократив тем самым на единицу количество параметров, определяющих траекторию движения каждого иона: Энерговклад вычисляется последовательно для каждой из частиц, для данного шага по времени. Пусть физическая область, в которой следует вычислить энерговклад от пучка ионов, покрыта набором криволинейных сеток {Л .УД , используемых для гидродинамического расчета. Сетка предполагается неподвижной в процессе вычисления энерговклада. Для каждой частицы и каждой сетки определяются координаты пересечения траектории частицы с границей сетки. Затем эти координаты упорядочиваются по координате х и потери энергии интегрируются от ячейки к ячейке, при этом вычисляется изменение энергии иона и удельной внутренней энергии вещества в ячейке по формулам: 7VR где ДЕ-изменение энергии иона, —- тормозная способность материала, Еп+і,Еп ді энергия иона на выходе из ячейки, энергия иона на входе в ячейку, еп+1,е„-пересчитанная внутренняя энергия вещества в ячейке, исходная внутренняя энергия вещества в ячейке, .-суммарная энергия ионов пучка пересекающих сечение с х=Х(ь за данный шаг по времени, Р-доля энергии приписанная иону с заданной траекторией, V -объем ячейки, Д/ -длина траектории иона при пробеге через ячейку, вычисляемая по формуле: Поскольку реальное количество ионов в пучке велико, то для нужд численного моделирования желательно описать пучок таким числом частиц и такими траекториями с весами, чтобы получить в переделах заданной точности правильное пространственное энерговыделение в мишени с возможно меньшими затратами машинного времени.
Так называемый многогрупповой подход, описываемый далее, предполагает разбиение всех траектории пучка на П дискретные группы: с К разбиениями по радиусу пучка, при x = x0f L(k) разбиениями по .у, и М(к) разбиениями nozj для каждого из разбиений по радиусу пучка. Рассмотрим далее так называемые «хорошо сфокусированные пучки», т.е. пучки в которых ионы достигают минимального расстояния до оси пучка г-( , )= , Vі г ПРИ г- - - г U) в окрестности точки с координатой xf- фокусе. В этом случае возможно провести такое разбиение на группы, для которых траектории ионов из разных групп по радиусу можно объединить в группы с подобными траекториями, для которых справедливы уравнения: На рис.1 показана схематическая иллюстрация представления пучка в виде групп с подобными траекториями. Представленные на рис.1 группы имеют один и тот же фокус у всех траекторий. Группы отличаются между собой величиной y"j и соответственно связанным с ней по второму из уравнений (7) z w Т.е. любая траектория в группе подобных траекторий, используя выражение (8), опуская несущественные индексы, может быть записана в виде: Обратное преобразование /" ( ) переводит эти траектории в прямые линии У =rk(x)f l(x) = Rk. Используя f \x) для каждой криволинейной сетки {ху,Ге\я, имеем: Энерговклад в мишень вычисляется на прямолинейной равномерной сетке {хі У к} образованную пересечением семейств прямых х = х,,у = ук. Для вычисления энерговклада необходимо также провести переинтерполяцию на сетку {х,,ук} параметров с криволинейных сеток: где Р и р - параметры мишени необходимые для вычисления энерговклада на криволинейной и прямолинейной сетке, Z - отображение криволинейной сетки на прямолинейную, которое легко определяется в силу прямолинейности последней. Расчет потери энергии иона вдоль каждой траектории производится по тем же формулам (4). Затем вычисленные значения энергии выделенной в ячейках сетки {х,,ук} добавляются к соответствующим на криволинейной сетке На практике реальное распределение ионов по траекториям, как правило, не известно. Из экспериментов по свечению газа определяют лишь огибающую траекторий ионов пучка. Ее можно использовать в качестве функции f(x) из уравнения (8). Экспериментально определяют так же распределение мощности по радиусу пучка, как правило, это функция Гаусса, обрезанная на некотором радиусе. Используя эту информацию можно построить одну группу с подобными траекториями, которую и использовать в дальнейших расчетах.
Вычисление критического диаметра ТНТ
Моделирование и предсказание критических условий возникновения и распространения детонации является важной проблемой связанной с использованием ВВ и безопасностью взрывоопасных объектов. Прямое экспериментальное исследование таких сложных проблем дорого и небезопасно. Одной из таких проблем является проблема нахождения критического диаметра заряда ВВ - минимального диаметра заряда ВВ, при котором еще может существовать устойчивая детонация, распространяющаяся вдоль заряда. Эта задача была поставлена и исследована качественно в [42, 43] далее в работах [44, 45, 46]. Численное моделирование и анализ детонации для зарядов конечного диаметра сделан в работах [47]. В хорошем согласии с экспериментом находится значение критического диаметра ТНТ, полученное в работе [32]. Наличие критического диаметра обусловлено мощной боковой разгрузкой для продуктов взрыва за фронтом волны детонации распространяющейся по цилиндрическому заряду. При использование стратегии подвижных сеток для данной задачи удобно использовать тот факт, что область влияния разлетающихся продуктов взрыва на скорость распространения фронта волны детонации ограничена так называемой звуковой линией, определяемой условием uf+c = D, где uf- скорость вещества в направлении распространения фронта волны, с-местная скорость звука, D-скорость фронта волны детонации. В проведенных расчетах счетная область ограничена фронтом волны детонации, осью симметрии заряда, свободной поверхностью продуктов взрыва, прямой линией перпендикулярной оси симметрии заряда, передвигающейся за фронтом волны детонации с такой скоростью, чтобы вдоль нее выполнялось условие где рх и єх- холодные составляющие давления и удельной внутренней энергии, константы В0=7.558 ГПа, n=3.575, m=2.5, V0=l/1.62 см3/г, Г(а)- коэффициент Грюнайзена, а- массовая для ВВ ( для непрореагировавшего ВВ а=1 и а=0 для продуктов взрыва). Кинетика превращения ТНТ описывается согласно [32]: Где Ef - прирост внутренней энергии во фронте ударной волны, к и у -постоянные материала к= І с Дж Па"1,7==0.3. Величина Ef изменяется скачком для каждой частицы прошедшей через фронт ударной волны и в дальнейшем остается постоянной величиной. где Pf и (V-Vf)- давление за фронтом ударной волны и изменение объема частицы при прохождении через ударную волну. Явное выделение фронта детонации позволяет значительно улучшить качество расчетов.
Для определения критического диаметра ТНТ рассчитывался заряд с диаметром, уменьшающимся от времени. фронтом волны детонации от времени. Радиус заряда уменьшается экспоненциально до 0.7 мм, когда стационарная детонация более не поддерживается. На рис. 2 сверху от оси симметрии изображены численные сетки на три различных момента времени. Снизу от оси Аналогично для тех же моментов времени на рис. 3 показаны изолинии массовой доли ТНТ в сверху и изолинии давления снизу от оси симметрии. Задача рассматривалась в плоской постановке для замкнутой емкости размером 1м х 0.5м. Емкость заполнена стехиометрической водородовоздушной смесью при давлении 1 атм. В начальный момент времени внутри области возник очаг с температурой ЗОООК и давлением 10 бар. Надо найти поля давлений, скоростей, концентраций и других параметров, возникающих при прохождении детонационной волны по газовой смеси и при взаимодействии со стенками емкости. Задача решалась для полной системы уравнений химических реакций такой же, как и для тестовой задачи описанной в пункте 2.5.3. Задача решалась двумя способами - с выделением фронта детонации подвижной сеткой (Д-метод), и на сетке, границы которой совпадают со стенками емкости (Ф-метод). На рис.1 показаны численные сетки в верхней половине каждого рисунка и изолинии давления в нижней. Верхние рисунки соответствуют начальному моменту времени, нижние - одному из более поздних моментов времени. Слева и справа показаны варианты насчитанные Д-методом, в середине Ф-методом.
Справа разрешение, т.е. средний размер ячейки в два раза больше, чем слева. Интересно, что сравнение результатов расчетов этих двух вариантов позволяет заметить значительное отклонение формы фронта от сферичности для варианта с более грубым разрешением. Это связано с тем, что размеры ячеек вдоль фронта различны и для грубого разрешения расчет дает заниженные величины для давления и плотности за фронтом детонационной волны, и, как следствие, заниженную скорость распространения фронта. Таким образом, искажение формы фронта для задачи сигнализирует о недостаточной точности расчета, что позволяет подобрать надлежащим образом разрешение сетки. Несмотря на то, что общие виды изолиний давления для Д и Ф методов показанных на рисунке 1. почти одинаковы, существует значительное расхождение в величине максимума давления, непосредственно за фронтом волны детонации.
Взаимодействие пучков тяжелых ионов с мишенями
В отличие от ускорителя KALIF в FZK (Карлсруэ, Германия) где пучки протонов разгоняются до энергии порядка 1 Мэв с пробегами несколько десятков микрон и фокусируются до размера 1 см, на ускорителе SIS в GSI (Дармштадт, Германия) ионные пучки имеют энергию 200-300 Мэв на нуклон с пробегом ионов порядка нескольких миллиметров и даже более сантиметра и размер пучка в фокусе около 0.5-5-1 мм в диаметре. Область энерговклада при взаимодействии со сплошной конденсированной мишенью в первом случае представляет собой тонкий слой микронного размера на поверхности мишени диаметром около сантиметра, что создает большие градиенты давления в направлении перпендикулярном поверхности мишени. Для плоской фольги градиенты вдоль поверхности мишени не велики, относительно градиентов в направлении нормали к поверхности, поэтому возможно проведение одномерных численных экспериментов при моделировании разгона металлических фольг протонным пучком [49]. Во втором случае область энерговклада имеет сравнимые по глубине и диаметру размеры, измеряемые миллиметрами. Соответственно поля градиентов давления и определяемых ими скоростей движения материала мишени имеют, как минимум, двумерную геометрию. На рис.1 показаны зависимости тормозной способности свинца для ионов аргона и урана в зависимости от энергии иона. Эти зависимости имеют ярко выраженный максимум, т.н. Брэгговский пик. Наличие Брэгговскоко пика приводит к возникновению зоны максимального энерговыделения на некотором расстоянии от места проникновения пучка в материал и является важной особенностью поглощения ионов с высокими энергиями. Для сфокусированных пучков тяжелых ионов, когда поперечные размеры пучка становятся сравнимы с длиной пробега отдельного иона в веществе, следует отметить такую особенность пучка при взаимодействии с веществом, как увеличение проникания пучка в материал на глубину превышающую пробег отдельного иона. Это явление связано с тем, что за время действия импульса вещество успевает разлететься из зоны энерговыделения пучка. Естественно, что чем быстрее падает плотность вещества на пути иона, тем дальше ион проникает в мишень.
Генерируемые при проникании пучка в вещество конфигурации ударных волн имеют некоторые аналогии с режимами обтекания тел при движении в среде. Можно выделить дозвуковой и сверхзвуковой режимы проникания пучка в вещество. Однако скорость проникания пучка даже в сверхзвуковом режиме напрямую связана с линейной плотностью вещества, которая определяются историей процесса, в отличие от обтекания тела, где режим обтекания определяется скоростью тела, не зависящей на прямую от процесса. Поскольку процесс проникновения пучка обладает рядом особенностей, здесь следует уточнить классификацию рассматриваемых режимов проникания. Определим дозвуковой режим проникания, как режим, когда зона энерговыделения пучка находится внутри области за формируемой ударной волной. Сверхзвуковым назовем режим, когда зона энерговыделения пучка опережает фронт ударной волны. Предполагается, что мощность пучка достаточна для генерации ударной волны. Различные режимы проникания наблюдались в расчетах взаимодействия пучка Аг+18 со свинцовой мишенью, с параметрами соответствующими условиям эксперимента [50], проведенного в GSI в 2000г. в В таблице 1 приведены параметры ионных пучков только для двух вариантов расчета, отличающихся не сильно, но при взаимодействии этих пучков с мишенью реализуются дозвуковой в первом случае и сверхзвуковой во втором режимы проникания пучка в мишень. В расчетах рассматривался осесимметричный пучок ионов движущихся параллельно, с распределением мощности по радиусу, описываемой функцией Гаусса, использовался экспериментально определенный профиль мощности пучка от времени рис. 1.
Профиль мощности имеет четыре выраженных пика. Мишень состояла из свинцовой пластины толщиной 2.5 мм и следующей за ней пластины из ПММА 5мм толщиной. На границе между пластинами на оси действия пучка был смонтирован датчик давления. Поскольку пробег ионов около 1.2 мм, ожидалось померить давление в ударной волне, генерируемой пучком. Однако сигнал на датчике показал отрицательное давление, чего не могло быть. Поэтому предположили, что сигнал на датчике отражал непосредственное взаимодействие ионного пучка с датчиком. В пользу этого достоверно, начало сигнала в точности соответствовало началу четвертого пика мощности рис.2. Результаты численного моделирования также подтвердили это предположение. На рис.3 показаны изолинии энергии и давления на момент выхода возмущения на границу свинец-ПММА для варианта 2 из таблицы 1. На рис. 36 хорошо виден конус Маха, свидетельствующий о сверхзвуковом проникновении пучка в свинец. Поскольку sin а = — , где а половина угла при вершине конуса Маха, М-число Маха, можно определить что проникание пучка в свинец происходит со скоростью примерно в 2 раза превышающей скорость звука.
