1. Введение

В настоящей работе моделируется воздействие ударной волны и продуктов сгорания ВВ на преграду, выполненную из жидкого свинца. Параметры взрыва подбирались таким образом, чтобы смоделировать ударно-волновую нагрузку на стенку камеры предполагаемого термоядерного реактора в концепции тяжелоионного синтеза (Fast Ignition Heavy Ion Fusion). Эта концепция была сформулирована в работах [1-4] и развивалась далее в [5-7].

Кроме тестирования расчетных параметров, характеризующих живучесть камеры при высокоинтенсивном импульсном воздействии, целью данной работы является описание вещества стенки с помощью широкодиапазонного уравнения состояния вещества на основе модели А.Б. Медведева [8-9]. Эта модель обобщает классическое уравнение состояния Ван-дер-Ваальса.

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

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

2. Исходные данные

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

Моделирование плоских УВ в веществе при воздействии РИ проводилось в [10] взрывом зарядов ТГ50/50 и ударом пластин из алюминия. Характерные параметры моделирующих устройств и нагрузки при воздействии по алюминию приведены в таблице 1, где: W – скорость ударника; ∆ - толщина ударника или заряда ВВ; I, E - удельные импульс и энергия удара взрыва;  τ = 2 ∆/a, P = I/τ – соответственно длительность и среднее давление при ударе; a - скорость звука в алюминии. Измерения проводились на относительных расстояниях:  при ударе -  0 <X/∆ <200; при взрыве - 0<X/∆< 10. Теплота взрыва Q = 4,87 МДж/кг, скорость детонации D = 7,65 км/с, скорость потока вещества за фронтом ударной волны U = 1,99 км/с, давление детонации P = 26,6 ГПа.

Таблица 1. Параметры моделирующих устройств

и создаваемых ими нагрузок

Таким образом, в тестовых расчетах можно считать, что в преграде формируется давление 20-40 ГПа. Длительность воздействия варьируется в широких пределах. Исходя из тех соображений, что длительность рентгеновского излучения от мишени весьма мала и не превосходит 1 мкс, в данном расчете она положена равной 500 мкс.

Теплофизические параметры свинца приведены в таблице 2.

Таблица 2. Теплофизические параметры свинца.

3. Математическая постановка задачи

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

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

Запишем постановку задачи сначала в эйлеровой форме (см. [11]). Обозначения:  – плотность,  – скорость,  – плотность внутренней энергии,  – давление,  – плотность энтальпии. Записываем уравнения законов сохранения массы, импульса и энергии в виде:

                          (3.1)

Начальные условия:

                                               (3.2)

Поставленная задача (3.1)-(3.2) решалась численно в лагранжевых координатах, следуя [12]. Кратко опишем ниже этот подход.

Вводится массовая координата

                                               ,                                          (3.3)

где n – размерность пространства. В данном случае рассматривается воздействие взрыва на плоскую преграду, так что n = 1. Скорость вещества определяется как . Системе (3.1) тогда отвечают уравнения

                                                 (3.4)

Здесь  – диссипативные члены с вязкостью. Начальные и граничные условия модифицируются аналогично.

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

4. Широкодиапазонное уравнение состояния

         Широкодиапазонное уравнение состояния представляет обобщение модели Ван-дер-Ваальса и имеет вид [9]:

      ГПа,  г/см³.                 (4.1)

Параметр  (давление отталкивания) определяется из трансцендентного уравнения как функция температуры и плотности:

                                                              (4.2)

Здесь  – газовая постоянная, а  – молярная масса свинца. Это уравнение имеет единственное решение относительно  при всех действительных значениях параметров . Перейдем к безразмерным величинам температуры, плотности и давления:

                                   .                         (4.3)

Тогда уравнение состояния (4.1) примет вид:

                 .                (4.4)

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

                                               .

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

                                          .                                   (4.5)

Равновесная кривая  находится из условия равенства температур, давлений и потенциалов Гиббса паровой () и жидкой () фаз:

                                                                  (4.6)

   .                 (4.7)

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

                          (4.8)

Используя систему (4.8), организуем итерационный процесс:

                                   (4.9)

Приращения находятся по формулам (температуру в аргументах опускаем):

                              (4.10) .         

