Численный метод SPH, использующий соотношения распада разрывов, и его применение в механике деформируемых гетерогенных сред тема автореферата и диссертации по механике, 01.02.04 ВАК РФ

Паршиков Анатолий Николаевич

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

Специальность 01.02.04 - механика деформируемого твёрдого тела

Автореферат диссертации на соискание ученой степени доктора физико-математических наук

6 ЛЕЗ 2014

Москва-2013

005544706

005544706

Работа выполнена в ФГБУН Объединенном институте высоких температур Российской академии наук.

Официальные оппоненты: доктор физико-математических наук, член-

корреспондент РАН,

зав. кафедрой вычислительной математики МФТИ Петров Игорь Борисович

доктор физико-математических наук, профессор, зав. кафедрой вычислительной математики и математической физики МГТУ им. Н.Э.Баумана Димитриенко Юрий Иванович

доктор физико-математических наук, профессор, зав. отделом ИПМ им. М.В.Келдыша РАН Гасилов Владимир Анатольевич

Ведущая организация: Центральный научно-исследовательский институт

машиностроения (ФГУП ЦНИИмаш)

Защита состоится 20 марта 2014 г. в 15 часов на заседании диссертационного совета Д 002.240.01 при Институте проблем механики им. А.Ю.Ишлинского РАН по адресу: 119526 Москва, проспект Вернадского, дом 101 корп.1, аудитория № 237.

С диссертацией можно ознакомиться в библиотеке Института проблем механики им. А.Ю.Ишлинского РАН.

Автореферат разослан « /?/» О/ 201 Цт.

Ученый секретарь

диссертационного совета Д 002.240.01

к. ф.-м.н. // Е.Я.Сысоева

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Численные методы решения уравнений динамики сплошных сред являются практически единственным инструментом для исследования процессов, происходящих в структуре гетерогенных сред при ударно-волновом нагружении, так как аналитического решения подобные задачи, как правило, не имеют. Перечень неоднородных сред, подвергаемых ударно-волновому воздействию, достаточно обширен - это керамики, материалы порошковой металлургии, пространственно-армированные композиционные материалы, компоненты конструкции перспективных энергетических реакторов, вспененные и волокнистые материалы, смеси взрывчатых веществ с инертными добавками. К неоднородным (пористым) средам можно отнести также и гомогенные материалы, находящиеся в режиме объемного вскипания.

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

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

Диссертация содержит описание разработанного автором численного метода «гладких частиц», или SPH (Smooth Particle Hydrodynamics), использующего соотношения распада разрывов, и результаты применения этого метода к решению ряда задач механики гетерогенных сред. Показано успешное применение разработанного метода SPH к моделированию ударного воздействия на следующие материалы:

— на пористый однокомпонентный материал;

— на двухкомпонентный материал без пор;

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

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

Актуальность темы. Частным случаем гетерогенных сред являются пористые металлы. Интерес к ударному сжатию пористых металлов был изначально вызван потребностью в данных о термодинамических свойствах вещества при высоких давлениях и температурах [1]. Различные металлы были протестированы в ГПа-ТПа диапазонах давлений [2] и получены пористые адиабаты Гюгонио. В результате были построены и верифицированы широкодиапазонные уравнения состояния. Значительное внимание к пористым материалам проявляется в связи с защитой от высокоскоростного удара [3]-[5]. При этом представляет интерес выявить основные физические механизмы поглощения энергии в пористых материалах. Зельдович и Райзер [1] показали, что схлопывание пор отвечает за повышенную диссипацию энергии при сжатии пористого материала. В первых экспериментальных исследованиях Боуда [6]-[7] было установлено, что при слабой динамической нагрузке уплотнение пористых металлов происходит в многоволновом режиме, с одной или несколькими волнами-предвестниками.

Одним из важных применений мезомеханического моделирования к ударно-волновому нагружению пористой среды является возможность построить расчётным путём адиабату Гюгонио в области неполного схлопывания пор, по известным ударным адиабатам сплошных веществ, которые измерены экспериментально для большого числа веществ и материалов [8]. В области высоких давлений, когда изначально пористый материал становится сплошным, расчёт ударной адиабаты достаточно проработан теоретически, но в области неполного схлопывания пор результат достигнут с помощью численного мезомеханического моделирования [9].

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

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

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

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

искажениям сетки.

Наиболее приемлемым для решения задач мезомеханики, в которых рассматриваются гидродинамические процессы на масштабе структуры среды, является метод «гладких частиц» SPH (Smooted Particle Hydrodynamics) [11]. Метод SPH интенсивно применяется для решения многомерных задач гидродинамики, и к его бесспорным техническим достоинствам следует отнести высокую простоту алгоритмической реализации при минимальном размере программного кода. Отсутствие расчетной сетки позволяет методу (в

рамках лагранжева формализма) естественным образом рассчитывать произвольные вращательные и сдвиговые течения, распад односвязных и слияние многосвязных расчетных областей. Более того, свободно-лагранжевы методы дают физически правильную картину эволюции течения в тех случаях, когда применение сеточных лагранжевых методов становится в принципе невозможным вследствие неприемлемых искажений расчётной сетки. Недостатком метода БРН, использующего искусственную вязкость (он называется ниже «стандартным методом»), является погрешность расчёта в окрестности контактных разрывов плотности и границ раздела компонент.

Таким образом, для задач мезомеханического моделирования ударно-волновых процессов в гетерогенных средах актуальной проблемой является повышение точности метода БРН в окрестности контактных границ.

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

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

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

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

Практическая значимость работы. Метод 8РН, использующий соотношения распада разрывов, реализован в виде свободно распространяемого комплекта компьютерных программ для ЭВМ, с интерфейсом пользователя, на языке ФОРТРАН-90. Этот комплект может применяться:

- в учебных и исследовательских целях, для изучения отклика гетерогенных сред на ударно-волновое нагружение;

- для решения инженерных задач по оценке стойкости броневой защиты и стёкол к ударному воздействию;

- при проектировании низкоплотных взрывчатых веществ с заданными параметрами, необходимыми для решения технологических задач штамповки взрывом.

Результаты исследований обобщены в виде следующих положений, выносимых на защиту:

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

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

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

4. Разработан мезомеханический подход к моделированию гетерогенных сред с помощью метода SPH, позволяющий путём вычислительного эксперимента рассчитать ударную адиабату пористого материала в области неполного схлопывания пор, если известна ударная адиабата сплошного вещества, составляющего пространственную мезоструктуру данного материала.

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

Достоверность результатов подтверждается их соответствием аналитическим решениям и экспериментальным данным.

Личный вклад автора. Автором впервые предложен и создан вариант численного метода SPH, основанный на решении задач распада разрывов. Все представленные результаты получены автором лично.

Апробация работы. Основные результаты исследований были доложены и обсуждались на:

- международной конференции по гиперскоростному удару. HyperVelocity Impact Symposium (Oct. 17-19, 1994, Santa Fe, New Mexico, USA);

- 24th International Symposium of Shock Waves (Beijing, China, July 11-16, 2004);

- международной конференции по гиперскоростному удару HyperVelocity Impact Symposium (Oct. 10-14,2005, Lake Tahoe, USA);

- международной конференции «Физика экстремальных состояний вещества» (Эльбрус, март 1-5, 2007);

- 2-й, 3-й, 4-й и 5-й Всероссийских школах-семинарах «Аэрофизика и физическая механика классических и квантовых систем», АФМ-2008-АФМ-2011 (Москва, 2008-2011);

- X Всероссийском съезде по фундаментальным проблемам теоретической и прикладной механики (авг. 24-30,2011, Нижний Новгород).

- международной конференции «Разностные схемы и их приложения», посвященной 90-летию профессора В.С.Рябенького, 27-31мая, 2013, Москва.

Публикации. По теме диссертации опубликовано 27 печатных работ, из них 14 в ведущих научных рецензируемых журналах, рекомендованных Перечнем ВАК РФ для публикации результатов докторских диссертаций.

Структура и объём работы. Диссертация состоит из введения, пяти глав и заключения; содержит 202 страницы, включая 85 рисунков, 13 таблиц, список литературы из 162 наименований и одно приложение.

СОДЕРЖАНИЕ РАБОТЫ

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

В первой главе рассмотрены основополагающие принципы метода БРН. Разностная сетка в методе БРИ отсутствует, сплошная среда заменяется дискретной системой расположенных в пространстве сглаженных частиц г'=1...ЛГ, допускающих произвольную связность друг с другом. Переход от континуума к дискретному представлению среды в виде сглаженных частиц предполагает, что вместо непрерывной функции /(г) , характеризующей какую-либо полевую величину (давление, энергию, плотность) вводится в рассмотрение её дискретный аналог, то есть дискретная функция /¡. Кусочно-постоянная величина /, определяется для каждой частицы г через сумму N величин £ из частиц окружения лежащих вокруг частицы I в соответствии с дистанцией сглаживания к, по уравнению:

N /,

(1)

>1 Р)

где IVу = ^(¡г, - г}|,й) есть весовая (сглаживающая) функция, или ядро; т] есть масса частицы у, ресть плотность вещества в частице у. На ядро накладывается ограничение

|ж(г,й)<й? = 1 (2)

При й—>0 сглаживающая функция превращается в 5-функцию. Согласно свойствам 5-функций, для вычисления производных от /(?) вместо непосредственного дифференцирования непрерывной функции /(?) нужно знать дискретное значение £ в частице и дифференцировать функцию Щ [12]:

яг N Л , х- - х5:

дх\; М } р] ' {П-гД

В этом заключается принцип бессеточного метода БРН. В практических применениях используют различные виды сглаживающих функций. В диссертации в качестве ядра был выбран кубический сплайн [13]:

/2 + Ъфг/а)/Ы, 0<ф<\

йф, - г;|, й) =■ ((2-фУ/ 4)/ N. \ <, ф <2 (4)

О, Ф*2

Здесь ф = \г,-^\/1г и Ы=1.5И для одномерной геометрии, N=0.7пк2 для двумерной (плоской и осесимметричной), Ы=пк3 для трёхмерной.

Производная от ядра запишется как:

IV

(-12ф + 9ф2)/^, 0<ф<1 (- 3(2 — )/ Л^!, 1<ф<2 О, ф-2:2

(5)

Здесь И^бк2 для одномерной геометрии, М]=28яЬ3 для двумерной (плоской и осесимметричной), N1=4-11}г4 для трёхмерной.

Сглаживающая функция (4) и её производная показаны на рисунке 1 а,б.

0,8-г

Рисунок 1 - Сглаживающая функция (а) и её производная (б) В главе 1 на примере уравнений для идеальной сжимаемой жидкости:

