CC BY
11
УДК 539.3
М. Е. Поварницын, А. В. Конюхов, П. Р. Левашов
МОДЕЛЬ МНОГОКОМПОНЕНТНОЙ УПРУГОПЛАСТИЧЕСКОЙ СРЕДЫ ДЛЯ ОПИСАНИЯ ВЫСОКОСКОРОСТНОЙ ДЕФОРМАЦИИ И РАЗРУШЕНИЯ, ВЫЗВАННЫХ ДЕЙСТВИЕМ ВЫСОКИХ ПОТОКОВ ЭНЕРГИИ
Ключевые слова: параллельный компьютерный код, упругопластика, многокомпонентные течения, лазерный нагрев.
Разработана модель многокомпонентной упругопластической среды, позволяющая исследовать в численном эксперименте динамику высокоэнергетического воздействия различной природы на многокомпонентные среды. На основе модели создан параллельный компьютерный код, с помощью которого проведено моделирование отклика вещества мишеней на фемтосекундное лазерное воздействие и высокоскоростной удар. Исследована динамика фазовых переходов и откольного разрушения среды в сильных волнах разгрузки.
Keywords: parallel code, elastoplastics, multi-material flows, laser heating.
We developed a multi-material elastoplastic model for simulation of dynamics of high-energy interaction with multi-component media. We developed a parallel code based on the presented model. Using the code, it has been performed simulation of matter response on ultrashort laser heating as well as high-velocity impact. Dynamics of phase transitions and spallation in intensive rarefaction wave has been investigated.
Введение
Активное изучениединамики высокоэнергетического воздействия на многокомпонентные конденсированные среды обусловлено многочисленными практическими приложениями. Разнообразие методов воздействия, таких как высокоскоростной удар, пучки заряженных и нейтральных частиц, лазерное излучение и др. дает возможность исследовать отклик вещества в широком диапазоне термодинамических параметров. В последнее время был достигнут значительный прогресс в изучении сверхбыстрых процессов, происходящих в веществе при интенсивном фемтосекундном лазерном воздействии, которое приводит к возникновению предельных режимов нагрева и деформацииматериалов в сильных ударных волнах и волнах разгрузки [1, 2]. В данной работе описана модель многокомпонентной упругопластической среды,учитывающая эффект поглощения фемтосекундного лазерного излучения. С помощью разработанного подхода можно моделировать задачи, связанные с ударным взаимодействием объектов на дозвуковых и сверхзвуковых скоростях, а также описывать процессы взаимодействия ультракоротких лазерных импульсов с веществом.
Модель многокомпонентной среды
Будем характеризовать состояние твердого тела скоростью и, энтропией 5 и градиентом деформаций Fj¡=дxjlдx0i, где дхи дх^обозначают пространственные и материальные (лагранжевы) координаты ненапряженного состояния, соответственно [3,4, 5]. Пусть в ячейке имеется одновременно несколько веществ, при этом сумма объемных долей fа в ячейке подчиняется очевидному условию:
^ а = 1. (1)
а
Тогда средняя плотность смеси р связана с плотностями отдельных веществ ра следующимсоот-ношением
P=YJ ара •
(2)
Уравнение сохранения массы для каждой компоненты имеет следующий вид
df ара 1 df apauk _ 0
dt
-+-
dxk
(3)
где ик - компоненты вектора скорости, а также подразумевается суммирование по повторяющимся нижним индексам.
Объемные доли не являются консервативными величинами, их изменение будем вычислять, следуя модели [6]
dfa df au k
- +
dt dx,
faKs £4 ка dxk'
(4)
где Ks _l ^ fajKaa I - эффективный модуль все-
Ка а л а ~
_ . _ _ .. , 5 = р X - истинный мо-
дуль сжатия вещества а, X® - квадрат скорости распространения продольной звуковойволны в веществе а. При этом член в правой части отвечает за изменение объемных долей вследствие сжатия или расширения вещества (обратно пропорционально их модулям всестороннего сжатия), что необходимо
для выполнения условия а = 1. Суммируя урав-
а
нение (1) по всем компонентам, получим тождество
дик =дик дхк дхк '
Уравнение сохранения импульса имеет следующий вид в «одножидкостном» (скорость компонент внутри смешанной ячейки одинакова) приближении:
dM + d(PuUk -yik) _ 0
dt
dxk
(5)
Здесь компоненты тензора напряжений с «смеси» вычисляются по модели [7]
'j.
а
С
=6 I
^а_а
с ва
Стт = КЭ - (6)
а 6 а КЭ
Эффективный модуль сдвига вычисляется по аналогии с модулем сжатия через модули сдвига компонент:
6 I 1 ва
(7)
6 а а 1 а ~ л ее
= р Л, — истинный модуль сдвига, а Л, — квадрат скорости распространения поперечной волны в веществе а.
Уравнение для полной энергии «смеси» записывается аналогично случаю однокомпонентной среды:
дрЕ д(ри к Е — ис)
д,
- + -
дхк
= Ql
(8)
- - и,и,
здесь Е = е +—^— полная энергия, а внутренняя
энергия смеси находится через внутренние энергии компонент е как
е = — У Уараеа , (9)
Р а
Источник объемного лазерного нагрева вычисляется по модели, описанной в работе [8].Для индивидуальных компонент уравнение внутренней энергии запишется в следующем виде
У араеа ду араеаик)
д,
дх
(
к
- А
са УаКэ + _а Уа6
скк--+ с к-
К« 'к
\ э У
ди, ,а а ' + У ара01
дх„
(10)
А для полной энергии компоненты «уравнение имеет вид:
дуараЕа д(УараЕаик — иС)
д,
дхк
= У ара01. (11)
Уравнения для компонент «усредненного» по ячейке тензора градиентадеформаций Р_ записываются следующим образом (без учета пластичности):
Р , д(р,ик-р1=,ки,)
__ + • ¡Г к У ¡к .. =
д дхк '
"' дхк
(12)
После нахождения средней деформации смеси определяем деформации компонент:
6 К
ра = 6 р Ра = 1 , (р — 1)
!*_ 6« »■ К тт
Это гарантирует, что IУаРа = Р .
а
Релаксация тензора напряжений для индивидуальных компонент проводится для каждого вещества в ячейке по итерационной однокомпонентной схеме для уравнений вида
У арар^а д
2Сага
с
Р р
(14)
Истинная плотность и истинный тензор напряжений с_ связаны между собой условиями для каждой компоненты а:
ра = р« Ра | (15)
деа
_а _а га
с_ = р р
дРа и1 _к
(16)
Для замыкания системы уравнений (1)-(16) используются уравнения состояния, которые задаются зависимостью внутренней энергии от тензора деформации и энтропии. Являясь скалярной функцией, внутренняя энергия зависит от инвариантов тензора деформаций. Симметрия тензора напряжений Коши гарантируется, если энергия выражена через инварианты симметричного тензора, в частности, тензора Фингераg=F-TF-1. При формулировке уравнения состояния внутренняя энергия рассматривается в виде суммы двух вкладов: гидродинамической составляющей е" и энергии деформации сдвига е55. Такое представление позволяет наиболее просто использовать уравнения состояния, разработанные для гидродинамических кодов [9, 10, 11]: е(/„/2,/з,Э) = е"(/з,Э) + е5(/„/2,/3,Э). ЗдесьЭ- энтропия, а инварианты тензора деформаций имеют вид:
/1= ^й), /2= (tr(g)2— tr(g2))2/2, /3= det|g|.
Следует отметить, что /3=(р/р)2 Например, уравнение состояния меди имеет вид [3]:
К
е" (/з' Э) = « («2 — 1)2 + ^2 [ехР (Э / Су) — 1],
В
е5 (/1,/2 ,/3 ,Э) = 32 (/1 / 3 — /2),
где К0, В0, Т0, су, а, у- материальные константы.
Уравнения решаются модифицированным методом типа Годунова [12]совместно с процедуройпо-строения многоуровневых пространственно-временных адаптивных сеток. В основе используемого алгоритма адаптивного измельчения сеток лежит метод Бергера и Олигера [13, 14] для задач, описываемых системой дифференциальных уравнений гиперболического типа.Под адаптивным построением вложенных сеток понимают рекурсивный процесс измельчения ячеек самой грубой сетки до тех пор, пока на уровне с самой мелкой сеткой решение не будет найдено с заданной точностью, причем измельчение происходит только в тех областях, где это необходимо (большие градиенты, разрывы, контактные границы). При моделировании нестационарных задач на сетках разных уровней дискретизации периодически оценивается ошибка вычислений, и если критерий ошибки в некоторых ячейках перестает удовлетворяться, то эти ячейки помечаются как «требующие измельчения» и добавляются в сетки более высокого уровня измельчения.
Результаты моделирования и их обсуждение
Разработанная модель и компьютерная программа использованы для описания воздействия интен-
1
сивных потоков энергии на многокомпонентные мишени. В качестве первого примера рассмотрим высокоскоростной удар стального ударника по двуслойной преграде из углерода и алюминия. Каждое из веществ описывается своим уравнением состояния. Диаметр сферического железного ударника -4мм,первая пластина 14 мм мишени сделана из углерода, вторая пластина 3 мм из алюминия. Скорость удара составляет 8 км/сек. Угол столкновения ударника с преградой 83 градуса (рис 1).
Ре-С-А1,4 тт, 83 с1ед, 8 кт/зес
4^ 1;
Рис. 1 - Высокоскоростной удар на скорости 8 км/сек стальным шариком под углом 83 градуса по двухслойной преграде из углерода и алюминия. Представлены трехмерные распределения плотности в моменты времени 0 мкс - верх-нийрисунок и 20мкс - нижний рисунок
Столкновение ударника с преградой приводит к формированию за несколько первых микросекунд кратера в первой пластине из углерода. На рис. 1Ьвидно формирование облака частиц углерода и железа, образовавшихся в результате взаимодействия. Возникший кратер имеет глубину примерно 1 см и поперечный диаметр порядка 2 см.
В другой задаче рассматривается воздействие импульса длительностью 120фс на алюминиевую фольгу толщиной 5мкм. Интенсивность излучения в центре пятна нагрева составляет 1018 Вт/см2. Радиус пятна фокусировки, на котором интенсивность излучения уменьшается в е2 раз, составляет 10 мкм. Импульс падает на мишень слева. Расчет длится 200 пикосекунд, использованы сетки четырех уровней 11, 12, 13, 14 с фактором разбиения 2, 2, 2, 4, соответственно. Уровень 10 состоит из 64 х 64 ячеек, эффективная сетка составляет 2048 х 2048 ячеек с ша-
гом по пространству ~ 50нм. Размер расчетной области 100 х 100 мкм.
Сверхбыстрый нагрев мишени лазерным импульсом длительностью всего 120 фс приводит к концентрации огромной плотности мощности в малом объеме вещества, что приводит к возникновению огромных давлений и температур. В процессе взаимодействия излучения с мишенью происходит изохорический нагрев вещества, и как видно на рисунке 2, плавление. Затем формируется ударная волна, за фронтом которой вещество сжимается более чем в 3 раза от нормальной плотности. С передней поверхности мишени происходит интенсивная разгрузка вещества в газо-плазменное состояние. В момент времени 120пс происходит выход ударной волны на тыльную поверхность мишени. Разгрузка с тыльной поверхности мишени также приводит к формированию газовой фазы при интенсивном испарении вещества мишени (рис. 2).
Пте=2.00786е-005х Юткэ
Рис. 2 - Распределение фазовых состояний (верх-нийрисунок) и плотности (нижнийрисунок). Распределение фазовых состояний: Красный - газ, оранжевый - жидкость-пар, голубой - жидкость, малиновый - плавление, зеленый -твердое. Соответствующие метастабильные состояния: синий - газ, желтый - метастабиль-ная жидкость, фиолетовый - метастабильное плавление, светло зеленый - метастабильное твердое состояние. Моменты времени 0 пс -верхний; 200 пс - нижний
Выводы
Нами разработана широкодиапазонная модель многокомпонентной упругопластической среды с учетом эффектов поглощения фемтосекундного лазерного излучения. Реализована параллельная ком-
пьютерная программа, позволяющая с помощью разработанного подхода моделировать задачи, связанные с ударным взаимодействием объектов на дозвуковых и сверхзвуковых скоростях, а так же описывать процессы взаимодействия ультаракорот-ких лазерных импульсов с веществом.
Работа выполнена при финансовой поддержке программфундаментальных научных исследований Президиума РАН № 111П и № 1.33П.
Литература
[1] С.И.Ашитков, М.Б.Агранат, Г.И.Канель, П.С.Комаров, В.Е.Фортов, Письма вЖЭТФ, 92, 8, 568-573 (2010).
[2] С.И.Ашитков, П.С.Комаров, М.Б.Агранат, Г.И.Канель, В.Е.Фортов Письма в ЖЭТФ, 98, 439-444 (2013).
[3] V. A. Titarev, E. Romenski, E. F. Toro, Int. J. Numer. Meth. Eng., 73, 897 (2008).
[4] P. T. Barton, D. Drikakis, J. Comput. Physics,229 5518-5540(2010).
[5] N. Favrie, S. Gavrilyuk, S. Ndanou, J. Comput. Phys., 270, 300-324 (2014).
[6] G. H. Miller, E. G. Puckett, J. Comput. Phys. 128, 134164 (1996).
[7] I. Lomov, R. Pember, J. Greenough, B. Liu, Patch-based adaptive mesh refinement for multimaterial hydrodynamics. In Five-Laboratory Conference on Computational Mathematics, Vienna, October 2005.
[8] M. E. Povarnitsyn, N. E. Andreev, P. R. Levashov, K. V. Khishchenko, O. N. Rosmej, Phys. Plasmas,.19, 023110 (2012).
[9] M. E. Povarnitsyn, K. V. Khishchenko, P. R. Levashov, Int. J. Impact Eng., 33,1, 625-633 (2006).
[10] K. V. Khishchenko, V. E. Fortov, I. V. Lomonosov, Int. J. Thermophysics, 26, 2, 479-491 (2005).
[11] V. E. Fortov, I. V. Lomonosov, Shock Waves 20, 1, 5371 (2010).
[12] С. К. Годунов, А. В. Забродин, М. Я. Иванов, А. Н. Крайко, Г. П. Прокопов, Численное решение многомерных задач газовой динамики, Наука, Москва 1976, с. 400.
[13] M. J. Berger, J. Oliger, J. Comput. Phys, 53, 484 (1984).
[14] М. E. Поварницын, А. С. Захаренков, П. Р. Левашов, К. В. Хищенко,Вычислительные методы и программирование, 13, 424-433 (2012).
© М. Е. Поварницын - с.н.с., ОИВТ РАН, Москва, [email protected];A. В. Конюхов - с.н.с., ОИВТ РАН, Москва, [email protected]; П. Р. Левашов - зав. лаб., ОИВТ РАН, Москва, [email protected].
© M. E. Povarnitsyn - senior researcher JIHT RAS, Moscow, [email protected]; A. V. Konukhov - senior researcher JIHT RAS, Moscow, [email protected]; P. R. Levashov - head of the laboratory, JIHT RAS, Moscow, [email protected]
