Динамика пространственных упругих и упругопластических тел Текст научной статьи по специальности «Физика»

НАУЧНЫЕ СТАТЬИ И ДОКЛАДЫ • МЕХАНИКА

УДК 539.3

А. Г. ХАКИМОВ

ДИНАМИКА ПРОСТРАНСТВЕННЫХ УПРУГИХ И УПРУГОПЛАСТИЧЕСКИХ ТЕЛ

Разработаны методы и алгоритмы расчета пространственных задач упругопластического деформирования твердых тел. Построены конечно-разностные уравнения типа Уилкинса для решения динамических трехмерных задач с учетом теплопроводности при различных граничных условиях и решен пример. Большие перемещения; упругие и упругопластические тела; пространственные

задачи

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

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

ПОСТАНОВКА ЗАДАЧИ

Основные уравнения. Основные уравнения, описывающие движение элемента твердого деформируемого тела, в прямоугольных декартовых координатах хг (* = 1,2,3) имеют вид:

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

ОТ. ■

—М- = рх\ (* = 1,2,3), уравнение неразрывности

V

V

дх1

дх1

• уравнение энергии (суммирование по повторяющимся индексам)

Е = —(Р + q)V + -V (Sjjijj + Sain),

• уравнения состояния: компоненты деви-атора напряжений

1 V

Sj_j = /4(1 + — 7^3ij у] + ''//•

(*, j = 1,2,3; не суммировать!) относительные деформации

дх^

ёп = ТГ-- , (не суммировать!)

ах’

£ij 2

1 ( дх' dxJ\ 2{1Ы + 0?)’

• гидростатическое давление (среднее напряжение)

Р = -k(InV - 3аТ),

V 5? охр + Зат) ,

Р =

3

• линейная и квадратичная искусственная вязкость

QL =

а Cl Ро

Чк =

V

/~<2 ?/~Л2

'-'ОРО V А

V

нормальные напряжения

Tij = Sij -(P + q), i

условие пластичности Мизеса

Sa S.;

—oz С 0,

'у .. • уравнение теплопроводности

сУТ

Ро с — = \\ У 2Т + Е.

Здесь Т^, £ц, £ц — компоненты тензора напряжений, девиатора напряжений, тензора

деформаций в сокращенной записи (

= 1,2,3), — символ Кронекера, ^ — по-

правки на поворот, — среднее напряжение, взятое со знаком минус, д = дь + дк — искусственная вязкость, — предел текучести материала, V — относительный объем, р, ро — плотность и ее начальное значение, — внутренняя энергия на единицу начального объема, , — модули объемного сжатия и сдви-

га, — коэффициент линейного расширения, — удельная теплоемкость, — коэффициент теплопроводности, — изменение температуры от стандартного значения, равного 288,15 К, £ — время, V — оператор Лапласа,

, — постоянные, — скорость звука, —

объем ячейки, точка над величинами означает производную по времени.

Граничные условия. На свободной поверхности , , — нормальное и ка-

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

