авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

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


Pages:     | 1 |   ...   | 4 | 5 || 7 |

«РОССИЙСКАЯ АКАДЕМИЯ НАУК СИБИРСКОЕ ОТДЕЛЕНИЕ ИНСТИТУТ ГИДРОДИНАМИКИ им. М. А. ЛАВРЕНТЬЕВА ИНСТИТУТ ВЫЧИСЛИТЕЛЬНОГО МОДЕЛИРОВАНИЯ ЧИСЛЕННОЕ РЕШЕНИЕ ...»

-- [ Страница 6 ] --

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

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

µn+1 µn, µn где µn = µ(1 + n /2), n значение на n-й итерации.

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

274 Глава 5. Динамика упругопластического деформирования 5.10.2. Устойчивость процесса вычислений Исследуем устойчивость процесса вычислений по начальным дан ным.

Используя энергетическое тождество (5.55) и условия непрерывно сти скоростей и напряжений на границах элементов, находим, что при отсутствии внешних сил в плоском случае имеет место равенство:

u1 u H1 H3 u1 + u3 + t t +11 e11 + 33 e33 + 213 e13 ] dd + Q} = 0, где означает суммирование по элементам, составляющим область S.

Замечая, что согласно (5.36) ij ij = ij +, 2 t и используя первую схему аппроксимации уравнений упругопластиче ского деформирования, находим, что для функций (5.36) и (5.39) имеют место следующие соотношения:

1 ij eij = + ij, e =.

2µ t 2µ K t Следовательно, 1 ij 1 H1 H3 ij eij dd = H1 H3 ij + dd+ 2µ t K t 1 + H1 H3 dd.

2µ ij ij В силу определения (5.2) величины 1 ij 1 H1 H3 ij eij dd H1 H3 ij + dd, 2µ t K t 1 следовательно, ij u1 u3 1 H1 H3 u1 + u3 + ij + dd + Q 0.

t t 2µ t K t 5.11. Сходимость приближенного решения к точному Таким образом, (EK + EP ) d + Q 0, t где EK кинетическая энергия;

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

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

5.11.1. Сходимость приближенного решения к точному в плоской задаче Пусть 11, 33, 13, u, u напряжения и скорости, соответствую 1 щие точному решению задачи (5.1), (5.29), (5.30) в области S, которая состоит из четырехугольных элементов, на промежутке времени [0, T ] при начальных условиях (5.34) и граничных условиях (5.33). Разбивая промежуток времени на конечное число интервалов длительностью и используя на каждом из интервалов изложенную выше процедуру, можно построить приближенное решение на всем промежутке време ни [0, T ]. Обозначим через E квадрат энергетической нормы разности между точным и приближенным решениями, соответствующими t = T :

T (u u1 ) (u u3 ) 1 H1 H3 (u u1 ) + (u u3 ) E= + 1 t t 0 S +(11 11 )(e e11 ) + (33 33 )(e e33 )+ 11 +2(13 13 )(e e13 )} dSdt. (5.99) 276 Глава 5. Динамика упругопластического деформирования Используя (5.29), (5.30), (5.37), (5.39), находим T (u u1 ) (u u3 ) 1 H1 H3 (u u1 ) + (u u3 ) dSdt = 1 t t 0 S T (11 11 ) (31 31 ) (u = u11 ) + (u3 u31 ) + x1 x 0 S (13 13 ) (33 33 ) +(u u13 ) + (u3 u33 ) + x3 x (11 11 ) (31 31 ) +(u11 u1 ) + (u31 u3 ) + x1 x (13 13 ) (33 33 ) +(u13 u1 ) + (u33 u3 ) + x3 x +(u u1 ) M13 H3 (13 13 ) M31 H1 (33 33 ) + 0 +(u u3 ) M31 H1 (13 13 ) M13 H3 (11 11 ) 0 dSd, (5.100) T H1 H3 [(11 11 )(e e11 ) + (33 33 )(e e33 )+ 11 0 S T 13 )(e +2(13 e13 )] dSdt = H1 H 0 S (u (u u31 ) u11 ) 1 (11 11 ) + (31 31 ) + x1 x (u u13 ) (u u33 ) 1 +(13 13 ) + (33 33 ) + x3 x (u u11 ) (u u31 ) 1 +(11 11 ) + (31 31 ) + x1 x (u u13 ) (u u33 ) 1 +(13 13 ) + (33 33 ) + x3 x +(u u0 ) [M31 H1 (33 33 ) M13 H3 (13 13 )] + 1 +(u u0 ) [M13 H3 (11 11 ) M31 H1 (13 13 )] dSdt. (5.101) 3 Из (5.99)–(5.101) следует E = E0 + E1 + E2, (5.102) 5.11. Сходимость приближенного решения к точному где T [(11 11 )(u u11 ) + (31 31 )(u u31 )] + E0 = 1 x 0 S [(13 13 )(u u13 ) + (33 33 )(u u33 )] dSdt, + 1 x T (u0 u1 )(M13 31 M31 33 ) + u M13 (31 H3 13 ) E1 = 1 0 S M31 (33 H1 33 ) + (u0 u3 )(M31 13 M13 11 )+ 0 +u M31 (13 H1 13 ) M13 (11 H3 11 ) 0 dSdt, T (11 11 ) (31 31 ) E2 = (u11 u1 ) + (u31 u3 ) + x1 x 0 S (13 13 ) (33 33 ) +(u13 u1 ) + (u33 u3 ) + x3 x (u u11 ) (u u31 ) 1 +(11 11 ) + (31 31 ) + x1 x (u u13 ) (u u33 ) 1 +(13 13 ) + (33 33 ) dSdt.

x3 x Из граничных условий (5.33) и неравенств (5.35) следует E0 0. (5.103) Разложим величины, соответствующие точному решению, в ряд:

u u 11 11 1 = + O( + h1 + h3 ), = + O( + h1 + h3 ), x1 x1 0 x1 x1 (5.104) M13 31 M31 33 = (M13 31 M31 33 )0 + O( + h1 + h3 ) и т. д., где обозначено: (...)0 = (...)|1 =3 ==0 ;

h1, h3 линейные размеры четырехугольного элемента.

Используя (5.104), находим T T (u0 (u0 u1 )O( + h1 + h3 )dSdt, u1 )(M13 31 M31 33 )dSdt = 1 0 S S 278 Глава 5. Динамика упругопластического деформирования T T M13 (31 H3 13 )u dSdt = 0 M13 (13 H3 13 )O( + h1 + h3 )dSdt.

0 S S Аналогично оцениваются и остальные слагаемые в выражении для E1.

Можно показать, что приближенное решение ограничено равномерно по, h1, h3 и, соответственно, ограничены величины u0 u1, u0 u3, 1 0 0 0 31 H3 13, 33 H1 33, 13 H1 13, 11 H3 11. Тогда E1 0 при, h1, h3 0. (5.105) Используя выражение (5.56), запишем E2 в (5.102) в виде T + (u31 u3 ) 31 + (u13 u1 ) 13 + E2 = (u11 u1 ) x1 x1 x 0 S u u + (11 11 ) 1 + (31 31 ) 3 + +(u33 u3 ) x3 x1 x u u +(13 13 ) 1 + (33 33 ) 3 dSdt Q, (5.106) x3 x3, где означает суммирование по элементам области S и интервалам, времени, составляющим промежуток [0, T ]. Из (5.102), (5.103), (5.106) находим T + (u31 u3 ) 31 + E+ Q E1 + (u11 u1 ) x1 x, 0S u 13 + (11 11 ) 1 + +(u13 u1 ) + (u33 u3 ) x1 x3 x u u u +(31 31 ) 3 + (13 13 ) 1 + (33 33 ) 3 dSdt. (5.107) x1 x3 x Используя разложения (5.104), находим T T (u11 u1 ) 11 dSdt = (u11 u0 ) 1 dSdt+ x1 x1 0 S S T + (u11 u1 )O( + h1 + h3 )dSdt, (5.108) 0 S 5.11. Сходимость приближенного решения к точному T T u u (11 11 ) 1 dSdt = 0 (11 11 ) dSdt+ x1 x1 0 S S T + (11 11 )O( + h1 + h3 )dSdt.

0 S Не приводя громоздких выкладок, заметим, что из (5.44), (5.56), (5.62), (5.63) следуют оценки:

1 (u 0 u0 )2 d 1 Q 1 Cd, 2 (5.109) 1 1 H 11 )2 d ( 1 Q 1 a11 d.

2 H В случае, когда выполнены неравенства (5.76), можно получить оценки, аналогичные (5.108), (5.109), и для остальных слагаемых в правой ча сти (5.107). Из этих оценок и ограниченности приближенного решения следует, что соответствующее любому фиксированному моменту време ни T решение по изложенной выше схеме сходится при, h, 0 к соответствующему этому же T точному решению:

E 0 при, h, 0, где h максимальный линейный размер составляющих область S четы рехугольных элементов;

максимальное из чисел,, = 1,..., 4.

5.11.2. Сходимость приближенного решения к точному в осесимметричной задаче Пусть, = 1, 2, 3;

13, u, u точное решение задачи (5.1), 1 (5.31), (5.32) в области S, состоящей из четырехугольных элементов, на промежутке времени [0, T ] при начальных условиях (5.47) и условиях на границе типа условий (5.46). Разбивая рассматриваемый промежуток времени на конечное число интервалов длительностью и используя на каждом из интервалов изложенную выше процедуру, можно построить приближенное решение на всем промежутке времени [0, T ]. Для раз ности между точным и приближенным решениями, соответствующими t = T, введем энергетическую норму, обозначив квадрат этой нормы 280 Глава 5. Динамика упругопластического деформирования через E:

T (u u1 ) (u u3 ) 1 (u E= H1 H2 u1 ) + (u3 u3 ) + t t 0S +(11 11 )(e e11 ) + (22 22 )(e e22 ) + (33 33 )(e e33 )+ 11 22 +2(13 13 )(e e13 )} dSdt. (5.110) В силу неравенств (5.35) имеем [(11 11 )(u u11 ) + (31 31 )(u u31 )] + 1 x S [(13 13 )(u u13 ) + (33 33 )(u u33 )] dS 0. (5.111) + 1 x Используя (5.31), (5.53) и (5.111), получим для величины E в (5.110) оценку T (11 11 ) (31 31 ) E (u11 u1 ) + (u31 u3 ) + x1 x 0 S (13 13 ) (33 33 ) +(u13 u1 ) + (u33 u3 ) + x3 x (u u11 ) (u u31 ) 1 +(11 11 ) + (31 31 ) + x1 x (u u13 ) (u u33 ) 1 +(13 13 ) + (33 33 ) + x3 x +(L L1 )(u u1 ) + (L L3 )(u u3 ) 1 1 3 (L L1 )(u u1 ) (L L3 )(u u3 )] dSdt, (5.112) 1 1 3 где L = M11 H2 13 M21 H1 22, L = M11 H2 11 M22 H1 22, 1 L1 = M11 H2 13 M21 H1 22, L3 = M11 H2 11 M22 H1 22.

