CC BY
43
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
8
УДК 517.9
Н. Г. Матюшев, И. Б. Петров
МОДЕЛИРОВАНИЕ ДЕФОРМАЦИОННЫХ И ВОЛНОВЫХ ПРОЦЕССОВ В МНОГОСЛОЙНЫХ КОНСТРУКЦИЯХ
Описываются численные методы и вычислительная система, позволяющие решать задачу моделирования деформационных и волновых процессов в многослойных конструкциях при умеренном динамическом нагружении.
Numerical methods and calculating system which enables to perform a modeling of deformation and wave processes inside multilayers constructions under moderate dynamical loading.
1. Введение
Численное исследование задач механики деформируемого твердого тела в случае тел сложной геометрии и неоднородного строения, который представляет практический интерес, предполагает не только исследование деформационных процессов, инициированных интенсивными импульсными нагрузками, но и моделирование волновых процессов, приводящих в конечном итоге к частичному или полному разрушению тела. Немаловажным аспектом моделирования является описание контактных границ, так как процессы многократного отражения и преломления волн на границах формируют волновую картину во всей области интегрирования. Цель моделирования — оптимизация структуры и состава сооружений, предотвращающая возникновение разрушений.
В работе рассматривается задача о соударении деформируемого ударника с многослойными преградами различной конфигурации. Для моделирования поведения материала преграды применялись реологические модели линейно-упругой (закон Гука [2]), упругопластической (модель Прандтля-Рейса с условиями пластичности Мизеса и Ми-зеса — Шлейхера [2; 3]) сред. Для численного решения рассматриваемых задач использовались сеточно-характеристический метод [1; 7], разработанный для исследования данного класса задач в [4], гибридная [5] и гибридизированная [6] сеточно-характеристические разностные схемы, хорошо зарекомендовавшие себя при решении задач с ярко выраженным волновым характером. Они учитывают распространение разрывов вдоль характеристических поверхностей, позволяют корректно строить численные алгоритмы на границах области интегрирования и поверхностях раздела сред, что реализовано в данной статье.
Вестник РГУ им. И. Канта. 2007. Вып. 10. Физико-математические науки. С. 8 — 14.
В ней представлены результаты расчетов для преград с большим количеством слоев и их различным относительным положением, без учета возможного разрушения.
2. Математическая модель твердого тела
Для математического моделирования волновых процессов в деформируемом твердом теле использовалась система динамических уравнений (1—4) в виде
Гр^ і =У ,
(1)
.. — // . - -с.* - - -і- .
і
Здесь Р — плотность среды, VI — компоненты скорости смещения, СТу, £у — компоненты тензоров напряжения и деформаций, V/ — ковари-антная производная по /-й координате, ¥/ — добавочная правая часть.
Вид компонент тензора 4-го порядка Цщ определяется реологией среды. Для линейно-упругого тела
Ц/ы = Х51]5ы + ц(5а 5 у + 8и8/к). (2)
Это соотношение является обобщением закона Гука, X и ц — параметры Ляме, 5у — символ Кронекера.
Плотность определяется из уравнения состояния
Р = Ро /К, (3)
где р = - -12стйс — давление, К = Х+ 2 ц — коэффициент всестороннего сжатия.
Уравнения (1) допускают запись в матричной форме:
д д д -*
—и + А1-----и + А2----и = /. (4)
И дх1 дх2
Здесь и = ^^^, v2, ст11, ст12, ст22, а33)т — вектор искомых функций, / — вектор правых частей той же размерности, А; — матрицы 6x6, явный вид которых приведен в [7], х1, х2 — независимые пространственные переменные, t — время.
3. Численный метод
9
3.1. Внутренние точки
Используемая сеточно-характеристическая разностная схема [1] во внутренних узлах регулярной сетки имеет вид
ип+1 = ип +т( / + й1 + й2). (5)
Здесь /, й1, й2 — составляющие приращения искомого вектора, соответствующие правой части и двум пространственным направлениям. Верхний индекс соответствует шагу по времени.
10
Величины di в формуле (5) аппроксимируют пространственные слагаемые в уравнении (4):
di ~ Ai-d-и.
dXi
Если умножить это выражение слева на левый собственный вектор Ю матрицы Ai, соответствующий значению X, получим
®Tdi ~ XroT и.
dxi
В случае явной схемы Куранта — Изаксона — Риса [1; 7], имеющей порядок 0(т, h), величины di определяются из соотношений
aTdi ~ |X|ST—i1——, (6)
hi
где и± i — значение в соседнем расчетном узле по данному направлению в зависимости от знака X (предыдущем при X > 0 и следующем при X < 0). Система уравнений (6) для различных ю и X используется для
определения di.
В настоящей работе для расчета внутренних узлов наряду с монотонной схемой первого порядка (6), не дающей нефизических осцилляций на разрывных решениях, использованы схема второго порядка, а также гибридная и гибридизированная схемы, сочетающие достоинства двух предыдущих схем [5; 6]. Гибридизированная схема является линейной комбинацией схем первого и второго порядка, коэффициенты которой постоянны и определяются экспериментально, в то время как в гибридной схеме переключение между первым и вторым порядком происходит локально в зависимости от свойств решения.
3.2. Граничные точки
Граничный узел сетки может не иметь соседей по одному или двум направлениям, поэтому для некоторых X соотношения (6) выписать нельзя, и вместо них привлекаются граничные условия. Здесь они формулируются в виде линейных уравнений для компонент и :
a(t, —n, X)T—n+1 = b(t, —n, X). (7)
Матрицы Ai для системы (1) имеют две пары ненулевых собственных значений, отличающихся знаком, отвечающие двум акустическим волнам (продольной и поперечной), которые распространяются в прямом и обратном направлении. Поэтому на каждой границе необходимо поставить по два условия.
Так, расчет узла на верхней границе может быть произведен в два
этапа. Сначала вычисляется d1 с привлечением формул (6) для всех X1,
r n +1 _ __
а затем решение и находится из системы линейных уравнений, со-
держащей два граничных условия (7) и четыре уравнения (5), умноженные слева на собственные векторы ю2/ соответствующие /-20.
Приведенные соображения не применимы для угловых узлов, так как в них оба значения й12 не определены. Поэтому они требуют формулировки специальных алгоритмов. Одним из способов расчета угловых узлов служит переход в систему координат, повернутую относительно исходной. При этом по одному направлению соседними станут два ближайших граничных узла, а по другому — внутренний диагональный. Так как границы реальных тел не имеют особенностей в виде абсолютно острых углов, такое приближение не вызывает физических противоречий.
3.3. Контактные границы
Условия на контактной границе формулируются в вице соотношений между величинами в двух соприкасающихся точках поверхностей, принадлежащих разным взаимодействующим телам. Например, в случаях полного слипания и свободного скольжения они имеют вид
VI = v'l, V2 = v'2, апп =апп, 7т=7т (слитнж);
Vn = Іїп , 7пп = 7пп , Стпт = Стпт = 0 (скольжение).
Здесь штрихованные величины относятся к одному телу, а не штрихованные — к другому. Индексы п и т обозначают соответственно нормальное и касательное направление.
Используемые контактные условия допускают запись в виде четырех соотношений, аналогичных (7):
аТи + (а')т и'= Ь, (8)
которые, как и в случае с граничными узлами, совместно с «внутренними» соотношениями (6) для каждого контактного узла дают 4 + 4 х 2 = 12 соотношений для определения всех неизвестных.
Расчет пары контактных узлов требует решения системы линейных уравнений вдвое большей размерности, однако, ввиду линейности соотношений (7, 8), задача может быть сведена к нескольким расчетам граничных узлов. Задав граничные условия в виде (стпт = 0, vn = 0} и (стпт = 0, vn = 1}, можно восстановить линейную зависимость стпп от vn для каждого контактного узла в отдельности, а затем найти такое значение ^ , при котором 7пп ^п ) = 7'пп ^п ).
При реальных расчетах часты случаи, когда взаимно однозначного соответствия между граничными узлами двух тел на контактной границе нет, например, когда их сетки имеют разный шаг по пространству, либо когда происходит взаимное движение тел. В этих случаях на границе тел вводятся фиктивные узлы. Значения внутренних величин в них определяются интерполяцией по истинным узлам. После расчета всех пар контактных узлов значения в фиктивных узлах снова интерполируются, давая решение в истинных узлах.
11
12
Отдельно должен решаться вопрос об определении момента касания или отрыва. Для этого вводятся критерии начала взаимодействия и отрыва, основанные на сравнении скорости сближения/ отрыва, расстояния между узлами и нормального напряжения с заданными пороговыми значениями. Скорости разлета при необходимости вычислялись путем предварительного расчета контактных узлов как узлов свободной границы. Использование пороговых значений обеспечивает устойчивое слипание / отрыв и отсутствие осцилляций.
4. Примеры расчетов
В работе приводятся примеры расчетов соударения деформируемого сферического ударника с многослойными преградами. На рисунках 1, 2 представлены результаты расчетов для ударов по пятислойным преградам с промежутками и без промежутков между слоями. На рисунке 3 можно видеть результаты расчетов для удара стального шарика по двадцатислойной преграде из тонких листов алюминия с плотно прилегающими слоями. В результате этих и других расчетов были выявлены характерные особенности распространения возмущений в периодических многослойных преградах — параболический или клиновидный фронт волны деформаций и наличие отраженных волн от границ слоев. Расчетная система корректно моделирует слипание, разъединение и соударение слоев. На рисунке 4 показано расслоение многослойной преграды после удара. В результате моделирования ударов по преграде из большого количества слоев получена волновая картина, соответствующая распространению вторичных волн, возникающих при соударениях слоев преграды и направленных против направления удара (рис. 5).
Рис. 1. Поля скоростей на разных этапах удара деформируемого ударника по многослойной преграде с промежутками между слоями: слева: начальная фаза соударения — характерный клиновидный фронт распространения волны деформаций; справа: отраженная волна, повторные соударения слоев
13
Рис. 2. Волновая картина на разных этапах удара деформируемого ударника по многослойной преграде с плотно прилегающими слоями: слева: начальная фаза соударения; справа: отраженная волна, повторные соударения слоев
Рис. 3. Нормальное напряжение при ударе стального шарика по преграде из 20 тонких листов алюминия
іррж
! -.--'Г.к ■% ■;*<і.-3№
V •/ • жШУі
ямнв
Рис. 4. Разделение слоев после отскока шарика
14
Рис. 5. Нормальное напряжение при ударе стального шарика по преграде из 20 тонких листов алюминия, разделенных воздушной прослойкой. Показаны значения в последовательные моменты времени
Список литературы
1. Магомедов К. М., Холодов А. С. Сеточно-характеристические численные методы. М.: Наука, 1988.
2. Новацкий В. К. Теория упругости. М.: Мир, 1975.
3. Новацкий В. К. Волновые задачи теории пластичности. М.: Мир, 1978.
4. Петров И. Б., Холодов А. С. Численное исследование некоторых динамических задач механики деформируемого твердого тела сеточно-характеристическим методом // Журн. вычисл. механ. и матем. физ. 1984. Т. 24. № 5. С. 722 — 739.
5. Петров И. Б., Холодов А. С. О регуляризации разрывных численных решений уравнений гиперболического типа // Там же. № 8. С. 1172 — 1188.
6. Петров И. Б., Тормасов А. Г., Холодов А. С. Об использовании гибридизированных сеточно-характеристических схем для численного решения трехмерных задач динамики деформируемого твердого тела // Там же. 1990. Т. 30. № 8. С. 1237—1244.
7. Магомедов К. М., Холодов А. С. О построении разностных схем для уравнений гиперболического типа на основе характеристический соотношений // Там же. 1969. № 9(2). С. 373 — 386.
Об авторах
Н. Г. Матюшев — асп., МФТИ.
И. Б. Петров — д-р физ.-мат. наук, проф., МФТИ.