(7 — (7ц,

т = \ а

'1'2

(7

13;

&Ы — Ту — (^к^’кгЬз-

Здесь — тензор напряжений, определенный в системе прямоугольных координат , — направляющие косинусы осей новой

системы координат; формулы перехода от одной системы координат к другой и обратно имеют вид

хк = 1ых\ х' = 1ЫХк,

Неподвижная произвольная граница, заданная уравнением

/(х')=0,

а) гладкая поверхность

т = О

^ г‘ = 0

о ? *

ах’

б) прилипание

а < О, т < /х|гг|, х’ = О,

где — коэффициент трения скольжения,

в) скольжение с трением

(7<0, г = /1Н, |^=0,

г) отрыв тела от границы

(7 = 0, т = 0.

Подвижная граница, параллельная плоскости : , — отрыв тела от

границы; , , , — прили-

пание; а < 0, т = Д|(7|, хг = Ц, — скольжение с трением, где — скорость перемещения границы вдоль оси .

Подвижная произвольная граница, заданная уравнениями

х^ —— х^ (ат, ч

где ( ) — криволинейные лагран-

жевы координаты на подвижной границе; а) прилипание

г = г.,

г = гг

б) скольжение с трением

дг,. дг,,

ч да1 да2 в) отрыв тела от границы (7 = 0, т = 0

= 0,

г = гг

где , ,

1»

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

Граничные условия для температурного поля. Температура на границе задана

Т(х\) = То,

где — заданная температура.

Теплообмен с окружающей средой описывается выражением

ОТ

= ас(Т — Тс),

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

Идеальный контакт двух тел

дТл

Тг ~ Т'2: _ ~Х2^1^ где индексы «1», «2» относятся к параметрам первого и второго тела.

ОТ-

2,2

А--’.. VI =

Выделение тепла на границе

тг2

дх' дх’

Ут = \г 1 - 1'2

(У кі Та = Ті •

Я = тУг,

тг2

дх’

Р'2^'2 ,, ; '

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

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

Тл — Т),

1] — 1 О у ;

Е = Е0,

Xі = (ж*)о,

Єу — Єо у;

V = Т'О;

Xі = Ж* (ДО, 0),

* = 0,

Т = То,

Р = Р(),

(І = 1,2,3),

где индекс «0» относится к параметрам твердого тела в начальный момент времени, ДО — лагранжевы координаты твердого тела.

Если ввести безразмерные величины

Xі = Xі/Ь, Г, і = Тп/ст-г. % = %/от,

Г = Р/ат, д = д/сгт, к = к/от,

ї = іа/Ь, а = аТ*,

ц = ц/ат,

Т = Т/Т*,

Хі = ХіТ*/(пЬат),

Р = Ра‘2 /ят, С = с/Г*/а2,

®т = ОсТ* / (аат),

то вид основных уравнений не изменяется (далее черточки над безразмерными величинами пропускаются). Здесь Ь — характерный размер, — характерная температура, — скорость звука.

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

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

$Р( п-ц)^ о Г 5

—- = І1Ш. -----------:-----------, п =

дх’ л-»-о А

с/,Б

СІ8

где — объем, ограниченный поверхностью 5, — вектор-элемент поверхности (вектор

площадки), п — орт внешней нормали к поверхности.

Рис.1

Применяя эти формулы к гексаэдру, объем которого равен А (рис. 1), для функции Р, определенной в вершинах 1 -г- 8, получим

где

№_1 Д

дх* .1 ^

т=1

^1 = (Р1 + ^2 + ^4)/3,

,ЧЬ; = 0,5(гі4 X Г12)/', . . .

А = — (Гі2 X П4 • Г15 О

+ І-:-, | X г32 + Г26 X Г23

г 1т = (4, - Х/)ІЬ

г7в X г78 • г73

■Г38

Г-28

Го8 X Г56 Г85 X Г84

Го2 Н Г82),

X Г£/ ‘ гтп —

хГ1 хт 2 хп 2 хгп г3 а _ гз т

X- — х\ 2 хз 2 — хі. X- _ гз .ь-

Х1 ~ хк 2 хі 2 — хк г3 х1 _ гз

Таким образом, эти величины дают производные в центре гексаэдра, с помощью которых можно получить выражения для производных |р-, |^, в заданной точ-

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

х

(х*‘)п+1/2 = о,5[(х*‘)п+1

Ап+1!2 = {щА»+1 +

где , , , — координаты

точек и объемы гексаэдров в моменты времени и £”+1 соответственно. Тогда конечноразностные соотношения приводят к точному равенству

дх'

дх1

которое равносильно уравнению неразрывности.

Рис. 2

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

1=.? + 1/2, А: + 1/2,1 + 1/2:...

1 = к, .

I = 1 + 1, к + 1,1 + 1;... ,

поэтому

ШЛ _ 1_

дх’) ■

8

т> X] В ~ X]

77=1 77 =1

<9/ = <9?4 + <9т8 + <‘’1112; • • • .

<912 = <91 + <92, ....

Конечно-разностные уравнения. Область, занятая телом, делится на гексаэдры сеткой - - , которая движется вместе с телом (рис. 2). Масса, соответствующая каждому гексаэдру в начальный момент времени, определяется умножением начальной плотности на объем гексаэдра. Например, масса в начальный момент времени для гексаэдра 1 вычисляется по формуле

Массы , , ... , вычисляются анало-

гично.

Сохранение массы

Р

п11 I р / 1

Уравнения движения

(л,г\п+1/2 — (л,г\п^1/2

Vх ~ \ Ьлл

А*п [дТг^п

дх,э

]Л:1

Р]ЛЛ - ^2 РтАт:

777 = 1

= (** )'и.1 + &%8/г дг”+1/2-

Деформации

4-1 / дт’ \ 11 ^2

(^1 2 = ( (не суммировать!),

77+а 1 (дх’ дх^\г' 2

^ ‘ ~2 ^ +

,4"+1/2 = 0,5(.4"+1 + .4”),

(Ае

(Ае;,)'Г1/2 = ) 1+1/2Д^+1/2, (* ф .?),

я+ 1/2

не суммировать!

__ / • \ 71+1/2

гг/1

•д,,у"+1/2 (?у+уг т \-;1+1-у;’ Т), = (V ) м‘+,=^тцг-

у*+1/2 = 0,5 (у-П+1 + Т.П).

Напряжения

/ о' У77+1 —

= (5,-,-)?

2//

(Де*.

1 АТ'" 3 ‘ Т~

77 + 1/2

не суммировать!

(в;.,);*1 = (5Ч)?+^(де,,)Г1/2«ч + (^-)

\П

'//л •

'012 = >• '013 = '023 = (|-

(*11)?

№1)1 (<922) \

2

2 (*912) 1^Ц^12;

'/22 ) 1

= -[V-

77

11/1;

уззя

= 0.

М”! = (виША -1&) -1]-

- [($п)? - (.Ь22)П/п4-

(4з)" = №»>г(в -1)- (здг<*

(4)'1 = №»>?(«- 1) + №*)?<»

/и = СОЙШз, 11 2 = МПО)3, Ш3 = -Ш3у,

8111=

( а" ул + 1 1°й/1

д^я + і/2 (д^2 ^ у^І х Я+1/2 дх1 дх2/1

(^)Г1 (не суммировать!).

( е \»7+1 1°й/1

і о уп+1 \°у)1

'023 \ гг+1

/ гґ N71+1

/'(Д^;/)| ' 1/2'/V/ + ('';/)I-Фі2 = ФіЗ = 0;

= 1,

'ірі

(Зд+1(« -1),

7"

'О

о о

» „

*22 *23

» „

ЬОО 60

/" -*22 — /" _ *33 —

77+1/2

ктаїїу = (#/)ї+1 (^.У" + 1

•013 =

УЇ7 + 1 _

'32 *33

СОЙШі, /23 = 8та;1;

''22; Ш1 — ш1у: д^7, + 1/2 /у^З ^ у^2Х

2 \ <У.т2 <У.т3 / 1

* І^и)ї+\ не суммировать!

р(Аєу),1+1/'2фї-

1; Фі2 = '023 = 0 / // ///

/ - -Г

*32 — *’

23;

С о"\п+1 10гг/1

/ е" уп + 1 1°у/1

У

/ -

'•кг —

Г -*11 —

*31 —

81110*2^ = р77+1 =

п+і/2 _ 41

'31 и *33

г'" г'"

С08 0>2 = /13 = 8111 Ш'2,

г'" ’

—‘13;ш2 — —<^2у!

д^77 + 1/2 /у^1 ^ ^зуя+1/2

2 уУз;3 дх1у1 Ы\^1+\

Сь + С'ЯЛ

4 77 + 1/2

2/3

Шї+1

Р

(^■)ї+1

°’Ї/А

77+1/2

77+1/2“

V

V

1

1-

(Рг‘ ' - + д"

Условие пластичности Мизеса

К”+1 = (ЗиУІ+1(ЗиУІ+1 - 2«тг/3.

Если , то каждое из напряжений

Бу умножается на (25^5^73) 1/,20т- Если же К",+1 < 0, то напряжения не изменяются.

Уравнение энергии элемента твердого деформируемого тела

(ЕУІ+1 = (ЕУ1 - [-

3 каТу

О, 5А;(1пТ/,,+1 4 (УГ1+1

+ (1.5V," ' 1/2 х (ЯуДєу + БцАєцУі ЗцАєц = ;Ь'| І Дг и + 3'2'2А22 + Д: .

ЧІ+1/2

1п V1'

- V*

п

ті + 1/2

Уравнение теплопроводности

■фі(Ег,+1 - ЕГІ)і + ( 8

грП=\ ___ грП

12

777 = 1 \?'=1

где ^ = рк ’ Хі = ^і(Аі)іДі”+1/2.

Устойчивость алгоритма численного расчета

1

Ді77+3/2 = і

сіп

з уГс

ь2

(тіп по jykyl)

Если Д^+3/2 > 1,1Д^”+1/2, то полагаем

Д^+З/2 = 1ДД^+1/2;

1

д^п+1 _ _ ^ д^тг+3/2 _|_ д^п+1/2^

Здесь — скорость звука, —

длина меньшей диагонали гексаэдра,

, , .

Граничные условия в конечных разностях. Свободная поверхность. Для узловой точки свободной поверхности на гра-

ни 1234 суперэлемента (рис. 3) все величины, относящиеся к воображаемым ячейкам 5, 6, 7, 8, принимаются равными нулю. Далее используются уравнения движения для обычной точки. Аналогично проводятся вычисления для точек свободной поверхности, находящихся на ребре 12 или вершине 1 суперэлемента.

Рис. 3

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

ІХ

-гуп+1/2/ ____

>].кЛ ~

■ г ч77+1/2

(X

.77+1/2 \

Скольжение относительно подвижной границы. Орт внешней нормали к поверхности

П = ±

дг,, дг,,

да1 да2

дг,, дг,,

да1 да2

Далее вводится новая система координат хг, причем ось ж1 направлена по внешней нормали к границе твердого тела, ось ж2 находится в плоскости х1Ох3. Напряжения в системе координат определяются по формулам

где 1ц = Щ, 1-21 = Іу.іНі/(1 - 1'п)1/2,

І22 = ^13 — Т”1^ Ьз = _ (1 _ ^1з) 1/ '2>

,,

.

Максимальное касательное напряжение на подвижной границе находится как

^ = \/"?2 - »Ь,.

Если , то напряжения и

с23 умножаются на множитель /і\ац\/т. Если , то напряжения остаются без

изменения. При изменении напряжений производится вычисление напряжений в системе координат хг по формулам Тц = ои^Іу. Далее проводится определение скоростей узловых точек на подвижной границе в соответствии с кинематическим граничным условием.

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

Начальные условия. Начальные условия в конечно-разностном виде записываются в виде

І = 0, (Ту) 0(1) = Т? 0;

{х'ЪзЫ = 0), (x’)ojkl = жг(6'-7, 0),...

Конечный элемент — гексаэдр. Рассматривается определение координат узловых точек конечного элемента в более общем виде (рис. 3). Пусть известны координаты вер-шии^ (т = 1,2,..., 8) конечного элемента — суперэлемента, радиусы кривизн ребер 12, 34, 56, 78; остальные ребра прямолинейные. Координаты узловых точек ребра 12 определяются

•''12/ = 4 + 12Р І1 < і < І2,

где — направляющие косинусы осей относительно осей , — прямоугольная си-

стема координат с началом в вершине 1 суперэлемента, причем ось проходит через вершину 2, ось направлена по внешней нормали к грани 1234, а ось направлена так, что система координат является правой, — координаты узловых точек ребра 12 в системе координат , , — значения на гранях

1584 и 2376 соответственно. Координаты узловых точек на ребрах 34,56,78 определяются аналогично:

■7-і — ! /43 £к і. _ г , /56{к

х8и = Х8 + Ф І1 < І < 32-

Причем

£г2; = £о + Кі2ІЇрі2 НІП/З^.

£г2г = £о + ІЇ12ІЇР12 ніп/З,,

£1% = 0, /Зі = + А «(І - л),

а = бі/2, с.] = \/(4 ^4) • (4-4):

а = 2 агсйіп ( —— )

\ІЇ12/

А а = ок/ (І2 ~ І1),

6 = \/я22 - я2, & = «, £2=6І\ГР12.

Здесь ^ — координаты центра окружности, — линейные размеры (рис. 4), ,

/б^ — углы, іїі2 — радиус кривизны ребра 12, — функция знака кривизны, причем

Npl2 = —1,0, +1. Если Ирі2 = 0, то ребро 12 является прямолинейным.

Рис. 4

Направляющие косинусы определяются

1ц = (ж^ — ж1)/|Г12|,

1-л; = (г 12 - Г14),:/|Г12 х п4|,

^21 = Ънкг — I 12-

^22 = ^33^11 — А",|/ I.'!-^23 = ^31^12 _ ^32^11;

где Г12, Г14 — радиусы-векторы, направленные от вершины 1 к вершинам 2 и 4 суперэлемента.

Координаты узловых точек на гранях 1234 и 5678 находятся по формулам

..

% % . Ж123(1 1 \

х12'Щк — ХГ2] ^ ^ \к — К1),

к\ ^ к ^ к-2,

где к\, к2 — значения к на гранях 1256, 3478 соответственно.

Координаты узловых точек суперэлемента определяются выражением

1^к1

'Г234 і к

хьв Щк х12Щк

к - к 1\ ^ І ^ І2'

где , — значения для узловых точек на

гранях 1234 и 5678 соответственно.

Рис. 5

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

Тестовый пример. Здесь приводятся результаты численного решения тестового примера. Кубик со стороной 10 мм выполнен из стали 4Х5МФС с плотностью 7850 кг/м, модулем упругости 0,215-1012 Па, коэффициентом Пуассона 0,3, пределом текучести 1,7 ГПа. Кубик находится в абсолютно гладких жестких направляющих и подвергается воздействию абсолютно жесткого штампа со скоростью 100 мм/с.

Анализ процесса динамического нагружения рассматриваемого образца показал, что бегущая волна давления монотонно движется вдоль выделенного элемента. Типичная картина распределения напряжений для момента времени 0,52 мкс приведена на рис. 5. Достоверность предложенной математической модели упругопластического деформирования твердого тела подтверждается сопоставлением полученных результатов с известным аналитическим решением [2]. За указанный интервал времени 0,52 мкс фронт волны прошел расстояние 2,7 мм. Напряжения за фронтом волны идентичны и равны МПа.

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

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

2. Соболев, С. Л. Уравнения математической физики / С. Л. Соболев. М.: Наука, 1992.432 с.

ОБ АВТОРЕ

Хакимов Аким Гайфуллино-вич, вед. науч. сотр. Ин-та механики УНЦ РАН, доц. каф. математики. Дипл. инж.-мех. (УАИ, 1970). Канд. физ.-мат. наук по мех. жидкости и газа (Казанск. гос. ун-т, 1977). Иссл. в обл. динамики взаимод. упругих и упру-гопластич. тел со средой.