(6)

Л Р

приведены стандартные БРН-уравнения, полученные с помощью (3): <И І Р]

Л ¿Г Р,Рі "\ г,-г\

Л І 2РіРі \7і-г\

(7)

(8)

(9) (10)

(П)

Система уравнений (9>—(11) приводит к погрешности при расчёте окрестности контактного разрыва, если частицы і и ] расположены по разные его стороны и вещество в одной из них расширяется, ускоряя и сжимая частицу по другую сторону контактного разрыва. Причину этой погрешности легко уяснить на примере формулы (9). Применим (9) к расчёту изменения

10

плотности, поочерёдно к паре частиц г и у, расположенных по разные стороны контактной поверхности. Для простоты исключим из рассмотрения все остальные частицы, чтобы оценить взаимовлияние частиц г и у друг на друга. Из выражения в круглых скобках у правой части (9) следует, что частица г и частица_/ будут или сжиматься обе (сближаясь), или расширяться обе (удаляясь друг от друга). Таким образом, ударную волну или волну разрежения уравнения стандартного метода будут описывать верно. Но процесс, при котором БРН-частица / расширяется, сжимая БРН-частицу ) (или наоборот) рассчитать по уравнениям (9Н11) не Удаётся, хотя это довольно распространённый случай взаимодействия двух сред через контактную границу. Типичными задачами такого рода являются задачи детонации взрывчатых веществ, содержащих инертные примеси, задачи взаимодействия продуктов взрыва с оболочками, задачи лазерной абляции. Например, если БРН-частица г представляет собой расширяющиеся продукты детонации, они сжимают инертный материал в БРН-частице / Если частица г поглощает лазерное излучение, расширяясь, она также будет сжимать частицу у,

прозрачную для излучения.

С помощью уравнений (9)410 нельзя обеспечить монотонность распределения расчетных величин в окрестности контактного разрыва. Немонотонность решения обычно сглаживается введением искусственной вязкости. Этому вопросу (корректному подбору искусственной вязкости) уделено внимание в ряде работ; особенно следует отметить [14], где сделана попытка построить искусственную вязкость на основе решения задачи о распаде произвольного разрыва.

Рассмотрим вывод модифицированной БРН-формы уравнений из уравнений (9)-{11). Пусть объем пространства У0 вместо непрерывной среды с плотностью ро заполнен сферическими (для определенности) БРН-частицами (рисунок 2). Расположение частиц в пространстве допускает как частичное пересечение или непересечение, так и точечный контакт частиц. Каждая частица г в такой БРН-среде может обмениваться импульсом и энергией с любой из у частиц окружения в пределах максимальной дистанции взаимодействия |гу-^|)/й = 2только вдоль оси взаимодействия Я, соединяющей центры частиц. Интенсивность этого обмена определяется видом сглаживающей функции. Значение /г в диссертации определялось как:

h = a{D,+DJ))/2 (12)

где величина а выбиралась в диапазоне 1 £ а < 1.4.

Введем теперь понятие точки касания частиц. Необходимо уточнить, что только в частном случае, когда -/¡|)/(/>у +£>)=! , будет иметь место

геометрическое касание сфер диаметрами £>, и Д в точке Ау на оси взаимодействия К.

ч

и

•я

Ь;

к к

Рисунок 2 - Схема взаимодействия частиц в вРН-жидкости

В общем случае под точкой касания частиц будем подразумеваеть точку Ау, делящую отрезок |г,-г(| на пропорциональные диаметрам взаимодействующих частиц части. Предположим далее, что указанная точка касания частиц в БРН-среде эквивалентна контактной поверхности (КП) в сплошной среде.

Тогда можно построить решение Римана для распада разрыва вдоль оси Л и вычислить скорость точки касания и\ и давление в ней Р*. Для простоты ограничимся акустическим приближением:

и*р/с'1+и*рр11-р/ +р, (13)

" " Рр\ + Рр\

р£>+РР1 ■

где С' есть скорость звука.

Скорости и1 , ¡7/ и давления Р, и Р} будут стремиться к распадным значениям, вычисленным из (1.20)-(1.21). Выполним в уравнениях (9Н11) замену

(с/,* + и«)/2 -» и;; ; (/> + Р, >2 -> я;

(15)

Тогда уравнения (9)—<11) преобразуются к иной SPH-форме и запишутся

как:

= zHUELrfr*-щ*) (16)

да J Pj

f = [ (I7)

л У PtPj Ь~Г>\

^ = (18) àt j ¿PiP]

Уравнение энергии предпочтительнее записать в форме уравнения сохранения полной энергии

£Î£+£ll = _Iv.(«?) (19)

dt{ 2) р

и преобразовать к виду:

d(Ei+\/2Uf) = 21mjp;u'f к (20)

dt j 2 p^j 4

Если перейти от уравнений (16),(17),(20), описывающих изотропную жидкость, к уравнениям прочной сжимаемой среды, то в характер взаимодействия частиц, показанный на рисунке 2, необходимо внести некоторые дополнения. Изотропные SPH-частицы обмениваются импульсом и энергией только в волнах сжатия (разрежения) вдоль оси R. Этот достаточно простой вид взаимодействия привёл к лаконичным по форме SPH-уравнениям.

В прочной сжимаемой среде такой обмен происходит как в волнах сжатия (разрежения), распространяющихся с продольной скоростью звука

d =V(3*+4G)/(3p) (21)

так и в волнах сдвига, распространяющихся с поперечной скоростью звука

С' =4ô7p (22)

Для разъяснения особенностей такого взаимодействия в прочной SPH-среде, позволяющих построить решения Римана для двух взаимодействующих частиц i и j, обратимся к рисунку 3, где и показаны эти частицы. Если через точку касания частиц i и j (обозначенную как точка А0) перпендикулярно оси взаимодействия O-R провести плоскость, она пересечет координатные оси системы XYZ в точках abc. Назовем плоскость abc плоскостью касания частиц. Участок этой плоскости в окрестности точки Aij и будет эквивалентен поверхности начального разрыва параметров в реальной среде. Разрыв распадается в общем случае на продольную волну и волну сдвига, входящие в

частицу i и две аналогичные волны, входящие в частицу у. В процессе распада разрыва напряженное состояние среды на площадке abc характеризуется вектором напряжений ë'R, имеющим нормальную к abc компоненту сг1у и касательные компоненты a*fR и <т*га, лежащие в плоскости abc. Системы координат RST и XYZ связаны через углы 9 и ср.

Рисунок 3 - Схема взаимодействия частиц в упругопластической вРН-среде

Матрица направляющих косинусов частицы./ относительно частицы і с использованием как традиционных, так и принятых на рисунке 1.3 обозначений записывается в виде:

"cos*1 cos"' cos*-"" cos ç sin в sin q> sin в COS0

cos& cos*' cosS: = l*> ls: = COSÇJCOS0 sin q> cos в -sinö (23)

cos7* cosTy cos7* p ¡ту -sin <p cosç 0

Матрица А обеспечивает переход из системы координат XYZ в систему RST. Обратный переход из системы RST в систему XYZ обеспечивает транспонированная матрица направляющих косинусов А .

Распадные значения компонент напряжений и скоростей легко вычисляются в системе координат ÄST в акустическом приближении: для поперечной волны

.5« аГрд+оГр.С'+рС'р.Сии3, -и?) ,25.

Р]С] +Р,Ч

,г и,р1С]+и1 р,с, -в, (26)

II ~ „ п> л. „ г' '

.„ <у™рР\ +*™р,с\. + Рр\Рр\<р*-и!)

для продольной волны

(27)

УЧрЛ+иГрд+оГ-аГ (28)

™ + +Р,с',р1с!(^-и?) (29,

= ■ ( }

Для прочной сжимаемой среды уравнение неразрывности (16) не изменится; в уравнениях (17) и (20) вместо давления Р* следует использовать тензор напряжений

а*<* =-Р'8а? (30)

где - девиатор напряжений в точке контакта и

8ч> =|1 пРи а = Р

