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

Мищенко, Александр Васильевич

Численное моделирование динамического деформирования и разрушения твердых тел в одномерном приближении методом разделения по физическим процессам тема автореферата и диссертации по механике, 01.02.04 ВАК РФ

Мищенко, Александр Васильевич АВТОР

кандидата физико-математических наук УЧЕНАЯ СТЕПЕНЬ

Москва МЕСТО ЗАЩИТЫ

2013 ГОД ЗАЩИТЫ

01.02.04 КОД ВАК РФ

Диссертация по механике на тему «Численное моделирование динамического деформирования и разрушения твердых тел в одномерном приближении методом разделения по физическим процессам»

Автореферат диссертации на тему "Численное моделирование динамического деформирования и разрушения твердых тел в одномерном приближении методом разделения по физическим процессам"

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ имени М.В. ЛОМОНОСОВА

МЕХАНИКО-МАТЕМАТИЧЕСКИЙ ФАКУЛЬТЕТ

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

УДК 539.3

Мищенко Александр Васильевич

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

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

Автореферат

диссертации на соискание ученой степени 2 8 НОЯ 2013 кандидата физико-математических наук

Научный руководитель: доктор физико-математических наук, профессор Киселев Алексей Борисович

Москва - 2013

005540167

Работа выполнена на кафедре газовой и волновой динамики механико-математического факультета Московского государственного университета имени М.В. Ломоносова и в Центре фундаментальных и прикладных исследований ВНИИ автоматики имени Н.Л. Духова.

Научный руководитель: Доктор физико-математических

наук, профессор Киселев Алексей Борисович

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

наук, старший научный сотрудник Пшеничное Сергей Геннадиевич

Доктор физико-математических наук, доцент Медведский Александр Леонидович

Ведущая организация: НИИ системных исследований РАН

Защита состоится «20» декабря 2013 г. в 16 часов на заседании диссертационного совета Д 501.001.91 при Московском государственном университете имени М.В. Ломоносова по адресу: 119992, Москва, Ленинские горы, Главное здание МГУ, механико-математический факультет, аудитория

12-66. 13--Q-4

С диссертацией можно ознакомиться в библиотеке механико-математического факультета МГУ (Главное здание, 14 этаж)

Автореферат разослан «14» ноября 2013 г. Ученый секретарь

диссертационного совета Д 501.001.91, доктор физико-математических наук, профессор

С.В. Шешенин

Общая характеристика работы

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

Цель диссертационной работы. Численное и аналитическое исследование одномерных упругопластических задач деформирования и разрушения твердых тел.

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

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

Научная и практическая значимость. Результаты работы могут быть использованы при тестировании новых численных методов и программных комплексов. Предложен новый численный метод, основанный на методе разделения по физическим процессам с использованием метода Годунова на подвижных эйлеровых сетках. Данный метод используется для решения широкого круга задач механики деформируемых сред. Метод положен в основу комплекса прикладных программ "ТИС". В его создании принимали участие И.С. Меньшов, А.Б. Киселев, П.П. Захаров, A.A. Серёжкин, М.И. Климов, A.B. Мищенко, являющиеся сотрудниками механико-математического факультета МГУ им. М.В. Ломоносова и Центра фундаментальных и прикладных исследований ВНИИ автоматики имени H.JI. Духова. Данный комплекс успешно применяется в ВНИИ автоматики имени H.JI. Духова для проведения расчетов динамики упругопластического деформирования сплошной среды.

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

Работа выполнена при поддержке РФФИ (гранты № 09-01-00144а и № 12-01-00425а).

Апробация работы. Основные результаты диссертации докладывались на следующих научных конференциях и семинарах:

1. Ломоносовские чтения МГУ. Москва (ноябрь 2011, апрель 2012, апрель 2013).

2. Международная научная конференция студентов, аспирантов и молодых учёных «Ломоносов-2013», Москва, апрель 2013.

3. XI Забабахинские научные чтения. Снежинск, 16-20 апреля 2012.

4. Advanced Problems in Mechanics. St. Petersburg, July 2-8,2012.

5. European Congress on Computational Methods in Applying Science and Engineering (ECCOMAS 2012). Vienna, Austria, September 10-14, 2012.

6. XII International Conference on Computational Plasticity. Fundamentals and Applications (COMPLAS XII). Barcelona, Spain. September 3-5,2013.

7. V-VII научно-технические конференции молодых ученых ВНИИ автоматики имени Н.Л. Духова (март 2011, март 2012, март 2013).

8. Научно-исследовательский семинар кафедры газовой и волновой динамики механико-математического факультета МГУ под руководством академика РАН Р.И. Нигматулина.

