Теплофизика высоких температур, 2021, T. 59, № 2, стр. 221-230

ВВЕДЕНИЕ

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

Мезоскопическое моделирование УВ [1–4] в системе координат неподвижного образца или неподвижного поршня, как схематически изображено на рис. 1, приводит к увеличению затрачиваемого вычислительного времени пропорционально квадрату длины пробега УВ (при постепенной подаче несжатого материала в расчетную область перед фронтом УВ), так как УВ распространяется по все большей области материала и расчетное время, затрачиваемое на шаг моделирования, растет пропорционально длине пробега УВ. При этом ударная волна в веществе с мезоструктурой в виде пор, гранул или смеси из разных веществ приводит к появлению зоны осцилляции напряжений за фронтом УВ. В работе [1] исследовано ударное сжатие пористого алюминия методом обращенного движения, где показана необходимость увеличения времени расчета для исследования стационарности УВ с ростом пористости и снижением амплитуды воздействия. Работа авторов [2] представляет первые попытки учитывать влияние неоднородностей, пор и ансамблей пор в трехмерном гидродинамическом моделировании. К настоящему времени мезоскопическое моделирование ударных волн в пористых материалах проводится в том числе методом молекулярной динамики [3] на нанометровом масштабе пор. Мезоскопическое моделирование распространения ударной волны по реагирующему материалу с переходом ударной волны в детонационную [4] может поставить отдельные задачи для крупномасштабного моделирования, где характерные времена термически активированных реакций разложения при схлопывании пор и длина пробега до детонации могут различаться на несколько порядков.

В организации входных и выходных границ в расчете методом сглаженных частиц в гидродинамике (SPH) имеется три доминирующих подхода. Первый из них представлен работами [5–8], где предлагается задавать границы слоями пограничных частиц, размещаемыми за проницаемыми границами. Положение этих частиц меняется в соответствии с их скоростью. Так, пограничная частица, переходя через границу входного условия, становится реальной расчетной. И наоборот, реальная частица, переходящая границу выходного условия, становится пограничной и затем удаляется по мере удаления от этой границы. В работе [5] введены неотражающие граничные условия входа и выхода, где несколько слоев частиц вблизи границ получают свое состояние на основе оценки значений инвариантов Римана для частиц из расчетной области и условия нулевых инвариантов Римана, соответствующих акустическим волнам, распространяющимся внутрь области от границ. Данные граничные условия применимы в том числе для дозвукового режима распространения течения и способны моделировать вытекание вещества из расчетной области.

Второй подход [9] является обобщением полуаналитических граничных условий с использованием искусственно введенных сегментов (отрезков), описывающих границу [10]. Так, новые частицы у входной границы появляются на месте точек-частиц контакта двух соседних сегментов в соответствии с дифференциальным уравнением изменения массы по времени, определяющим плотность потока массы. Удаление частиц на выходе происходит с поглощением массы точками-частицами. Третий подход реализует зеркальные отражения частиц за границы [11–13], при этом отражение частиц происходит лишь по координате. В настоящей работе используется подход, наиболее близкий к представленным в [5, 6], так как он является самым простым для реализации из трех рассмотренных.

Моделирование в подвижном окне наблюдения позволяет достигнуть истинно стационарного состояния ударной волны в пористом материале. В работе [14] был предложен метод молекулярно-динамического моделирования в подвижной системе координат, привязанной к фронту УВ. В этом методе фронт УВ покоится в расчетной области, а несжатое вещество втекает в него с постоянной скоростью, равной по модулю скорости движения УВ. Для получения эквивалентной ударной волны в системе координат неподвижной жесткой стенки (рис. 1) необходимо произвести моделирование столкновения образца при начальной скорости −u_p о неподвижную жесткую стенку. В результате столкновения со стенкой вещество останавливается и в нем формируется УВ, распространяющаяся со скоростью u_s – u_p.

В настоящей работе этот метод модифицирован для целей гидродинамического моделирования стационарных УВ с помощью параллельного кода [15], использующего бессеточный контактный метод сглаженных частиц (contact smoothed particle hydrodynamics – CSPH) [16, 17]. Метод CSPH является модификацией метода SPH, где скорости и напряжения в точке контакта SPH-частиц определяются из решения задачи Римана.

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

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ СРЕДЫ

Модель формулируется в приближении сплошной среды, в котором используются законы сохранения массы, импульса и энергии, дополненные материальными уравнениями, описывающими механический отклик среды на деформацию [18]:

Для определения девиатора напряжений в упругой среде используется закон Гука