За исходную точку итерационного процесса берем критические значения .

5. Описание расчетного гидродинамического кода

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

                                                                         (5.1)

Система (5.1) аппроксимировалась на двух взаимосвязанных сетках: на основной сетке  с шагом (вообще говоря, неравномерным)  (т.е. , m=1,…, M-1),  и промежуточной (, m=0,…, M), причем основная сетка служит для аппроксимации величин , а промежуточная – для .

Следуя [12], для аппроксимации системы (5.1) была использована двухслойная (с шагом t), чисто неявная, полностью консервативная схема. Избегая громоздкости формул, в величинах на промежуточной сетке будем использовать целочисленные индексы, так что, например, g_m означает g_m+1/2 . Вводим оператор

                                                             (5.2)

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

Система разностных уравнений (5.3) задает общую схему численного решения системы. Общим методом ее решения служит метод раздельных прогонок [12], когда система (5.3) делится на две группы: «динамическую» и «тепловую». Каждая из них решается отдельно некоторым итерационным методом (причем величины из другой группы считаются замороженными), с последующими итерациями между группами.

Рассмотрим, например, «динамическую» группу для k=0. Величина , определяющая псевдовязкое давление, имеет вид

      ,             (5.4)

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

         (5.5)

решается с помощью прогонки. Тогда находятся значения пары величин на следующем слое:

                                            (5.6)

По известным значениям  давление  находится из уравнения состояния с помощью итерационного метода Ньютона. Если при этом  пересекает кривую равновесия (использовалась табулированная кривая на 10000 точек, полученная по описанной выше методике – см. п. 4), то выбирается значение давления в точке пересечения.

Для решения «тепловой» части запишем формулу для энергии (в приближении однократной ионизации) в виде:

                              ;                         (5.7)

,        ,

                                               .       

Теперь, опуская для краткости индекс m в правой части, получаем:

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

6. Моделирование РИ механическим воздействием

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

В работе [10] предложены расчетные и экспериментальные методы моделирования воздействия импульсов мощного рентгеновского излучения на жидкометаллические теплоносители и жесткие стенки камер реактора ядерных энергетических установок путем воздействия на них плоской ударной волной или ударом тонких пластин, метаемых продуктами взрыва химических взрывчатых веществ. Там же приведены результаты экспериментальных и расчетных исследований поведения ударных волн при взрывном и ударном нагружении различных материалов. С использованием разработанных методик проведены исследования поведения жидкой теплозащитной пленки из Li₁₇Pb₈₃ на стенке взрывной камеры, когда вся энергия облучения идет на испарение и возникает максимальный импульс отдачи.

Моделирование воздействия РИ проводилось в [10] с использованием химических зарядов ВВ специальной конструкции и осуществлялось двумя способами:

 - воздействием плоской ударной волны (УВ), генерируемой зарядом заданной толщины; 

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

Подчеркнем, что мягкое РИ поглощается в поверхностном слое облучаемого объекта, имеющего толщину 10^-6 см, за время порядка 10^-15 с, что на много порядков меньше, чем в случае механического воздействия. Воздействие РИ приводит к испарению поверхностного слоя вещества и возникновению УВ под действием импульса отдачи. Повреждения взрывной камеры при таком облучении связано с действием механической нагрузки генерируемой УВ, вызывающей деформацию корпуса, разрушение и расслоение теплозащитного покрытия, а также откольное разрушение внутренней поверхности несущей оболочки.

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

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

В качестве уравнения состояния во всей области течения в [10] использовались ударные адиабаты. Расчет затухания УВ проводился в приближении центрированной простой волны с прямыми характеристиками, исходящими из одной точки на границе мишени [4].

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

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

Численное моделирование откольного эффекта при импульсном воздействии излучений различных типов на жидкометаллические и, возможно, слоистые преграды является центральным пунктом нашей программы. Особенность предлагаемого метода состоит в том, что как таковых откольных напряжений в расчетную модель не вводится. В результате численных расчетов наблюдается движение вещества в рамках гидродинамической модели, в которой, еще раз подчеркнем, используется реальное уравнение состояния вещества. Последнее означает, что в отсутствие внешней среды (воздуха) как составной части модельной системы вакуум между отдельными частями среды (расчетными ячейками) не образуется. Напомним, что моделируется камера реактора, где кроме паров теплоносителя с давлением около 0,01 Па ничего нет. Разрыв (откол) тогда интерпретируется как почти свободное движение двух достаточно плотных ячеек, одна из которых (убегающая) имеет более высокую скорость, чем вторая, и эти две ячейки разделены третьей, плотность которой по сравнению с ними очень мала – на 3-5 порядков меньше.

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

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

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