[0 при а ф р

есть символ Кронекера и а=х,у,2 и (3= х,у,г.

Тогда уравнения (17) и (20) с учётом записи векторов в системе Л5Г примут вид:

(з»

л /АЛ 4^+1/2^)^¿-«и;; [г, (з2)

Л "7 Ару

Уравнения (16),(31),(32) есть модифицированные уравнения БРИ, позволяющие построить монотонную схему расчета без искусственной вязкости. В диссертации приведен полный алгоритм для расчёта прочной сжимаемой среды, включая пересчёт напряжений из системы координат ХН. в систему РБТ и обратно.

Если уравнение энергии ограничить рассмотрением только механизма теплопроводности, то оно примет вид

— = -сІічд (33)

Л

где есть вектор теплового потока. БРН-аппроксимация (1.41) запишется как:

(34)

Л і РІР]

БРН-аппроксимация (34), полученная из уравнения энергии (33), соответствует закону Фурье. Для вычисления потоков ¿7, и в [15] используется температура контакта частиц Ту.

(34)

Лг,

. Т, — Т*,

ъ

где X есть теплопроводность.

Температура контакта Т'ц и Лг, и Лг] определяются в [15] как

^--ХІі-^і. (35)

' Аг,

„я _ „я Я І =ЧІ

(36)

(37)

В этом случае Т* можно найти как

, _ Я,Г, + Я/у {38)

9~ А,+ЛУ

В соответствии с приведенными выше соотношениями в [15] получена следующая 8РН-аппроксимация уравнения энергии для случая только

механизма теплопроводности:

^ = (39)

а ! р{р, Ц + у}-г,\

В работе [15] проведено тестирование уравнения (39) и получено соответствие с аналитическим решением для контакта различных пар материалов. Интегрирование уравнения (39) должно производиться с шагом

Л1п = /ЗрСук1 /Я (40)

где Су есть теплоёмкость и р<0.15.

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

Начальные условия следующие: Г(х,0)=Г„ х<0 и Т(х,0)=Тр х>0. Распределение температуры в этом случае определится как [16]

х

т-т„ =

(т'-ф/-^- при Х<0

2 д/а,/

при х > О

(41)

где Т* есть температура в точке контакта х=0

п _

рСу

есть температуропроводность. Комбинируя (1.50) с (1.42Н1-44)» получаем следующую форму уравнения энергии:

(42)

(43)

Л у

(44)

Следует заметить, что в случае теплового контакта частиц г и у из одного материала уравнения (39) и (44) полностью эквивалентны. Уравнение (44) точнее, как показано ниже на примере тестовых расчётов, определяет профиль температуры при контакте частиц из различных материалов.

Для описания упругопластического течения используется процедура Уилкинса [17]

= (45)

где корректирующий множитель Кр в соответствии с критерием Мизеса есть

(46)

Кр = \г0^27?

Го

/>2702

и / = 35^5^ . Уравнение состояния, замыкающее систему (16),(31),(32), принимается в форме Ми-Грюнайзена

Р-Рг=УзР(Е-Ег) (47)

здесь Рг(р) И Ег( Р) - опорные кривые.

Для инертных материалов опорными кривыми являются ударные адиабаты и, =Са+5аир при плотности выше начальной и упругие кривые при плотности ниже начальной

\Рн'р > ро \Е"'Р > р0

где = о - »)/[»„ - sa(v0 - v)]2 , Ен = P„(v0 - v)/2 ,

Pc=K{v0-v)/v0> Ec=Pc{v0-v)/2, v = l/p

Уравнение состояния в форме (48), использующее в качестве опорной кривой ударную адиабату для нагружаемого материала, выбрано из тех соображений, что коэффициенты Са и Sa ударной адиабаты вида Us=Ca+SaUp измерены достаточно точно для большинства материалов. Они многократно проверены независимыми исследователями для различных материалов при давлениях до ¡=20-500ГПа, что с избытком перекрывает диапазон давлений, реализуемых при ударном нагружения алюминия с неполным схлопыванием пор, а также диапазон давлений, реализуемых при детонации смеси насыпной плотности из взрывчатого вещества с инертной добавкой (решения этих задач представлены в главах 3 и 5 диссертации). Не менее важным аргументом в пользу уравнения (48) является и то, что решение задачи Римана с двучленным уравнением состояния позволяет учитывать свойства среды достаточно общего вида.

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

Во второй главе приведены тестовые расчеты и сравнение

разработанного метода со стандартным методом SPH, использующим

искусственную вязкость. Для расчёта распада разрыва в упругопластической

среде уравнение состояния принималось в виде

р = к(р-р0)/р (49)

где изотермический модуль объёмного сжатия K=const. Расчёт проводился в одномерном приближении. Начальный разрыв параметров напряженного состояния среды в точке x/L=0.5 есть о? = 4 ГПа и =0, слева и справа от разрыва соответственно. Во всей области р,=р2=2700кг/м3. Это соответствует алюминию, и в расчёте принимались следующие его свойства: К=73 ГПа, G=23 ГПа, У0=0.3 ГПа. Длина расчётного интервала 1=0.1 м была разбита на 200 SPH-частиц. На рисунке 4 показаны результаты расчётов для момента времени 5 мкс. Сплошной линией показано аналитическое решение. Скорости и амплитуды

упругих и пластических волн, полученные в расчете, практически совпадают с аналитическим решением [1].

0,0 0,2 0,4 0,6 х/Ь

Рисунок 4 - Распад разрыва в упругопластической среде. Сплошной линией показано

аналитическое решение

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

0,20-г

Рисунок 5 - Распад разрыва в газе по модифицированному методу вРН. Скорость (а), давление (б), плотность (в) и внутренняя энергия (г) относительно расстояния. Сплошная линия показывает аналитическое решение

Начальные параметры газа задаются на момент времени i=0. Исходное положение контактной поверхности x/L=0.5. Расчетная область размером Z,=0.1m содержит 200 расчетных частиц равных размеров. Слева от контактной поверхности: Л =3x104 Па, р,=1500 кг/м3, £/,= 0. Справа от контактной поверхности: Р2= 104Па, р2=1200 кг/м3, U2=0. Показатель адиабаты газа yi=r2=3. Уравнения (1.24)-(1.26) были дополнены уравнением состояния идеального газа На рисунке 5 проведено сравнение расчетных профилей U(x), Е(х), Р(х), р(х) с

точным решением.

На рисунке 6 показаны для сравнения результаты решения этой же задачи стандартным методом SPH с помощью уравнений (1.15)-(1.17), и с применением искусственной вязкости

q, = 20 Dfp,s\s\ + 0.5 D^Cje (50)

Рисунок 6 - Распад разрыва в газе по стандартному методу вРН. Скорость (а),

давление (б), плотность (в) и внутренняя энергия (г) относительно расстояния.

Сплошная линия показывает аналитическое решение

При решении задачи стандартным методом в окрестности контактного разрыва наблюдается значительная погрешность вычисления параметров течения газа, природа которой объяснялась ранее.

Для сравнения разработанного метода (44) со стандартным (39) была решена задача теплового контакта фарфора с горячим газом высокой

20

температуропроводности. Начальные условия: Г(х,0)=Г1=300К, х/Ь< 0.1 (фарфор) и 71(л:,0)=Г2=10000К, х/Ь>0.\ (горячий газ). На рисунке 7а показан профиль температуры в момент времени 4х10"9с для области 1=0.1мм, рассчитанный обоими методами.

На рисунке 76 сравниваются относительные ошибки обоих методов с профилем температуры Та, полученным из аналитического решения. Относительная ошибка в расчете температуры горячего газа вблизи стенки на порядок ниже для разработанного метода, чем для стандартного метода].

Рисунок 7 - Профили температуры и относительная ошибка вычисления температуры с использованием уравнения (44) (о)и (39) (+)

В разделах 2.3-2.7 диссертации представлены результаты тестирования кода БРН на ряде задач, а именно: расчёт взрывной волны, расчёт сдвигового течения в жидкости, расчёт соударения резиновых цилиндров, расчёт вращения упругой пластины, расчёт разрушения стёкол в волне разрушения. Целью последнего теста являлась проверка, насколько правильно расчётная схема описывает процесс разрушения предварительно сжатого ударной волной материала (стекла). Разрушение материала происходило в волне разрушения, распространяющейся от поверхности нагруженного ударом образца.

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

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

самоподдерживающееся разрушение при сжатии. Степень разрушения материала описывается параметром разрушения £) (0<0<\) . Для неразрушенного материала />=0, для полностью разрушенного материала П= 1. Следуя феноменологическим подходам работ [18-20], предложенным для определения эволюции параметра разрушения в нагруженном материале, считалось, что разрушение материала представляет собой самораспространяющийся волновой процесс в теле, в котором до приложения нагрузки уже имелись локализованные зоны или поверхности поврежденного материала с заданным начальным распределением параметра И. Волны разрушения распространяются от поврежденных участков тела после того, как в них достигаются пороговые значения критериев нагружения. Скорость волны разрушения С/ при этом полагается заданной характеристикой материала, определяемой в эксперименте. Скорость накопления повреждения в лагранжевой БРН-частице принимается пропорциональной модулю градиента параметра разрушения О. Нелинейное волновое уравнения для параметра разрушения И записывается в виде

^сШШ (51)

л Ч^,

Предельное состояние материала в процессе разрушения описывается параметрической связью между эквивалентным напряжением и давлением Р. Была решена одномерная задача о нагружении упругого полупространства, граница которого в начальный момент времени начинает двигаться с постоянной скоростью к/. На границе для параметра разрушения выполняется условие £>=1, и при достаточно высокой интенсивности упругой волны сжатия за ней будет распространяться волна разрушения. Построено аналитическое решение для такой двухволновой конфигурации. На фронте головной упругой волны, движущейся со скоростью Си выполняются условия сохранения массы и импульса

р(С,-и) = р0С, (52)

стх + ри(С/ - и) = 0 (53)

На фронте волны разрушения выполняются уравнения сохранения

р/(С/ + и-и/) = рС/ (54)

(55)

условие Друкера-Прагера для стекла и уравнение состояния (49)

-0}гл*е+рг=а + ър,, а= 1.1-109Па, 6=1.2 (56)

Р, =К^-Ра/р,) = РоС2ь (1 -ро /Р/) (57)

Индекс/в уравнениях (54)-(57) относится к параметрам разрушенного материала. Уравнения (56)-(57) указывают на то, что в рассматриваемом случае разрушенный материал ведет себя как предварительно напряженное упругое тело с модулем объемного сжатия, равным модулю объемного сжатия неповрежденного материала К{ - К, и модулем сдвига = 3К(Ь -1) / 4. Выразим <ух и ст^-ег, через р и рг\

ах = -р0С,и = -РоС}(\ -р0/р) (58)

ах/-ах= -рС/(иг -и) = -рС}{1 - р/р}) (59)

Исключая р/р{ из уравнений (58)-(59), приходим к кубическому уравнению для определения р/ рй

1 - (Р0С/ + РоСуС, - а - ЪрйС] [-I Рй) Ро\Ро)

ЬЛ{С/-и/+С1)-С1

Р , ьс2ьс, =0

(60)

Ро

Уравнение (60) указывает на то, что значение плотности р за фронтом упругой волны является функцией скорости м/, а также свойств материала и параметров критерия разрушения. Упругое состояние за фронтом головного скачка определяется в результате решения полной задачи в отличие от задачи для упругопластической волны, в которой состояние за фронтом упругого предвестника фиксировано заданием предела упругости. По вычисленному значению р/р0 определяются и, ах , р, и а4 . На рисунке 8 показаны распределения продольного (<тх , кривая 1) и поперечного ( ау, кривая 2) напряжений за фронтом ударной волны и в волне разрушения при скорости границы и/ = 1000м/с. Прямые линии - аналитический расчет, кривые -численное решение, полученное в результате расчетов по двумерному коду БРН с граничными условиями, имитирующими одномерное приближение.

В разделе 2.8 диссертации представлены результаты моделирования разрушения стеклянных пластин в плоском двумерном приближении. Решалась задача удара стеклянной пластины размерами /гх/ = 40.4мм><3.2мм о жесткую стенку с начальной скоростью м0=1000м/с (пластина движется справа налево). Стенка располагается в х=0. Условия на вертикальной жёсткой стенке задаются следующим образом: для каждой расчётной БРН-частицы г с параметрами {т,р,0^,их,иу,Р!,Т1,Зх*,Зуу,$ху}1 на каждом расчётном шаге создаётся виртуальная частица] с параметрами {т,р,й,-х,у,- их, иу,Р^ Г,,^,^,-^}

о

с£ -8

ьм

-12

0,00 0,01 0,02 0,03 0,04 х, см

Рисунок 8 - Аналитическое и численное решения задачи о волновом разрушении

стекла

Взаимодействие расчётной г и виртуальной } частиц имитирует жёсткую стенку с проскальзыванием. Расчётная область вмещала 12800 8РН-частиц размерами 0.1ммх0.1мм каждая. Все частицы, расположенные на границах расчетной области, в момент времени ¿=0 полагались полностью разрушенными (£>=1) и тем самым имитировались поверхностные дефекты стекла, являющиеся источником зарождения волны разрушения. Результаты расчетов иллюстрируются на рисунке 9, где построены двумерные распределения (уровни значений) параметра разрушения ¿>.и эквивалентного напряжения .