где S – девиаторная часть тензора напряжений; σ – тензор напряжений; P(ρ, e) – давление (шаровая часть тензора напряжений, взятая с обратным знаком, определяемая из уравнения состояния); ρ – плотность; U – вектор массовой скорости; E = E_int + (U²/2) – сумма внутренней и кинетической энергий единицы массы вещества; I – единичный тензор; G – модуль сдвига упругого материала; ε – тензор девиаторной части скорости деформации; Ω – матрица угловой скорости частицы. В случае жидкости девиаторная часть напряжений равна нулю.

Тензор девиаторной части скорости деформации

где полная скорость деформации и угловая скорость в частице определяются из уравнений ⊗ – знак тензорного произведения. Использованный в работе метод сглаженных частиц, в отличие от сеточных схем, вместо континуума рассматривает дискретное представление среды в виде системы частиц, допускающих произвольную связность друг с другом. Использование решения задачи о распаде произвольного разрыва (задачи Римана) позволяет обойтись без добавления искусственной вязкости и обеспечить монотонность решения вблизи контактных границ. Используется консервативная формулировка закона сохранения энергии. При расчете тензора скорости деформации учитывается его вращение аппроксимацией из [16, 18].

Уравнение состояния меди задано в форме Ми–Грюнайзена:

где e – удельная внутренняя энергия; Γ – параметр Грюнайзена; в качестве опорных кривых давления и внутренней энергии используются их величины на ударной адиабате:

При плотности вещества ρ выше начальной ρ₀ (x = ρ₀/ρ < 1) состояние вещества описывается ударной адиабатой вида u_s = c_a + s_au_p, где u_s – скорость ударной волны, u_p – скорость вещества, c_a, s_a – коэффициенты. При плотности вещества ниже начальной (x > 1) использовалась линейно-упругая кривая давления с согласованным модулем объемного сжатия:

При моделировании ударных волн в стандартной постановке с неподвижным поршнем материал рассматривался как сжимаемая жидкость, уравнение состояния для которой задано в форме Ми–Грюнайзена. Константы вещества приведены в таблице (модель 1).

Для моделирования в подвижной системе координат используется модель 2, константы которой соответствуют результатам молекулярно-динамического моделирования меди. Аппроксимация ударной адиабаты получена в моделировании монокристаллической меди методом [14]. Модель упругопластического течения использует условия на цилиндре текучести Мизеса, ограничивая величину интенсивности напряжений

соответствующим масштабированием компонент девиатора напряжений где Y – предел текучести.

МЕТОД ПОЛУЧЕНИЯ СТАЦИОНАРНОЙ УДАРНОЙ ВОЛНЫ

Моделирование в подвижной системе координат имеет целью получение стационарного режима распространения УВ, при котором ударный фронт остается неподвижным относительно расчетной области (окна наблюдения). Расчетная область моделирования, схематически изображенная на рис. 2, состоит из зон подачи вещества через правую границу окна, ударного сжатия и выхода вещества из расчетной области с последующим удалением вещества на левой границе окна. Детали алгоритма обсуждаются ниже на примере моделирования методом сглаженных частиц. При моделировании CSPH методом стационарной УВ в окне наблюдения происходит создание потока новых частиц на входе в окно и удаление такого же потока частиц при выходе из окна. Таким образом, расчет производится с конечным количеством частиц и может производиться неограниченно долго, в отличие от расчетов в неподвижной системе координат, где требуется включать в расчет все новые и новые частицы, переходящие через ударный фронт в связи с пробегом УВ по образцу. Помимо этого очевидного преимущества, моделирование в подвижном окне дает возможность расчета с более мелким размером SPH-частиц, что особенно важно для УВ в пористых материалах, так как ширина фронта волны в пористом материале значительно больше, чем в сплошном веществе, а из-за колебаний давления при коллапсе пор необходима постановка зоны выхода вещества и границы удаления на значительных расстояниях. Относительная величина колебаний давления во фронте УВ растет с увеличением пористости образца и уменьшается с ростом скорости УВ, что необходимо учитывать при моделировании. Поэтому моделирование первоначально проводится для сильной УВ, а затем полученное течение используется для получения УВ с меньшим давлением путем перезапуска расчета с преобразованными скоростями частиц. Программная реализация метода подвижного окна выполнена в виде модуля для параллельного программного комплекса CSPH&VD³, разрабатываемого во ФГУП “ВНИИА им. Н.Л. Духова” авторами [15].