Используя (5.59), запишем (5.112) в виде T + (u31 u3 ) 31 + E+ Q (u11 u1 ) x1 x, 0S u + (u33 u3 ) 33 + (11 11 ) 1 + +(u13 u1 ) x1 x3 x 5.12. Алгоритм решения динамической упругопластической задачи... u u u + (13 13 ) 1 + (33 33 ) 3 + +(31 31 ) x1 x3 x +(1 u1 )L1 + (L1 L1 )u1 + (3 u3 )L3 + (L3 L3 )u3 dSdt. (5.113) u u Из (5.53), (5.65), (5.67), (5.70), используя (5.66), (5.71) и (5.80), можно получить оценки 1 1 d (u 0 u0 )2 d 1 Q 1, 2 H1 H (5.114) 1 1 H H2 11 )2 d ( 1 a11 Q2 d.

2 H Оценки, аналогичные (5.114), имеют место и для остальных слагаемых в правой части (5.113).

Из (5.113), (5.108), (5.114) и оценок, аналогичных (5.108), (5.113) для остальных слагаемых в правой части (5.113), следует, что E 0 при, h, 0, где h максимальный линейный размер четырехугольных элемен тов, составляющих область;

максимальное из чисел,,, ( = 1, 2, 3), 1.

5.12. Алгоритм решения динамической упругопластической задачи для тел вращения при неосесимметричном нагружении 5.12.1. Построение решения для элемента тела вращения Рассмотрим в плоскости Orz произвольный четырехугольник с вер шинами 1, 2, 3, 4 (рис. 5.4). Повернем его вокруг оси Oz на угол.

Получим трехмерный элемент тела вращения, показанный на рис. 5.5.

282 Глава 5. Динамика упругопластического деформирования z x x - - r Рис. 5.4. Элементарный четырехугольник Введем в элементе систему координат 1, 2, 3 :

x = r( 1, 3 ) cos, y = r( 1, 3 ) sin, z = z( 1, 3 ), 4 r= Ni ri, z= Ni z i, i=1 i= 1 N1 = (1 + 1 )(1 3 ), N2 = (1 + 1 )(1 + 3 ), 4 1 N3 = (1 )(1 + ), N4 = (1 1 )(1 3 ), 1 4 2 i + i+ = +, = i i1, 2 ri, zi (i = 1, 2, 3, 4) координаты вершин четырехугольника. Граням элемента соответствуют значения координат i = ±1.

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

1i,i 22 = 1, 2i,i + 12 = 2, 3i,i = 3, (5.115) 2 где g1i = G(r z,j rz r,j ), 12 = gr, 5.12. Алгоритм решения динамической упругопластической задачи... t v t x v t x x v Dj Рис. 5.5. Элемент тела вращения 284 Глава 5. Динамика упругопластического деформирования g2i = G(r z,j z r,j ), 22 = g, g3i = G(rz z,j z r,j ), 32 = gz, i, j = 1, 2, 3;

= 1 при i = 1, j = 3;

= 1 при i = 3, j = 1;

vi i = G Fi, G=r g, g = r,1 z,3 r,3 z,1 ;

t vi, Fi, r, z,, r, z, rz компоненты вектора скорости, вектора массовых сил, тензора напряжений в цилиндрической системе коорди нат. Здесь и в дальнейшем выражение, содержащее два и более двух немых индексов означает, если нет специальной оговорки, сумму по этому индексу от 1 до 3. Немым считается индекс, который отсутствует хотя бы в одном из слагаемых уравнения, (),i обозначает дифференци рование по i.

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

z,3 v1,1 z,1 v1,3 z,3 v2,1 z,1 v2,3 v1,2 v +, er =, 2er = g g r G v1 g z,1 v2,3 z,3 v2,1 v3, +, e = + v2,2, 2ez = (5.116) r G g G r,1 v3,3 r,3 v3,1 r,1 v1,3 r,3 v1,1 + z,3 v3,1 z,1 v3, ez =, 2erz =.

g g Уравнения упругопластического течения примем в виде 1 er = [r ( + z )] + (r ),..., erz = rz + rz, (5.117) E 2µ где E, µ, модуль упругости, модуль сдвига, коэффициент Пуассона соответственно;

функция скоростей деформаций и напряжений.

Задача о построении решения в элементе на отрезке времени [tk, tk+1 ] ставится следующим образом.

Найти функции vi, r, z,, r, z, rz в области :

{ i, [1, 1], = 2[t (tk + tk+1 )/2]/, = tk+1 tk }, удовлетворяющие уравнениям (5.115)–(5.117), условиям на гранях элемента:

(a± vi + b± i1 )1 =±1 = ±, i i i (ci vi + d± i2 )2 =±1 = i, ± ± (5.118) i (e± vi + fi± i3 )3 =±1 = ± (i = 1, 2, 3) i i и начальным условиям {vi, r, z,, r, z, rz }=1 = {vi, r, z,, r, z, rz }. (5.119) 5.12. Алгоритм решения динамической упругопластической задачи... В (5.118) величины a±, b±,..., ± заданные постоянные, причем i i i a± b± 0, c± d± 0, e± fi± 0, (5.120) ii ii i |a± | + |b± | = ± ± ± ± 0, |ci | + |di | = 0, |ei | + |fi | = 0 (i = 1, 2, 3).

i i При построении решения в элементе принимается линейная аппрок симация по времени для скоростей и напряжений 1 1 vi = vi + vi, r = r + r,..., rz = rz + rz. (5.121) На среднем слое по времени напряжения и скорости вычисляются по формулам:

r = [aer + b(e + ez )] + r,..., z = [aez + b(e + er )] + z, 2 (5.122) r = 2µer + r,..., rz = [2µerz + b(e + er )] + rz, 2 vi vi = vi +, 2 t где 4 2 µ a = K + µ, b = K µ, µ =, 3 3 1 + / r z r = r +, = +, z = z +, 1 + /2 1 + /2 1 + / r z zr r = r +, z = z +, zr = zr +.

1 + /2 1 + /2 1 + / В случае упругого деформирования величины r,..., rz значе ния соответствующих напряжений на нижнем слое по времени, а коэф фициенты a, b, µ не зависят от скоростей деформаций. При пластиче ском деформировании величины r,..., rz, a, b, µ зависят от скоростей деформаций и для их вычисления применяется итерационная процеду ра.

В (5.115), (5.122) ускорения и скорости деформаций вычисляются по формулам r v1 1 v2 = 1 + F1, = 2 + F2, (5.123) t g r t g r + v3 1 Fi = 3 + F2, i = ij,j, = Fi d, (5.124) t g 286 Глава 5. Динамика упругопластического деформирования z,3 v11,1 z,1 v13,3 z,3 v21,1 z,1 v23,3 v12,2 v +, er =, 2er = g g r G v g z,1 v23,3 z,3 v21,1 v32, e = 1 + +, v22,2, 2ez = r G g G r,1 v33,3 r,3 v31,1 r,1 v13,3 r,3 v11,1 + z,3 v31,1 z,1 v33, ez =, 2erz =.

g g В (5.123), (5.124) величины ij, vij линейные полиномы, аппрок симирующие напряжения и скорости на среднем слое по времени ij = ij + ij j, vij = vij + vij j, 1 (5.125) удовлетворяющие на гранях элемента условиям (5.118) так, что (a± vi1 + b± i1 )1 =±1 = ±,..., (e± vi3 + fi± i3 )3 =±1 = ±.

i (5.126) i i i i Скорости и напряжения на верхнем слое по времени находятся по формулам vi (t + ) = 2vi vi, r (t + ) = 2r r,..., rz (t + ) = 2rz rz.

(5.127) 5.12.2. Энергетическое тождество Согласно приведенным выше формулам, для построения решения в элементе нужно определить коэффициенты полиномов (5.125). Усло вий (5.126) для этого недостаточно. Для функций (5.121) и полиномов (5.125) выполняется следующее энергетическое тождество:

vi vi d + W d + q1 + q2 = Fi vi d + ij,j vij d, (5.128) t где W = r er + e + z ez + 2(r er + z ez + rz erz ), q1 = (vij vi )ij,j d, q2 = (ij ij )vij,j d, v1 = v1 + 1 S1, d = Gd1 d2 d3 d, r g1 4r 2µ v2 = v2 + 2 2er, v3 = v3 +, r g r g2 4r 1 S1 = ae + b(er + ez ), v1 = (v1 + F1 ), 1 2 2r 1 v2 = (v2 + F2 + r ), v3 = v3, 2 2 2r 5.12. Алгоритм решения динамической упругопластической задачи... 2a 2µ 1 = 1 +, 2 = 1 +, 4r2 4r r 11 = 11 + E1 + Q1 e, Q1 = z,3 b, 2 2 r E1 = [z,3 (aer + be + bez ) 2r,3 µerz ], 22 = + E2 + Q2 e, E2 = g(ber + ae + bez ), Q2 = ga.

2 Структура выражений для 33,..., аналогична выше приведенным.

Подставляя выражения для скоростей и напряжений на среднем слое по времени в q1 и q2, получим:

q = q1 + q2 = (v11 v1 )+ r g1 1 ) + 13 (v13 v + (v12 v1 )+ r g1 2 r g1 + 21 (v21 v2 2 ) + 22 (v22 v )+ r g2 1 r g2 + 23 (v23 v2 2 ) + 31 (v31 v )+ r g 3 r g2 + 32 (v32 v3 3 ) + 33 (v33 v )+ r g 2 r g v11 v12 µ + (11 11 E bz,3 v1 ) + (12 12 E6 + g v2 )+ 1 2 22 2 2 2 r v13 v21 + (13 13 E4 + bz,1 v1 )+ (21 21 E7 + µz,3 v2 )+ 3 2 22 1 2 v22 a v23 + (22 12 E g v1 ) + (23 23 E µz,1 v2 )+ 2 2 2 r 3 2 v31 v33 + (31 31 E5 + br,3 v1 )+ (33 31 E br,1 v1 )+ 1 2 22 3 2 v32 + (32 32 E9 ) d.

2 288 Глава 5. Динамика упругопластического деформирования 5.12.3. Дополнительные уравнения Недостающие уравнения для определения коэффициентов полино мов (5.125) выберем так, чтобы выполнялись следующие условия:

1) вычисление решения сводится к решению системы уравнений, расщепленных на ряд независимых систем;