Рисунок 9 указывает на то, что в начальный момент времени в ударяющейся пластине формируются три волны разрушения: прямая фронтальная, идущая от жесткой стенки, и две косых боковых, идущих от свободных поверхностей пластины внутрь материала, нагруженного ударной волной. Косые волны разрушения образуются за счет нагружения свободных поверхностей стекла скользящей ударной волной, и свободная поверхность становится источником волны разрушения после прохождения ударной волны. Одновременно от боковых поверхностей распространяются волны разгрузки. При этом происходит боковой разлет материала, начинающийся также после прохождения ударной волны. Вблизи жесткой стенки материал растекается вдоль ее поверхности, его внешняя поверхность имеет форму, характерную для течения разупрочненной среды. Тестирование показало преимущества разработанного метода вРН по сравнению со стандартным методом в точности расчётов и не обнаружило алгоритмических или схемных дефектов нового метода.

Параметр разрушения

х, см см х, см

а б в

Рисунок 9 - Распределение параметра разрушения в стеклянной пластине в моменты времени 1 мкс (а), 1,6 мкс (б) и 2мкс (в) при скорости удара о жесткую стенку

110= 1000м/с

Дальнейшая валидация разработанного метода БРН была проведена с помощью исследований по моделированию разрушения хрупких материалов при ударном нагружении (по модели разрушения Ш-2 [21]) и по моделированию взаимодействия ударников с преградами конечной толщины при средних и высоких скоростях соударения. Результаты представлены в разделе 2.10 диссертации. Во всех случаях достигнуто удовлетворительное согласие с экспериментами.