В начальный момент времени расчетная область подготавливается в виде прямоугольного параллелепипеда. Область плотно заполняется SPH-частицами по узлам гранецентрированной кубической (ГЦК) решетки. После этого запускается алгоритм “встряхивания” частиц [19] для достижения упаковки, близкой к жидкой, в которой дальний кристаллический порядок отсутствует. Подготовка “жидкой” упаковки состоит в воздействии на частицы случайных объемных сил, способствующих встряхиванию частиц, с постепенным затуханием их амплитуды, а также во введении сил трения для полной остановки частиц. Затем в образце вырезаются поры требуемого размера.

РАСЧЕТНАЯ ОБЛАСТЬ И ГРАНИЧНЫЕ УСЛОВИЯ В ПОДВИЖНОМ ОКНЕ

Поперечные размеры образца составляют L_y × × L_z = 20 × 20 мкм², вдоль этих направлений заданы периодические граничные условия. Длина образца вдоль оси X составляла L_x ≈ 0.57 мм. Начальный размер SPH-частиц равен 1/$\sqrt 2 $ мкм, что соответствует 20 периодам ГЦК-решетки в поперечном направлении расчетной области. Сферические поры создаются в узлах другой ГЦК-решетки, причем ось X является направлением [100] решетки. В одной кубической ячейке размерами 20 × 20 × 20 мкм³ располагаются четыре поры с координатами центров (–5, 5, 5), (–5, –5, 5), (‒5, 5, –5), (–5, –5, –5) в микрометрах, которые отсчитываются от центра ячейки при выборе осей вдоль ее ребер. Решетка пор имеет период, равный поперечному размеру расчетной области, а радиус пор составляет r_p = 4.9854 мкм и соответствует пористости $m = {{L_{z}^{3}} \mathord{\left/ {\vphantom {{L_{z}^{3}} {\left( {L_{z}^{3} - 4 \times {{4\pi r_{p}^{3}} \mathord{\left/ {\vphantom {{4\pi r_{p}^{3}} 3}} \right. \kern-0em} 3}} \right)}}} \right. \kern-0em} {\left( {L_{z}^{3} - 4 \times {{4\pi r_{p}^{3}} \mathord{\left/ {\vphantom {{4\pi r_{p}^{3}} 3}} \right. \kern-0em} 3}} \right)}} = 1.35.$

SPH-частицы удаляются после пересечения плоскости с фиксированной координатой x = x_l левой границы окна. Приграничными считаются частицы внутри расчетной области вблизи плоскости удаления, а именно, частицы с координатами x < x_l + d₀, где d₀ – ширина зоны выхода. В этой зоне вместо системы уравнений гидродинамики решается лишь уравнение движения для координат приграничных частиц, при этом их скорости одинаковы и равны отрицательной скорости выхода u_out → u_p – u_s, а поперечные компоненты скорости полагаются равными нулю. Все остальные, а именно, плотность, энергия и другие термодинамические величины, остаются неизменными и равными последним своим значениям перед попаданием в зону выхода. Здесь и далее сохраняются для скорости поршня u_p и скорости УВ u_s их обычные положительные значения, измеренные в системе координат неподвижного несжатого вещества. Для создания ударно-волнового сжатия в веществе необходим механизм передачи давления от частиц изнутри приграничной зоны к частицам вне ее, который реализуется посредством их взаимодействия, что естественным образом создает УВ при соответствующей разнице массовых скоростей, давления и плотности между этими частицами. При этом межчастичным взаимодействиям запрещено менять состояние приграничных частиц.

Граница подачи вещества имеет переменное положение, которое может варьироваться в ходе расчета на величину в несколько слоев SPH-частиц. Это связано с программной реализацией CSPH&VD³, где вставка нового слоя возможна лишь в момент обновления списка соседних частиц. Обновление списка соседей происходит через небольшое число шагов интегрирования Δt_b (около 10), но переменный шаг по времени в гидродинамическом моделировании приводит к тому, что правая граница, двигаясь с постоянной скоростью u_s, проходит разное расстояние за одинаковое число шагов интегрирования. Положение границы равномерно смещается в отрицательном направлении оси X, но скачком изменяет свое положение в положительном направлении при вставке новых слоев SPH-частиц. Поэтому граница вблизи подачи вещества характеризуется дискретностью своего положения. Для такой послойной подачи подготовка малого образца с порами производится до начала моделирования. При этом используется образец с поперечными размерами L_y × L_z, равными поперечному размеру исходного образца, и такого же размера вдоль оси X. В данном образце укладывается ровно один период вышеописанной ГЦК-решетки сферических пор, а подготовка образца проводится подобно основному образцу. Из данного образца формируется структура из набора связных списков. Каждый связный список хранит в себе положения SPH-частиц, находящихся в одном слое ГЦК-решетки вдоль оси X. Таким образом, в момент подачи вещества определяются номера связных списков из данного набора и происходит последовательная вставка всех частиц из этих слоев. Вдоль оси X задано периодическое граничное условие, т.е. после последнего слоя подается первый слой.