2) q = q1 + q2 0 (это обеспечивает устойчивость процесса вычис лений, величина q при этом имеет смысл искусственной диссипации);

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

Покажем, что этим условиям удовлетворяют уравнения + 1 v11 v1 d = 0, r g1 (5.129) + 1 2 v 11 (az,3 + µr,3 ) d = 0;

2g2 + 2 v21 v2 d = 0, r g2 (5.130) + 2 µ(z 2 + r2 ) v21 d = 0;

21,3, 2g2 + 3 v31 v3 d = 0, r g (5.131) + 3 (ar2 + µz 2 ) v31 d = 0;

31,3, 2g2 + 1 v12 v1 d = 0, r g1 (5.132) ( + 1 ) g v12 d = 0;

µ 12 5.12. Алгоритм решения динамической упругопластической задачи... + 2 v22 v2 d = 0, r g2 (5.133) ( + 2 ) ga v 22 d = 0;

+ 3 v32 v3 d = 0, r g (5.134) ( + 3 ) g v32 d = 0;

µ 32 + 1 v13 v1 d = 0, r g1 (5.135) + 1 2 v 13 (az,1 + µr,1 ) d = 0;

2g2 + 2 v23 v2 d = 0, r g2 (5.136) + 2 2 v 23 µ(z,1 + r,1 ) d = 0;

2g2 + 3 v33 v3 d = 0, r g1 (5.137) + 3 2 v 33 (ar,1 + µz,1 ) d = 0;

2g2 где 11 = 11 + bz,3 v1, 21 = 21 + µz,3 v2,....

22 Выполнение первого условия следует непосредственно из вида уравне ний (5.129)–(5.137).

290 Глава 5. Динамика упругопластического деформирования Используя соотношения (5.129)–(5.137), выражение для q можно представить в виде 2 1 2 12 11 + q= 1 r g1 1 2 1 13 2 13 11 + 1 2 +...+ 3 1 3 3 2 a g v22 1 1 r 2 v + az,3 +4 4 g 1 r v11 v 4 z,3 b +... d = 0. (5.138) 1 Остальные слагаемые в (5.138) имеют аналогичный вид и не выписы ваются ввиду их громоздкости. Рассмотрим выражение, стоящее в пер вой квадратной скобке подынтегрального выражения. Если 2 1, 2 1, 2 1 то 2 12 13 1 2 + 1 2 + 1 2 1 2 3 13 + 1 2 3 3 2 2 11 11 12 + + 0.

1 2 1 3 2 Так как для реальных материалов b a, то при выполнении неравенств 4 1, 4 2 и вторая квадратная скобка в подынтегральном выра жении (5.138) также неотрицательна:

2 2 a g 1 1 r 2 v11 v22 v11 v az,3 +4 4 z,3 b 4 g 1 r 2 1 g v 1 r 2 v b z,3 2 0.

4 g 1 r Аналогичным образом оцениваются и остальные слагаемые в подынтегральном выражении в (5.138). В результате этих оценок полу чим, что величина q будет неотрицательной (условие 2) при выполнении 5.12. Алгоритм решения динамической упругопластической задачи... неравенств:

1 4, 1 2, 1 4, 2 2, 2 4, 2 2, (5.139) 3 4, 3 2, 3 4.

Запишем уравнения (5.129)–(5.137) в виде:

[vij vi ( + ij )Cij ij,j ]d = 0, (5.140) [ij ij ( + ij )Dij vij,j ]d = 0.

В (5.140) vi, ij зависят от скоростей и напряжений на нижнем слое по времени;

Cij, Dij выражаются через функции формы четырех угольника и их производные по 1, 3.

Уравнений (5.129)–(5.137) и условий (5.126) достаточно для опреде ления коэффициенты полиномов (5.125) и, следовательно, для постро ения решения в элементе тела вращения, представленного на рис. 5.4.

Заметим, что система (5.140) система алгебраических уравнений.

Интегралы по оставлены для удобства дальнейших выкладок.

5.12.4. Сходимость численного решения к решению дифференциальной задачи в энергетической норме Утверждение 1. Каждая из систем (5.140) однозначно разреши ма при любых правых частях.

Доказательство. Умножим первое из уравнений (5.140) на ij,j, второе на vij,j. В результате получим:

{[ij ij ( + ij )Dij vij,j ]vij,j }d+ + {[vij vi ( + ij )Cij ij,j ]ij,j }d = или ( + ij ) (vij,j Dij + ij,j Cij )d = [(ij,j ),j ij vij,j vi ij,j ]d.

Так как в силу диссипативности краевых условий (ij,j ),j d 0, то 292 Глава 5. Динамика упругопластического деформирования для соответствующей однородной системы имеем (vij,j Dij + ij,j Cij )d и, следовательно, однородная система имеет только тривиальное реше ние, а неоднородная система (5.140) имеет единственное решение при любых правых частях.

Утверждение 2. Сумма кинетической и потенциальной энергий при переходе с одного временного слоя на другой не возрастает при отсутствии массовых и поверхностных сил.

Доказательство. Пусть S область, которая состоит из четы рехугольных элементов типа, приведенных на рис. 5.4. Используя энер гетическое тождество (5.128) и условия непрерывности скоростей и на пряжений на гранях элементов, при отсутствии массовых и поверхност ных сил получим:

vi vi + W + q d = 0, (5.141) t где означает суммирование по всем элементам, составляющим об ласть S. Используя аппроксимацию уравнений упругопластического те чения (разделы 5.10, 5.11), и учитывая, что ij = ij + ( /2)( ij /t), можно записать равенство ij 1 1 W d = ij + d + ij ij d. (5.142) 2µ t K t 2µ В силу неотрицательности величины из (5.141) и (5.142) следует нера венство ij 1 W d ij + d, 2µ t K t из которого, в свою очередь, следует неравенство ij vi 1 vi + ij + d + qd 0.

t 2µ t K t Таким образом, (EK + EP ) d + Q 0 (Q 0). (5.143) t 5.12. Алгоритм решения динамической упругопластической задачи... Неравенство (5.143) обеспечивает устойчивость вычислений и огра ниченность приближенного решения.

Утверждение 3. Пусть ij, vi напряжения и скорости, со ответствующие точному решению задачи (5.115)–(5.120) в момент времени T в области S. Разбивая промежуток времени [0, T ] на ко нечное число интервалов длительностью, можно на этом интерва ле построить приближенное решение по изложенной выше процедуре, которое сходится в энергетической норме T (ij ij ) (vi vi ) E= (vi vi ) + (ij ij ) dSdt (5.144) t j 0 S к точному решению, соответствующему моменту времени T.

Доказательство. В силу непрерывности напряжений и скоро стей на гранях элементов и диссипативности краевых условий (5.118) следует неравенство (ij ij )(vi vij )dS 0.

(5.145) S Из (5.115), (5.123) и (5.145) следует T (ij ij ) (vi vij ) E (vij vi ) + (ij ij ) dSdt. (5.146) j t 0 S Используя выражение для q = q1 +q2 и неравенство (5.145), неравенство (5.146) можно привести к виду T T ij vi E+ qdsdt (vij vi ) + (ij ij ) dSdt. (5.147) j j 0 S S Предполагая точное решение гладким, запишем ij ij vi vi = +O( +h1 +h2 ), = +O( +h1 +h2 ), (5.148) j j j j 0 где индекс 0 означает, что соответствующая величина вычисляется при 1 = 2 = 3 = = 0. Тогда выражение в правой части неравенства 294 Глава 5. Динамика упругопластического деформирования (5.147) можно записать в виде:

T ij v + (ij ij ) i dSdt = (vij vi ) j j 0 S T ij vi 0 = (vij vi ) + (ij ij ) dSdt+ j j 0 0 S + (vij vi )O( + h1 + h2 ) + (ij ij )O( + h1 + h2 ). (5.149) Используя уравнения одномерных задач (5.129)–(5.137) и ограничения на шаг по времени (5.139), можно показать, что T T vi )2 dSdt (vij C1 qdSdt, 0 S S T T 0 ij )2 dSdt (ij C2 qdSdt, 0 S S где C1, C2 постоянные, зависящие от констант диссипации и констант, входящих в уравнения упругопластического деформирования. Следова тельно, T 1/ T ij 1 ij dSdt dSdt (vij vi ) j j 0 S S 1/ T T C1 qdSdt + (vij vi )O( + h1 + h2 )dSdt, (5.150) 0 S S 1/ T T vi vi dSdt dSdt (ij ij ) j j 0 S S 1/ T T C2 qdSdt + (ij ij )O( + h1 + h2 )dSdt. (5.151) 0 S S Из последних неравенств и ограниченности решения (утверждение 2) следует утверждение 3.

5.12. Алгоритм решения динамической упругопластической задачи... 5.12.5. Построение решения для тела вращения Пусть область, занятая телом вращения, разбита на элементы типа, представленного на рис. 5.5. Пусть M N число четырехугольников, на которое разбито меридианальное сечение, L количество элементов в окружном направлении при фиксированных значениях r и z. Тогда построение решения во всей области с учетом условий непрерывности полиномов (5.125) на гранях элементов и краевых условий сводится к решению 3 (M N + L N + L M ) систем вида:

pi pi1 i1/2 (qi qi1 ) 2pi1/2 = 0, qi qi1 i1/2 (pi pi1 ) 2qi1/2 = 0, (5.152) a0 q0 + b0 p0 = 0, ak qk + bk pk = k.

Эти системы могут быть решены прогонкой. Если для каждого из элементов выполняется условие i1/2 i1/2 = 1, то решение системы (5.152) вычисляется по явным формулам. Последнее условие наряду с неравенствами (5.139) приводит к ограничению на шаг по времени при вычислении решения по явной схеме:

2 g 2 g min min,, 25(az 2 + µr2 ) ga/r 2 16(µ(z,3 + r,3 ) gµ/4r,3, 2 g r r r,,,, 3µ 9µ a(/2)2 25a µ(/2) 2 5 µz,3 + ar, 2 g 2 g 2 g,,.

16(µ(z 2 + r2 ) gµ/4r 5 µz 2 + ar2 25(az 2 + µr2 ) ga/r,1,2,1,1,1, (5.153) 5.12.6. Определение скоростей на оси вращения При отсутствии сосредоточенных сил на оси вращения должно вы полняться равенство 0 i1 vi1 d = 0. Выражая в этом равенстве ком поненты вектора скорости vi через компоненты vx, vy, vz в декартовой системе координат и учитывая, что оно должно выполняться при лю бых значениях vx, vy, vz получим 2 (11 cos 21 sin ) d = 0, (11 sin 21 cos ) d = 0, (5.154) 0 296 Глава 5. Динамика упругопластического деформирования 31 d = 0.

Равенства (5.154) используются для определения скоростей на оси вра щения. Поясним процедуру вычисления скоростей для случая, когда плоскость x = 0 является плоскостью симметрии. Если ввести обозна чения {k1, vk1 }1 =1,=i1/2 = {pk, q0,i }, k 0,i то 1 2 q0,i = C2 sin i1/2, q0,1 = C2 cos i1/2, q0,i = C и из соотношений (5.154) в случае явной схемы следует p1 = 1/2,i C2 sin i1/2 + fi1, p2 = 1/2,i C2 cos i1/2 + fi2, (5.155) 1 0,i 0,i p3 = 1/2,i C3, fik = 1/2,i q1/2,i, k = 1, 2, 3.

1 k k 0,i Заменяя в (5.154) интегрирование суммированием по формуле прямо угольников, с учетом (5.155) получим L L (fi1 sin i1/2 + fi2 cos i1/2 ) fi i=1 i= C2 =, C3 =.

L L 1 2 cos (1/2,i sin i1/2 + 1/2,i i1/2 ) 1/2,i i=1 i= Если решение находится по явной схеме, то для вычисления ско ростей на оси вращения используется аналогичная процедура. Но для получения соотношений типа (5.155) нужно во всех сечениях = const при фиксированных значениях z выполнить прямую прогонку, начиная с внешней границы по направлению к оси вращения.

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

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

Рассмотрим два соседних элемента. Величины, которые относятся к элементу, расположенному с правой стороны от границы, отметим знаком +, величины, относящиеся к элементу, который расположен с левой стороны, знаком. На рис. 5.6 приведены три типа кра евых условий на границе между элементами (через и обозначены нормальные и касательные напряжения на границе между элементами, через u и v компоненты вектора скорости по нормали и касательной к границе между элементами). Условия 1 соответствуют неразрушенной границе, условия 2 и 3 моделируют разрушение по границам элемен тов. Если имеют место условия 2, то будем говорить, что по границе между элементами произошел разрыв, если же имеют место условия произошло расслоение.

Вычисление решения с учетом разрушения на каждом шаге по вре мени состоит из двух этапов:

1) вычисления соответствующей этому слою времени системы рас положения разрывов и расслоений с учетом возможного захлопывания трещин и образования новых;

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

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

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

298 Глава 5. Динамика упругопластического деформирования - + - + s =s =0, u =u - + - + - + t =t =0, v =v - + s =s = + - - + t =t = - + s =s - + - + u =u * s s, t t* - + + d Dt + -u ) - + (u * ss - + - + * tt, s0 (u --u + - + )Dt d d s0 - + - + * tt, s0 + s + Рис. 5.6. Моделирование процессов хрупкого разрушения 5.14. Примеры численных решений 5.14. Примеры численных решений Одним из способов проверки качества численных схем решения ди намических задач может быть сравнение численных результатов с ана литическими, соответствующими установившемуся режиму.

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

5.14.1. Упругопластическое деформирование кольца Рассматривается плоская деформация тонкого кольца под дей ствием внутреннего давления. Радиус кольца R = 94,1 см, толщина H = 1,03 см. Зависимость давления от времени p(t) = p0 sin(t/t0 ), p0 = 6,89 МПа, t0 = 3,133 мс, 0 = 0,63 ГПа. Принимается кусочно линейная диаграмма одноосного растяжения с модулем упругости E = 200 ГПа и касательным модулем E /E = 0,01;

предел текучести ма териала s = 689 МПа, плотность = 2,86 г/см3. На волновой стадии, когда напряжения и деформации практически постоянны по сечению кольца, можно выписать аналитическое решение как на стадии упруго го, так и на стадии пластического деформирования.

На волновой стадии упругого деформирования (0 t t ) движе ние кольца описывается уравнениями d2 w 1+ w + a2 w = ( 33 ), 22 = E22, 22 =, (5.156) e h dt R где w прогиб кольца;

22 окружное напряжение;

плотность;

+ a2 = E/(R2 );

33 = 0, 33 = p0 sin(t/t0 ). Решение уравнения (5.156) e записывается в виде p0 E w= 2 (sin t sin ae t), 22 = w. (5.157) h[a2 ( t0 ) ] t0 ae t0 R e Момент времени t, при котором сечение кольца переходит в пла стическое состояние, определяется условием 22 (t ) = s. Обозначим 22 (t ) = 22, w(t ) = w.

300 Глава 5. Динамика упругопластического деформирования На стадии упругопластического деформирования (t t t ) уравнения движения кольца имеют вид:

d2 w E p0 sin t + + a2 w = 1, (5.158) p dt2 h t0 R E d22 E dw e, 22 = 22 + (w w ), = dt R dt R E где ap = R2. Решение уравнения (5.158) записывается в виде:

p0 E sin t+ 22 R w = c1 sin ap t+c2 cos ap t+ 1, (5.159) 2 ( )2 ] h[ap t0 E E t константы c1, c2 определяются из условий непрерывности dw dw w(t + 0) = w(t 0), (t + 0) = (t 0).

dt dt Момент t определяется из условия d22 /dt = 0 и при t t снова начинается упругое деформирование:

d2 w p0 d E dw + a2 w = sin t, =, (5.160) e dt h t0 dt R dt E 22 = 22 + (w w ), 22 = 22 (t ), w = w(t ), R p0 R sin t + w, (5.161) w = c1 sin ae t + c2 cos ae t + 2 ( )2 ] h[ae t0 E t константы c1, c2 также определяются из условий непрерывности dw dw w(t + 0) = w(t 0), (t + 0) = (t 0).

dt dt Для приведенных выше параметров задачи t = 1,634 мс, t = 2,945 мс.

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

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

5.14. Примеры численных решений Таблица 5. Сравнение результатов численного решения с аналитическим t, 22 /0 wE/R 103 с аналит. реш. числ. реш. аналит. реш. числ. реш.

0,5 0,0493 0,0478 0,0492 0, 1,0 0,3454 0,3369 0,3452 0, 1,5 0,9233 0,9004 0,9234 0, 1,6 1,0551 1,0288 1,0551 1, 1,7 1,0951 1,0930 1,1798 1, 1,8 1,0963 1,0943 1,3074 1, 1,9 1,0976 1,0956 1,4318 1, 2,0 1,0988 1,0969 1,5523 1, 2,1 1,1000 1,0981 1,6685 1, 2,2 1,1011 1,0994 1,7783 1, 2,3 1,1021 1,1005 1,8811 1, 2,4 1,1030 1,1016 1,9750 1, 2,5 1,1039 1,1026 2,0581 2, 2,6 1,1046 1,1035 2,1289 2, 2,7 1,1051 1,1042 2,1851 2, 2,8 1,1055 1,1048 2,2244 2, 2,9 1,1057 1,1051 2,2447 2, 3,0 1,1023 1,1000 2,2431 2, 5.14.2. Деформирование цилиндрической оболочки под действием синусоидальной нагрузки Ниже приводятся результаты численного решения задачи о де формировании цилиндрической оболочки под действием равномерного внутреннего давления с синусоидальной зависимостью от времени.

Константы материала оболочки и зависимость нагрузки от време ни те же, что и в предыдущем примере. Радиус оболочки R = 91,1 см, длина L = 468,4 см, толщина H = 1,03 см. Параметры оболочки вы браны для сравнения результатов счета такие же, как и в [161, 184].

При численном решении поперечное сечение оболочки разбивалось на 49 элементов по длине и 5 по толщине. Шарнирное закрепление торцов оболочки моделировалось условиями: 11 = 0, v = 0 при = 0 и = L.

302 Глава 5. Динамика упругопластического деформирования w E/(Rs0 ) t, 0 3 1 2 Рис. 5.7. Деформирование цилиндрической оболочки под действием синусо идальной нагрузки Вычисления проводились по неявной схеме с различными шагами интегрирования по времени. Результаты решения с шагами = 0,05 мс и = 0,01 мс практически совпадают.

Результаты счета сравнивались с результатами работ [161, 184].

В [161] решение получено в виде рядов по собственным формам колеба ний с использованием уравнений оболочек типа Доннелла – Муштари.

В [161] максимальный прогиб w E/(R0 ) = 1,79 достигается в момент времени = 2,5 мс. В наших расчетах w E/(R0 ) = 1,72 при = 2,5 мс.

В [184] решение этой же задачи получено численным интегрирова нием уравнений оболочек методом, предложенным в [176]. Результаты численного счета в этой работе представлены в виде графиков. В мас штабе графиков они совпадают с нашими расчетами (рис. 5.7).

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

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

Зависимость нагрузки от времени и ее распределение по внеш ней поверхности арки показаны на рис. 5.8. Были приняты следую щие характеристики материала арки: модуль упругости E = 200 ГПа, касательный модуль E /E = 0,01, предел текучести материала s = = 689 МПа, плотность = 2,86 г/см3. Края арки считались защемлен ными.

Задача решалась по явной схеме. Сечение арки разбивалось на N M криволинейных четырехугольников (N число элементов в радиальном направлении, M в окружном). Шаг интегрирования по времени полагался равным максимальному значению, допускаемому неравенствами (5.97), что обеспечивало минимальную величину искус ственной диссипации при заданном разбиении области. На рис. 5.8, 5. представлены результаты счета при N = 36, M = 18 и N = 54, M = 18.

В первом случае шаг = 1,48107 с, во втором случае = 9,87108 с.

На рис. 5.8 приведено распределение радиального напряжения r в центральном сечении арки для нескольких моментов времени. Сплош ные кривые соответствуют разбиению сечения арки на 54 18 элемен тов, штриховые разбиению на 36 18 элементов. На рис. 5.9 показано развитие зон пластических деформаций.

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

Были приняты следующие характеристики для однородного коль ца: модуль Юнга E = 7 105 кг/см2, коэффициент Пуассона = 0,33, плотность = 2,7 106 кг·с2 /см4, предел прочности на растяжение = 3,45 103 кг/см2, предел прочности на сдвиг = 3 103 кг/см2.

304 Глава 5. Динамика упругопластического деформирования P P sr /P0 t t P 0, R1 R 60o 0, r 0 3 Рис. 5.8. Распределение радиальных напряжений в центральном сечении арки 5.14. Примеры численных решений Рис. 5.9. Зоны пластических деформаций в сечении круговой арки 306 Глава 5. Динамика упругопластического деформирования P P t t R R Рис. 5.10. Схема нагружения кольца Слоистое кольцо состоит из трех слоев. Наружный и внутренний слои из материала с выше приведенными свойствами, свойства мате риала среднего слоя следующие: E = 3104 кг/см2, = 0,3, = кг/см2, = 6 102 кг/см2, = 1,15 106 кг·с2 /см4.

Геометрические параметры колец следующие. Для однородного и слоистого колец внутренний радиус R1 = 100 см, внешний радиус R2 = = 105 см. Для слоистого кольца толщина внутреннего и наружного слоев h1 = h2 = 1 см, среднего слоя h2 = 3 см. К внешней поверхности кольца прикладывалось боковое давление, распределенное по закону p = p(t) cos (рис. 5.10). Зависимость p(t) показана на рис. 5.10.

Для построения численного решения использовалась явная схема с шагом по времени = 8 108 с. В силу симметрии задачи решение вычислялось для половины кольца 0. При = 0 и = стави лись условия симметрии напряженного и деформированного состояния.

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

5.14. Примеры численных решений ) ) Рис. 5.11. Схема разрушения однородного и слоистого колец 308 Глава 5. Динамика упругопластического деформирования На рис. 5.11 показаны трещины нормального разрыва и расслоений в поперечном сечении однородного и слоистого колец для двух момен тов времени. Двойная черта соответствует трещине разрыва, одинарная трещине расслоения.

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