В третьей главе представлены результаты по мезомеханическому моделированию ударного нагружения пористого алюминия. Эта задача решена для термо-упруго-пластической среды в плоской двумерной геометрии. Периодическая структура пористого материала задавалась явно и обладала свойствами сплошного алюминия. Алюминий выбран для моделирования, поскольку его физические свойства хорошо известны и адиабаты Гюгонио определены экспериментально как для сплошных, так и для пористых веществ.Набор параметров вещества для моделирования следующий: начальная локальная плотность р0=2700кг/м3, изотермический модуль объемного сжатия Х=73ГПа, модуль сдвига С=23ГПа, предел текучести У0=0.4 ГПа, удельная теплоемкость Су=880Дж/(кг-К), коэффициент Грюнайзена у=2.17, константы адиабаты Гюгонио Са=5350 м/с и 5а=1.35, теплопроводность А,=200 Вт/(м К), объемный коэффициент теплового расширения (3=23 хЮ'6 К'1.Пористая пластина считается перфорированной квадратными отверстиями по нормали к плоскости

расчета. Поры расположены в этой плоскости регулярно, с равными интервалами вдоль горизонтального и вертикального направлений. На рисунке 10 представлен один горизонтальный слой, вырезанный из бесконечной в вертикальном направлении пластины. Этот слой представляет собой расчетную область. Один структурный элемент, с помощью которого создаётся периодическая мезоструктура пористого вещества, выполнен как квадратная рамка с внешними размером 21, заполненная 8РН-частицами. Структурный элемент содержит квадратную пору площадью /х/. Это определяет коэффициент пористости, как т = р0/роо = 4/3, где роо - это средняя плотность материала. Длина перфорированного образца равна I и он содержит N пор. На

рисунке 9 показана начальная укладка 8РН частиц. у

21

О

ро эо

бот со»

элемент поры

ц

9«

пора N

й

оо:

ш

линия симметрии _

оооооооооооо

оо оо, оо оо

.-соодоооооосо оооо

О ООО ОО " ОО

пора 1

о

ообоооообб

ПО.---ОООППООООО'

линия симметрии

Рисунок 10 - Начальная расчетная область и расположение 8РН частиц

Восемь частиц располагаются вертикально. Общее число 8РН частиц составляет 1232 на 25 пор в ряду. В некоторых расчётах использовалось 4928 8РН частиц на 25 пор. Размеры и массы частиц были различными в зависимости от различных размеров пор / = 0.04, 0.4, 4, 40 и 400 мкм. В частицах а и Ъ, выбранных контрольными точками, расчетные параметры сохранялись в отдельных файлах для каждого шага интегрирования, чтобы выявить характер релаксации материала после прохождения ударной волны. Эти частицы расположены на расстоянии Ь/4 от жесткой стенки. На всех границах расчетной области (кроме свободной поверхности), задаются условия для жесткой стенки с проскальзыванием. Вертикальная жесткая стенка в х=Ь означает, что для каждой расчётной БРН-частицы / с параметрами {т,Р,0^,у,их,иу,РьТь31а,8уу,8ху}1 создаётся на каждом расчётном шаге виртуальная частица с параметрами {т,р,0,2Ь-х,у,-их,иу,РьТь5т5уу,-5ху}1. Для горизонтальной стенки в у=0 и у=21 создаются виртуальные частицы с

І,>8XXі$ууі'} і

и

{т,р, 0,Х,41-Уі, и„-иу,Р І9^Іі^ХХУ^ууі'^Ху} І

соответственно. Взаимодействие БРН-части ц из расчётной области и виртуальных частиц, расположенных вне расчётной области, имитирует жёсткие стенки с проскальзыванием. В начальный момент времени при / = 0 все БРН-частицы расчётной области приобретают скорость [/х=Е/0.Когда пористый материал начинает взаимодействовать с жесткой стенкой, в нём формируется ударная волна, распространяющаяся от жесткой стенки к свободной поверхности. На некотором расстоянии за фронтом ударной волны прекращаются колебательные процессы, обусловленные структурной неоднородностью, и параметры вещества приходят к значениям, определяемым адиабатой Гюгонио для пористого материала. Рассматривала динамическая и термическая релаксация Динамическая релаксация в основном проявляется в эволюции колебаний давления, которые наблюдаются за ударной волной. Термическая релаксация вызывается теплопроводностью и приводит к выравниванию температурных неоднородностей в материале за ударной волной. Для алюминия термическая релаксации ожидается значительной при достаточно малых размерах структурных поровых элементов. Расчеты при различных значениях размера элемента свидетельствуют о том, что при /=0.4мкм теплопроводность существенно влияет на процесс динамичной релаксации. Эволюция структуры ударной волны в пористом алюминии для различных скоростей удара пластины о жесткую стенку, иллюстрируется на рисунке 11, где профили давления для всех БРН-частиц построены в разные моменты времени. Рисунок 11а соответствует скорости удара 100м/с. Видно, что волна сжатия характеризуется сильной неоднородностью давления по направлению ее распространения. Интенсивность волны уменьшается, а ширина фронта увеличивается. В головной части волны возникают две пульсации. Уровень нагрузки материала достаточно низок и упругая деформация преобладает в течении материала. Когда скорость удара о стенку возрастает, становится существенной пластическая деформация пористой структуры. Это проявляется на рисунке 116, где при скорости удара 300м/с ясно видна двухволновая структура. Структура упругого предвестника и его эволюция, аналогичны тем, которые показаны на рисунке 11а. Амплитуда упругого предвестника равна 0.2 ГПа. Фронт пластической волны еще неоднороден, но его ширина намного меньше, чем для упругих волн. На рисунках 11в и 11г, при скорости £/о=500м/с и 1000м/с соответственно, интенсивность пластической волны увеличивается, а ее ширина продолжает сокращаться. Упругий предвестник постепенно сокращается, тогда как скорость пластической волны растет. На скоростях удара 2000 м/с и 3500 м/с (рисунок11д и рисунок Не соответственно) образуется ударная волна. За фронтом ударной волны возникают колебания в результате схлопывания пор.

Рисунок 11 - Профили давления при скорости удара t/^ЮОм/с (а), 300м/с (б), 500м/с (в), 1000м/с (г), 2000м/с (д), 3500м/с (е). Размер поры /=0.04 мкм

Профили давления на рисунке 11а-е были использованы для определения фазовых скоростей: скорости ударной волны и скорости перемещения страт

давления. В системе координат, привязанной к жесткой стенке, скорость ударной волны определяется как

^о=(*2-*.И'2-'.) (61)

Координаты х\ и х2 приняты для полуамплитуды скачка давления во время t\ и f2 соответственно. Ударная волна распространяется от жесткой стенки в отрицательном направлении оси X, поэтому значение Us0 отрицательно. Скорость ударной волны Us и скорость частиц Up в системе координат, движущейся со скоростью удара пластины U0, равняются

Us = U0-Us 0, ир=и0 (62)

Адиабата Гюгонио в координатах Up, Us представлена на рисунке 12 сплошными треугольниками (А).

Экспериментальные данные [22] для т=\АЗ и т=1.25 также показаны значками (®) и (а), соответственно. Расчетная кривая Гюгонио 1 располагается между экспериментальными кривыми. При низких значениях скорости частицы, Up <1000 м/с, вычисляемая Гюгонио отклоняется от линейной интерполяции £/¿=1850+1.976^ (штрихпунктирная линия 2). Это отражает преобладание пластичного поведения материала в процессе уплотнения. На скорости частиц менее 200м/с, когда при уплотнении материала упругая деформация преобладает, волновая скорость возрастает до звуковой скорости (»5200м/с).

10000 8000

„ 6000

—-к

г

£5* 4000 2000 0

0 1000 2000 3000 4000

и„, м/с

Рисунок 12 - Вычисленные значения Up, Us для адиабаты Гюгонио (А, кривая 1) в алюминии для пористости т=1.333, линейная экстраполяция t7s=1850+1.9t/p для т=1.333 (штрихпунктирная линия 2) и {/¿=5350+1.35 Up для т=1 (пунктирная линия 3).

Экспериментальные данные для /и=1.43 (®) и для m=1.25 (s) взяты из [22]

3 * •

D

-г -У « 1

У с

Точки А,В соответствуют чисто упругому отклику материала на слабое ударное нагружение. Структура материала не претерпевает существенных деформаций, и слабое возмущение распространяется со скоростью звука. Участок В-С соответствует процессу нагружения с разрушением перемычек между порами. На участке адиабаты C-D происходит частичное компактирование материала в пластических деформациях. Участок адиабаты D-Е соответствует достаточно сильному ударноволному сжатию, при котором процесс схлопывания пор уже не оказывает существенного влияния на линейный характер зависимости U, =AUP). Таким образом, участок A-B-C-D отображает режимы нагружения пористого материала с неполным закрытием пор.

В работе [23] при помощи интерферометра VISAR и датчиков давления были получены экспериментальные профили волн сжатия в образце из пористого алюминия с т= 1.1 при ударе по нему алюминиевым ударником со скоростью Uo=314м/с. Экспериментальные профили из [23] хорошо согласуются с представленными на рисунке 11.

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

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

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

плоской двумерной постановке в термо-упруго-пластическом приближении. Неоднородность среды задавалась включением в несущую фазу (скелет) цепочки восьмиугольных включений примесной фазы. Каждый восьмиугольник примесной фазы содержал 120 8РН-частиц и был окружён 136 частицами несущей фазы, что составляло элементарную квадратную ячейку размером 16x16 частиц. Отношение объемных долей несущей фазы и примеси составляет в этом случае 136:120 = 1.13, т.е. близкую к единице величину. Элементарная квадратная ячейка имеет площадь /х/. Расчетная область ограничена единичным слоем высотой / с жесткими горизонтальными и вертикальной левой стенками. Расчетная область имеет начальную протяженность от соударяемой жесткой стенки до свободной поверхности 1=321, то есть содержит 32 элементарных ячейки. На рисунке 13 показан фрагмент расчётной области в начальный момент времени.

Число БРН-частиц для области протяжённостью I составляет 8192. Размер элементарной ячейки среды принимался различным с целью исследовать влияние на характер релаксации масштабного фактора: /*/ = 0.16x0.16 мкм (X = 5.12 мкм), /х/=1.6х1.6 мкм (I = 51.2 мкм) и 1x1 = 16x16 мкм (Ь = 512 мкм). В расчёте на каждом шаге по времени запоминаются значения термодинамических и кинематических параметров для каждой 8РН-частицы, что позволяет построить профиль давления, температуры и скорости после прохождения ударноволнового фронта по материалу.

элементарная ячейка среды

Жесткая стенка

Рисунок 13 - Схема расчётной области

Моделируемый образец материала в момент времени ?=0 приобретает скорость ир = 2000м/с в направлении начала координат. Ударно-волновой

31

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

Термодинамические свойства компонент описываются уравнением состояния (48). Проведенные расчеты показали, что течение за фронтом ударной волны характеризуется пульсациями в полях давления, температуры и скорости. Проанализирован физический эффект скоростной неравновесности, ответственный за релаксацию полей гидродинамических параметров к установившимся значениям на фронте ударного импульса, без возбуждения колебаний за фронтом, как показано на рисунке 14, где представлены результаты расчёта удара о жесткую стенку слоя лития, содержащего изолированные включения твердого вольфрама. В таком материале ширина фронта увеличивается в 3-4 раза по сравнению с образцом с твёрдым скелетом и жидкими включениями, пульсации давления за фронтом импульса практически отсутствуют. Как видно из рисунка 14, где показан профиль давления для двух моментов времени, профиль близок к стационарному.

100

Рисунок 14 - Профиль давления и скорости в литии с изолированными включения вольфрама (/ =0.16 мкм) на моменты времени ?=0.5нс и ¿=1нс (•- вольфрам, о - литий)

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

процессах релаксации для обоих материалов. В первом случае ширина ударного фронта, как показано на рисунке 15а, составляет примерно 0.2 мкм, а зона релаксации давления примерно 2 мкм. Зона релаксации скорости (рисунок 156) несколько меньше и составляет 1.5 мкм.

-2000 ---1-1--.--•----!--

0,0 0,5 1,0 1,5 2,0 2,5 3,0 ДС.МКМ

Рисунок 15 - Профиль давления и скорости в вольфраме, содержащем изолированные литиевые включения (I =0.16 мкм) на момент времени г=1 не (• - вольфрам, о - литий)

В пятой главе представлены результаты численного моделирования мезоструктуры течения за фронтом детонационной волны, распространяющейся по гетерогенным взрывчатым веществам. Решалась плоская двумерная задача о движении детонационной волны от жесткой стенки по прямоугольной области, заполненной взрывчатым веществом, содержащем вакуумные пустоты. Верхняя и нижняя границы области являлись жесткими стенками. Инициирование детонации осуществляется коротким ударом по левой границе взрывчатого вещества пластиной из вольфрама со скоростью 1600м/с, для чего всем БРН-частицам ударника через некоторое время после удара (обычно это время составляло в расчете 0.05мкс) приписывалась скорость, равная нулю. Расчетная область, вертикальный период которой представлен на рисунке 16, составляла х*у=1000x40 с}2 (где ¿-диаметр БРН-частицы), а толщина пластины равнялась 5 0й?. В расчетной области располагалось от 10000 до 20000 БРН-частиц. Изменение средней плотности ВВ при таком задании расчетной области осуществляется изменением диаметра пор. Детонационная волна распространялась вправо.

Поры моделировалась регулярной по координате X цепочкой окружностей. Располагались поры в шахматном порядке, что достигалось сдвигом каждой чётной (по координате У) цепочки пор на 'Л периода по координате X. В расчетах изменялся размер пор 5.

x=1000d

пластина

Рисунок 16 - Расчетная область (d - диаметр SPH-частицы)

При моделировании учитывались структурные неоднородности порядка 10"4-Н0'5м. Форму и количество пор необходимо задавать на этапе подготовки исходных данных, исходя из представлений о мезоструктуре пористого ВВ. Допускается наличие трёх типов SPH-частиц:

- частицы, содержащие взрывчатое вещество;

- частицы, содержащие продукты детонации;

- «смешанные» частицы, содержащие одновременно ВВ и ПД.

В последнем случае обе фазы (ВВ и ПД) перемешаны по всему объему SPH-частицы и находятся в равновесии по давлению, скорости и температуре. Для расчета взрывчатого разложения ВВ к «смешанным» SPH-частицам применимы макрокинетические уравнения горения, а также уравнения состояния, позволяющие корректно описать обе фазы.

Для конденсированного состояния и продуктов детонации использовалось известное уравнение состояния Джонса-Уилкинса-Ли (JWL) [24] с константами, взятыми для взрывчатого вещества нормальной плотности. Разложение взрывчатого вещества моделировалось с помощью макрокинетического уравнения «ignition and growth» [25], согласованного с уравнениями состояния JWL. Для валидации метода было проведено численное моделирование эффекта кумуляции продуктов взрыва в микроканале взрывчатого вещества (тэн). Получено согласие с экспериментами по вхождению детонационной волны в микроканал взрывчатого вещества.

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

образующихся на поверхности поры при ударе кумулятивной струи продуктов детонации или струи непрореагировавшего взрывчатого вещества (рисунок 17). Ы),2901мкс

Рисунок 17 - Процесс движения детонационной волны по взрывчатому веществу (тэн) с начальной плотностью р0о=Ю06кг/м3. Размер пор 8=0.1мм

Для определения зависимости скорости детонации РЕТЫ (тэн) от средней плотности были проведены расчеты для ВВ с размерами пор 8=0.1мм, 5=0.4мм и 5=1мм при неизменной расчетной области. Результаты расчетов представлены на рисунке 18. Там же нанесены экспериментальные значения из [26] и [27] для ВВ (тэн). В работе [26] приведена эмпирическая зависимость скорости детонации от плотности ВВ (тэн), полученная обработкой большого числа опытов:

|4780 + 3.7(Ро-800), р„<1650 (7920+3.05 (р0-1650); /30>1650

Из рисунка 18 можно видеть, что результат двумерного моделирования по коду БРИ (расчет при диаметре пор 0.4мм) практически совпадает с (63) и другими экспериментальными данными.

Была решена задача о скользящей детонации в слоях гетерогенного взрывчатого вещества насыпной плотности, состоящем из частиц гексогена (ШЗХ) и соли (№С1). Полученные результаты согласуются с экспериментальными данными. При численном моделировании процесса детонации смеси гексоген/КаС1 выделяется расчетная область, в которой располагаются БРН-частицы, содержащие гексоген или КаС1. Верхняя граница области полагается свободной, нижняя граница является жесткой стенкой.

800 1000 1200 1400 1600 1800

Плотность р, кг/м3

Рисунок 18 - Скорость детонации пористого PETN (тэн) в зависимости от плотности ВВ; о - расчет по модели ЗНД с константами из [25], А - расчет 8РН для пор диаметром 0.1мм, ♦ - расчет БРИ для пор диаметром 0.4мм, ▼ - расчет 8РН для пор диаметром 1мм, ® - эксперимент [27], в - эксперимент [26], расчет по соотношению

(63)

Расчетная область представляется в виде периодического повторения квадратной подобласти, называемой шаблоном (рисунок 19). Шаблон в свою очередь разбивается на пространственные ячейки, каждая из которых заполняется (полностью или частично) БРН-частицами. При частичном заполнении ячеек поры в смеси формируются естественным образом как пустые подобласти, в которых БРН-частицы отсутствуют.

Форму и количество пор, а также форму и количество гранул компонент необходимо задавать на этапе подготовки шаблона, исходя из представлений о мезоструктуре смеси. Алгоритм заполнения расчетной области БРН-частицами в описанном подходе сводится к периодическому повторению шаблона вдоль

каждой из координат.

Шаблон представляет собой подобласть пространства из 10x10 квадратных ячеек (рисунок 19), в каждой из которых размещается одна БРН-частица, соответствующая ВВ или ИаС! (либо ячейка остаётся пустой и

соответствует вакууму). Геометрические размеры БРН-частицы, содержащей ВВ или №С1, равны геометрическим размерам ячейки. Расчетная область содержит, таким образом, 45x5 =■ 225 шаблонов. Все 8РН-частицы (как содержащие ВВ, так и содержащие №С1) имеют одинаковые размеры.

зШШ* Ї&з'ЇО: ІШЙ 1 ш щш

Рча ■ ■ 2\

■й • в шк пора ШШ

і У к 11 і

0Й ш||||

т

ш II

Г ■

1

50(1

свободная поверхность

шмшщш ~

11ТТІТ ГГГ П 4Ч-ГІІТГГТ жесткая стенка

поршень

500(1

Рисунок 19 - Расчетная область (<=0) при использовании шаблона, задающего мезоструктуру смеси из 18 вРН-частиц гексогена и 28 вРН-частицШСІ (¿-размер

ячейки шаблона)

При начале расчета пластина-ударник имеет заданную скорость 2500м/с. Через некоторое время { после удара (обычно в расчете £ принималось от 0.5 до 1мкс) всем 8РН-частицам пластины-ударника приписывается скорость, равная нулю. Таким образом, металлический поршень инициирует детонацию и останавливается в момент і. На рисунке 19 толщина пластины-ударника равняется 50<і. Детонационная волна распространяется вправо.

Из рисунка 19 можно предположить, что практически все БРН-частицы инертной примеси будут взаимодействовать с БРН-частицами продуктов детонации, то есть реализуется случай, описанный в главе 1, когда одна из

частиц расширяется, сжимая вещество другой. Моделировать смесь насыпной плотности, представленную на рисунке 19, не представляется возможным с помощью стандартного метода БРН или конечно-разностным способом.

В работе [28] экспериментально были измерены скорости детонации в смеси ВВ (гексоген) и соды для различных толщин слоя смеси и различных размерах гранул ВВ, входящих в смесь. В эксперименте использовались порошки гексогена с размером частиц 1+3 мкм или с размером частиц =70мкм. Первый случай достаточно сложен для численного моделирования на уровне мезомасштаба, так как требует чрезмерных вычислительных затрат. В показанных выше расчётах минимальный размер 8РН-частицы при моделировании полагался 40 мкм. Результаты численного моделирования оказались близки к данным эксперимента (рисунок 20). Расхождение расчётных и экспериментальных результатов составило ± 150м/с, или примерно 5% от установившегося значения скорости детонации. При этом достаточно точно определены критическая толщина слоя смеси (=2мм) и ход изменения скорости детонации в зависимости от роста толщины слоя Я.

3000

2500

2000

■Н 1500 £

а

1000

500 0

0 5 10 15 20 25 30

Н, мм

Рисунок 20 - Скорость детонации в смеси ВВ (гексоген) и инертной добавки рч'аСІ) в

зависимости от толщины слоя смеси Я; в - расчет с использованием при размере частиц ВВ 40+320 мкм, ■ - эксперимент (смесь гексоген/сода при размерах частиц ВВ

1+3 мкм).

Расчеты детонации смеси ВВ с инертным материалом проводились с использованием единого шаблона для описания мезоструктуры смеси (рисунок

19), в различных расчётах изменялись только абсолютные размеры шаблона. Таким образом, каждый из слоёв смеси толщиной Я содержал одинаковое число частиц ~13000 независимо от величины Я и с увеличением Я изменялись исходные размеры частиц.

В Приложении представлен алгоритм для расчета разрушения хрупких материалов по модели Джонсона-Холмквиста (Ш-2), использованный в главе 2 диссертации для валидации кода БРН.

ЗАКЛЮЧЕНИЕ

В диссертации рассмотрено применение решения задач распада разрывов (гидродинамического и температурного) к среде из «гладких частиц» в численном методе БРН, для описания взаимодействия БРН-частиц. Впервые были получены новые уравнения численного метода БРН, основанные на решении задач распада разрывов. Показано, что разработанный метод обладает более высокой точностью в окрестностях контактных границ, чем стандартный метод БРН, использующий искусственную вязкость. Разработанный метод обеспечивает монотонность решения в окрестности контактных границ и показал свою эффективность при мезомеханическом моделировании ударно-волновых процессов в гетерогенных средах с большим числом контактных разрывов плотности. Рассчитанные интегральные характеристики отклика среды на ударное воздействие находятся в хорошем соответствии с результатами экспериментов. При этом получены следующие результаты:

1. При решении тестовых задач проведено моделирование волн разрушения в стеклянных пластинах и обнаружены двухволновые и трёхволновые конфигурации волн разрушения.

2. Проведено исследование влияния масштаба мезоструктуры пористого алюминия в диапазоне изменения размера пор 40нм-400мкм на динамическую и тепловую релаксацию за фронтом ударной волны. Определены характерные значения размера пор, при которых динамическая и тепловая релаксация происходит независимо, а также диапазон размера пор, в котором оба типа релаксации реализуются одновременно. В последнем случае обнаружено взаимное влияние динамической и тепловой релаксации на установление равновесного состояния пористого материала, сжатого ударной волной. Предложена методика определения ударных адиабат пористого материала в переменных «скорость частиц - скорость ударной волны» по расчётным распределениям параметров течения. Построена ударная адиабата для пористого алюминия в килобарном диапазоне давлений, для которого отсутствуют экспериментальные данные. Получены картины расщепления

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

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

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

5. Проведено численное моделирование мезоструктуры течения в детонационной волне, распространяющейся в пористом взрывчатом веществе. Проведен анализ эволюции структуры течения при распространении детонационной волны по пористой среде в зависимости от размера пор при одинаковой средней плотности. Установлено, что при достаточно большом размере пор распространение детонационной волны определяется ее дифракцией на перегородках пористой структуры. При сопоставимых размерах пор и зоны разложения ВВ наблюдается инициирование детонации в горячих пятнах, образующихся на поверхности поры при ударе кумулятивной струи. Достигнуто хорошее согласование данных по скорости детонационной волны в пористом взрывчатом веществе РЕТО (тэн), полученных в двумерных расчетах, с данными экспериментов и расчетов по модели Зельдовича-Неймана-Дёринга.

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Зельдович Я.Б., Райзер ЮЛ. Физика ударных волн и высокотемпературных гидродинамических явлений. -М.: Физматлит, 2008. 656с.

2. Трунин Р.Ф., Крупников К.К., Симаков Г.В., Фунтиков АЛ. Ударно-волновое сжатие пористых металлов // Ударные волны и экстремальные состояния вещества / Под ред. В.Е.Фортова, Л.В.Альтшулера, Р.Ф. Трунина, А.И.Фунтикова. -М.: Наука, 2000. С. 121.

3. Seitz M.W., Skews B.W. Effect of compressible foam properties on pressure amplification during shock wave impact II Shock Waves. 2006. V. 15. P. 177.

4. Zhao H. et al. Perforation of aluminium foam core sandwich panels under impact loading // Int. J. Impact Eng. 2007. V. 34. P. 1147.

5. Ioilev A.G. et al. Numerical model of ductile fracture kinetics: comparison of results of 2-D simulations to experimental data // Int. J. Impact Eng. 2003. V.29. P. 369.

6. Boade R.R. Compression of porous copper by shock waves // J. Appl. Phys. 1968. V. 39. P. 5693

7. Boade R.R. Dynamic Compression of Porous Tungsten // J. Appl. Phys. 1969. V. 40. P. 3781

8. Экспериментальные данные по ударно-волновому сжатию и адиабатическому расширению конденсированных веществ. / Под ред. Трунина Р.Ф. - Саров: РФЯЦ-ВНИИЭФ, 2001. 446 с.

9. Thoma К, Riedel W., Hiermailer S. Mesomechanical Modeling of Concrete Shock Response Experiments and Linking to Macromechanics by Numerical Analysis // ECCM'99. Munich. Germany. 1999. URL: http://hsrlab.gatech.edu/AUTODYN/papers/paperl04.pdf

10. Menikqff R. Interfaces and Reactive Flow // Los Alamos National Laboratory Report. LA-UR-06-7005. 2006.

11. Monaghan J.J. Particle methods for hydrodynamics // Comput. Phys. Rep. 1985. V.3. No. 2. P.71.

12. БомД. Квантовая теория. - M.: ГИФМЛ. 1961. 728с. С.260.

13. Wingate С.А., Fisher H.N. Strength Modeling in SPHC / Los Alamos National Laboratory Report. LA-UR-93-3942. 1993.

14. Monaghan J.J. SPH and Riemann Solvers // J. Сотр. Phys.1997. V.136. P. 298.

15. Cleary P.W., Monaghan J.J. Conduction Modeling Using Smoothed Particle Hydrodynamics // J. Сотр. Phys. 1999. V.148, P. 227.

16. Карслоу Г., ЕгерД. Теплопроводность твердых тел. - М.: Наука, 1964.

17. Уилкинс МЛ., Расчёт упруго-пластических течений.// Вычислительные методы в гидродинамике. / Под ред. Б.Олдера, С.Фернбаха, МРотенберга-М.: Мир, 1967. С.212-263.

18. Капель Г.И., Разоренов С.В., Фортов В.Е., Абахезов М.М. Влияние волны разрушения на динамику импульса сжатия в стекле // 4-е Всесоюзное совещание по детонации - Черноголовка: ИФХ АН СССР,1988. Т.2. С.104-110.

19. Kanel G.I., Rasorenov S.V., Fortov V.E. The failure waves and spallation in homogeneous brittle materials // Shock Compression of Condensed Matter. / Eds. Shmidt S.C. et al- Amsterdam: Elsevier, 1992. P. 451-454.

20. Partom Y. Modeling failure waves in glass // Int. J. Impact Engng. 1998. V.21. No.9. P. 791-799.

21. Johnson G.R., Holmquist T.J. Response of boron carbide subjected to large strain, high strain rates, and high pressures // J. Appl. Phys. 1999. V.85. No.12. P. 8060-8073.

22. Мак-Куш P. и dp. Уравнение состояния твердых тел по результатам исследований ударных волн // В сб. «Высокоскоростные ударные явления» / Под ред. Кинслоу Р. -М.: Мир, 1973. 533 с.

23. Воппап S. et al. Experimental characterization of quasi static and shock wave behavior of porous aluminum // J. Appl. Phys. 1998. V. 83. P. 5741

24. Kapila A.K., Schwendeman D.W., Bdzil J. В., Henshaw W.D. A study of detonation diffraction in the ignition-and-growth model // Combust. Theory and Modeling. 2007. V. 11. P. 781.

25. Lee E.L., Tarver C.M. Phenomenological model of shock initiation in heterogeneous explosives //Phys. Fluids. 1980. V.23(12). P. 2362.

26. Физика взрыва / Под ред. Л.П.Орленко - М.: ФИЗМАТЛИТ, 2002 - Т.1 -832с.

27. Куропатенко В.Ф. Модели механики сплошных сред - Челябинск: Челяб. Гос. Ун-т, 2007. 302с

28. Andreevskikh L.A., Deribas А.А., Drennov О.В., Mikhailov A.L. et al. Mixed Explosives for Explosive Welding of Thin Materials // X International Symposium on EPNM-2010 - Sep. 7-11, Bechichi, Montenegro. 2010.

ОСНОВНЫЕ ПУБЛИКАЦИИ АВТОРА

1. Паршиков А.Н. Применение решения задачи Римана в методе частиц // ЖВМиМФ. 1999. Т.39. №7. С. 1216-1225.

2. Parshikov Anatoly N., Medin Stanislav A., Loukashenko Igor I., Milekhin Valery A. Improvements in SPH method by means of interparticle contact algorithm and analysis of perforation tests at moderate projectile velocities // Int. J. Impact Engng. 2000. v.24. No.8. P. 779-796.

3. Parshikov A.N., Medin S.A. Smoothed Particle Hydrodynamics Using Interparticle Contact Algorithms //J. Сотр. Phys. 2002. v.180. P. 353-382.

4. Медин C.A., Паршиков А.Н. Развитие метода SPH и его применение в задачах гидродинамики конденсированных сред // ТВТ. 2010. Т.48. №6. С. 973-980.

5. Иванов М.Ф, Паршиков А.Н. Численное моделирование динамики ударных волн в композиционном материале // ТВТ. 1993. Т.31. №1. С.92-96.

6. Kanel G.I, Ivanov M.F, Parshikov A.N. Computer simulation of the heterogeneous materials response to the impact // Int. J. Impact Engng. 1995. v.17. No.1-3. P. 455-464.

7. Basko M., Churazov M., Ivanov P., Koshkarev D., Medin S.A.. Orlov Yu., Parshikov A., Sharkov В., Suslin V. Power plant conceptual design for fast ignition heavy-ion fusion // Nuclear Instruments and Methods in Physics Research A. 2005. v.544.

8. Orlov Yu.N., Basko M.M., Churazov M.D., Ivanov P.P., Koshkarev D.G., Medin S.A., Parshikov A.N., Sharkov B.Yu. and Suslin V.M. Energy conversion in a reactor chamber for fast в ignition heavy ion fusion // Nucl. Fusion. 2005.V.45.

9. Fortov V.E., Lebedev E.F., Luzganov S.N., Kozlov A.V., Medin S.A., Parshikov A.N., Polistchook V.P. and Shurupov A.V. Railgun experiment and computer simulation of hyper-velocity impact of lexan projectile on aluminum target // Int. J. Impact Engng. 2006. v.33. P. 254-263.

10. Медин C.A., Паршиков A.H., Орлов Ю.М, Лозицкий И.М. Термомеханические процессы в бланкете реактора ИТС при циклическом воздействии нейтронного флюенса // Атомная энергия. 2011. Т.110, вы п.2. С.92-100.

11.Минеев В.Н, Набоко И.М., Паршиков А.Н., Петухов В.А., Фортов В.Е, Гостинцев ЮЛ. Гусев П.А. Горение и взрыв в замкнутой конической полости. Физический эксперимент // ТВТ. 1999. Т.37. №2. С.313-318.

12. Минеев В.Н, Набоко И.М., Паршиков А.Н., Петухов В.А., Фортов В.Е, Гостинцев ЮЛ. Гусев ПЛ. Горение и взрыв в замкнутой конической полости. Численный эксперимент // ТВТ. 1999. Т.37. №3. С.457-463.

13. Паршиков А.Н., Медин С.А. Численное моделирование волн разрушения при ударном сжатии стекол // Вестник Нижегородского университета им. Н.И. Лобачевского. 2011. №4 Т.5. С.2417-2418.

14. Медин СЛ., Паршиков А.Н. Моделирование распространения волн разрушения при ударном сжатии хрупких материалов (стекол) // Механика Твердого Тела, №2,2012. С.102-113.

15.Паршиков А.Н., Медин С.А. Применение решений распада разрывов в методе SPH / Математическое моделирование: проблемы и результаты / Под ред. О.М.Белоцерковского и В.А.Гущина - М.: Наука, 2003. 478 с. С.320-358.

16.Иванов М.Ф, Паршиков А.Н. Численное моделирование распространения ударных волн в композиционных материалах при импульсном нагружении / Воздействие мощных потоков энергии на вещество / Под ред. Фортова В.Е., Кузьменкова Е.А. - М.: Научное объединение ИВТАН (РАН), 1992,263 с, С.210.

11.Иванов М.Ф., Паршиков А.Н. Влияние микроструктуры композиционного материала на дисперсию ударных волн // Препринт ФИАН им. П.Н.Лебедева №68. Москва, 1992. С.22.

18 .Иванов М.Ф., Паршиков А.Н. Моделирование микромеханики композиционного материала при импульсном нагружении // Препринт ФИАН им. П.Н.Лебедева №69. Москва, 1992. С.31.

19. Паршиков А.Н. Метод SPH на основе решения задачи Римана // Препринт ОИВТ РАН № 2-414. Москва, 1998. С. 17.

20. Медин С.А., Орлов Ю.Н., Паршиков А.Н., Суслин В.М. Моделирование отклика первой стенки камеры и бланкета реактора ИТС на микровзрыв // Препринт ИПМ им. М.В.Келдыша № 41. Москва, 2004. С. 32.

21. Golub V.V., Lu F.K., Medin S.A., Mirova O.A.Parshikov A.N.,Petukhov VA., Volodin V. V. Blast wave attenuation by lightly destructible granular materials // 24ft International Symposium of Shock Waves // Beijing, China, July 11-16, 2004. v. II. part X. P.989-994.

22. Паршиков A.H., Медин С.А. Релаксационные процессы при ударно-волновом нагружении пористых материалов / Физика экстремальных состояний вещества-2007 // Под ред. Фортов В.Е. и др., Черноголовка: ИПХФ РАН. 2007.340 с.

23. Медин С.А, Паршиков А.Н. Численное моделирование структуры ударных волн в гетерогенных двухкомпонентных средах // Физико-

химическая кинетика в газовой динамике. 2008. Т.7. http://www.chemphys.edu.ru/pdf/2008-09-01-015.pdf

24. Медин С.А., Паршиков А.Н. Моделирование мезоструктуры течения при распространении детонации в гетерогенных ВВ // Физико-химическая кинетика в газовой динамике. 2010. Т.9. http://www.chemphvs.edu.ru/pdf/2010-01-12-008 .pdf

25. Паршиков А.Н., Лозицкий И.М. Численное моделирование кумулятивного эффекта в микроканале взрывчатого вещества // Физико-химическая кинетика в газовой динамике. 2011. Т.11. http://www.chemphys.edu.ru/media/files/2011-02-01-

019 Parshikov Lozitskii.pdf

26. Медин С.А, Паршиков А.Н. Использование уравнения состояния JWL и макроскопического уравнения разложения ВВ в методе SPH // 4-я Всероссийская школа семинар «Аэрофизика и физическая механика классических и квантовых систем»: Сборник научных трудов. - М.: ИПМех РАН. 2011. 172с. С.97-102.

27. Медин С.А., Паршиков А.Н. Моделирование скользящей детонации в мелкодисперсной смеси взрывчатых и инертных веществ. // Физико-химическая кинетика в газовой динамике. 2013. Т.15. http://www.chemphvs.edu.ru/pdf72013-04-29-024.pdf.

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

Работа выполнена в Объединенном институте высоких температур РАН (ОИВТ РАН) в период с 1998 по 2013 г.г. при поддержке:

- программы фундаментальных исследований Президиума РАН «Вещество при высоких плотностях энергии» (координатор - академик В.Е. Фортов);

- Межсекционной программы фундаментальных исследований Отделения энергетики, машиностроения, механики и процессов управления РАН «Интегрированные модели физической механики» (координатор - академик Д.М. Климов);

- программы фундаментальных исследований Президиума РАН «Информационные, управляющие и интеллектуальные технологии и системы» (координатор - академик C.B. Емельянов);

- U. S. Civilian Research and Development Foundation, Grant No. RE2-2481-MO-02.

Паршиков Анатолий Николаевич

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

Специальность 01.02.04 - механика деформируемого твёрдого тела

Автореферат диссертации на соискание ученой степени доктора физико-математических наук

Подписано к печати 19.12.2013. Заказ N0 42-2013. Тираж 100 экз.

Отпечатано на ризографе Института проблем механики им. А.Ю.Ишлинского РАН 119526, Москва, пр-т Вернадского, 101,1

Федеральное государственное бюджетное учреждение науки Объединенный институт высоких температур РАН (ОИВТ РАН)

На правах рукописи

/

05201450773

Паршиков Анатолий Николаевич

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

Специальность 01.02.04 - механика деформируемого твёрдого тела

Диссертация на соискание ученой степени доктора физико-математических наук

Москва-2013

ОГЛАВЛЕНИЕ

ВВЕДЕНИЕ................................................................................................................4

ГЛАВА 1. ПРИМЕНЕНИЕ СООТНОШЕНИЙ РАСПАДА РАЗРЫВОВ В ЧИСЛЕННОМ МЕТОДЕ БРН.............................................................................19

1.1. Стандартная формулировка метода БРН.........................................................22

1.2. Модифицированные уравнения БРН................................................................24

1.3. Уравнения 8РН для упругой среды..................................................................27

1.4. Теплопроводность в ЭРН...................................................................................30

1.5. Алгоритм решения трёхмерных упругопластических задач.........................32

ГЛАВА 2. ТЕСТИРОВАНИЕ МОДИФИЦИРОВАННОГО МЕТОДА 8РН39

2.1. Расчёт распада разрыва в упругопластической среде....................................39

2.2. Расчёт распада разрыва в газе...........................................................................40

2.3. Расчёт взрывной волны......................................................................................43

2.4. Расчёт сдвигового течения в жидкости............................................................45

2.5. Расчёт распада температурного разрыва..........................................................47

2.6. Расчёт соударения резиновых цилиндров........................................................49

2.7. Расчёт вращения упругой пластины.................................................................51

2.8. Расчёт разрушения хрупких материалов (стёкол) по волновой модели при ударном сжатии..........................................................................................................52

2.9. Расчёт разрушения хрупких материалов по модели Джонсона-Холмквиста (Ш-2)...........................................................................................................................69

2.10. Сравнение натурных экспериментов и результатов моделирования, проведенного разработанным методом 8РЫ..........................................................82

ГЛАВА 3. МОДЕЛИРОВАНИЕ УДАРНОВОЛНОВОГО НАГРУЖЕНИЯ ПОРИСТЫХ МАТЕРИАЛОВ..............................................................................98

3.1. Формулировка задачи и исходные данные......................................................99

3.2. Динамическая релаксация...............................................................................101

3.3. Термическая релаксация..................................................................................104

3.4. Эволюция структуры ударной волны.............................................................107

3.5. Расчетное построение адиабаты Гюгонио.....................................................111

ГЛАВА 4. МОДЕЛИРОВАНИЕ УДАРНОВОЛНОВОГО НАГРУЖЕНИЯ ГЕТЕРОГЕННЫХ ДВУХКОМПОНЕНТНЫХ СРЕД..................................118

4.1. Постановка задачи............................................................................................122

4.2. Твёрдая несущая фаза с жидкими включениями..........................................125

4.3 Жидкая несущая фаза с твёрдыми включениями...........................................130

4.4. Масштабный фактор.........................................................................................135

ГЛАВА 5. МОДЕЛИРОВАНИЕ ДЕТОНАЦИИ СМЕСЕВЫХ И ПОРИСТЫХ ВЗРЫВЧАТЫХ ВЕЩЕСТВ......................................................139

5.1. Моделирование детонации пористого взрывчатого вещества....................139

5.2. Детонация взрывчатого вещества с включениями парафина......................159

5.3. Моделирование скользящей детонации в мелкодисперсной смеси взрывчатых и инертных веществ...........................................................................166

ЗАКЛЮЧЕНИЕ....................................................................................................179

СПИСОК УСЛОВНЫХ ОБОЗНАЧЕНИЙ......................................................182

СПИСОК ЛИТЕРАТУРЫ...................................................................................184

ПРИЛОЖЕНИЕ....................................................................................................200

ВВЕДЕНИЕ

Численные методы решения уравнений динамики сплошных сред являются практически единственным инструментом для исследования процессов, происходящих в структуре гетерогенных сред при ударно-волновом нагружении, так как аналитического решения подобные задачи, как правило, не имеют. Перечень неоднородных сред, подвергаемых ударно-волновому воздействию, достаточно обширен - это керамики, материалы порошковой металлургии, пространственно-армированные композиционные материалы, компоненты конструкции перспективных энергетических реакторов, вспененные и волокнистые материалы, смеси взрывчатых веществ с инертными добавками. К неоднородным (пористым) средам можно отнести также и гомогенные материалы, находящиеся в режиме объемного вскипания.

Ударное воздействие на гетерогенный материал приводит к сложному силовому взаимодействию между компонентами, составляющими материал, к процессу волнообменов на масштабе структуры материала (или в мезоструктуре) и вызывает интегральный отклик, доступный для экспериментальной регистрации в виде профилей давления и перемещений свободных границ испытуемого образца. Адекватная трактовка результатов таких экспериментов затруднительна без математического моделирования процессов, происходящих в мезоструктуре гетерогенной среды. Альтернативными мезомеханическому рассмотрению гетерогенной среды являются два подхода: первый предполагает, что гетерогенная среда эквивалентна (в смысле физико-механических свойств) некоторой однородной среде [1,2], для которой вводятся эффективные упругие характеристики, связывающие усреднённые компоненты напряжений и деформаций через матрицу жёсткости. Компоненты матрицы жёсткости вычисляются через технические постоянные материалов, составляющих гетерогенную среду (через модули упругости, сдвига, коэффициенты Пуассона) в зависимости от объёмного коэффициента и схемы армирования. После вычисления

эффективных упругих характеристик расчёт напряжённо-деформированного состояния этой однородной анизотропной среды проводится в рамках обобщённого закона Гука. Основную трудность при данном подходе представляет проблема осреднения (то есть определение эффективных характеристик), которая решаетя различными способами: например, методом интегральных сечений [3], с использованием симметрии трансляционного элемента [4], дифференциальным и самосогласованным методами [5]. Подобный подход к моделированию гетерогенных сред справедлив только при низких уровнях ударной нагрузки, когда применимо односкоростное приближение.

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

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

Диссертация содержит описание разработанного автором численного метода «гладких частиц», или SPH (Smooth Particle Hydrodynamics), использующего соотношения распада разрывов, и результаты применения этого метода к решению задач механики гетерогенных сред. Показано успешное применение разработанного метода SPH к моделированию ударного воздействия на следующие материалы:

- на пористый однокомпонентный материал;

- на двухкомпонентный материал без пор;

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

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

Актуальность темы.

Поведение гетерогенных сред при динамическом нагружении описывается различными физическими моделями. Пористый материал, в частности, в области высоких давлений хорошо моделируется осредненными уравнениями сохранения с эффективным уравнением состояния [1], при этом прочность, а также теплофизические свойства должны определяться экспериментально. На атомистическом уровне молекулярная динамика является мощным инструментом для моделирования поведения пустот при нагрузке с большими напряжениями в материале [8]. Но из-за ограниченной производительности компьютера молекулярная динамика встречается с трудностями при моделировании течений на таких пространственных и временных масштабах, какие реализуются в эксперименте. Альтернативная методика, как отмечалось выше, состоит в использовании мезомеханического подхода к численному моделированию [9]-[11], то есть в рамках механики сплошных сред, но с учётом структуры материала. Структура гетерогенной среды задаётся (при лагранжевом формализме описания) такой конфигурацией

расчетных подобластей на поле мезочастиц (для метода частиц), какая соответствует реальной структуре моделируемого материала. Каждая из мезочастиц в этом случае является однокомпонентной, и границы раздела компонент отслеживаются автоматически. Для расчёта передачи импульса и энергии через контактные границы не требуется дополнительного описания. Если необходимо учитывать адгезионные явления и поверхностное натяжение на границах раздела компонент, то для описания взаимодействия мезочастиц необходимы соответствующие модели. Интегральное проявление всех упомянутых физико-механических процессов взаимодействия мезочастиц формирует специфику, характерную для отклика гетерогенной среды на ударно-волновое воздействие. Таким образом, при мезомеханическом описании среды реализуется многоуровневый подход [12] к моделированию процессов в гетерогенной среде.

Частным случаем гетерогенных сред, как уже сказано выше, являются пористые металлы. Интерес к ударному сжатию пористых металлов был изначально вызван потребностью в данных о термодинамических свойствах вещества при высоких давлениях и температурах [13]. Различные металлы были протестированы в ГПа-ТПа диапазонах давлений [14] и получены пористые адиабаты Гюгонио. В результате были построены и верифицированы широкодиапазонные уравнения состояния. Значительное внимание к пористым материалам проявляется в связи с защитой от высокоскоростного удара [15]-[17]. При этом представляет интерес выявить основные физические механизмы поглощения энергии в пористых материалах. Зельдович и Райзер [13] показали, что схлопывание пор отвечает за повышенную диссипацию энергии при сжатии пористого материала. В первых экспериментальных исследованиях Боуда [18]-[19] было установлено, что при слабой динамической нагрузке уплотнение пористых металлов происходит в многоволновом режиме, с одной или несколькими волнами-предвестниками. Эти предвестники являются упругими волнами и волнами разрушения. Адиабаты Гюгонио были определены в области неполного сжатия. Структура течения вещества за ударной волной в

алюминии наблюдалась экспериментально в работах [20]—[21]. Временные профили напряжений, измеренные за ударной волной, характеризуются затухающими осцилляциями, порожденными схлопыванием пор. Мезомеханическое моделирование в двумерном приближении, осуществленное Шуваловым [11] для смеси гранитных зерен и воздуха, показывает аналогичный эффект. Кумулятивные явления, развивающиеся в ходе схлопывания пор, вызывают сильные локальные неоднородности распределения температуры в уплотненном материале [22].

Одним из важных применений мезомеханического моделирования к ударно-волновому нагружению пористой среды является возможность построить расчётным путём адиабату Гюгонио в области неполного схлопывания пор, по известным ударным адиабатам сплошных веществ, которые измерены экспериментально для большого числа веществ и материалов [23]. В области высоких давлений, когда изначально пористый материал становится сплошным, расчёт ударной адиабаты достаточно проработан теоретически [24,25], но в области неполного схлопывания пор результат достигнут с помощью численного мезомеханического моделирования [26].

Моделирование прохождения ударных волн через многокомпонентные среды представляет собою другой обширный круг задач, решение которых важно для практического применения. В их число входит исследование процессов дисперсии ударных волн и диссипации энергии ударного воздействия в композиционных материалах [27]-[31], являющихся основными конструкционными материалами в аэрокосмической технике. Эта задача своим формализмом постановки близка к задаче по моделированию прохождения ударной волны сквозь компоненты устройств перспективных термоядерных реакторов взрывного типа [33]-[37], где элементом конструкции является пористая стенка, насыщенная жидким металлическим теплоносителем.

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

является крайне важным для технологии штамповки взрывом изделий из материалов порошковой металлургии, которые обрабатывать иным способом не представляется пока возможным. При изучении детонации неоднородных взрывчатых веществ (в частности - пористых и гранулированных) накоплен и обобщён значительный экспериментальный материал [38]. Получены эмпирические зависимости параметров детонации от плотности и состава взрывчатого вещества, которые пригодны в большинстве инженерных приложений. Для численного моделирования детонации гетерогенных взрывчатых веществ разработаны многоскоростные (многожидкостные) модели [39]-[41], теоретические основы которых хорошо развиты [42]. Но в последнее время возрос интерес к изучению явления детонации на мезомасштабе взрывчатого вещества [43,44], что связано с необходимостью интерпретации экспериментов по детонации смесевых, насыпных, пористых, флегматизированных, агатированных и содержащих тяжёлые инертные добавки взрывчатых веществ [45].

Интерпретация результатов экспериментов по ударному воздействию на гетерогенные среды требует реалистичного моделирования с учётом мезоструктуры этих сред. Для успешного решения многих прикладных задач достаточно вычислять эффективные характеристики материала [46-48], но при изучении интенсивных процессов, соизмеримых по масштабу со структурными неоднородностями гетерогенной среды, необходимо использовать мезомеханический подход. Это позволяет исследовать непосредственно из результатов численного моделирования те физические особенности отклика гетерогенной среды на воздействие, какие не могут быть получены с помощью смесевых моделей, заменяющих структурно-неоднородную среду однородной средой с эффективными параметрами. Необходимо заметить также, что если ударное воздействие на неоднородную среду приводит к течениям с большими локальными относительными перемещениями компонент [27]-[30], это делает затруднительным даже применение лагранжевых конечно-разностных методов на треугольных адаптирующихся сетках. Аналогично, если рассматривается

детонация низкоплотной среды или газа в области со сложной геометрией [49,50], это также приводит к неприемлемым для сеточных методов локальным искажениям сетки.

Наиболее приемлемым для решения задач мезомеханики, в которых рассматриваются гидродинамические процессы на масштабе структуры среды, является метод «гладких частиц» SPIi (Smooted Particle Hydrodynamics) [51,52]. Метод SPH интенсивно применяется для решения многомерных задач гидродинамики, и к его бесспорным техническим достоинствам следует отнести высокую простоту алгоритмической реализации при минимальном размере программного кода. Отсутствие расчетной сетки позволяет методу (в рамках лагранжева формализма) естественным образом рассчитывать произвольные вращательные и сдвиговые течения, распад односвязных и слияние многосвязных расчетных областей. Более того, свободно-лагранжевы методы дают физически правильную картину эволюции течения в тех случаях, когда применение сеточных лагранжевых методов становится в принципе невозможным вследствие неприемлемых искажений расчётной сетки. Недостатком метода SPH, использующего искусственную вязкость (он называется ниже «стандартным методо