На схеме рис. 3 показана зона входа частиц в расчетную область. Новые наборы слоев частиц изображены темно-красным цветом, синим – частицы, уже участвующие в расчете, а оставшиеся красные частицы будут вставляться по мере освобождения места вблизи правой границы образца. Положение вставляемого слоя учитывает положение последнего вставленного слоя. Когда вставка из вещества доходит до конца красного образца, то слои начинают вставляться, начиная с первого, т.е. процесс периодически повторяется. Частицы в зоне входа вещества имеют координаты x > x_r – d_i, где d_i – ширина зоны входа, x_r – правая граница расчетной области. Границы зоны входа показаны синими пунктирными линиями на рис. 3. Частицы внутри зоны входа движутся равномерно и прямолинейно вплоть до попадания в зону расчета межчастичного взаимодействия, при этом массовая скорость всех частиц направлена вдоль оси X и равна скорости входа u_in = −u_s. Отметим еще раз, что скорость входа является фиксированной на протяжении всего времени моделирования и равна скорости УВ в стационарном режиме, а скорость выхода является переменной контролируемой величиной, которая требует подстройки для достижения стационарности ударного фронта. Значение скорости выхода изменяется таким образом, чтобы стабилизировать ударный фронт в желаемом положении внутри подвижного окна.

ПОЛУЧЕНИЕ СТАЦИОНАРНОГО УДАРНОГО ФРОНТА

Пусть в начальный момент времени в окне задан образец с определенной мезоскопической структурой, имеющий среднюю начальную плотность ρ₀₀ и состоящий из N₀ частиц. Если в системе распространяется УВ и сжатое вещество обладает большей плотностью, то при достижении стационарного режима следует ожидать увеличения полного числа SPH-частиц N_f относительно начального числа, т.е. N_f > N₀. Условие стационарности требует такого положения ударного фронта, при котором количество частиц равно N_f и ударный фронт неподвижен. Однако в случае высокой пористости возможен аномальный ход ударной адиабаты, когда плотность за ударным фронтом уменьшается, и тогда N_f < N₀. В этом случае вместо числа частиц предпочтительнее иметь дело со средним давлением в окне наблюдения

которое в УВ всегда больше начального давления P₀ в несжатом образце. Таким образом, можно определить целевую функцию, которая однозначно задает положение ударного фронта в зависимости от числа частиц N:

Если целевой величиной является среднее давление, то

Опишем итерационный алгоритм обратной связи для подбора массовой скорости u_out на выходе из окна, а также технику проведения такого расчета на примере сплошного образца. Цель этого алгоритма состоит в нахождении неизвестной скорости u_p, при которой ударный фронт стабилизируется при u_out = u_p – u_s.

Для этого будем управлять скоростью выхода в зависимости от числа частиц или среднего давления с помощью простой управляющей функции

Выбранная зависимость является убывающей функцией, что необходимо, чтобы гарантировать сходимость алгоритма к u_out = u_p – u_s при некотором ω *. Чтобы в стационарном режиме получить ω = 1, заранее неизвестный подгоночный параметр скорости β должен стать равным по модулю скорости выхода. С целью достижения β → –(u_p – u_s) введем корректирующее уравнение для изменения параметра β за одну итерацию:

где i – номер итерации, α – управляющий параметр, а также используются скорости наблюдаемого изменения целевой функции Δω_i = (ω_i − ω_i – 1)/Δt за время Δt и “желаемой” скорости сходимости $\Delta \omega _{i}^{{\text{w}}}$ = (1 − ω_i)/τ. Таким образом, если скорость сходимости, наблюдаемая в моделировании, отличается от скорости, желательной для того, чтобы решение сходилось к ω = 1 за характерное время τ, параметр β будет скорректирован в нужном направлении.