9. Научно-исследовательский семинар кафедры механики композитов механико-математического факультета МГУ под руководством профессора Б.Е. Победри.

10. Научно-исследовательский семинар кафедры теории пластичности механико-математического факультета МГУ под руководством член-корр. РАН Е.В. Ломакина.

11. Научно-исследовательский семинар кафедры теории упругости механико-математического факультета МГУ под руководством профессора И.А. Кийко.

На защиту выносятся:

- аналитическое исследование волн нагружения и разгрузки в твердом упругопластическом теле при одноосном деформировании в полной постановке;

- численное исследование задачи о плоском соударении тонких пластин с откольным разрушением методом разделения по физическим процессам.

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

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

Содержание работы

Диссертационная работа состоит из введения, двух глав, заключения и списка литературы. В работе содержится 43 рисунка, 81 библиографическая ссылка. Общий объем диссертационной работы составляет 104 страницы.

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

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

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

Математическая модель упругопластической среды Прандтля-Рейса без

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

р+р{ег+кёд) = 0 (1)

■ уравнение движения:

- уравнение энергии:

8(тг , , сг— <7Й

= -в- (2)

or г

PeJJ^l+k^ (3)

or г

- определяющие уравнения модели Прандтля-Рейса:

- критерий пластичности Мизеса:

(4)

Se+XSe=\nH{k)-p-k)£e-Er) (5)

% = 6+5 4 Зк2 Sr2+2k(2-k)SrSe + 2H(k)S^<jY2 (6)

Здесь и далее используются следующие обозначения: р - плотность материала; и - скорость вдоль оси г ; р — давление; <Jr=—p+Sr, <7g=—p + S0 — радиальная и кольцевая компоненты тензора напряжений; Sr, Se — радиальная и кольцевая компоненты девиатора тензора

- с, . 8и

напряжении; Su - интенсивность девиатора тензора напряжении; ег= —,

8г

£д=--радиальная и кольцевая компоненты тензора скоростей деформаций;

с , И2 ^

е = д+— — полная удельная энергия на единицу массы; д — удельная

внутренняя энергия на единицу массы; Н(х) - единичная функция Хевисайда; к - коэффициент симметрии:

к = 0 - для случая плоской симметрии, к = 1 - для случая цилиндрической симметрии, к = 2 — для случая сферической симметрии.

Предел текучести при простом растяжении Y и модуль сдвига р. в данной главе считаются постоянными ( Y = Y0 = const, р — pQ = const).

Среда предполагается термодинамически двухпараметрической, т.е., три термодинамических параметра (давление, плотность и удельная

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

Р = Р{Р,4) (7)

В качестве УРС используются следующие уравнения:

1) УРС твердого тела ("логарифмический закон"):

^НУтА

Здесь Кв - объемный модуль, ау - коэффициент объемного расширения, cv - теплоемкость при постоянном объеме, р0 - плотность недеформированного материала.

2) двучленное УРС:

р = {у-\)рй+с1(р-р0)

Здесь у - показатель адиабаты, с0 - скорость недеформированном материале.

3) УРС Ми-Грюнайзена:

Р = Ро4/{Я)+Ро г0£»

(9)

звука в

где

Ро

(10) (П)

Здесь Г0 - коэффициент Грюнайзена, 5 - константа, связывающая скорость ударной волны Б и скорость частицы среды и: £> = с0 + ,?м.

В материальных соотношениях (4), (5) положим параметр Л = 0 (гипоупругое приближение) и запишем уравнения для нахождения пластических деформаций. Тогда систему определяющих дифференциальных уравнений (1) - (5) в дивергентном векторном виде можно записать следующим образом:

+ - н

где 0 =

/ \ р ри

ри ри2+р-