Численные решения, полученные с шагом интегрирования по вре мени = 4 108 с, практически совпадают с приведенными на рис. 5.10, 5.11.

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

Глава Динамика упругопластического деформирования при больших деформациях В ряде твердых тел, например в металлических телах, при произ вольной величине поворотов и деформаций элементов тела компоненты девиатора упругих деформаций величины порядка отношения пре дела текучести на сдвиг к модулю Юнга и, следовательно, малы по сравнению с единицей. Ниже на основе результатов работы [64] фор мулируются уравнения изотропного упругого и идеального упругопла стического деформирования таких тел. Проведено сопоставление полу ченных уравнений с известными [64, 72, 152, 153]. Для простоты записи уравнений ниже рассматриваются только адиабатические процессы де формирования.

6.1. Уравнения упругопластического деформирования при произвольной величине поворотов и деформаций Обозначим через g, g (, = 1, 2, 3) базисные векторы лагранже вой системы координат, порождаемой декартовой системой координат xi с базисными векторами ei = e i (i = 1, 2, 3).

Пусть g g = g g = ij e i e j какой-либо симметричный тензор. Дифференцируя по времени t формулы связи компонент, с компонентами ij находим (d /dt) g g = (Dij /Dt + sj esi + si esj ) e i e j, (6.1) d /dt g g = (Dij /Dt sj esi si esj ) e i e j, 310 Глава 6. Динамика упругопластического деформирования...

где eij = (1/2)(ui /xj + uj /xi ), ui компоненты вектора скорости;

Dij /Dt производная Яуманна [149] Dij /Dt = dij + ki kj + kj ki, ij = (1/2)(ui /xj uj /xi ).

Из (6.1), в частности, следует Dij /Dt + ki ekj + jk eki = eij, (6.2) где ij e i e j = g g, = (1/2)( ), g = g ·g.

g Очевидно, уравнения (6.2) можно записать в виде Dij /Dt + (2/3)eij = aij, = ij ij, (6.3) aij = eij ik ekj jk eki, ij = ij ij.

Так как в системе координат xi, образованной главными осями дефор маций, a11 = (1 211 )e11, a11 = (1 + 233 )e12,..., то при 1 + ij 1 (6.4) тензор с компонентами aij эквивалентен тензору скоростей деформаций eij. Поэтому в случае (6.4) можно вместо (6.3) использовать уравнения Dij /Dt = (1 (2/3))eij. (6.5) Из (6.5) и уравнения неразрывности d/dt + e = следует 2 Dij /Dt = 1 eij, d/dt = 1 e, 3 (6.6) 3/ 1 eij = eij ij e, e = ij eij, = 0 1, 3 где 0 плотность в недеформированном состоянии.

Аналогично находим из (6.1), (6.6), что в случае (6.4) для записи уравнений (6.6) в лагранжевой системе координат можно использовать формулы d /dt g g = d /dt g g = Dij /Dt e i e j, (6.7) = g, = g s g k sk, g s = g ·g s.

6.1. Уравнения упругопластического деформирования... Полагаем, что при изотропном упругом деформировании внутрен няя энергия E функция энтропии S и инвариантов, 2, 3 тензора деформаций 1 E = E(, 2, 3, S), 2 = ij ij, 3 = ik kj ij. (6.8) 2 Из (6.6), (6.8) находим dE 2 E E dS E = 1 e + ik kj eij + T, T=. (6.9) dt 3 2 3 dt S Так как упругое деформирование обратимый процесс, то из закона сохранения энергии ij eij = dE/dt (6.10) и (6.6), (6.9) следуют уравнения = E/, ij = ij + ik kj 2 ij, dS/dt = 0, (6.11) 5/ = 0 1, = E/2, = E/3.

Для того, чтобы уравнения (6.11) были уравнениями упругого де формирования, необходимо, чтобы они однозначно определяли величи ны деформаций по заданным величинам напряжений. Из (6.11) находим ik kj = Aij + Bik kj + Cij, 2 1 B = 2 2 2, C = 2 3 + 2, A= 2 + 3, 3 1 J2 = ij ij = 2 2 + 2 2 + 33, 2 1 2 1 J3 = ik kj ij = 2 2 2 2 2 + 3 + ( 3 + 3 3 )3.

3 3 9 Отсюда и из (6.11) следует, что в случае (E/)/ = 0, 3 2 (2 + 3 ) = 0, (J2 /2 )(J3 /3 ) (J2 /3 )(J3 /2 ) = уравнения (6.11) разрешимы относительно, ij и их можно записать в виде ij = aij + bQij, Qij = ik kj J2 ij, = (, J2, J3, S), (6.12) где a, b можно рассматривать как функции, J2, J3, S.

312 Глава 6. Динамика упругопластического деформирования...

Для того, чтобы уравнения (6.6), (6.12) с заданными функциями a, b, были уравнениями упругого деформирования, достаточно, что бы a, b, удовлетворяли условиям разрешимости уравнений (6.12) от носительно ij, a3 b2 (aJ2 + bJ3 ) = 0, / = 0, (2 /J2 )(3 /J3 ) (2 /J3 )(3 /J2 ) = 0, (6.13) 2 = a J2 + b2 J2 + 3abJ3, 2 2 1 3 = bJ2 a2 J2 b2 J2 + abJ3 + J3 (a3 + b3 J3 ) 3 9 и закону сохранения энергии (6.10) с внутренней энергией E, рассмат риваемой как функция, J2, J3, S. Так как по (6.6), (6.12) имеет место 1 eij = D(aij + bQij )/Dt, (6.14) 1 ij eij = adJ2 /dt + 2bdJ3 /dt + 2J2 da/dt + 3J3 db/dt, то закон сохранения энергии (6.10) при E = E(, J2, J3, S) будет выпол нен, если E/J2 = a + 2J2 a/J2 + 3J3 b/J2, E/J3 = 2(b + J2 a/J3 ) + 3J3 b/J3, E/ = + 2J2 a/ + 3J3 b/ и, следовательно, a, b, должны удовлетворять условиям a/J3 = b/J2, 3 a + 2J2 a/J2 + 3J3 b/J2 = 1 (/J2 + a/), (6.15) 5 3 2(b + J2 a/J3 ) + 3J3 b/J3 = 1 (/J3 + b/).

5 Таким образом, уравнения упругого деформирования в случае (6.4) можно строить, задавая a, b, как функции, J2, J3, S так, что бы были выполнены условия (6.13), (6.15). В частности, эти условия будут выполнены при a=, µ = µ(, S), b = 0, = / J2 d(1/2µ)/d, 2µ (6.16) = (, S), E = J2 /2µ2 + (, S).

6.1. Уравнения упругопластического деформирования... Если относительное изменение плотности мало по сравнению с еди ницей, то величину 1 (2/3) в (6.5), (6.6) и во всех последующих урав нениях можно заменить на единицу.

Если угловые скорости элементов среды и скорости деформаций сдвига величины одного порядка ij eij ui /xj (i = j), то тензор с компонентами eij ik uk /xj jk uk /xi эквивалентен тензору скоростей деформаций. В этом случае производ ную Яуманна в (6.5)–(6.7), (6.14) можно заменить на производную по времени.

6.1.1. Уравнения упругопластического деформирования при малых компонентах девиатора упругих деформаций При упругопластическом деформировании наряду с обратимыми (упругими) деформациями имеют место необратимые (пластические) деформации. Обозначим eij = ee + ep, ij = e + p, (6.17) ij ij ij ij p скорости упругих и пластических деформаций;

ij, p e e где eij, eij ij упругие и пластические деформации. Из (6.1), (6.2), (6.17) следует De /Dt + e esj + e esi = ee = eij ep. (6.18) ij is js ij ij Если компоненты девиатора упругих деформаций малы по сравне нию с единицей 1 + ij 1, ij = e ij ks e, e e (6.19) ij ks то тензор с компонентами ij e e ij = eij is esj ij esi эквивалентен тензору скоростей деформаций eij и уравнения (6.18) можно записать в виде De /Dt + eij ks e = ee = eij ep. (6.20) ij ks ij ij Полагаем, что объем элемента среды изменяется упруго e = ij eij = ks e. (6.21) ks Из (6.20), (6.21) и уравнения неразрывности находим 1 eij = Dij /Dt + p, e (6.22) ij 314 Глава 6. Динамика упругопластического деформирования...

3/ 2 d/dt = 1 e, = 0 1.

3 e Если в качестве зависимости ij от ij использовать уравнения (6.12), (6.16), а в качестве зависимости ep от ij использовать уравнения иде ij ального пластического течения с условием пластичности Мизеса [88], то получим ij = ij /2µ, ep = ij, e (6.23) ij 5/ = / J2 d(1/2µ)/d, µ = µ(, S), = 0 1, = (, S), E = (, S) + J2 /2µ2, T dS/dt = 2 2, T = E/S;

= 0, если f 0 или f = 0, df /dt 0;

(6.24) f = J2 2, 0, если f = 0, df /dt = 0, = (, S) ( предел текучести при сдвиге).