В дополнение к описанному алгоритму для ускорения выхода ударного фронта на стационарное положение используется объемная сила для всех частиц с давлением выше некоторого порогового P > P *, пропорциональная отклонению массовой скорости SPH-частицы от u_out. Данная сила отключается автоматически, если система приближается к стационарному положению |ω − 1| < ε (ε = 0.03), однако при отклонении от стационарного значения включается режим “грубой” настройки скорости выхода каждый раз, когда система оказывается достаточно отклонена от стационарного состояния. На рис. 4 представлен пример сходимости скорости выхода на стационарное значение u_out = −2.34 км/с. При этом уровень шума скорости выхода в стационарном режиме составляет приблизительно 0.05%. Поэтому ошибка расчета давления УВ также равна 0.05% при фиксированной начальной плотности и скорости УВ, так как $P = {{\rho }_{{00}}}{{u}_{s}}{{u}_{p}}.$

СРАВНЕНИЕ ВОЛНОВЫХ ПРОФИЛЕЙ, ПОЛУЧЕННЫХ В ПОДВИЖНОМ ОКНЕ И В СИСТЕМЕ НЕПОДВИЖНОГО ЖЕСТКОГО ПОРШНЯ

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

На рис. 5 красными линиями представлены профили плотности, массовой скорости, давления и девиаторов напряжений, полученные в подвижной системе координат. Скорость ударной волны u_s = 3.5 км/с равна (с противоположным знаком) скорости входа вещества u_s = –u_in в расчетную область. В результате моделирования получена стационарная скорость выхода u_out = = −2.34 км/с, откуда вычисляется массовая скорость за фронтом ударной волны u_p = u_s + u_out = = 1.16 км/с относительно покоящегося несжатого образца.

Для получения эквивалентной ударной волны в системе координат неподвижной жесткой стенки (рис. 1) необходимо произвести моделирование столкновения образца меди с начальной скоростью вдоль оси X U_x = −u_p с неподвижной жесткой стенкой. В результате столкновения вещество останавливается и образуется УВ, распространяющаяся со скоростью U_s = u_s – u_p. На рис. 6 черные пунктирные линии показывают соответствующие величины плотности и напряжений, а также профиль преобразованной массовой скорости, полученный из исходного профиля скорости вычитанием скорости УВ в данной системе координат ${{\dot {U}}_{x}}$ = U_x – U_s. Наблюдается согласие в полученных профилях, что означает эквивалентность в возможности задания граничного условия для создания ударных волн в системе координат неподвижного ударного фронта.

Отметим, что плотность SPH-частиц имеет периодический характер, что связано с явным заданием мезоскопической структуры и периодичностью геометрии пор. При этом сохраняющиеся осцилляции плотности за областью релаксации давления являются результатом неоднородной плотности SPH-частиц, вызванной неоднородным прогревом частиц, вовлеченных в процессе схлопывания пор. Показанные на рис. 5 профили подобраны таким образом, чтобы распространяющаяся УВ была в одинаковой фазе в обеих системах отсчета, что позволило сопоставить профили явно.

РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ УДАРНЫХ ВОЛН В ПОРИСТОЙ МЕДИ

Гидродинамическая постановка. Моделирование ударного нагружения пористой меди производится в приближении сжимаемой жидкости для SPH-частиц меди с использованием модели 1 и метода обращенного движения в системе неподвижной жесткой стенки. Исходный образец готовится с явным заданием мезоскопической структуры пор с размерами, соответствующими заданной пористости. Для моделирования распространения ударной волны в пористой меди использовалась структура, состоящая из последовательности кубов с парой пор, имеющих центр на главной диагонали, как описано в [22]. Размеры образцов составляли 1036 × 80 × 80 мкм³.

Результаты проведенного моделирования для меди с различной пористостью хорошо согласуются с экспериментальными ударными адиабатами в имеющемся диапазоне экспериментальных данных. Также проведено моделирование ударного нагружения пористой структуры при низких скоростях воздействия и начальной плотности образца 5.7 г/см³.

Ударная адиабата сплошного вещества и ударная адиабата, полученная в CSPH-расчетах, представлены на рис. 6a. Результаты масштабирования при пропорциональном увеличении линейных размеров системы без изменения числа SPH-частиц показаны на рис. 6б. На этом же рисунке показан результат моделирования с более мелким размером SPH-частиц при фиксированной геометрии пор, что показывает сходимость полученных результатов по размеру SPH-частиц. На рис. 6а представлено сравнение с экспериментальными данными из [21], где черными линиями показана линейная аппроксимация полученных в моделировании ударных адиабат при различной пористости. Представленные данные находятся в пределах погрешности экспериментальных ударных адиабат в рассматриваемом диапазоне. Однако для моделирования сильных УВ в материалах с большей пористостью необходимо применять модель сплошного вещества для более высоких плотностей энергии [23].

