Разработка алгоритма расчета условий на свободной и контактной границах для моделирования деформирования материалов методом Sph Текст научной статьи по специальности «Физика»

УДК 539.3

Разработка алгоритма расчета условий на свободной и контактной границах для моделирования деформирования материалов методом SPH

А.В. Герасимов, P.O. Черепанов

Научно-исследовательский институт прикладной математики и механики ТГУ, Томск, 634050, Россия

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

Ключевые слова: гидродинамика сглаженных частиц, метод SPH, математическое моделирование, контактные границы

Computational algorithm for free and contact surface conditions in SPH simulation of deformation of materials

A.V. Gerasimov and R.O. Cherepanov

Research Institute of Applied Mathematics and Mechanics of Tomsk State University, Tomsk, 634050, Russia

A new algorithm is proposed for realization of free and contact surface conditions in numerical simulation of a dynamically deformed material with internal structure by the smoothed particle hydrodynamics (SPH) method. The peculiarity of the proposed algorithm is in the retention of the approximation order at the boundaries and in the possibility to combine the SPH method either with finite difference or finite element methods for no change of numerical schemes. The algorithm allows consideration of an arbitrary number of contacts of arbitrary solids calculated by different methods.

Keywords: smothed particle hydrodynamics, SPH method, mathematical modeling, contact boundaries

1. Введение

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

Бесспорным достоинством метода является отсутствие численной сетки и использование весьма простых

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

© Герасимов А.В., Черепанов P.O., 2010

2. Краткое описание метода

В оригинальном методе SPH [1-4] ядерная аппроксимация произвольной функции в точке j может быть получена путем умножения этой функции на функцию сглаживания и интегрирования полученного произведения по всему рассматриваемому объему:

f -J f (x)W(x - xi, h)dx, (1)

где h — параметр сглаживания, выбираемый достаточно произвольно; x — пространственная координата. Производная этой функции может быть найдена как:

fa=f -1 f (x) Wa(x - x'’ h)dx- (2)

dxa

Соответствующая (2) узловая аппроксимация имеет вид:

f!a=f - Е fW (Xk - X, Й)Д/, (3)

dxa k

где xk, fk, Дvk — радиус-вектор, значение аппроксимируемой функции и некоторый ассоциированный объем, соответствующие k-й точке.

На функцию сглаживания W(xk - X, h) накладываются следующие ограничения:

J W (x, h)dx = 1, J Wa (x, h)dx = 0, (4)

W(x, h) = 0, IX > kh, lim W(x, h) = S(x).

1 1 h^0

Достаточно подробный анализ разработок в области метода SPH приведен в [4], однако из-за весьма бурного развития метода в последнее время данные быстро устаревают.

Ранее было показано [1], что аппроксимация (3) имеет первый порядок точности во внутренних областях при однородном распределении частиц и нулевой — вблизи границ расчетной области и внутри области при неоднородном распределении частиц. Эта проблема в англоязычной литературе получила название particle inconsistency, что можно перевести как узловая несогласованность, т.е., несмотря на то что ядерная аппроксимация имеет высокий порядок точности, соответствующая ей узловая аппроксимация имеет более низкий порядок точности, вплоть до нулевого (говорят, что узловая аппроксимация не согласована с ядерной). Существуют различные подходы к решению этой проблемы [5], но в данной работе мы будем основываться на методе, предложенном в [6].

Идея метода заключается в использовании разложения функции в ряд Тейлора и применения к этому разложению формулы (1):

J f (x)W (x - x )dx = f (xt) J W (x - x )dx +

+ f,a (xi)J (xa- Xa)W(x - x )dx> (5)

J f (x) Wp (x - xi- )dx = f (x) J We (x - xi- )dx +

+ f,a (xi)J (xa- xa)We (x - xi )dx- (6)

Соответствующие узловые аппроксимации будут иметь вид:

X f (Xп )Ж (хп - Xі )уп = Г (хп - Xі )уп +

N

+ & 1 (4 - 4 Ж(Xп - Xі у, (7)

XГ(Xп)Гр(Xп -Xі)уп = /ІЖв(Xп -Xі)уп +

N

+ Га 1 (ха - 4 Же (Xп - Xі )уп. (8)

В (5)-(8) и далее принято правило суммирования по повторяющемуся нижнему индексу. Для верхних индексов правило суммирования по повторяющемуся индексу не применяется, если специально не оговорено иное.

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

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

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

3. Основные идеи предлагаемого подхода

Рассмотрим сначала следующую задачу: некоторая SPH-точка а получила приращение скорости и приращение координат в результате какого-либо внешнего воздействия. Требуется выяснить, каково полное изменение импульса соответствующего тела и какое влияние полученные приращения оказали на состояние материала в окрестности этой точки. В данном случае речь идет именно о внешнем воздействии. Изменения координат и скоростей, вызванные интегрированием уравнений движений по времени, являются частью самого метода и здесь не рассматриваются.

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

вестно поле скоростей

V = ?(х), (9)

тогда

I =| у(х)р(х^К (х). (10)

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

| / (х)Ш (х - X )ёх = /1Ш (х - X )ёх +

+ /% I (X*- 4 )Ш (х - хг' )ёх, (11)

I / (х)Шр (х - хг' )ёх = /1(х - хг' )ёх +

+ I(ха- х‘а )^р (х - х*^ (12)

а = 0, ..., й - 1. (13)

Здесь и далее dx = dV(x) — соответствующий элемент объема; индекс й — размерность задачи (2 — для двумерной, 3 — для трехмерной). Верхний индекс соответствует номеру точки, к которой относится рассматриваемая величина, нижние индексы показывают номер компоненты для векторов и тензоров; если перед нижним индексом стоит запятая, это обозначает производную по соответствующей координате. Нумерация координат начинается с 0 (возможно, читателю это покажется необычным, однако полученные в такой записи формулы легче преобразовывать и они проще для программирования).

Для того чтобы уравнения (11), (12) привести к виду, более удобному для дальнейших выкладок, введем следующие обозначения:

Г(х)

,а*-1,

Г (хг'), а = -1,

А(х, У)а =

ха- У а

(14)

[1, а = -1, а = -1, ..., й - 1.

Теперь систему (11), (12) можно записать в виде одного уравнения:

IГ (х)Шр (х - хг' )ёх =

= х* )а (х - х* ^ (15)

а, Р = -1, ..., й.

Введем новые обозначения:

^ (Г, х) = 1Г (У )Щр (У - х)ау,

7р«(х) = 1А(у; х)а (У - х^У> (16)

Вав (х) = [ëР(х)] 1 а, Р = -1, ..., й.

Используем их для записи произвольной функции /:

Г(х),а = (Г, х) • Вав (х) =

= [1 ДУ^У - x)dy]• в«р(х)- (17)

Если в (17) вместо / подставить удельный импульс ур и использовать получившуюся аппроксимацию в (10), то с учетом того, что значению функции соответствует значение а = -1, получим следующую формулу:

1 = I[Р(х)У(хМх] = 1[^ (ур, х) • (x)dx] =

= 1[ (х)ёх! р(У)у(У )Шр (У - x)dy ] =

= II Р(У) • (х) • у(у ) • (у - х^х, (18)

а = -1.

Формула (18) является ядерной аппроксимацией полного импульса тела. Переходя к узловым аппроксимациям, имеем:

1 = IIР(У) • Вав (х) • У(У) •(У - x)dxdУ =

= Х

N

А • ВфЕр"Vм • А

м

(19)

а=-1

где Ат — площадь SPH-точки в двумерных задачах и объем — в 3D.

Изменим порядок суммирования:

I = £| утрт • Ат •£ А” •

1

Введем

пт

ар • п,в

м I N /а=-1

Мт =

•X А”

N

о” Ш”т • Бав •

(20)

(21)

а=-1

В таких обозначениях уравнение (20) можно записать в виде:

I = Х (Vм ^Мт). (22)

м

Выражение (21) представляет собой эффективную массу SPH-точки.

Таким образом, мы получили очень простую формулу для вычисления импульса SPH-точки. Заметим, что получившаяся эффективная масса вовсе не равна произведению объема точки на ее плотность и может изменяться со временем.

В случае решения осесимметричной задачи в выражении (21) член А” должен быть заменен наАпуп = = А^п.

Изменение полного импульса тела, вызванное изменением скорости одной его SPH-точки вычисляется по формулам (20), (21).

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

ращений скоростей и координат в методе Джонсона [7], считая SPH-точки «ведомыми», а сегменты «сеточной» границы — «ведущими». Однако этого недостаточно для использования метода Джонсона, так как необходимо еще точно вычислить ускорения свободной поверхности.

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

= 1/2 ^. +SRj.i.), (23)

где , 8^ — компоненты радиус-вектора текущего

положения точки и компоненты вектора смещения соответственно.

Однако в предположении малости приращений, подставляя в (23) выражение для производных (17), можно получить выражение, связывающее приращения координат узла т с приращениями деформаций в узле п:

8^.= (8^. + 8^., )/2 =

= ^ (8Д-, R) • Б в (К) +

+ ^ (8Л,-, R) • Вв (R) +

5Я1

Я

= 1(У Ж,р (у - ^ • В—р(R) dy + +1 Щ (У) Жр (у - R) • Врр (R) dy + ^Я|(У)

+ 1

Яі(у)

Жр (у - R) • В-1в (^у.

(24)

В (24) и далее члены, помеченные *, присутствуют только в осесимметричной двухмерной постановке.

Формула (24) есть ядерная аппроксимация приращения удельного объема точки. Переходя от нее к узловой аппроксимации аналогично (19), получим следующую формулу:

8е” = X (8^п • Б% + 8Rm • Б”) • ж”” • Ат +

м

(

I

м

Я

ттт-тп гуп лп

• В-1;Р •А

(25)

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

приращения компонент вектора координат в чистом виде. Воспользуемся для этого преобразованиями вида:

&£■ = 8гу • Впр + 8ят8Л • Вр) • )т +

м

м

Я

тп п т

Жр • В-1;Р • А

Х5ят ( В”р+5 вр) • Ат +

м

Я

тп п т

Шр • В-1;р^ А

(26)

Отсюда имеем:

5е” =!\5Ят • АтЖт м I

+ 5іу • В1р + 5Л • Вір

51у51і51 ] • В-п1;р

Я1т

(27)

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

(28)

Подставляя в (28) полученные ранее выражения для приращений деформаций и меняя порядок суммирования, получим:

5Е = 12 а” А Vn 115^т • АтЖт х

м

51у51і51 — • В-п1;р

Я1т

+ 5іу * В]р+5-У • Вір

= І]5Я; • Аті2а”АпЖ,ртпх

м

51у51і51 — • В-1;р

Я1т

+ 5iy• В”р + 5—у • Вір

= 2{{т^т|,

(29)

где

F™ = Л”AVnW”n х

SiASij '

+ SiY

' B^/P + 8 jY

гр

(30)

Теперь рассмотрим ситуацию соударения двух тел, одно из которых аппроксимируется конечно-разностной сеткой, а другое — методом SPH.

В наиболее известных алгоритмах расчета условий на контактных границах для метода сглаженных частиц [8-10] граница SPH-облака (а следовательно, и моделируемого тела) представляет собой переходный слой ненулевой толщины. Мы же можем считать, что граница тела проходит точно через центры узлов и представляет собой «строгую» границу.

На рис. 1 приведен пример взаимодействия облака частиц (5PH) с разностной сеткой (FDM). Узел а взаимодействует с сегментом Ы1-Ы1. Для определения приращения координат и скоростей необходимо вычислить эффективную массу узла а и использовать ее для вычисления приращений скоростей и координат по методу Джонсона [7]. Вычислив приращение координаты узла а, необходимо по формуле (27) вычислить приращения деформаций всех его соседних узлов (на рис. 1 лежат внутри круга радиуса 2Н) и по новым значениям скорректировать напряжения в этих точках.

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

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

граничного слоя ячеек для конечно-разностной схемы и кинетической энергии слоя толщиной Н для метода SPH. Даже на достаточно редких сетках (20x20), они составляют порядка 2-3 %, а на более густых эти потери уже можно считать пренебрежимо малыми. Если такая точность не удовлетворительна, можно провести вычисления, подобные (23)-(27), для того чтобы определить, какой вклад в изменение внутренней энергии вносит смещение точек и какой вклад в изменение внутренней энергии вносят потери при жестком соударении. Впрочем, последнее легко оценить:

ДЕ = М •

(31)

Формула (31) верна как для метода SPH, так и для конечно-разностной схемы, только в последнем случае это изменение энергии необходимо каким-либо образом распределять между соседними с точкой ячейками. Это вызвано тем, что одна точка может быть узлом сразу нескольких ячеек, а внутренняя энергия вычисляется именно в ячейках. Таким образом, навязанное извне изменение скорости точки должно вызывать приращение внутренней энергии всех связанных с ней ячеек, и для определения этих приращений необходимо проводить дополнительный анализ.

Из (27) следует, что зависимость между приращениями деформаций и приращениями координат узлов линейная, и может быть выражена в виде:

(32)

jmn

yij •

Se” =£SRm •Y”

M

Это позволяет вычислить работу внутренних сил на перемещении узла m:

nnnVn =

SEm =£SRm Y”

n

5 V ITfmn^njrn 5d” 17”

= SRy Е YYj j = SRY FY •

(33)

Рис. 1. Взаимодействие SPH-облака с FDM-сеткой

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

ат = рт1 (Гт рт). (34)

Главной особенностью формулы (34) является то, что она позволяет вычислить ускорение любой, как внутренней точки расчетной области, так и граничной. Данный факт заслуживает особенного внимания.

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

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

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

Описанный алгоритм можно представить в виде последовательности шагов:

1. Во всех граничных SPH-узлах вычисляются их эффективные массы по формуле (21).

2. Во всех граничных SPH-узлах вычисляются ускорения по формулам (30), (34).

Рис. 2. Взаимодействие двух БРН-облаков

3. Вычисляются новые координаты и скорости всех SPH-узлов.

4. Производится расчет движения всех остальных тел, рассматриваемых в задаче, тем методом, который для них выбран.

5. Проводится поиск точек и узлов, пересекших границы других тел.

6. Проводится поиск таких приращений координат и скоростей точек, при которых будут выполнены условия непротекания на контактных границах и законы сохранения импульса и момента импульса по формулам Джонсона [7], если рассматривается контакт SPH-FDM или SPH-SPH. Если для моделирования используется какой-либо другой метод, то для него необходимо провести анализ, подобный предлагаемому в данной работе, и определить эффективные массы точек, поправки к деформациям и т.д.

7. Для SPH-точек, получивших приращения скорости или координаты, и всех их соседей проводится корректировка деформаций (27) и вычисляются соответствующие приращения напряжений.

8. Пункты 6-8 повторяются до тех пор, пока не останется ни одной расчетной точки, для которой нарушено условие непротекания.

9. Осуществляется переход на новый временной слой.

4. Результаты расчетов

В качестве модельной задачи был выбран процесс высокоскоростного соударения алюминиевых цилиндров. Результаты расчетов приведены на рис. 3. Начальная скорость цилиндров составляла 400 м/с, радиус г — 20 см, длина цилиндров Ь — 20 см. Для моделирования использовались совместно метод SPH (481 частица, расположенная в узлах квадратной сетки размером 20x20) и метод Уилкинса [11] с квадратной сеткой 20x20 ячеек. Расчет граничных условий осуществлялся методом Джонсона [7]. Использовалась модель идеально упругопластического тела с критерием текучести Мизе-са. В методе SPH параметр сглаживания выбирался автоматически таким образом, чтобы у каждого из узлов было не менее 4 и не более 16 соседей.

Проведенные расчеты с использованием предложенных в данной работе подходов показали, что закон сохранения импульса при аппроксимации граничных условий выполняется точно. Использование формул (37), (38) для вычисления ускорений узлов приводит к нарушению консервативности схемы по импульсу, в отличие от классического метода SPH [2] или метода CSPM [5]. Полные потери импульса при моделировании соударения двух одинаковых цилиндрических тел составляют порядка 1.5-2.0 %, что, на наш взгляд, является вполне приемлемым результатом. Полные потери энергии в ходе расчетов составили порядка 2.0 % даже при дли-

Рис. 3. Результаты расчета деформаций и давлений при соударении двух алюминиевых цилиндров (г = 20 см, Ь = 20 см) методом БРН (слева) и методом Уилкинса (справа). Начальная скорость соударения — 400 м/с. t = 50 (а), 200 мкс (б)

тельных (20-50 пробегов упругой волны) расчетах, что тоже является хорошим результатом.

Результаты сравнений расчетов методом БРН и методом Уилкинса показали, что качественно и количественно решения совпадают с хорошей точностью, что косвенно подтверждается данными, приведенными на рис. 3. Давление в ударной волне, рассчитанное методом Уилкинса, составляет 62 кбар, давление в ударной волне, рассчитанное методом БРН, составляет 60 кбар, частично разница в давлениях может быть объяснена влиянием искусственной вязкости [11, 12].

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

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

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

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

Работа выполнена при частичной финансовой поддержке грантов РФФИ №№ 10-08-00633-а и 10-08-00453-а.

Литература

1. Bonet J., Kulasegaram S. Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations // Int. J. Numer. Meth. Engng. - 2000. - V. 47. - No. 6. -P. 1189-1214.

2. Lucy L.B. A numerical approach to the testing of fusion hypothesis // Astron. J. - 1977. - V. 82. - P. 1013-1024.

3. Gingold R.A., Monaghan J.J. Smoothed particle hydro-dynamics: theory and applications to non-spherical stars // Mon. Not. Roy. Astron. Soc. - 1977. - P. 181:375-389.

4. Libersky L.D., Randles P. W., Carney T.C., Dickinson D.L. Recent improvements in SPH modeling of hypervelocity impact // Int. J. Impact Eng. - 1997. - V. 20. - No. 6-10. - P. 525-532.

5. Chen J.K., Beraun J.E., Jih CJ. A corrective smoothed particle method for transient elastoplastic dynamics // Comput. Mech. - 2001. - V 27. -No. 3. - P. 177-187.

6. Liu M.B., Liu G.R. Restoring particle consistency in smoothed particle hydrodynamics // Appl. Numer. Math. - 200б. - V. 5б. - No. 1. -P. 19-36.

7. Johnson G.R., Stryk R.A. Symmetric contact and sliding interface algorithms for intense impulsive loading computations // Comput. Meth. Appl. Mech. Eng. - 2001. - V. 190. - No. 35-36. - P. 4531-4549.

8. Campbell J., Vignjevic R., Libersky L. A contact algorithm for smoothed particle hydrodynamics // Comp. Meth. Appl. Mech. Eng. -2000. - V. 184. - No. 1. - P. 49-65.

9. Campbell P.M. Some new algorithms for boundary value problems in SPH // DNA-TR-88-286, Mission Research Corporation, Albequerque, 1989.

10. Belytschko T., Neal M.O. Contact-impact by the pinball algorithm with penalty and Lagrangian methods // Int. J. Num. Methods Engng. -1991. - V. 31. - No. 3. - P. 547-572.

11. Уилкинс М.Л. Расчет упругопластических течений // Bычис-лительные методы в гидродинамике. - M.: Mир, 1967. - С. 212-2б3.

12. Nejad-AsgharM., Khesali A.R., Soltani J. Artificial viscosity in simulation of shock waves by smoothed particle hydrodynamics // Astro-phys. Space Sci. - 2008. - V. 313. - No. 4. - P. 425^30.

Поступила в редакцию 20.04.2008 г., после переработки 2б.03.2010 г.

Сведения об авторах

Герасимов Александр Владимирович, д.ф.-м.н., снс, зав. отд. НИИ ПММ ТГУ, [email protected] Черепанов Роман Олегович, мнс НИИ ПММ ТГУ, [email protected]