ре {ре+р-

рїг - вектор консервативных р= риБг

р$в переменных, риБд

РЕРГ риєрТ

,РЕВ, риєрв

(12)

вектор потока,

н,

ри

ри2-Бг+8в

(ре + р-Б,.)и

риБг

риБд

рие?

риерв

- вектор симметрии в левой части,

Н м = р

Г/"о

(ди_

I дг '

к и

г"

ди

дг

2 (ди__к_ П 51, з[аг 2 г) 2ца

ди

дг) 2//0

вектор правой части

Для численного решения системы используется метод разделения по физическим процессам. Система (12) расщепляется на две подсистемы

от дг г

л ~

(13)

(14)

и, соответственно, расчетный цикл временного шага разбивается на 2 этапа, условно называемые "гидродинамический" (эйлеров) и "упругопластический" (лагранжев).

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

Следует отметить, что точное решение существует только для двучленного УРС. При использовании в расчете другого уравнения состояния, его необходимо аппроксимировать двучленным. Пусть р, р, <; -значения плотности, давления и внутренней энергии в данной ячейке на некотором временном шаге, £;0 - выражение внутренней энергии через те же давление и плотность из двучленного УРС. Тогда неизвестные параметры аппроксимации у, р0 и с0 находятся путем приравнивания внутренних энергий £,п и их производных по плотности и давлению. Для повышения точности схемы используется метод подсеточного кусочно-линейного восполнения МиЯСЬ. С явным интегрированием по времени это приводит к монотонной схеме второго порядка по времени и координате. Устойчивость схемы обеспечивается выполнением условия Куранта-Фридрихса-Леви.

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

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

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

др ( д(ри)_0

5? дх

д(ри) д[ри1 + р-8)

дх

= 0

(15)

д(ре) 8((ре+р-8)и) ~дГ~ дх

л 4 ди ЪИо~дх

иф

Здесь и далее х — продольная координата, $ = ЯХ - компонента девиатора напряжений вдоль оси х, которая вследствие интегрирования определяющего соотношения модели (4-е уравнение в системе (15)) с использованием уравнения неразрывности (1-е уравнение в системе (15)), является функцией только плотности:

Т^о,

—±у 1 л0>

при р<ру =р0 ехр , при ру < р < ру при р> ру =р0 ехр

-2К+3&

4 А,

(16)

2Ур + 3£0

4 Мо

Здесь ру, ру - плотности материала при переходе в состояние текучести при сжатии и растяжении, соответственно. Индексом "0" обозначено начальное (невозмущенное) состояние, |50| < 2У0/3.

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

два соотношения, определяющие состояние среды за волной. Это - адиабата Гюгонио (АГ)

(17)

и линия Михельсона-Релея (МР)

сг-о-й=(т)2(у-у0), (18)

где У = \/р — удельный объем, т = р[и-П) = р0[и0-О) - интенсивность сжатия (£> — скорость волны).

В работе приводится подробный анализ зависимости взаимного расположения кривых АГ и МР от интенсивности сжатия т. Показано, что существуют три возможных волновых режима нагружения в упругопластическом материале: одноволновой упругий режим, двухволновой режим с упругим предвестником и одноволновой пластический режим. Схематически все эти случаи изображены на рисунке 1. Черная кривая соответствует адиабате Гюгонио, цветные пунктирные линии — линиям Релея-Михельсона для различных волновых режимов.

---(ш) <(тг)

ОДНОВОЛНОВОЙ упругий рвЖИМ

двухволновой упруголластмческий режим

О») >ю

одноволновой пластический режим

Рис. 1. Взаимное расположение кривых АГ и РМ при различных значениях массового расхода. Возможные волновые режимы нагружения.

Волна разгрузки в твердом теле - решение системы уравнений (15), зависящее только от одной автомодельной переменной Л = ^.

Напряжение <т и удельный объем у связаны обыкновенным дифференциальным уравнением:

с1о а1

с/у

ей"

где

2 , л \2 др с/у 8 V а = (и-А) -щ—

В результате интегрирования этого уравнения с начальным условием =Сг0=— р0 + Я0 (у>у0), получается уравнение кривой сг = сх(У),

являющейся аналогом изэнтропы в газовой динамике. В виду слабого разрыва данной кривой в точке у = Уу возможны два волновых режима разгрузки: одноволновой упругий режим и двухволновой упругопластический режим. Схематически это отображено на рисунке 2.

ОДНОВОЛНОВОЙ двухволновой

упругий режим упругопластическии .

режим 1

0 ¿у

уг У

Рис.2. График кривой а = а(у) в волне разгрузки.

Возможные волновые режимы разгрузки.

В порядке иллюстрации данного анализа приводятся решения задач об ударе по жесткой стенке и об ударном растяжении пластины. В первой задаче полубесконечная пластина налетает на абсолютно жесткую стенку со скоростью и0, во второй задаче скорость и0 направлена от стенки. Материал в начальный момент находится в ненапряженном состоянии. В качестве УРС выбран УРС Ми-Грюнайзена. Таким образом, постановки обеих задач выглядят следующим образом:

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

р\,=0 = Ро

"1=0 = +и

4=0 = 0

= 0

»1=0 = 0

= 0

Граничные условия:

Ввиду нелинейной зависимости давления от плотности в УРС Ми-Грюнайзена, полностью аналитически разрешить обе задачи невозможно. В связи с этим, в задаче об ударе для решения алгебраических уравнений используется метод Ньютона, а в задаче о растяжении для решения

обыкновенных дифференциальных уравнений используется метод Рунге-Кутта. В таблице 1 приведены константы и критические скорости для трех различных материалов. Здесь (и0)у и (и0), - скорости удара, при которых происходит переход в двухволновой упругопластический и одноволновой пластический режимы, соответственно, (щ)у - скорость растяжения, при которой происходит переход в двухволновой упругопластический режим.

Таблица 1. Константы материала и рассчитанные критические скорости для задач об ударе по жесткой стенке и ударном растяжении

алюминий медь бериллий

константы материала

р0, кг/м3 2780 8930 1845

Го 2 2 2

1,338 1,49 1,124

а0, м/с 5330 3970 12870

/л, ГПа 27,6 45 151

У, ГПа 0,29 0,09 0,33

критические скорости

34,03 4,75 18,13

(и0)., м/с 839,38 517,59 3297,2

(м0);, м/с 33,79 4,74 18,1

На данных задачах была проведена верификация численного алгоритма. Результаты приведены на рисунках 3 и 4.

а) 0.12-

га 0.08

•0.04

- аналитика

- численное решение

б>1.4

и()= 10 м/с 2 4

0.2 0

••• аналитика І I— численное решение

и„= 100 м/с

X. см

в) 20_ 16

