CC BY
21
УДК 539.3
С. А. Мурзов1'2'3, А. Н. Паршиков1'3, С. А. Дьячков1'3, М. С. Егорова1'3,
В. В. Жаховский1'3
1 Всероссийский научно-исследовательский институт автоматики им. Н. Л. Духова, Росатом
2 Московский физико-технический институт (национальный исследовательский университет)
Объединенный институт высоких температур Российской академии наук
Моделирование стационарных ударных волн в пористом веществе методом БРН
Предлагаемый в работе метод мезоскопического расчета отклика пористых материалов на ударное сжатие в подвижном окне открывает возможность прямого расчёта ударных адиабат пористых материалов. Мезоскопическое моделирование описывает пористую среду явным заданием структуры пор и использует информацию только об уравнении состояния и прочностных характеристиках сплошного материала. Для достижения истинно стационарного режима распространения ударной волны разработан метод подвижного окна наблюдения, суть которого состоит в следующем: несжатое вещество втекает в расчетную область с постоянной скоростью, а скорость вытекания из области подбирается итерациями так, чтобы фронт волны установился неподвижным в окне наблюдения. Моделирование ударных волн было проведено как в стандартной постановке с неподвижным поршнем (методом обращенного движения), так и в системе подвижного окна наблюдения. Продемонстрировано, что после достижения стационарного режима профили, полученные обоими методами, идентичны.
Ключевые слова: ударные волны, гидродинамическое моделирование, упруго-пластика, метод сглаженных частиц.
S. A. Murzov1'2'3, A. N. Parshikov1'3, S. A. Dyachkov1'3, M. S. Egorova1'3,
V. V. Zhakhovsky1'3
1Dukhov Research Institute of Automatics 2 Moscow Institute of Physics and Technology 3Joint Institute for High Temperatures
Stationary shock waves in porous copper by SPH
simulation
The mesoscopic simulation of stationary shock waves simulation in porous materials by a moving window method is proposed and developed for smoothed particle hydrodynamics. The mesoscopic simulation describes the material by the explicit definition of porous structure and uses an equation of state and strength properties of the bulk material only. To produce stationary shocks in the simulation we apply the constant inflow velocity of porous material to the computational domain but an iterative method for outflow velocity. This iterative method stabilizes the shock front at the prescribed position. The simulation of the shock produced by the standard solid wall boundary condition is consistent with the stationary shock produced in a moving-window.
Key words: shock waves, hydrodynamic simulation, elastoplasticity, smoothed particle hydrodynamics.
© Мурзов С. А., А. Н. Паршиков, Дьячков С. А., Егорова М. С., Жаховский В. В., 2020
(с) Федеральное государственное автономное образовательное учреждение высшего образования
«Московский физико-технический институт (национальный исследовательский университет)», 2020
1. Введение
Моделирование стационарных быстро распространяющихся процессов, например ударная волна (УВ) или волна детонации (ДВ), требует высоких вычислительных затрат. Моделирование УВ в системе координат неподвижного образца или неподвижного поршня, где УВ распространяется но все большей и большей области вещества, приводит к увеличению затрачиваемых на расчет вычислительных ресурсов квадратично физическому времени счета, так как расчетная область растет пропорционально длине пробега ударной волны в веществе. При моделировании в таких системах координат в работе [1] исследовано ударное сжатие пористого алюминия методом обращенного движения, где показана необходимость увеличения времени расчета для исследования стационарности УВ с ростом пористости и снижением амплитуды воздействия.
В работе [2| был предложен метод молекулярно-динами ческохх) моделирования в подвижной системе координат (метод подвижного окна MW-MD). Данный метод позволяет получить покоящийся фронт УВ, а в расчетную область втекает постоянный ноток вещества со скоростью, равной по модулю скорости УВ us в системе координат неподвижного образца. Условие стационарности фронта определяет скорость потока на выходе из расчетной области равной ир — us. Аналогично, переходя из одной системы координат в другую, можно получить такую же УВ в системе координат покоящейся жесткой стенки, распространяющуюся со скоростью us — ир, если начальная скорость образца равна —ир (см. рис. 1). В том числе данная система координат дает возможность расчета с более мелким размером частиц, особенно для УВ в пористых образцах, где волны давления за фронтом могут осциллировать на больших длинах, например, при коллапсе пор, что ведет к необходимости постановки зоны выхода вещества и плоскости удаления на значительных расстояниях. Организация входных и выходных границ в методе сглаженных частиц (SPH) имеет три основных подхода: задание слоев пограничных частиц за проницаемыми границами [3,4], полуаналитические (semi-analytical) граничные условия с искуственным введением храничных сегментов (отрезков) [5], зеркальные отражения частиц за границы только по координате в несжимаемой жидкости [6].
В настоящей работе метод подвижного окна наблюдения модифицирован для целей гидродинамического моделирования стационарных У В с помощью параллельного кода [7], использующих) CSPH (contact SPH) [8,9]. Параллельный программный комплекс показывает высокую эффективность обратная линейная зависимость времени расчета от числа процессорных ядер, задействованных в расчете, в том числе для динамических задач.
Po> Ро,Т0
I
Shock
p,p,T p,p0,T0
Moving «Piston»
«Piston at rest»
«Shock at rest»
H up us \ i up=0 Us-Up -up up-us -us
V
Рис. 1. Различимо системы координат в моделировании ударной волны
1.1. Математическая модель
Модель формулируется в приближении сплошной среды, в котором используются законы сохранения массы, импульса и энергии, дополненные материальными уравнениями, описывающими механическое поведение среды:
р + р V-U = 0,
(1)
рЬ - V • а = 0, (2)
рЁ - V • (стИ) = 0. (3)
В случае жидкости девиаториая часть напряжений равна нулю [10]. Для определения девиатора напряжений упругой среды используется закон Гука:
а = Б - Р(р,е)1, (4)
& - ОБ + БО = 2Се, (5)
Где § — девиаториая часть тензора напряжений, а — тензор напряжений, Р(р, е) — давление, определяемое из уравнения состояния, р — плотность, и — вектор скорости, Е = е + И2/2 — сумма внутренней и кинетической энергий единицы массы вещества, I — единичный тензор, С — модуль сдвига упругого материала, е — тензор девиаторной части скорости деформации, О — матрица угловой скорости частицы. Девиаториая часть скорости деформации
е = е - ¿г(ё)/3, (6)
где полная скорость деформации и угловая скорость в частице определяется из
е = ^® И)т + V® И, (7)
О = 0.5^® И - (V® И)т). (8)
Тензорное произведение (V ® И)т = VIгде V = д/дгг — частная производная вдоль оси п. Мы используем метод сглаженных частиц, который, в отличие от сеточных схем, вместо континуума рассматривает дискретное представление среды в виде системы частиц, допускающих произвольную связность друг с другом. Использование решения задачи задачи Римана позволяет обойтись без добавления искусственной вязкости [11] и обеспечить монотонность решения вблизи контактных границ. В расчетах, представленных в данной статье, используется СБРН, где взаимодействие частиц обеспечивает межчастичное решение задачи Римана [12] на контактной границе между частицами. При этом мы используем консервативную формулировку закона сохранения энергии. При расчете тензора скорости деформации учитывается его вращение аппроксимацией из [8,10].
1.2. Модель среды
Уравнение состояния на основе ударной адиабаты для меди задано в форме Ми— Грюнайзена:
Р - Рге/ =Тр(е - еге4), (9)
где е — удельная внутренняя энергия, Pref (р) и eref (р) — опорная кривая, Г — параметр Грюнайзена. При плотности вещества р выше тачальной ро (х = ро/р < 1) состояние вещества описывается ударной адиабатой вида и3 = са + ир, где и3 — скорость ударной волны, ир — скорость вещества, са,8а — коэффициенты. Опорные ударные адиабаты для давления и энергии:
_ 2 1 Ж _ 1 Ж у _ ч
= Р^а [1 - ^(1 - Х)]2 , ^ = ^ ■ (10)
При плотности вещества ниже начальной (х > 1) использовалась линейно-упругая кривая давления с согласованным модулем объемного сжатия В = ро(?а'-
Рте/ = р0С2а(1 - X), еге/ = Са(1 - х)2/2. (11)
Модель упругопластического течения материала определяется на основе предела текучести и условия на цилиндре Мизеса, ограничивая величину интенсивности напряжений [10] тед = \!г(Б2) соответствующим масштабированием компонент девиаторов напряжений
^ = («ау (12)
в^У/Тед, если (Тед > У,
где У — предел текучести.
Таблица1
Параметры упруго-пластической модели меди
Ро [г/см3] са [км/с] 8а О [ГПа] У [ГПа]
8.96 4.148 1.4205 43 0.35
Задача распада упругопластического разрыва имеет аналитическое решение в случае линейной модели объемного сжатия. Модель одномерной линейной упругопластической среды описывается уравнениями движения без учета уравнения энергии и включает в себя систему уравнений сохранения массы и импульса, уравнения состояния, закона Гука и условия текучести Мизеса. Зададим уравнение состояния меди с начальной плотностью ро = 8.96 г/см3 с модулем объемного сжатия В = 153.86 ГПа, модулем сдвига О = 43 ГПа и пределом текучести У = 0.23 ГПа. Рассмотрим задачу распада разрыва, где состояние
( = 5
( = 0
Уравнение состояния определяется дифференциальным соотношением
Р = Ву, (13)
где V — объем, а интегрирование уравнения проводится при условии Р = 0 при V = V), у0 — объем, соответствующий плотности ро. Следовательно, уравнение состояния имеет вид
Р = В 1п( V). (14)
Отметим, что уравнение отличается от привычной формы записи закона Гука Р = В(У — V))/Уо, который предполагает малые изменения объема, но сжатие до нулевого объема происходит при конечном внешнем давлении.
В упруго-пластической среде напряжение может быть разложено на шаровую и девиа-
4
закона необходимо лишь для девиаторной части напряжений, то есть для уравнения (5). Шаровая часть напряжения определяется уравнением состояния.
С использованием перечисленных предположений решение распада произвольного разрыва с начальными условиями слева и справа соответственно:
Т ТЮ Т 0
р = рю р = ро
и 0 и 0
где Тю = — 5 ГПа, рю = ро(1 — (2/3У — тю)/В), то есть вещество задается в состоянии одноосного сжатия и находится в пластическом состоянии. При условии контакта в начальный момент времени в точке = 0, для случая пластического течения слева и справа
т* = 0.5(тю — (1 — тсгсь(р1 )/сг (р))), р* = рю (1 — (2/3У — т10 )/В),
— 2асг (Гг Gcr
Pi = р* с1(р*)/(с1 (р*)+ щ) Рг = Ро ci(po)/(ci (ро) — ur
2асг/р*/а(р*) ur <?cr/ро/С1(ро)
* (2р1 щ сь (р*) + (аш ) - (0 + асг))
р1С1 (р* ) + рг С1 (ро) '
= -ч(р*^, хег = с1(ро)г, ХРг = -(сь(р*) - ЩХрг = (сЬ(ро) - иг где (Тег = У (В + 40/3)/(20) — критическое одноосное напряжение, при котором течение
е е Р Р
становится пластическим; х\,хг,х1 ,Хг — положения левых и правых упругих и пластических фронтов волн; р,и,а — плотность, скорость и напряжение. Сравнение решения задачи распада упругоилаетичеекого разрыва методом СБРН с аналитическим решением, представленное на рис. 2, показывает, что с увеличением числа БРН-частиц решение сходится к аналитическому. Вопрос о порядке аппроксимации метода БРН рассмотрен в [13]. Для интегрирования уравнений по времени используется метод Эйлера первого порядка аппроксимации. Продольная и объемная скорость звука определяются соотношениями
а (р) =
/
В + 4G/3
сь(Р) = \ I ^т-
г
р
5x109 4xl09
£ 3x109
И 9
е 2xio9
I
l09 0
o 4x10-8 8x10-8 1.2x10-7
axis x (m)
Рис. 2. Распад упругопластического разрыва в меди: слова ударно-сжатое вещество при 5 ГПа, справа при пулевом давлении. Результат сравнения с аналитическим решением двух разрешений сетки при одномерной постановке с использованием двумерного контактного метода сглаженных частиц. При этом синему графику соответствует 120 SPH частиц в расчетной области на длину образца, а красному 1920 SPH частиц. Результат моделирования показан на момент времени t = 6 пс
V
1
1
I
_i_L
2. Метод получения стационарной ударной волны 2.1. Начальная расчетная область и граничные условия
Расчетная область моделирования, схематически изображенная на рис. 3, состоит из зоны подачи вещества с правой границы расчетной области, зоны ударного сжатия и зоны выхода вещества из расчетной области с последующим удалением. При моделировании методом сглаженных частиц происходит создание частиц при входе в расчетную область и их удаление при выходе из расчетной области. Таким образом, расчет производится с некоторым ограниченным количеством частиц и может производиться на больших временах счета по сравнению с расчетом в неподвижной системе координат.
В образце вырезаются сферические поры, которые находятся в узлах решетки кубического гранецентрированпого кристалла (ГЦК), причем ось X является направлением [100] решетки. Поперечные размеры образца составляют Ьу х = 20 х 20 рм2, вдоль этих направлений заданы периодические граничные условия. Начальный размер БРН-частиц равен равен 1/\/2 рм, что соответствует 20 периодам ГЦК-решетки в поперечном направлении расчетной области.
После этого запускается алгоритм упаковки частиц [14], близкой к жидкой. Метод подготовки жидкой упаковки суть временная замена вещества на другое достаточно жесткое вещество и включение случайных объемных сил, способствующих перемешиванию частиц с постепенным затуханием амплитуды объемной силы, а также введение сил скоростного трения для остановки частиц и последующей подменой вещества на исходное вещество с новыми положениями БРН-частиц. После данной процедуры в образце вырезаются поры необходимого для моделирования размера.
Решетка пор имеет период, равный поперечному размеру расчетной области, а радиус пор составляет 4.9854 рм и соответствует пористости т = р/ро = 1.35. Длина образца вдоль оси X составляла в пашем случае приблизительно 0.57 мм. Размер образца определяется положением правой границы подачи вещества и левой фиксированной границей удаления.
БРН-частицы удаляются после пересечения плоскости с фиксированной координатой х = хд[. Граничными считаются частицы внутри расчетной области, вблизи плоскости удаления, а именно для частиц с координатами х ^ ха1+йо, где йо — ширина зоны выхода, ха1 — левая граница расчетной области. В этой зоне решается лишь уравнение движения для координат БРН-частиц, при этом скорость всех БРН-частиц равна скорости выхода и,х = иоиг, а поперечные компоненты скорости полагаются равными нулю. Все остальные величины — плотность, энергия и другие термодинамические величины — остаются неизменными и равны последним своим значениям перед попаданием в область действия условия выхода. Частицы в граничной зоне являются видимыми для частиц внутри зоны расчета уравнений среды, создают и поддерживают ударную волну, соответствующую разнице массовой скорости, давления и плотности. При этом частицы внутри граничной зоны двигаются со скоростью выхода и фиксированным термодинамическим состоянием.
Граница подачи вещества имеет переменное положение, которое может варьироваться в ходе расчета на величину в несколько слоев БРН-частиц. Это связано с тем, что вставка нового слоя возможна лишь в момент обновления списка соседствующих БРН-частиц. Обновление списка соседей происходит через одинаковое число шагов интегрирования (примерно 5-15), что соответствует времени А^, но переменный шаг по времени в гидродинамическом моделировании приводит к тому, что правая граница, двигаясь с постоянной скоростью п3, проходит разное расстояние за одинаковое число шагов интегрирования. Положение границы равномерно двигается в отрицательном направлении оси X, но скачком изменяет свое положение в положительном направлении при вставке новых слоев БРН-частиц. Поэтому граница вблизи подачи вещества характеризуется дискретностью своего положения. На схеме рис. 4 изображена зона входа частиц в расчетную область, где появление нового набора слоев частиц выделено темно-красным цветом, синим показаны частицы, для которых решается система уравнений (1) - (5), а оставшиеся красные частицы будут вставляться по мере освобождающегося места вблизи правой границы образца. Положение вставляемого слоя внутри подаваемого учитывает положение последнего вставленного слоя. Когда вставка вещества доходит до конца красного образца, то слои начинают вставляться с первого, то есть периодически повторяются. Частицы в зоне входа вещества имеют координаты х ^ хдг — йг, йг — ширина зоны входа, хдг — правая граница расчетной области. Границы зоны входа показаны синими пунктирными линиями на рис. 4. Частицы внутри зоны входа движутся равномерно и прямолинейно вплоть до попадания в зону расчета уравнений гидродинамики, при этом массовая скорость всех частиц направлена вдоль оси X и равна скорости входа пх = пгп-
Скорость входа является фиксированной на всем времени моделирования и равна скорости ударной волны в стационарном режиме. Скорость выхода является переменной величиной, которая требует определения из условия стационарности ударного фронта. Значение скорости выхода изменяется, стремясь остановить ударный фронт в данной системе координат в желаемом положении ударного фронта внутри расчетной области.
outflow
а
м
"3 •а
ир=иои()
Рис. 3. Схема расчетной области в моделировании ударной волны в подвижной системе координат методом сглаженных частиц. Скорость входа является постоянной величиной, равной скорости ударной волны и81 при этом скорость выхода изменяется во времени Ио-^^)
Программная реализация метода подвижного окна произведена в рамках архитектуры параллельного программного комплекса СЯРН-УБЭ, разрабатываемого авторами [7,14].
Рис. 4. Схема границы входа вещества в расчетную область. Хотя вставка нового слоя происходит через одинаковое число шагов интегрирования, но, вообще говоря, разное время ДЬь, что приводит к переменной толщине вставляемого слоя новых частиц, кратной ширине одного слоя ГЦК-решетки начальной упаковки частиц. Вставка новых слоев частиц показана темно-красным цветом, при этом последующие слои будут учитывать последнее положение слоя внутри потенциально вставляемого объема (красные частицы)
2.2. Получение стационарного ударного фронта
Пусть в начальный момент времени в окне задан образец с определенной мезоскопиче-ской структурой, состоящий из Ио частиц и со средней начальной плотностью роо- Если в системе распространяется УВ и вещество сжатое будет при большей плотности, то при достижении стационарного режима следует ожидать увеличения полного числа ЯРН-частиц N1 относительно начального, то есть Nf > N0. Условие стационарности требует такого положения ударного фронта, когда количество частиц равно Nf и ударный фронт неподвижен. Однако в случае большой пористости возможен аномальный ход ударной адиабаты, когда плотность за ударным фронтом уменьшается, и тогда Nf < N0. В этом случае предпочтительнее иметь дело не с числом частиц, а со средним давлением в окне наблюдения Р — / р (х)йх/Ь ж, которое в У В всегда больше начально го давления Р^ > Ро — Ро-
Опишем итерационный алгоритм обратной связи для подбора массовой скорости на выходе из окна, а также технику проведения такого расчета на примере сплошного образца. Цель этого алгоритма состоит в схождении к такой иои1, при которой ударный фронт стабилизируется в заданном положении в окне.
Задание точного решения — ир — и8 в начале сформирует ударный фронт у левой границы. Для распространения начальной УВ к центру окна ее скорость должна быть больше, чем и8, чтобы ударный фронт достиг или После этого амплитуда У В должна уменьшаться, ударный фронт будет замедляться, затем начнет сноситься втекающим пото-
ком к желаемому положению. Алгоритм обратной связи должен быть сконструирован так, чтобы, управляя и^, привести положение ударного фронта к стационарному положению.
Искомый алгоритм должен автоматически адаптироваться к УВ разной амплитуды, чтобы корректировка скорости выхода вещества всегда приводила к формированию стационарной УВ. Для этого вводится управление скоростью выхода с помощью явной зависимости от скорости выхода щп = — и8, безразмерной наблюдаемой в моделировании целевой величины ш и некоторого подгоночного параметра
Uout = uin( 1 - exp(-^w)),
(15)
где измеряемая в процессе моделирования целевая функция от числа частиц ш = ^ 1 при N ^ Nf. Если целевой величиной управления является среднее давле-
ние Р, то ш = (Рг — Ро)/(Р/ — Ро) ^ 1 при Р ^ Р/. График функции, задаваемой уравнением (15) при начальном значении параметра @ = 1, показан зеленой пунктирной кривой на рис. 5. Если бы параметр @ не менялся, то процесс выхода на стационарный режим шел бы по этой кривой до пересечения с горизонтальной линией и^ = ир — и8. При этом расчет начинается со значения ш = 0 (А = Ао ми Р = Ро ), в образце возникает сильная УВ, что приводит к росту ш и соответствующему изменению скорости выхода согласно уравнению (15). В результате положение ударного фронта сходится к некоторому наперед неизвестному положению в окне при некотором ш* = 1 в общем случае.
Чтобы установить ударный фронт в желаемом положении, определяемом значением целевой функции ш = 1, мы запишем параметр в виде
Р * =
ln(1 - Uout/Uin)
ш*
(16)
Можно было заметить, что в стационарном режиме и^ = ир — и8,щп = —и8, поэтому произведение = сош^ тогда искомое значение параметра @* = ш*. Теперь для получения стационарного фронта будет использоваться функция, заданная уравнением (15) при р = 0*, которое приведет УВ в желаемое положение при ш = 1. Для более гладкого воздействия на образец скорость выхода ограничена значением в данном расчете величиной 2 км/с. График функции скорости выхода от числа частиц показан сплошной зеленой линией на рис. 5.
time (цс)
scaled particles number w
Рис. 5
Можно использовать другую зависимость, отличную от заданной уравнением (15), например линейный закон, однако функция должна быть обязательно убывающей, чтобы гарантировать сходимость целевой функции:
и^а = —¡3/2(и + 1), (17)
где подгоночный параметр @ теперь является размерным и равен по модулю скорости выхода в стационарном режиме данная функция изображена синей линией на рис. 5. Значение
параметра @ было выбрано так, чтобы кривая, заданная уравнением (17), пересекала кривую стационарной скорости выхода при ш — 1.
На левом графике рис. 5 показаны истории скорости выхода при использовании двух различных зависимостей и^, заданных уравнениями (17) и (15). При этом переход к стационарному режиму и установление ударного фронта в желаемом положении проводилось в две стадии:
1) вывод па стационарное значение ш* при некотором пробном значении параметра @* (мы полагали — 1);
2) корректировка параметра @ таким образом, чтобы провести расчет с выводом в стационарный режим при желаемом значении ш — 1.
Для исключения двухстадийности расчета требуется управление параметром @ в процессе сквозного моделирования с целью ш — 1. Рассмотрим уравнение для скорости выхода от числа частиц, заданное линейным убывающим законом (17). Чтобы прийти к @ ^ и^ и ш ^ 1, введем корректирующее уравнение на изменение параметра @ за одну итерацию:
ДА — аДи(Дшг — Ди?), (18)
где г - номер итерации, а - управляющий параметр, а также используются скорости наблюдаемого изменения целевой функции Дшг — — Шг-к)/Д£ за время Д£ за к шагов итерации и «желаемой» скорости сходимости Дш'^ — (1 — ш£)/т. Таким образом, если скорость сходимости, наблюдаемая в моделировании, отличается от скорости желательной для того, чтобы решение сходилось к ш — 1 за характерное время т, параметр @ скорректируется в соответствии с уравнением (18).
По мере приближения к значению ш — 1 выгодно применять более топкую подстройку скорости выхода и^, так как желаемая скорость сходимости за время т может становиться фактором, тормозящим сходимость, если для достижения нужного положения УВ требуется время, меньшее т. Предлагается применять более топкую подстройку, ослабляя воздействие уравнения (18) при условии \ш —1|<е, где полагается е — 0.03. С этой целью модифицируем управляющее уравнение для скорости выхода (17) и зависимость управления параметром [3 (18) следующим образом:
и<шг — —Р/2[(ш + 1) — (ш — 1) ■ С И], (19)
ДД — аД1ь( 1 — ((иг))(Диг — Ди?), (20)
где функция ((ш) представляет собой купол единичной высоты с некоторой конечной шириной:
<И — Г+ТоЬм • (21)
где ширина а — 50 функции ((ш) была выбрана эмпирически. Измененная таким образом зависимость скорости выхода и^ показана пунктирной красной линией на рис. 5. Качественный смысл приведённого алгоритма управления скорости выходного течения состоит в следующем:
1) при — 1| ^ 1/а, управление параметром @ осуществляется с использованием корректирующего уравнения (20) и управляющего уравнения (19);
2) при приближении ш ^ 1 контроль скорости выхода ослабляется и корректирующее уравнение слабо воздействует на параметр что снижает уровень шума вносимый в профиль У В при его изменении. Куполообразная форма функции позволяет достигать слабого воздействия на управляющий параметр, пока система стационарном режиме;
3) если система выходит из положения близкого к равновесию механизм управления снова становится более жестким;
3. У пру гоп л астическое моделирование
Получены стационарные ударные волны в пористой меди с использованием подвижного окна наблюдения и вышеописанной методики установления стационарного ударного фронта. Моделирование проводилось при скорости ударной волны us = —щп = 3.5 км/с.
-0.2 -0.1 о
axis x (mm)
-0.2 -0.1 0 axis x (mm)
-0.2 -0.1
axis x (mm)
Рис. 6. Профили плотности, давления, девиатора напряжений и скорости при моделировании стационарной ударной волны со скоростью и3 = 3.5 км/с сравниваются с результатом моделирования ударной волны методом обращенного движения. Двумерный срез середины образца с явным изображением ЯРН-частиц точками имеет раскраску, соответствующую плотности частиц
Скорость выхода иои1 позволяет изменить характер выхода ударной волны в стационарное значение. При этом в начальный период времени по среде должна пройти УВ большей амплитуды, чем ожидается в стационарном положении, чтобы иметь и большую скорость распространения, чем скорость подачи вещества. Начальное значение скорости выхода и0иг= 0) оценено из известных ударных адиабат [15] для пористой меди. Использовалась ударная адиабата меди с начальной плотностью 6.3 г/см3. Плотность кристаллической меди используемого уравнения состояния 8.96 г/см3, для пористости т = 1.35 плотность об-
3
нение массовой скорости ир в ударной волне при малом изменений пористости для фиксированной скорости ударной волны и3, можно положить ир(ро = 6.637) ^ 6.3/6.637^р(р0 = 6.3).
Профили скорости и давления имеют осциллирующий характер за фронтом стационарной УВ (см рис. 6), причем колебания находятся практически в противофазе, что объясняется многочисленными волнами сжатия и разрежения за фронтом, а также вследствие перехода кинетической энергии в потенциальную энергию. Стоит отметить, что амплитуда колебаний затухает и этому способствуют различные механизмы диссипации механической энергии в тепло - процесс многократного сжатия и разгрузки, в том числе схемные эффекты метода ЭРН.
На рис. 6 величина плотности ЭРН-частиц до прихода УВ имеет периодический характер, что связано с явным заданием мезоскопической структуры и периодичностью геометрии пор. При этом осцилляции плотности за областью релаксации давления являются результатом неоднородной плотности ЭРН-частиц, вызванной неоднородным прогревом в процессе схлопывания пор. Показанные профили были подобраны таким образом, чтобы распространяющаяся УВ была в фазе в обеих системах отсчета, что позволило получить явное сопоставление профилей. Представлено сравнение профилей, полученных методом подвижного окна и методом обращенного движения. Моделирование проводилось в следующей последовательности: во-первых, была получена стационарная УВ методом подвижного окна наблюдения. Стационарная У В из метода подвижного окна наблюдения определяет значение скорости вещества ир из соотношения со скоростью выхода и0ы = —2.372, откуда ир = 1.128 км/с. Во-вторых, в методе обращенного движения жесткая стенка [16] неподвижна (см. рис. 1), при этом образец имеет скорость их = —ир. В результате столкновения образца со стенкой образуется УВ, которая имеет одинаковые характеристики с УВ, полученной из моделирования методом подвижного окна наблюдения.
4. Заключение
Метод моделирования в подвижном окне наблюдения позволил провести моделирование стационарных УВ в пористой меди методом SPH. Мезоскопическое моделирование УВ позволяет наблюдать процесс схлопывания пор в пористой среде, используя данные для сплошного вещества — уравнение состояния и механические свойства. Полученные стационарные У В имеют протяженную зону осцилляций напряжений за фронтом. Моделирование методом подвижного окна согласуется с моделированием УВ, полученной столкновением пористого образца с неподвижной жесткой стенкой. Данная методика моделирования применима для получения ударных адиабат пористых материалов, если известны лишь свойства каркаса сплошного вещества.
Литература
1. Medin S.A., Parshikov A.N. Mesomechanical Simulation of Shock Compaction of Porous Aluminum // Math. Models Comput. Simul. 2014. V. 6, N 5. P. 435-444.
2. Zhakhovskii V. V, Zybin S. V, Nishihara K., Anisimov S.I. Shock Wave Structure in Lennard-Jones Crystal via Molecular Dynamics // Phvs. Rev. Lett. 1999. V. 83. P. 1175 1178.
3. Lastiwka M., Basa M., Quinlan N.J. Permeable and non-reflecting boundary conditions in SPH // Int. J. Numer. Meth. Fluids. 2009. V. 61. P. 709-724.
4. Sun P.-N., Colagrossi A., Zhang A. Numerical simulation of the self-propulsive motion of a fishlike swimming foil using the ¿+-SPH model // Theor. App. Mech. Lett. 2018. V. 8. P. 115-125.
5. Ferrand M., Laurence D. R., Rogers B. D., Violeau D., and Kassiotis C. Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method // Int. J. Numer. Meth. Fluids. 2012. V. 71. P. 446-472.
6. Kunz P., Hirschler M., Huber M., Nieken U. Inflow/outflow with Dirichlet boundary conditions for pressure in ISPH //J. Сотр. Phvs. 2016. V. 326. P. 171-187.
7. Egorova M.S., Dyachkov S.A., Parshikov A.N., Zhakhovsky V.V. SPH modeling using dynamic domain decomposition and load balancing displacement of Voronoi subdomains // Comput. Phvs. Commun. 2019. V. 234. P. 112-125.
8. Parshikov A.N., Medin S.A. Smoothed particle hydrodynamics using interparticle contact algorithms // J. Comput. Phvs. 2002. V. 180. P. 358-382.
9. Medin S.A., Parshikov A.N. Development of smoothed particle hydrodynamics method and its application in the hydrodynamics of condensed matter // High Temp. 2010. V. 46. P. 973-980.
10. Олдер В., Фернбах С., Ротенберг M. (ред.) Вычислительные методы в гидродинамике. Москва : Мир, 1967.
11. Monaghan J.J. SPH and Riemann Solvers // J. Comput. Phvs. 1997. V. 136, N 2. P. 298307.
12. Dukowicz J.K. A general, non-iterative Riemann solver for Godunov's method //J. Comput. Phvs. 1985. V. 61. P. 119-137.
13. Quinlan N. J., Basa M., Lastiwka M. Truncation error in mesh-free particle methods // Int. J. Numer. Methods. Eng. 2006. V. 66, N 13. P. 2064-2085.
14. Dyachkov S.A., Parshikov A.N., Zhakhovsky V.V. SPH simulation of boron carbide under shock compression with different failure models //J. Phvs.: Conf. Ser. 2017. V. 815. P. 012012.
15. S.P Marsh LASL shock Hugoniot data. Los Alamos National Laboratory, Los Alamos, NM, USA, 1980.
16. Takeda H., Miyama S.M., Sekiya M. Numerical simulation of viscous flow by smoothed particle hydrodynamics // Prog. Theor. Phvs. 1994. V. 92, N 5. P. 939-960.
References
1. Medin S.A., Parshikov A.N. Mesomechanical Simulation of Shock Compaction of Porous Aluminum. Math. Models Comput. Simul. 2014. V. 6, N 5. P. 435-444.
2. Zhakhovskii V. V, Zybin S. V, Nishihara K., Anisimov S.I. Shock Wave Structure in Lennard-Jones Crystal via Molecular Dynamics. Phvs. Rev. Lett. 1999. V. 83. P. 11751178.
3. Lastiwka M., Basa M., Quinlan N.J. Permeable and non-reflecting boundary conditions in SPH. Int. J. Numer. Meth. Fluids. 2009. V. 61. P. 709-724.
4. Sun P.-N., Colagrossi A., Zhang A. Numerical simulation of the self-propulsive motion of a fishlike swimming foil using the ¿+-SPH model. Theor. App. Mech. Lett. 2018. V. 8. P. 115-125.
5. Ferrand M., Laurence D. R., Rogers B. D., Violeau D., and Kassiotis C. Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method. Int. J. Numer. Meth. Fluids. 2012. V. 71. P. 446-472.
6. Kunz P., Hirschler M., Huber M., Nieken U. Inflow/outflow with Dirichlet boundary conditions for pressure in ISPH. J. Comp. Phvs. 2016. V. 326. P. 171-187.
7. Egorova M.S., Dyachkov S.A., Parshikov A.N., Zhakhovsky V.V. SPH modeling using dynamic domain decomposition and load balancing displacement of Voronoi subdomains. Comput. Phvs. Commun. 2019. V. 234. P. 112-125.
8. Parshikov A.N., Medin S.A. Smoothed particle hydrodynamics using interparticle contact algorithms. J. Comput. Phvs. 2002. V. 180. P. 358-382."
9. Medin S.A., Parshikov A.N. Development of smoothed particle hydrodynamics method and its application in the hydrodynamics of condensed matter. High Temp. 2010. V. 46. P. 973980.
10. Older B., Fernbach S., Rotenberrg M. (ed.) Vvchislitelnie metodi v gvdrodynamvke. Moscow : Mir, 1967.
11. Monaghan J.J. SPH and Riemann Solvers. J. Comput. Phvs. 1997. V. 136, N 2. P. 298-307.
12. Dukowicz J.K. A general, non-iterative Riemann solver for Godunov's method. J. Comput. Phvs. 1985. V. 61. P. 119-137.
13. Quinlan N. J., Basa M., Lastiwka M. Truncation error in mesh-free particle methods. Int. J. Numer. Methods. Eng. 2006. V. 66, N 13. P. 2064-2085.
14. Dyachkov S.A., Parshikov A.N., Zhakhovsky V.V. SPH simulation of boron carbide under shock compression with different failure models. J. Phvs.: Conf. Ser. 2017. V. 815. P. 012012.
15. S.P Marsh LASL shock Hugoniot data. Los Alamos National Laboratory, Los Alamos, NM, USA, 1980.
16. Takeda H., Miyama S.M., Sekiya M. Numerical simulation of viscous flow by smoothed particle hydrodynamics. Prog. Theor. Phvs. 1994. V. 92, N 5. P. 939-960.
Поступим в редакцию 30.12.2019