В моделировании УВ с начальной плотностью 5.7 г/см³ в области малых значений амплитуды УВ наблюдается отклонение от линейной зависимости u_s от u_p, которая характерна для высоких скоростей удара u_p > 750 м/с. Данная часть кривой может быть приближенно описана процессом полного схлопывания вакуумированных пор при пренебрежимо малом скачке давления. Если допустить, что коллапс пор происходит без сопротивления и нагрева, то плотность вещества за ударной волной можно оценить величиной плотности сплошной меди ρ₀. Тогда закон сохранения массы на ударном скачке записывается в виде

где ρ₀₀ – начальная плотность пористого материала, ρ₁ – плотность вещества в ударной волне. В приближении ρ₁ ≈ ρ₀ [24, 25]:

Оценка с использованием уравнения (2) показана пунктирной линией на рис. 5a для исходной плотности 5.7 г/см³. При скорости удара u_p < 200 м/с скорость ударной волны выходит на постоянное значение, при этом из-за низких амплитуд давления сжатия повышается относительная амплитуда осцилляций давления на фронте, в результате чего снижается точность определения скорости фронта и давления на фронте.

Упругопластическая постановка. При экспериментальном определении отклика материала на ударное сжатие измеряется напряжение в сжатом материале, которое зависит от режима упругопластического деформирования. В мезоскопическом моделировании пористого материала предлагается учесть прочностные свойства каркаса сплошного материала, которые дают отдельный вклад в наблюдаемое напряженное состояние. Следует отметить, что данный эффект особенно важен для определения ударной адиабаты при такой амплитуде воздействия на материал, когда давление в УВ сравнимо с девиаторной частью напряжений. В используемой упругопластической модели величина девиаторов напряжений ограничена условием на цилиндре Мизеса, радиус которого по порядку величины равен пределу текучести Y, что может служить оценкой амплитуд ударных волн, где ударная адиабата пористого материала определяется прочностными эффектами в сплошной меди. Рассмотренный диапазон напряжений ударного сжатия имеет пластический характер. Получены стационарные ударные волны в пористой меди с использованием упругопластической модели вещества для SPH-частиц. Использовалось подвижное окно с вышеописанной методикой установления стационарного ударного фронта. Моделирование проводилось, начиная со скорости ударной волны u_s = –u_in = 5 км/с. Выбор функции скорости выхода u_out частиц из окна позволяет изменить темп перехода УВ на стационарный режим. В начальный период времени по среде должна пойти УВ большей амплитуды, чем ожидается в стационарном режиме, чтобы иметь и большую скорость распространения, чем скорость подачи вещества. Начальное значение скорости выхода u_out для старта расчета (при t = 0) оценено из известных ударных адиабат для пористой меди. В частности, использовалась ударная адиабата меди с начальной плотностью ρ₀₀ = = 6.3 г/см³. Плотность сплошной меди для используемого уравнения состояния ρ₀ = 8.96 г/см³, что для пористости m = 1.35 составляет ρ₀₀ ≈ ≈ 6.637 г/см³. Полагая линейной поправку на изменение массовой скорости u_p в УВ при малом изменении пористости для фиксированной скорости ударной волны u_s, можно положить u_p(ρ₀₀ = = 6.637) ≈ 6.3/6.637u_p(ρ₀₀ = 6.3). Полученные профили скорости и давления имеют осциллирующий характер за фронтом стационарной УВ (рис. 7), причем колебания находятся практически в противофазе, что объясняется адиабатическим изменением внутренней энергии волнами сжатия и разрежения, возникших после коллапса пор за фронтом УВ. Приход ударного фронта к поре приводит к высокой скорости схлопывания пор, что может вызвать кумулятивный эффект и ускорение струй вещества, которые тормозятся при ударе о внутреннюю поверхность поры, вызывая пики давления в профиле. Стоит отметить, что амплитуда колебаний затухает с расстоянием от фронта, и этому способствуют различные механизмы диссипации механической энергии в тепло: процессы многократного сжатия и разгрузки и схемные эффекты метода SPH. Ширина окна выбиралась достаточно большой, чтобы исключить воздействие выходного граничного условия на процесс затухания колебаний давления, и гарантирует установление стационарного режима распространения УВ.

ЗАКЛЮЧЕНИЕ

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

Исследование выполнено за счет гранта Российского научного фонда (проект № 19-19-00697).