га р С

— аналитика

— численное решение

и„= 1000 м/с

0° 2 4ХСИ 6 8

Рис. 3. Распределение давления по координате для трех характерных скоростей удара для алюминия.

СО

0,16 0,12

-.0.08

0.04

0 О

— аналитика

численное решение

— аналитика \ 1 — численное решение !

Рис. 4. Распределение напряжения по координате для двух характерных скоростей растяжения для алюминия.

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

р\,=о =Ро

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

Граничные условия:

I R,

4=о =-у

°>L,=0

СГ„

= 0

В качестве УРС выбран УРС Ми-Грюнайзена, в качестве материала -бериллий. Данная задача решена аналитически (Howell В. P., Ball G. J. A Free-Lagrange Augmented Godunov Method for the Simulation of Elastic-Plastic Solids//! Сотр. Phys.-2002.-Vol. 175.-Pp. 128-167).

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

1.5»

0.5:

°0

ПОЛИЛИ IHCprим

кинетическая энергия внутренняя энергия

60 80 t, мкс

100 120

••• пмеиший радиус (аналшикп) j в внешний pajlllVC (Ч1ГСЛСННО) Н внутренний рииус (шюлитика) j — wiyipomiuft pa. щус (числении) 1

60. SO t, МКС

[00 120

Рис. 5. Зависимость энергии (а) и радиусов (б) цилиндрической оболочки от времени.

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

£ев=Н{к)

\(ди+кгЛ+^ Ъ\дг г) 2р

£ =

3{дг

Sг_ 2 rj

2 р

( 2 •Я

, 3 " .3 V /

(20)

1- —— V3SL

■Н\ S,

Здесь г] — динамическая вязкость, которую в данном случае считаем постоянной (г/ = т]0= const).

Для расчетов использовался алюминий со следующим набором параметров материала:

р0 =2780 кг/м3; сс=< 6,72-Ю"5 1/АГ; cv =924,3 Дж/(кг-£) К0 = 78,06 ГПа; //0 = 27,6 ГПа; У0 = 0,29 ГПа; % = 700 Па • с

На рисунке 6а показаны зависимости энергии от времени. На рисунке 66 представлены зависимости радиусов оболочки от времени. Видно, что численное и аналитическое (Киселев А.Б. К исследованию процесса нестационарного расширения толстостенных сферических и цилиндрических вязкопластических оболочек // Вестн. Моск. ун-та. Сер. 1. Матем. Механ. -2012. — № 6. - С. 20-25.) согласуются с высокой точностью (погрешность порядка 0,2 %).

а) 4

б) 52|

501

полная энергия кинетическая энергия внутренняя энергия

К 48

0/,

0 50 100 150 200 250 300 1. мке

0 50 100 150 200 250 300 I, мке

Рис. 6. Зависимость энергии (а) и радиусов (б) сферической оболочки от времени.

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

Рис. 7. Зависимость скоростей внутренней (а) и внешней (б) границ оболочки от времени.

Во второй главе представлена используемая математическая модель повреждаемой упруговязкопластической среды (Киселев А.Б. Математическое моделирование динамического деформирования и комбинированного микроразрушения термоупруговязкопластической среды // Вестн. МГУ. Сер. 1. Матем. Механ. - 1998. - № 6. - С. 32-40), корректировка вычислительной схемы при учете повреждаемости и разрушения, а также подробно рассмотрена задача о плоском соударении тонких пластин с откольным разрушением.

Математическая модель повреждаемой упруговязкопластической среды помимо законов сохранения (1) - (3) включает в себя определяющие уравнения модели Соколовского-Пэжины (20). В модели возможен учет

о *

о 50 100 150 200 250 300 I. МКС

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

У0=¥т+/3^, (21)

где еЦ = ~ интенсивность тензора пластических деформаций, ¡3 -

параметр упрочнения. В случае одноосного деформирования е? =\£р\-

В общем случае микроразрушение материала описывается с помощью двух скалярных параметров поврежденности со и а. Они характеризуют наличие микроразрушений типа пор сферической формы и типа полос адиабатического сдвига соответственно. В случае одноосного деформирования можно положить а = 0. Кинетическое уравнение для нахождения параметра объемной поврежденности со ^0<®<1 ^ выглядит следующим образом:

о> = в(^-аЫ^-а.]+Ф?^Н{а-<тЛ+со?^Н{а--ст), (22) (Л-® ) и-© ) 4% * > Ащ ^ >

где

а+ = ——¥0\па>, а = -сг+; В, сг — константы материала.

Параметры материала при наличии повреждений меняются следующим образом:

К = К0(1-а),М = М0{1-со), ц = т}0{1-а>), Г = У0(1-ю)

В качестве критерия начала макроразрушения материала используется критерий предельной удельной диссипации:

£> = /-(</м+^)Л = А, (23)

о р

где с1и = сг^ЁЦ - механическая диссипация, (1Г = Аа>2 - диссипация континуального разрушения; и - время начала разрушения; А — предельная удельная диссипация; Л — константа материала.

В качестве уравнения состояния используется УРС твердого тела, который с введением поврежденности преобразуется следующим образом:

р=к

(24)

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

др д(ри)_

а &

д(ри) д(ри2 + р-8)_

сЧ дх

д(ре) | д((ре+р-8)и) 0Г дх

Ее решение строится описанным выше способом. Разница лишь в изменении компонент расчетных векторов (2, Г и Н„.

Пусть в некоторой ячейке области выполнился критерий разрушения: £>>£).. В этом случае считается, что через данную ячейку проходит поверхность разрушения. Для определенности полагаем, что поверхность разрушения проходит через центр ячейки. Тогда исключаем разрушенную ячейку из дальнейшего расчета и производим разбиение данной области на две новых области по центру данной ячейки. На образовавшихся новых границах ставим условие свободной поверхности. Пересчет параметров в соседних ячейках производится с выполнением закона сохранения массы.

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

-/г, О К

Рис. 8. К задаче о соударении пластин. Экспериментально данная задача была подробнейшим образом изучена в работах Г.И. Канеля с соавторами для множества различных материалов и

скоростей соударения.

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

Ч=о = А/

Ударник:

Мишень:

4=о=о

4=0= Ах

4=0=0

4=0 = 0

Здесь и далее параметры с нижним индексом " г" относятся к ударнику, а с индексом " I" - к мишени. Верхними индексами " Ь " и " Я " обозначены, соответственно, левая и правая границы каждой из пластин.

Начальные координаты пластин зададим следующим образом:

4\ =0

X* (=0 = хЛ =0 ' 1<=0

X? 1=0

На левой границе ударника и правой границе мишени ставится условие свободной поверхности:

г \х=х<-

сг,1* =0

Условие на контактной границе между пластинами (х = 0) имеет вид:

«и* = =и" ^Ц" =^Ц£ ПРИ М <0 ^и^^Ц^0 ПРИ К1>0

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

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

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

На рисунке 9 приведен график сравнения численного расчета и экспериментальных данных для теста по соударению алюминиевых пластин, в котором толщина ударника равна /?, = 2 мм, толщина мишени И, = 4,1 мм, скорость соударения и, = 690 м/с. Параметры материала выбраны следующими:

= 2610 кг/м3; с„ =924,3 Дж/кг-Л"; ау = 6,72-Ю 5 1/К;

К0 = 71,94 ГПа; /и0 =26,22 ГПа; У00 =0,18 ГПа

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

г/0 = 78,7 Па• с; /3 = 0,5 ГПа; сг.=0,65ГПа;

Л=9,36 1 04 Па с; В = 9,ЗЗ-Ю"3 1/(#а-с); а=30кДж/кг

600

"I 400

(Л

200 0

Рис. 9. Скорость свободной поверхности мишени при соударении алюминиевых пластин.

В данном расчете упругий предвестник выходит на свободную поверхность в момент времени te = 0,65 мкс. Разрушение происходит в момент времени ?,= 1,13 мкс, относительная толщина отколотой части пластины-мишени (откольной тарелки) составляет /г, = 0,47 (h, = \-X,, где

X, — относительная координата разрушения, X = XR Х' L — текущая

■

vAíJ

эксперимент q расчет

0.65................1.................................l'X..........................2........................2*5

t (ИКС

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

Также был рассмотрен тест с соударением алюминиевой (ударник толщиной А, =2 мм) и титановой (мишень толщиной /г, =10 мм) пластин со скоростью = 700 м/с.

Параметры материала для титана выбраны следующим

/>0 =4450 кг/м3; с„ = 520,7 Дж/кг-^; «„ = 2,52-Ю 5 1/К; ^„ = 116,65 ГПа; /лй =38,74 ГПа; 700 = 1,08 ГПа

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

7/0 = 688 Па с; /? = 0,1 ГПа; ег, =3,85 ГПа;

Л=5,33-104 Па-с; 5 = 2,74-Ю"2 1 /(Яа-с); А = 50 кДж/кг

На рисунке 10 приведен график скорости свободной поверхности пластины-мишени. Здесь также можно отметить близкое совпадение экспериментальных данных и численного расчета. В данном расчете упругий предвестник выходит на свободную поверхность в момент времени 4 = 1,6 мкс, а разрушение происходит в момент времени Г, = 2,4 мкс.

Толщина откольной тарелки составляет =0,17 от толщины пластины-

мишени.

600

о 400 2

^ 200

Рис. 10. Скорость свободной поверхности мишени при соударении алюминиевой и титановой пластин.

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

Рис. 11. Распределение удельной диссипации (а) и деформации (б) по толщине пластины-мишени момент разрушения.

На рисунках 12а и 126 изображены зависимости времени разрушения

и относительной толщины откольной тарелки пластины-мишени от скорости

удара. Как можно видеть по данным графикам, время разрушения имеет

тенденцию к снижению при увеличении скорости удара. Толщина откольной

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

составляет приблизительно пятую часть от толщины пластины-мишени.

«) Зг

2,8

2 2,6 s i 2.4

2,2

650 700

900

1000

650 700

u0, М/С

1000

Рис. 12. Зависимость времени разрушения (а) и относительной толщины откольной тарелки (б) от скорости удара.

В заключении приводятся основные результаты работы и выводы:

1) Были впервые подробно и в полной постановке аналитически исследованы волны нагружения и разгрузки в твердом упругопластическом теле при одноосном деформировании. Рассчитаны параметры перехода из одного волнового режима в другой для различных материалов.

2) Предложен новый численный метод расчета упругопластических задач - метод разделения по физическим процессам. Метод успешно верифицирован на одномерных задачах об ударе по жесткой стенке, ударном растяжении пластины, сжатии цилиндрической оболочки и расширении сферической оболочки.

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

Список публикаций

1. Киселев А.Б., Мищенко A.B. Одномерные упругопластические задачи в плоской постановке. Аналитические и численные решения // Вестн. Моск. ун-та. Сер.1. Матем. Механ. -2014. -№ 2.

2. Меньшов И.С., Мищенко A.B., Серёжкин A.A. Численное моделирование упругопластических течений методом Годунова на подвижных эйлеровых сетках // Математическое моделирование. - 2013. - Т. 25. - № 8. - С. 89-108.

3. Menshov I., Mischenko A., Serezhkin A An eulerian Godunov-type scheme for calculation of the elastic-plastic flow equations with moving grids // Europ. Congress on Comput. Methods in Appl. Sc. and Eng. (ECCOMAS 2012). J. Eberhardsteiner et. al. (eds.). Vienna, Austria, September 10-14, 2012. CD format, 2012. Article 2164. 20 p.

4. Мищенко A.B., Серёжкин A.A., Меньшов И.С., Киселев А.Б. Метод разделения по физическим процессам для моделирования деформирования и разрушения твердых тел // Забабахинские научные чтения: сб. материалов XI Межд. конф. 16-20 апреля 2012. - Снежинск: Изд-во РФЯЦ-ВНИИТФ, 2012. -С. 306.

5. Киселев А.Б., Меньшов И.С., Мищенко A.B. Программный комплекс «ТИС»: тестирование на задачах динамики твердого тела // Ломоносовские чтения. Научная конф. Секция механики. Апрель 2012 года. Тезисы докладов. - М.: Изд-во Моск. ун-та, 2012. - С. 90-91.

Текст научной работы диссертации и автореферата по механике, кандидата физико-математических наук, Мищенко, Александр Васильевич, Москва

московским государственный университет

имени М.В. ЛОМОНОСОВА

механико-математическии факультет

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

04201365740

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

Мищенко Александр Васильевич

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

Научный руководитель: доктор физико-математических наук, профессор Киселев Алексей Борисович

Москва-2013

Содержание

Введение..................................................................................................4

Глава 1. Математическое моделирование упругопластического деформирования твердых тел в одномерном приближении без учета микроповреждений и разрушения.................................12

1.1. Описание математической модели упругопластической среды.... 12

1.2. Численный метод и основные характеристики вычислительной схемы....................................................................................................20

1.3. Аналитические решения упру го пластических задач в случае одноосного деформированного состояния.......................................33

1.3.1. Волна нагружения в твердом теле....................................34

1.3.2. Волна разгрузки в твердом теле............................................41

1.4. Аналитическое решение задачи о расширении толстостенной сферической оболочки.......................................................................46

1.5. Верификация численного алгоритма........................................51

1.5.1. Задача об ударе пластины по жесткой стенке......................51

1.5.2. Задача об ударном растяжении пластины......................55

1.5.3. Задача о сжатии толстостенной цилиндрической оболочки..............................................................57

1.5.4. Задача о расширении толстостенной сферической оболочки..............................................................61

Глава 2. Математическое моделирование упругопластического деформирования и разрушения повреждаемых твердых тел в одномерном приближении................................................................66

2.1. Математическая модель повреждаемой упруговязкопластической среды.....................................................................................................66

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

2.3. Задача о плоском соударении тонких пластин. Постановка

и валидационные расчеты.................................................................80

Заключение..........................................................................................94

Список литературы............................................................................96

Введение

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

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

К настоящему времени разработаны несколько численных подходов [65] для моделирования упругопластических волновых процессов, которые обладают определенными индивидуальными преимуществами и недостатками. Например, в лагранжевых методах [48, 64, 78] расчет изменения параметров среды происходит в каждой конкретной частице, что упрощает постановку граничных условий и позволяет отслеживать положение поверхности материала в процессе соударения. Однако при больших деформациях может происходить значительное искажение расчетной сетки, что приводит к потере точности результатов. Для эйлерова подхода [44, 52, 53, 54, 57, 58, 72, 76], когда изменение параметров рассматривается в неподвижной точке пространства, менее актуальны трудности, связанные с большими деформациями. Но, например, отслеживание изменения положения контактных границ и свободной поверхности является более сложной задачей, поскольку возникают счетные ячейки, частично заполненные различными средами. Решение данной проблемы, к примеру, методом концентраций [36] приводит к размытию границы, и, как следствие, потере точности.

Кроме того, часто используются методы конечных элементов [49, 73], на основе которых создан ряд коммерческих вычислительных комплексов, например, ANSYS, LS DYNA [56] и др. Также для решения данных задач используются сеточно-характеристические методы, например [13, 38, 39]. Как альтернатива сеточным, активно исследуются и нередко используются бессеточные методы, например, SPH (Smooth Particle Hydrodynamics) метод -метод сглаженных частиц [62, 70].

При решении упругопластических задач в последнее время применяются также гибридные методы [63], которые используют преимущества как лагранжевых, так и эйлеровых схем. К этому классу методов можно отнести рассматриваемый в настоящей работе метод [37, 66], который основан на принципе разделения по физическим процессам [32, 50, 80] и использует подвижные эйлеровы сетки. Решение задачи на каждом временном шаге ищется в два этапа. На первом этапе решается система уравнений так называемого в литературе гидродинамического приближения в предположении постоянства в каждой лагранжевой частице среды ее упругопластических параметров (девиаторные компоненты тензора напряжений, пластические деформации и поврежденности). Решение строится на подвижной эйлеровой сетке. На втором этапе данное решение корректируется с учетом выбранной модели упругопластического деформирования. Сетка на втором этапе остается неподвижной.

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

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

2 Под валидацией численного решения понимается сравнение данного решения с экспериментальными данными

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

Кроме верификации на аналитических решениях в данной работе приводится валидация численного алгоритма на хорошо изученной экспериментально задаче о плоском соударении тонких пластин [17, 18, 19]. Необходимо отметить, что толщины пластин малы по сравнению с их размерами, благодаря чему, данная задача рассматривается в постановке одноосной деформации. Таким образом, при решении данной задачи

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

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

Научная и практическая значимость. Результаты работы могут быть использованы при тестировании новых численных методов и программных комплексов. Предложен новый численный метод, основанный на методе разделения по физическим процессам с использованием метода Годунова на подвижных эйлеровых сетках. Данный метод используется для решения широкого круга задач механики деформируемых сред. Метод положен в основу комплекса прикладных программ "ТИС" [37, 66]. В его создании принимали участие И.С. Меньшов, А.Б. Киселев, П.П. Захаров, A.A. Серёжкин, М.И. Климов, A.B. Мищенко, являющиеся сотрудниками механико-математического факультета МГУ им. М.В. Ломоносова и Центра фундаментальных и прикладных исследований ВНИИ автоматики имени H.JI. Духова. Данный комплекс успешно применяется в ВНИИ автоматики имени H.JI. Духова для проведения расчетов динамики упругопластического деформирования сплошной среды.

На защиту выносятся:

Основные результаты работы представлены в следующих публикациях:

1. Киселев А.Б., Мищенко А.В. Одномерные упругопластические задачи в плоской постановке. Аналитические и численные решения // Вестн. Моск. ун-та. Сер.1. Матем. Механ. - 2014, №2.

2. Меньшов И.С., Мищенко А.В., Серёжкин А.А. Численное моделирование упругопластических течений методом Годунова на подвижных эйлеровых сетках // Математическое моделирование. — 2013. - Т. 25 - №8. - с. 89-108.

3. Menshov I., Mischenko A., Serezhkin A. An eulerian Godunov-type scheme for calculation of the elastic-plastic flow equations with moving grids // Europ. Congress on Comput. Methods in Appl. Sc. and Eng. (ECCOMAS 2012). J. Eberhardsteiner et. al. (eds.). Vienna, Austria, September 10-14, 2012. CD format, 2012. Article 2164. 20 p.

5. Киселев А.Б., Меньшов И.С., Мищенко A.B. Программный комплекс

«ТИС»: тестирование на задачах динамики твердого тела //

Ломоносовские чтения. Научная конф. Секция механики. Апрель 2012

года. Тезисы докладов. - М.: Изд-во Моск. ун-та, 2012. - с. 90-91.

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

- Ломоносовские чтения МГУ. Москва (ноябрь 2011, апрель 2012, апрель 2013).

- Международная научная конференция студентов, аспирантов и молодых учёных «Ломоносов-2013», Москва, апрель 2013.

- XI Забабахинские научные чтения. Снежинск, 16-20 апреля 2012.

- Advanced Problems in Mechanics. St. Petersburg, July 2-8, 2012.

- European Congress on Computational Methods in Applying Science and Engineering (ECCOMAS 2012). Vienna, Austria, September 10-14, 2012.

- XII International Conference on Computational Plasticity. Fundamentals and Applications (COMPLAS XII). Barcelona, Spain. September 3-5, 2013.

- V-VII научно-технические конференции молодых ученых ВНИИ автоматики имени Н.Л. Духова (март 2011, март 2012, март 2013).

Научно-исследовательский семинар кафедры газовой и волновой динамики механико-математического факультета МГУ под руководством академика РАН Р.И. Нигматулина.

Научно-исследовательский семинар кафедры механики композитов механико-математического факультета под руководством профессора Б.Е. Победри.

- Научно-исследовательский семинар кафедры теории пластичности механико-математического факультета МГУ под руководством член-корр. РАН Е.В. Ломакина.

- Научно-исследовательский семинар кафедры теории упругости механико-математического факультета МГУ под руководством профессора И.А. Кийко.

Работа выполнена при поддержке РФФИ (гранты № 09-01-00144а и № 12-01-00425а).

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

1.1. Описание математической модели упругопластической среды

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

упругую Б? и пластическую Б? составляющие:

е,=ч1+е: (1.1)

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

¿„=¿¡4' (1-2)

Пластическое течение предполагается несжимаемым:

¿¿=о (1.3)

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

аи=ойа+8у, а = ^ = (1.4)

где р - давление, 8 - символы Кронекера. Из определения следует:

3*=0 (1.5)

В данной главе процессы упругопластического деформирования твердого тела рассматриваются без учета накопления микроповреждений и разрушения. В качестве упругопластической модели выбрана классическая модель Прандтля-Рейса. В общем (трехмерном) случае ее материальные соотношения имеют следующий вид [14, 15, 16, 20, 48]:

2/4 (1.6) Условие текучести в форме Мизеса [15, 16, 48]:

«<§Г С1-7)

Здесь = • Зу - интенсивность девиатора напряжений, ё^ — —^-Зу -

девиатор скоростей деформаций. Предел текучести У и модуль сдвига ¡л считаются постоянными: У = У()= СОШ(, // = //0 = СОШ(. В данных предположениях эта модель является моделью идеальной пластичности.

3МА4

В материальных соотношениях Я = 0 в упругой области и X

У 2

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