Для среды, удовлетворяющей условию (6.19), уравнения (6.22)–(6.24) образуют систему уравнений упругопластического деформирования, корректную при произвольной величине поворотов и деформаций элементов среды. Уравнения (6.22)–(6.24) совпадают с известными уравнениями [72], [152], [153] в случае, когда µ = const, = const, а относительное изменение плотности мало по сравнению с единицей и, следовательно, 1 (2/3) = (/0 )( 2/3) 1, 0, d/dt e.

При записи уравнений (6.22)–(6.24) в лагранжевой системе коорди нат зависимость компонент девиатора скоростей деформаций от напря жений можно записать, используя (6.6) в виде ee e 1 e = d /dt +, e = /2, остальные уравнения (6.22)–(6.24) остаются без изменений.

Так как по (6.22), (6.24) 2 1 ij eij = dJ2 /µdt + J2 d(1/µ)/dt + 2J2, 3 то условия (6.24), определяющие функцию в (6.23), можно заменить условиями 1 = c/ 2, = 1 ij eij d( /µ)/dt, 2 6.2. Алгоритм решения динамических задач упругопластического... c = 0, если f 0 или f = 0, 0;

c = 1, если f = 0, 0.

Для вычисления по скоростям деформаций малых приращений напряжений за малый интервал времени можно вместо уравнений (6.22)–(6.24) использовать предложенную в [153] процедуру корректи ровки девиатора напряжений. При этом приращения напряжений до корректировки вычисляются по уравнениям (6.11), (6.16).

Используя (6.12), (6.13), (6.15), (6.22), можно для среды с услови ем (6.19) сформулировать уравнения идеального упругопластического деформирования с более общим, чем в (6.22)–(6.24), законом упругого деформирования.

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

Локальная система координат В качестве системы координат Эйлера примем декартову систему координат (x1, x2, x3 ) (рис. 6.1), e1, e2, e3 векторы базиса этой системы координат. Пусть (1, 2, 3 ) система координат Лагранжа, связанная с материальными частицами среды, g1, g2, g3 базисные векторы.

Если r радиус-вектор материальной точки, то справедливы сле дующие соотношения:

r r g1 = r,1 = xi,1 ei, g2 = r,2 = xi,2 ei, g3 = e3, 1 g2 e3 g1 = = (x2,2 e1 x1,2 e2 ), D D e3 g1 g2 = = (x2,1 e1 + x1,1 e2 ), D D D = (g1 g2 ) · e3 = x1,1 x2,2 x1,2 x2,1, (6.25) где g 1, g 2 векторы контравариантного базиса лагранжевой системы координат;

D якобиан перехода от эйлеровых координат к лагран жевым.

316 Глава 6. Динамика упругопластического деформирования...

x x - g g1 - 1 x r e e x Рис. 6.1. Эйлерова и лагранжева системы координат В четырехугольнике с вершинами 1, 2, 3, 4 (рис. 6.1) введем ло кальную систему координат (1, 2 ):


x 1 = Nk x k, x 2 = Nk x k, (6.26) 1 где 1 N1 = (1 + 1 )(1 2 ), N2 = (1 + 1 )(1 + 2 ), 4 1 N3 = (1 1 )(1 + 2 ), N4 = (1 1 )(1 2 ).

4 6.2.2. Уравнения динамической упругопластической задачи с учетом больших перемещений В лагранжевой системе координат (1, 2 ) уравнение неразрывности записывается в виде:

i,i Dv = 0, (6.27) где i = D g ii ti = D ij gj, точка обозначает полную производную по времени t;

i векторы усилий на гранях i = const, ti векторы напряжений на этих же гранях;

ij компоненты тензора напряжений в лагранжевом базисе g1, g2 ;

v вектор скорости точки среды. Обозначим 6.2. Алгоритм решения динамических задач упругопластического... через ij компоненты тензора напряжений в декартовом базисе ei :

ij ei ej = g g. (6.28) Пусть ij проекции векторов j на оси xi. Из равенства i = ij ej = i = D g с учетом (6.25) получим:

ij = Djk i,k.

(6.29) Из (6.27) и (6.28) получим связь между ij и : 11 = 11 x2,2 12 x2,2, 21 = 12 x2,2 22 x1,2, 12 = 12 x1,1 11 x2,1, 22 = 22 x1,1 12 x2,1, или 1 11 = (11 x1,1 + 12 x1,2 ), 22 = (21 x2,1 + 22 x2,2 ), D D (6.30) 1 21 = (21 x1,1 + 22 x1,2 ), 12 = (11 x2,1 + 12 x2,2 ).

D D Умножая (6.29)) на v и интегрируя по ячейке и промежутку времени [t1, t2 ], получим уравнение, выражающее закон сохранения энергии:

1 (v · v + ij eij )Dd1 d2 d = (i · v),i d1 d2 d, (6.31) 1 где 2 tn + tn+ = t, = tn tn+1, 1 e11 =(v1,1 x2,2 v1,2 x2,1 ), e22 = (v2,2 x1,1 v2,1 x1,2 ), D D e12 = (v1,2 x1,1 v1,1 x1,2 + v2,1 x2,2 v2,2 x2,1 ), 2D eij декартовы компоненты тензора скоростей деформаций;

vi де картовы компоненты скорости вектора v.

Уравнение неразрывности в лагранжевых переменных можно запи сать в виде:

(dx1 dx2 ) = (Dd1 d2 ) = 0 или D = 0 D0, (6.32) t t где нижним индексом 0 обозначены значения величин при t = 0.

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

1 + ij 1. (6.33) 318 Глава 6. Динамика упругопластического деформирования...

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

ij = eij jk ki ik kj ij, = e, = 0 3/2, = 1, 3 (6.34) eij = eу + ij, ij ij = 0, eij = eij e ij, ij ij = ij ij, e = eij ij, = ij ij, где ij компоненты тензора упругих деформаций;

ij компоненты тензора скоростей необратимых (пластических) деформаций, 1 vi vj ij = 2 xj xi компоненты тензора скоростей поворотов. В плоском случае ненуле выми компонентами будут:

v1 v 12 = 21 = = x2 x = (v1,1 x1,2 + v1,2 x1,1 v2,1 x2,2 + v2,2 x1,1 ). (6.35) 2D Полные деформации ij удовлетворяют уравнению:

vk vk ij = ki + kj = eij. (6.36) xj xi При построении определяющих соотношений между напряжениями, де формациями и температурой считаем, что внутренняя энергия U явля ется функцией первого и второго инвариантов тензора упругих дефор маций и энтропии S:

U = U (, J2, S), J2 = ij ij. (6.37) Используя закон сохранения энергии ij eij = U (6.38) и соотношения (6.34), получаем:

e + ij eij = eU + ij (eij ki kj kj ki ij )UJ2 + US S.

Из последнего соотношения следует:

= U, ij = ij UJ2, US S = ij ij UJ2. (6.39) 6.2. Алгоритм решения динамических задач упругопластического... Обозначив US = T (T температура), последнее из соотношений (6.39) перепишем в виде:

TS =. (6.40) ij ij Через свободную энергию = U ST соотношения (6.39), (6.40) запи сываются в виде:

=, ij = ij J2, (6.41) + T T + T J2 T J2.

T T T T = ij ij Величина T T T имеет смысл теплоемкости при постоянных деформа циях и обозначается через cV.

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

Если принять, что теплоемкость cV не зависит от и J2, среднее напряжение не зависит от J2, а зависимость ij линейна по ij, то функцию можно представить в виде:

= [() + T () + 2µJ2 + f (T )], где µ модуль сдвига. При этом уравнения (6.41) записываются в виде:

= 5/2 ( + T ), ij = 2µ2/5 ij, (6.42) = 1 ij ij + eT.

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

= K5/2 cV (T T0 ), (6.43) где K = E/(3(1 )) коэффициент объемного расширения;

E мо дуль Юнга;

коэффициент Пуассона;

коэффициент Грюнайзена;

T абсолютная температура;

T0 температура, соответствующая нор мальным условиям (обычно 300 K).

Уравнение (6.43) соответствует зависимостям:

1 3 cV = K 2 + cV T0 0 ln, = 0.

2 2 Следовательно, третье из соотношений (6.42) принимает вид:

T= T e. (6.44) cV ij ij 320 Глава 6. Динамика упругопластического деформирования...

Скорость необратимых деформаций определим соотношениями:

ij = ij. (6.45) Уравнения (6.43)–(6.45) определяют модель упругопластического тела.

При построении модели сжимаемой жидкости считаем, что внут ренняя энергия зависит только от и S:

U = U (, S).

Пренебрегая в (6.34) компонентами девиатора упругих деформаций, по лучим eij = ij, = e, = 0 3/2, = 1.

(6.46) Принимая зависимость (6.43) и полагая ij = 0 (случай идеальной жид кости), уравнение для температуры запишем в виде:

T = T e. (6.47) Соотношения (6.46), (6.47) определяют модель сжимаемой идеаль ной жидкости.

При построении модели политропного газа считаем, что внутренняя энергия зависит только от температуры:

U = U (T ).

Зависимость давления p = от плотности и температуры задается уравнением Клапейрона p = RT, откуда следует S S p = A(S), A(S) = R exp, S = R ln + cV ln T + S0, (6.48) cV где = 1 + R/cV показатель политропы, R универсальная газовая постоянная. При этом в (6.42) = 0, = (0 /)R, и уравнение для T принимает вид:

cV T = RT e. (6.49) Уравнения (6.48), (6.49) определяют модель политропного газа.

Таким образом, задача сводится к нахождению функций v, ij, ij, ij,, T по их начальным значениям при T = 0 и граничным условиям:

B11 v + B21 1 = ± при 1 = ± ± ± (6.50) и B12 v + B22 1 = ± при 2 = ±1, ± ± (6.51) ± где матрицы Bij невырожденные диагональные матрицы с постоян ными элементами.

6.2. Алгоритм решения динамических задач упругопластического... 6.2.4. Алгоритм численного решения Значения величин на нижнем слое tn интервала [tn, tn+1 ] будем обо значать верхним индексом, на среднем слое tn+1/2 = tn + / индексом 0.

Для кинематических соотношений (6.34) на интервале [tn, tn+1 ] при нимаем следующую аппроксимацию:

0 = 0 e0, 0 = 0 (0 )3/2, 0 = 1 0, 2 3 (6.52) 0 00 00 ij ij = ( eij ik kj jk ki ij ).

Разрешая (6.52) относительно компонент девиатора упругих напря жений на среднем слое, получаем:

ij ij =, (6.53) где = (1 + )[(1 + )2 + ( 0 )2 ], 2 0 2 11 = [(1 + )2 + (12 )2 ]11 + (1 + ) 12 12 + 0 0 ( ), 2 12 2 0 2 0 ( ) (1 + ) 12 12 + [(1 + )2 + (12 )2 ] 22, 0 22 = 2 12 0 0 0 12 = 21 = 12 (1 + )(22 11 ) + (1 + )2 12, ij = ij + 0 eij, = 0.

0 2 Уравнения для полных деформаций (6.36) аппроксимируем уравнения ми:

vk vk 0 ij = e0 0 + 0. (6.54) ij ij ki kj 2 xj xi Определяющие соотношения упругопластического тела аппрокси мируем соотношениями:

0 = K(0 )5/2 0 cV 0 (T 0 T0 ), ij = 2µ(0 )5/2 ij, 0 (6.55) T 0 T = 0 0 T 0 e0.

cV 0 0 ij ij Соотношениями 0 = K(0 )5/2 0 cV 0 (T 0 T0 ), T0 T = T 0 e0 (6.56) аппроксимируются определяющие уравнения идеальной сжимаемой 322 Глава 6. Динамика упругопластического деформирования...

жидкости, соотношениями p0 = p A(S 0 )(0 ) e0, T 0 T = T 0 e0 (6.57) 2 уравнения политропного газа.

Векторы скорости и напряжений представляются в виде:

v = v 0 + v 1 = (v1 + v1 )e1 + (v2 + v2 )e2, 0 1 0 (6.58) 0 0 1 0 i = i + i = (1i + 1i )e1 + (2i + 2i )e2.

Уравнение движения аппроксимируем соотношением v 0 = v + 0 (p1,1 + p2,2 ), (6.59) 2 D где p1, p2 полиномы, аппроксимирующие векторы напряжений:

p1 = p11 e1 + p21 e2 = (p0 + p1 1 )e1 + (p0 + p1 1 )e2, 11 11 21 (6.60) p2 = p12 e1 + p22 e2 = (p0 + p1 2 )e1 + (p0 + p1 2 )e2.

12 12 22 Скорости деформаций аппроксимируются выражениями:

1 e11 = (u11,1 x2,2 u12,2 x2,1 ), e22 = (u22,2 x1,1 u21,1 x1,2 ), D D e12 = (u12,2 x1,1 u11,1 x1,2 + u21,1 x2,2 u22,2 x2,1 ), (6.61) 2D 12 = (u12,2 x1,1 u11,1 x1,2 u21,1 x2,2 + u22,2 x2,1 ).

2D Аппроксимацию определяющих соотношений запишем в матричной форме:

i = si + (A1i u1,1 + A2i u2,2 ), i = 1, 2, (6.62) 2D где u1 = u11 e1 + u21 e2 = (u0 + u1 1 )e1 + (u0 + u1 1 )e2, 11 11 21 u2 = u12 e1 + u22 e2 = (u0 + u1 2 )e1 + (u0 + u1 2 )e2.

12 12 22 Для упругопластического тела:

K1 x2 + K2 x2 K3 x1,2 x2,2 2,2 1, A11 = +R, K1 x2 + K2 x K3 x1,2 x2,2 1 1,2 2, K1 x2,1 x2,2 K2 x1,1 x1,2 +RD K3 x1,1 x2,2 K2 D A21 = + K3 x2,1 x1,2 + K2 D K1 x1,1 x1,2 K2 x2,1 x2,2 +RD 0 + Rg12, 6.2. Алгоритм решения динамических задач упругопластического... K1 x2,1 x2,2 K2 x1,1 x1,2 RD K3 x1,2 x2,1 + K2 D A12 = + K3 x1,1 x2,2 K2 D K1 x1,1 x1,2 K2 x2,1 x2,2 RD 0 + Rg12, K 1 x2 + K 2 x2 K3 x1,1 x2,1 2,1 1, A22 = + Rg11, K 1 x2 + K 2 x K3 x1,1 x2,1 1 1,1 2, где 2 4µ K1 = (0 )7/2 K + (1 + )2 + (12 )2, (6.63) 3 µ K2 = (0 )7/2 (1 + )2, µ K3 = (0 )7/2 K + (1 + )2 + 2 (12 ), µ R = (0 )7/2 (1 + ) 12, 11 x2,2 12 x1,2 11 x2,1 + 12 x1, s = s =,, 1 12 x2,2 22 x1,2 12 x2,1 + 22 x1, 2µ 0 5/ 11 = K(0 )5/2 cV 0 (T 0 T0 ) + ( ) 2 0 2 2 (1 + )2 + 11 + (12 )2 22 + (1 + ) 12 12, 0 ( ) 2 12 2µ 0 5/ 22 = K(0 )5/2 cV 0 (T 0 T0 ) + ( ) 2 0 2 2 (12 ) 11 + (1 + )2 + (12 )2 22 (1 + ) 12 12, 0 2 2µ 0 5/2 0 (1 + ) 12 (22 11 ) + (1 + )2 12.


12 = ( ) Формулы для идеальной несжимаемой жидкости получаются, если в (6.63) положить µ = 0.

Полиномы (6.58) должны удовлетворять начальным условиям, а полиномы (6.60) краевым условиям:

B1i ui + B2i pi = i±, ± ± i = 1, 2.

324 Глава 6. Динамика упругопластического деформирования...

6.2.5. Дополнительные уравнения.

Ограничение на шаг по времени при вычислении решения по явной схеме Для определения коэффициентов полиномов (6.58) и (6.60) урав нений (6.59) и (6.62) недостаточно. Необходимо сформулировать еще четыре дополнительных уравнения. Для полиномов (6.58) и (6.60) име ет место разностный аналог закона сохранения (6.31):