Рис. 1. Профили давления в жидкой преграде в начале воздействия.

Рис. 2. Профиль давления в жидкой преграде

после отражения УВ от ее внешней границы.

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

Рис. 3. Распределение плотности в преграде через 4,5 мкс после импульса.

Рис. 4. Распределение плотности в преграде через 7 мкс после импульса.

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

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

Сравнивая графики плотности и температуры на Рис. 4 и Рис. 5 видим, что низкой плотности отвечает и низкая температура. Однако состояние разреженного вещества при этом является газообразным.

В «оторвавшихся» каплях, напротив, температура возрастает вследствие роста давления.

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

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

Рис. 6. Распределение давления в преграде через 7 мкс после импульса.

Рис. 7. Распределение скорости среды через 7 мкс после импульса.

Рис. 8. Распределение температуры в среде через 20 мкс после импульса.

Рис. 9. Распределение давления в среде через 20 мкс после импульса.

Рис. 10. Распределение плотности в среде через 20 мкс после импульса.

Рис. 11. Распределение скорости в среде через 20 мкс после импульса.

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

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

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

Благодарности

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

Работа выполнена при частичной финансовой поддержке гранта РФФИ, проект № 07-08-00389, и при частичной поддержке Программы Президиума РАН № 14.

Литература

1. Koshkarev D.G. Charge-Symmetric Driver for Heavy-Ion Fusion. // IL Nuovo Chimento, Vol.106 A, No.11, p.1567-1573. 1993.

2. Koshkarev D.G., Korenev I.L., Yudin L.A. Conceptual Design of Linac for Power HIF Driver. /  CERN 96-05, VI, p.423-426. 1996.

3. Чуразов M.Д., Аксенов A.Г., Забродина E.A. Зажигание термоядерных мишеней пучком тяжелых ионов. // ВАНТ, Сер. Математические модели физических процессов, Вып. 1, №.20. 2001.

4. Medin S.A. et al. Evaluation of a power plant concept for fast ignition heavy ion fusion // Laser and Particle Beams, 2002. V.20. P.419-423.

5. Medin S.A. et al. Reactor Chamber and Balance-of-Plant Characteristics for Fast-Ignition Heavy-Ion Fusion Power Plant // Fusion Science and Technology, 2003. V.43. No.3. P.437-446.

6. Medin S.A., et al. Conceptual Analysis of Energy Conversion in Power Plant for Fast Ignition Heavy Ion Fusion / 30-th EPS Conference on Controlled Fusion and Plasma Physics. Russia, S-Petersburg, July 7-11, 2003.

7. Basko M. M., Churazov M. D. and  Aksenov A. G. Prospects of heavy ion fusion in cylindrical geometry. // Laser and Particle Beams, 2002. V.20.P.411-414.

8. Копышев В.П., Медведев А.Б. Термодинамическая модель сжимаемого коволюма. Саров: РФЯЦ-ВНИИЭФ, 1995.

9. Медведев А.Б. Модификация модели Ван-дер-Ваальса для плотных состояний. / В сб.: Ударные волны и экстремальные состояния вещества. Под ред. В.Е. Фортова, Л.В. Альтшулера, Р.Ф. Трунина и А.И. Фунтикова. М.: Наука, 2000.

10. Соловьев В.О., Христофоров Б.Д. Моделирование механического воздействия импульсного рентгеновского излучения на стенки камер ядерных реакторов. /

11. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика, Т.VI. Гидродинамика. М.: Наука, 1988.

12. Самарский А.А., Попов Ю.П. Разностные методы решения задач газовой динамики. М.: Наука, 1980.

13. Физика взрыва. Ред. Орленко Л.П. Том 2. М. Физматлит. 2004. 656 с.