1 i0 ui,i )Dd1 d2 d ( v v + +Q= (pi ui ),i d1 d2 d, (6.64) 1 где Q = Q1 + Q2, 1 (ui v 0 )pi,i d1 d2 d, (pi 0 )ui,i d1 d2 d.

Q1 = Q2 = 1 Дополнительные уравнения строим так, чтобы расщепить исход ную двумерную задачу на две одномерные векторные задачи и чтобы величина Q была неотрицательной.

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

(u1 v p1,1 )d1 d2 = M1 p1,1 d1 d2, 2D 2D (p10 s A11 u1,1 )d1 d2 = C1 u1,1 d1 d2, 2D 2D (6.65) (u2 v p2,2 )d1 d2 = M2 p2,2 d1 d2, 2D 2D (p20 s A22 u2,2 )d1 d2 = C2 u2,2 d1 d2, 2D 2D где матрицы Mi, Ci выбираются так, чтобы квадратичные формы Q1 = [(M1 p1,1 · p1,1 ) 2 (p1,1 · p2,2 ) + (M2 p2,2 · p2,2 )]d1 d2 d, 2D Q2 = [(C1 u1,1 · u1,1 ) 2 (A21 u1,1 · u2,2 ) + (C2 u2,2 · u2,2 )]d1 d2 d 2D были неотрицательны.

6.2. Алгоритм решения динамических задач упругопластического... Матрицы A11, A22 представим в виде Aii = Aii ± Rgii.

1 Собственными векторами матрицы A11 являются векторы g 1 и g2, а матрицы A22 векторы g1 и g 2 :

A11 g 1 = 1 g22 g 1, A22 g1 = 2 g11 g1, A22 g 2 = 1 g11 g 2, A11 g2 = 2 g22 g2, где 1 = K1, 2 = K2.

Следовательно, матрицы A11 и A22 приводятся к диагональному виду невырожденными преобразованиями T1 = g 1, g2 и T2 = g1, g 2.

Если в качестве матриц Mi взять матрицы 1i 1 M1 = T1 1 T1, M2 = T2 2 T2, i =, 0 2i а в качестве матриц Ci матрицы 1 C1 = T1 1 T1, C2 = T2 2 T2, 1 11 g22 0 2 12 g11 1 =, 2 =, 0 2 21 g22 0 1 22 g то одномерные задачи (6.65) приводятся к виду (u1 v X1 p1,1 )d1 d2 = 0, 2D (p10 s Y1 u1,1 )d1 d2 = 0, 2D (6.66) 0 (u2 v X2 p2,2 )d1 d2 = 0, 2D (p20 s Y2 u2,2 )d1 d2 = 0, 2D 326 Глава 6. Динамика упругопластического деформирования...

где + 11 0 + 12 T 1, T 1, X 1 = T1 X2 = T + 21 1 + 22 0 ( + 11 )1 g22 T 1, Y1 = T ( + 21 )2 g22 ( + 12 )2 g11 T 1.

Y2 = T ( + 22 )1 g11 При таком выборе матриц для неотрицательности квадратичных форм Q1 и Q2 достаточно выполнения условий:

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

X1 Y1 = I, X2 Y2 = I. (6.68) Равенства (5.197) эквивалентны равенствам:

( + ij )( + ij )k gll = 4D2, i, j, l = 1, 2. (6.69) Если положить ij =, равенства (6.68) записываются в виде:

( + ij )k gll = D, (6.70) а матрицы Xi, Yi становятся диагональными:

+ 1i T 1 = 2 I, Xi = Ti + 2i i ( + 1i )1 gjj T1 = D2 I.

Yi = T i 0 ( + 2i )2 gjj Неравенства (6.67) и равенства (6.70) дают ограничение на шаг по времени при вычислении решения по явной схеме:

D min min. (6.71) k gii i,k 6.3. Квазиодномерная модель процессов высокоскоростного соударения 6.3. Квазиодномерная модель процессов высокоскоростного соударения Задачи, излагаемые в этом разделе, можно рассматривать как те стовые. С другой стороны, они представляют определенный и практи ческий интерес.

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

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

Одной из первых одномерных моделей высокоскоростного соударе ния является гидродинамическая модель Лаврентьева [113] для оцен ки глубины проникания цилиндрического стержня в полубесконечную преграду. Согласно этой модели (l/l0 )/(пр /ст )1/2 = const, где l0 ис ходная длина стержня;

l глубина проникания;

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

328 Глава 6. Динамика упругопластического деформирования...

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

Наиболее полно существующие одномерные и двумерные модели высокоскоростного соударения изложены в [144].

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

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

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

6.3.1. Квазиодномерная модель для оценки предельной глубины проникания ударника в преграду Ударник рассматривается как абсолютно жесткий цилиндр радиуса R, высоты H и массы M. Преграда изотропная упругопластическая среда.

Рассмотрим в преграде область D (расчетную область), находя щуюся под ударником и ограниченную боковой поверхностью A1 A (рис. 6.2).

6.3. Квазиодномерная модель процессов высокоскоростного соударения z R H A D r A Рис. 6.2. Расчетная область На боковой поверхности A1 A2 в процессе соударения действуют ка сательные и нормальные напряжения rz и r, изменяющиеся во време ни. Для определения напряженно-деформированного состояния в обла сти D достаточно знать закон изменения этих напряжений во времени и по пространственной координате. Для определения этих напряжений необходимо решать двумерную задачу. При построении квазиодномер ной модели воздействие на область D через поверхность A1 A2 со сторо ны остальной части преграды будем моделировать краевыми условиями определенного вида.

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

Радиальное напряжение r полагаем равным нулю. Очевидно, что среди всех возможных вариантов краевых условий для r это условие максимально ослабляет сопротивление преграды прониканию ударни ка. При отсутствии радиальных напряжений на боковой поверхности A1 A2 появляются скорости в радиальном направлении, которые приво дят к истечению части материала из расчетной области. Часть материа ла, вытекающая через поверхность A1 A2, исключается из дальнейших расчетов. Поскольку эта масса обладает скоростью в осевом направле нии, то исключение ее из расчета означает, что из системы изымает 330 Глава 6. Динамика упругопластического деформирования...

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

Краевые условия на поверхности A1 A2 для uz скорости в осевом направлении и rz касательного напряжения формулируются следу ющим образом.

В начале принимается, что uz на боковой поверхности равно ну лю. Если при этом условии |rz | ( предел текучести на сдвиг или сопротивление срезу), то решение вычисляется с условием uz = на A1 A2. Если же |rz |, то решение вычисляется при условии |rz | = ± ;

при этом знак rz принимается тем же, какой был у rz при нулевом значении uz. Заметим, что при вычислении решения не требу ется итераций для реализации таких краевых условий (раздел 5.12).

6.3.2. Квазиодномерная модель для оценки предельной глубины зоны тыльных отколов Одним из основных факторов, определяющих тыльные отколы, яв ляется амплитуда давления на фронте ударной волны. Очевидно, что амплитуда давления будет максимальной, если ударник в процессе про никания считать недеформируемым абсолютно жестким телом, а среди вариантов краевых условий для радиальных составляющих скорости и напряжения на поверхности A1 A2 принять условие ur = 0.

Краевое условие для касательного напряжения rz принимается та ким же, как и 5.13.

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

Исходные данные задачи были следующими: начальная скорость V0 = 2 км/с, радиус ударника R = 0,5 см, высота ударника H = 1,0 см, /0 = 0,4, 0,6, 0,8, 1,0 (0 плотность материала ударника), модуль Юнга E = 200 ГПа, касательный модуль E = 0, коэффициент Пуассо на = 0,3, предел текучести при сдвиге s = 500 МПа, коэффициент 6.3. Квазиодномерная модель процессов высокоскоростного соударения 0, l/H 0, 0, r0 / r 0,4 0,8 1, Рис. 6.3. Зависимость глубины проникания от соотношения плотностей мате риалов ударника и преграды Грюнайзена = 2,0, теплоемкость cv = 8,96 106 см2 /(с2 К), начальная температура T0 = 300 К, толщина преграды H = 5 см, число конечных элементов в столбце N = 50.

На рис. 6.3 по оси ординат отложена величина l/H (l глубина проникания ударника), по оси абсцисс величина 0 /. Крестиками отмечены точки, полученные в предположении, что на границе расчет ной области радиальная скорость равна нулю, кружочками точки, полученные в предположении, что на границе расчетной области ради альное напряжение равно нулю. Сплошная кривая соответствует гид родинамической модели Лаврентьева. При указанных входных данных оба предельных варианта условий на границе расчетной области при водят к близким результатам.

332 Глава 6. Динамика упругопластического деформирования...

6.3.4. Зависимость глубины проникания цилиндрического стержня в преграду от скорости соударения В [150] приведены экспериментальные данные по скоростному уда ру цилиндрических стержней из мягкой стали по мишени из того же материала. Там же приведено сравнение экспериментальных данных с расчетами по модифицированной гидродинамической модели [150]. Ди намический предел текучести материала Yp = 300 МПа. В модифициро ванной гидродинамической модели присутствуют два параметра R и Y, характеризующие прочностные свойства ударника и преграды, которые определяются при некоторых предположениях по адиабате Гюгонио и пределу прочности. Их определение весьма условно и отношение R/Y в [150] предлагается выбирать так, чтобы расчетные данные соответ 1. D/L R/Y = 0. R/Y = R/Y = 0. V* 0 2 Рис. 6.4. Зависимость глубины проникания от скорости ударника 6.4. Задача теплопроводности ствовали экспериментальным. В [150] расчетные зависимости глубины проникания от скорости ударника приведены для ряда значений R/Y.

На рис. 6.4 представлены зависимости глубины проникания от ско рости ударника. Треугольниками, квадратами и кружочками нанесены экспериментальные точки [150]. Треугольники соответствуют отноше нию длины ударника к его диаметру l/d = 5, квадраты l/d = 7, кружочки l/d = 10. Поскольку кривые на рис. 6.4 (за исключением штриховой) взяты из работы [150]), мы сочли целесообразным оставить на нем единицы измерения, принятые в [150]. Сплошные кривые соот ветствуют расчетам по модифицированной гидродинамической модели [150], штриховая расчетам по схеме пункта 2 настоящего раздела.

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

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

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

Решению задач теплопроводности в областях, составленных из нерегулярных элементов посвящены, например, работы [66, 103]. В [103] решение строится вариационно-разностным методом, в [66] методом факторизованных смещений, аналогичным методу переменных на правлений. В этих работах схема вычислений неявная. В этом раз деле излагается алгоритм решения плоских задач теплопроводности в областях из произвольных выпуклых четырехугольных элементов, ос нованный на аппроксимации температуры несколькими линейными по линомами. Введение нескольких аппроксимаций позволяет обеспечить непрерывность средних величин потоков тепла и температуры на гра нях элементов и так же, как и в динамических задачах упругопласти ческого деформирования, произвести расщепление исходной задачи на одномерные, решения которых можно вычислить скалярными прогон ками, либо по явной схеме. При этом ограничение на шаг по времени, 334 Глава 6. Динамика упругопластического деформирования...

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

6.4.1. Уравнения теплопроводности для четырехугольного элемента Задачу теплопроводности на промежутке времени t [t0, t0 + ] в четырехугольном элементе можно сформулировать в следующем виде:

q T T q i + cik a + = 0, = 0, t (6.72) (k T + k q k )k =±1 = ±, ± ± T |t=t0 = T0, k = 1, 2, k где g = |e1 e2 |, q i = gq · e i, cik = ge i · e k ;

a = c g, плотность;

c теплоемкость;

коэффициент теплопровод i ности;

e, e i базисные векторы косоугольной системы координат, в которой четырехугольному элементу соответствует квадрат i i [1, 1];

, T0, ±, k, k ± ± заданные функции, причем k + ± ± k k 0, k k 0, k 0, |k + k | = 0. (6.73) В (6.72) и ниже слагаемое, содержащее два и более одинаковых немых индексов, означает сумму по этому индексу от 1 до 2. Под немым ин дексом понимается индекс, который отсутствует хотя бы в одном из слагаемых уравнения.

При построении приближенного решения задачи (6.72) примем ап проксимацию температуры и потоков тепла в виде:

qk T k T + k = 0, q i +ik k = 0, a c t 0 2 0 T =T + T, = (t t0 ), q =q + q i i,i i (6.74) 0 1 qk T i =T i + T i i, T T i = ik k, T |=1 = T0.

( ± T k + ± q ± )k =±1 = ±, i, k = 1, 2, (6.75) k k k k k k k где T, T i, q (k = 0, 1;

i = 1, 2) i постоянные в пределах элемента ве cik, T осредненные по = { 1, 2 [1, 1]} величины личины;

a,, 6.4. Задача теплопроводности ± ± k a,, cik, T0 ;

k, k, ± осредненные по соответствующей грани эле мента и промежутку времени [1, 1] величины k, k, ± ;

ik ± ± k параметры, значения которых удовлетворяют условиям:

2 11 22, 11 22, 11 0, 22 0, 11 + 22 = 0. (6.76) 12 В (6.74) полином T характеризует осредненную по элементу темпера туру, полиномы T i, q i осредненные по времени [1, 1] средние i значения в сечениях = const температуры и потоков тепла.

Из (6.74) следует зависимость между величинами T i |i =±1, qi =± i средними значениями температуры и потоков тепла на гранях элемента.

Условимся называть эту зависимость уравнениями теплопроводности элемента.

Условия (6.73), (6.76) обеспечивают возможность однозначно опре делить из уравнений теплопроводности элемента и граничных условий (6.75) коэффициенты полиномов T i, q i. Температура элемента в момент времени t0 + вычисляется по формуле k 0 + ( q ).

T (t0 + ) = T (6.77) k 6.4.2. Уравнения теплопроводности для системы элементов Уравнения теплопроводности для системы элементов состоят из уравнений теплопроводности каждого из элементов, условий непрерыв ности T k, q k на общих гранях соседних элементов и условий (6.75) на гранях, составляющих границу системы элементов.

Нетрудно установить, что решение уравнений теплопроводности си стемы элементов должно обладать энергетическим свойством T i T k qi qk T + cik i + ik i + k + i (T i q i ) d = 0, (6.78) a k t где знак означает сумму по всем элементам, входящим в состав си стемы. Из (6.78) следует, что в случае (6.73), (6.76) можно однозначно определить из уравнений теплопроводности для системы элементов ко эффициенты полиномов T i, q i во всех элементах, вычислить по формуле (6.77) температуру элементов в момент времени T0 + и последующими шагами по времени вычислить температуру в любой момент времени.

Можно показать, что в случае, когда точное решение достаточно гладкое, получаемое по изложенной методике численное решение при 336 Глава 6. Динамика упругопластического деформирования...

h, 0 сходится к точному решению (h характерный размер эле мента).

6.4.3. Итерационное решение уравнений теплопроводности для системы элементов Полагая (6.74) 12 = 21 = (6.79) 2a и считая величины c12 (T 2 / 2 ), c21 (T 1 / 1 ) равными их значениям 12 2 2 21 1 [ (T / )], [ (T / )] на предыдущей итерации (на первой ите c c рации [ (T / 2 )], [21 (T 1 / 1 )] считаются равными нулю), урав 12 c c нения теплопроводности элемента можно записать в виде:



Pages:     | 1 |   ...   | 4 | 5 || 7 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

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