CC BY
47
Вычислительные технологии
Том 9, № 4, 2004
ИНТЕГРИРОВАНИЕ ЗАДАЧИ ДИНАМИЧЕСКОГО УПРУГОПЛАСТИЧЕСКОГО ИЗГИБА АРМИРОВАННЫХ СТЕРЖНЕЙ ПЕРЕМЕННОГО ПОПЕРЕЧНОГО СЕЧЕНИЯ ОБОБЩЕННЫМИ МЕТОДАМИ РУНГЕ - КУТТЫ*
Ю. В. Немировский, А. П. Янковский Институт теоретической и прикладной механики СО РАН,
Новосибирск, Россия e-mail: [email protected]
A problem on the elasto-plastic dynamic core bending is formulated. A method of numerical integration of the problem based on a generalization of the Runge — Kutta methods is developed. Comparison of the obtained numerical solution with the known analytical solution in the elastic case has been done. The efficiency of the suggested numerical procedure is shown on the calculations of the non-elastic dynamics of isotropic and reinforced cores of sections with constant or variable cross-section.
Введение
Стержни составляют основу силовых элементов многих защитных ограждений и ответственных элементов судостроительной, машиностроительной и авиационной техники, объектов стройиндустрии. При воздействии взрывных нагрузок высокой интенсивности их повреждаемость определяет возможность дальнейшего функционирования указанных объектов. Поэтому проблема динамического расчета таких элементов конструкций является одной из важнейших в механике деформируемого твердого тела. Большинство существующих на данный момент решений основано на модели идеального жесткопластического тела. Их обзор и анализ приведены в [1]. Проблема упругопластического изгиба при оценке повреждаемости в основном рассматривалась применительно к однородным изотропным стержням [2, 3]. Настоящая работа посвящена разработке численного метода интегрирования задач динамики неупругих армированных стержней переменного сечения и в этом смысле является продолжением исследований, опубликованных авторами в [4], где разработан и апробирован на решениях задач математической физики один подход к обобщению методов Рунге — Кутты на многомерные случаи. Авторам представляется, что предлагаемый в настоящем исследовании обобщенный метод Рунге — Кутты может быть эффективно использован в дальнейшем при решении задач динамики неупругого деформирования
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 02-01-00115).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.
неоднородных композитных пластин и оболочек со сложными структурами армирования, изучение которых находится пока в зачаточном состоянии [5].
1. Постановка задачи неупругой динамики стержней
Рассмотрим прямолинейный стержень длиной I. Поперечное сечение стержня может быть переменным и имеет две ортогональные оси симметрии. Свяжем со стержнем декартову систему координат (х, у, г): ось х совместим с осью центров тяжести поперечных сечений, оси у, г направим по осям симметрии поперечных сечений. В продольном направлении стержень усилен N семействами армирующих элементов. В поперечном сечении структура армирования является регулярной и симметричной относительно направлений у, г, в частности квазиоднородной. (Под усиливающими элементами можно понимать армирующие волокна или тонкие, регулярно чередующиеся слои, т. е. стержень может иметь волокнистую или слоистую структуру.)
Поперечная распределенная нагрузка действует в направлении г, поэтому при сделанных выше предположениях о геометрии и структуре армирования стержня в случае малых прогибов будет реализоваться поперечный изгиб стержня в плоскости (х, г). Уравнение динамического изгиба такого стержня имеет вид [2]
м (х,г),хх = -р(х,г) + я(х),ш(х, ¿),й , (1)
где
М(x,t) = JJ агйЯ, а(х,у,г,Ь) = ^2 а^к, Я(х) = ^ Рк // Шк^Я; (2)
р к к р
М — изгибающий момент; а — осредненное продольное напряжение в стержне; Я(х) — площадь поперечного сечения стержня; р — погонная поперечная нагрузка, действующая в направлении г; ,ш — прогиб в направлении г; рк — объемная плотность материала к-й фазы композиции (при к = 0 здесь и далее будем понимать связующее, а при к =
1, 2,..., N — к-е семейство арматуры); а к — продольное напряжение в к-й фазе композиции;
Шк(х,у,г) — интенсивность (плотность) армирования к-й фазой композиции ^^ Шк = 1^ ;
Е(х) — погонная масса стержня, в случае квазиоднородного армирования (шк = Шк(х)) имеющая выражение
Я(х) = Я (х)(х), р(х) = £ ркШк (х), (3)
к
р — приведенная объемная плотность композиции; суммирование производится по указанному индексу от 0 до N; нижний индекс после запятой означает частное дифференцирование по пространственной переменной х или времени t.
Стержень подчиняется гипотезе Кирхгофа, поэтому при малых прогибах связь между продольной деформацией е и прогибом ,ш определяется так [6]:
е(х,г^) = -гт,хх (х^), |г| < Н(х)/2, (4)
где Н(х) — высота поперечного сечения.
Связь между деформацией е и напряжениями ак определяется диаграммой растяжения — сжатия материала к-й фазы композиции. Будем предполагать, что диаграммы растяжения — сжатия всех фазовых материалов имеют линейное упрочнение (при необходимости могут быть использованы и другие диаграммы а к ~ е, в рамках настоящего
исследования это непринципиально). Тогда связь между Ok, £ до начала разгрузки определяется равенствами [6]
ok Ek£} |£| < £sk 1 |z| < hk!2,
(5)
Ok = Sgn(£)Osk + Esk (£ - sgn(£)£sk), £sk < |£ < £yk, hk/2 < |z| < H/2,
где
Esk (oyk Osk) / (£yk £sk) 1 £sk Osk/ Ek 1 £yk 8k + Oyk/ Ek (k 0, 1 1 N) 1 (6)
Ek 1 Esk — модули упругости и упрочнения материала k-й фазы композиции соответственно; osk 1 oyk — пределы текучести и временного сопротивления k-й фазы; £sk — предельная упругая деформация k-й фазы; £yk — деформация разрушения, соответствующая значению временного сопротивления oyk материала k-й фазы; 8k — остаточное относительное удлинение при разрыве материала k-й фазы; ±hk/2 — аппликаты границ между упругим и неупругими слоями материала k-й фазы (hk(x} t) > 0 — толщина упругого слоя материала k-й фазы), определяемые из условия |£ = £sk, т. е. с учетом (4), (6) — из равенства
hk/2 - £sk/ | W}XX | - O sk/ \Ekw }XX | 1 k - °) 11 ■■■} N- (7)
(Если из (7) получаем hk > H, то в (5) следует формально положить hk — H, так как в сечении x в k-й фазе композиции пластическое состояние не достигается. Если материал k-й фазы является упругохрупким, то в (6) нужно формально принять osk — oyk и 8k — 0.)
Если в сечении x имеет место первичная разгрузка (разгрузка после первичной пластичности (5)), то напряжения в этом сечении определяются равенствами [6]
(8)
ак = Еке, \г\ < Кк/2, а к = а*к + Ек (е — е^), Н*к/2 < \г \ < Н /2, к = 0,1,..., N,
где, согласно (4) и второму равенству (5),
ек = —г'ш*хх , ак = 8§а(ек)а.3к + Е,3к(ек — sgn(ек)еsk), Кк/2 < \г\ < Н/2; (9)
/ш*хх , ек, ак — кривизна стержня, деформация и напряжение в пластическом слое к-й фазы композиции в момент времени tк начала разгрузки ('ш*хх = /ш,хх (х, tк), ек = е(х, г, tк), ак = ак(х,г^к)); Ккк(х) — толщина упругого слоя в к-й фазе композиции в сечении х в момент tк, определяемая из (7) при замене Кк ^ Ккк и гш,хх ^ ш*х:х (Ккк(х) = Кк(х^к)).
Первичная разгрузка в сечении х начинается в момент времени tк = tк(x), характеризуемый соотношениями е,1 (х,г^к) = 0, е,и (х,г^к) = 0 (\г\ < Н/2) или с учетом (4)
Ы,ххЬ (х^к) = 0, -Ш,ххЫ (х^к) = 0, (10)
где 1ш,ххц (х^к) можно определить, дважды продифференцировав по х уравнение (1). Если в процессе разгрузки в сечении х выполняются соотношения
-Ш,хх = ,Ш*хх , ^п('Ш*хх )'Ш,ххЬ > 0, е = ек, sgn(ек)е,t > 0, (11)
то в этом сечении вновь будет развиваться первичная пластичность и от равенств (8), (9) следует перейти к (5) — (7). Если же в процессе разгрузки в сечении х выполняются соотношения
\е — ек\ = 2е8к, \е\ < \ек\ при \г\ = Н/2, (12)
то в крайних “волокнах” стержня (\г\ = Н/2) в к-й фазе композиции с учетом идеального эффекта Баушингера [6] достигается вторичная пластичность (противоположная по знаку первичной пластичности). При выполнении неравенств
Щп(е к )(е — е к ) < 2евк, Щп(Ш*хх )(/Ш,хх —/^*хх )Н/2 < 2евк = 2а вк / Ек (13)
в сечении х в к-й фазе композиции развивается вторичная пластичность. Напряжения в таком сечении определяются равенствами
ак = Еке (\г\ < Ккк/2), ак = ак + Ек(е — ек) (Ккк/2 < \г\ < Нк/2), (14)
ак Евке й^(ек')(авк Евкевк) (Нк/2 < \г\ < Н/2')1
где ±Нк /2 — аппликаты границ между слоями первичной и вторичной пластичности в к-й фазе композиции (или между упругим слоем и слоями с вторичной пластичностью, если в сечении х слои с вторичной пластичностью полностью перекрывают слои с первичной пластичностью, т. е. при 0 < Нк(х,Ь) < Ккк(х)). Функция Нк(х^) > 0 определяется из предельного равенства первого неравенства (13), откуда с учетом (4) имеем
Нк (хА) 2евк 2авк Н . .
^ ; 7-----вк-----7-----------------вк-< —, к = 0,1,...^. (15)
2 \т,хх (х, ^ — т*хх (х) \ Ек \ш,хх (х, ^ — т£х (х) \ 2
Если в сечении х имеет место “вторичная” разгрузка (разгрузка после вторичной пластичности (14)), то напряжения в этом сечении определяются равенствами
ак = Еке (\г\ < Ккк/2), ак = ак, + Ек (е — е к) (Кк к/2 < \г\ < Нк к/2), (16)
ак = акк + Ек (е — ек к) (Нк к/2 < \г\ < Н/2),
где, согласно (4) и третьему равенству (14),
е к к = —г-ш*х*х , Нк к (х) = Нк (х^ к к), ак = Евке к к — й^(е к )(авк — Евк евк) (Н кк/2 < \г\ < Н/2)1
(17)
ш,х*х , екк, акк — кривизна стержня, деформация и напряжение в слоях вторичной пластичности к-й фазы композиции в момент времени tк к начала вторичной разгрузки (/ш*хх = ы,хх (х^ к к), е к к = е(х,г^ к к), акк = ак (х,г^ к к)); функция Нк к (х) определяется из (15) при замене Нк ^ Н кк и гш,хх ^ "ш*хх.
Вторичная разгрузка в сечении х начинается в момент времени tкк = tкк(х), характеризуемый соотношениями (10), в которых следует заменить tк на tкк ^ к к > tк), и в к-й фазе композиции продолжается до момента возникновения “третичной” пластичности (противоположной по знаку вторичной пластичности), определяемого с учетом эффекта Баушингера равенством
^п(е к )(е — е к к) = 2евк при \г\ = Н/2 или с учетом (4) — равенством
sgn(w*x*x )(,ш,хх —ы*х*х )Н/2 = 2евк = 2авк/Ек, к = 0,1,..., N.
(18)
Если в процессе вторичной разгрузки в сечении х выполняются соотношения
'Ш,хх = т,*х*х , sgn(w*x*x -т,*хх )т,ххг > 0 (е = е**, sgn(е**- е*)е,г > 0), (19)
то в нем вновь будет развиваться вторичная пластичность и от равенств (16), (17) следует перейти к (14), (15).
Если из (7), (15) следует Нк < к * к или Н* к < к * к, то в (14) или (16), (17) нужно соответственно формально принять к * к = Нк или к * к = Н* к.
Можно выписать связь между напряжениями а к и деформацией е в сечении х стержня и при развитии третичной пластичности и последующей “третичной” разгрузке этого сечения, а также и при последующих циклах знакопеременной пластичности. Не будем этого делать, так как это не является целью настоящего исследования, тем более что с увеличением числа циклов знакопеременной пластичности такая связь становится все более громоздкой.
Подставим выражения для напряжений (5), (8), (14), (16) во второе равенство (2), а результат — в первое равенство (2). Тогда после вычисления интеграла с учетом соотношения (4) получим следующее представление для момента М в сечении х:
М = А(х, w,xx ^*хх V,хх +Б(х^,хх ^*хх ^,*хх ), (20)
где А, В — определенные функции указанных аргументов.
Если поперечное сечение стержня является прямоугольным и имеет высоту Н(х) и ширину Ь(х) (в направлении у), то при квазиоднородном армировании поперечного сечения (шк = Шк(х)) функции А, В имеют следующие выражения: при упругом деформировании всех фаз композиции или развитии в них первичной пластичности
B(x, w,xx ) = -sgn(w,xx ) Y “k(x) (ask - Esk£sk) (H2 - k\) ,
k
где kk = kk(w,xx) определены в (7); при первичной разгрузке и вторичной пластичности A(x, w,xx , w,*xx )= ^ Е “k(x) \HlEk + Esk (H3 - Hf)] ,
(x) rh3 E E H3 ^
12
k (21)
A(x,w,xx ) = uk(x) HEk + Esk (H3 - kk)] ,
В(x, w,xx , w,Xx ) = Uk(x) [(Ek — Esk) (H3 - k3k) w*xx +
k
+3sgn(w,xx ) (<Jsk — Esk£sk) (H2 — 2Hk + klk)] , (22)
где k*k = kk(w*xx ), Hk = Hk(w,xx ,w*xx ) определены в (7), (15); при вторичной разгрузке
A(x) = H 3(x) E uk (x)Ek, В (x,w,lx ,w" ) = E “k (x) {[(H3 — Hi„ ) w" +
12 k 12 k
+ (H*k — k* k ) w*xx] (Ek — Esk) + 3sgn(w,xx ) (ask — Esk £sk) (H2 — 2Hîk + k2k )} , (23)
где k * k = kk (w*xx ), H* k = Hk (w*xx ,w,x*x ) определены в (7), (15).
Если из (7) следует kk > H, то в (21) нужно формально принять kk = H (упругий изгиб k-й фазы композиции). Если из (15) получаем Hk > H, то в (22) нужно положить
Нк = Н (первичная разгрузка к-й фазы). Если из (7), (15) следует к * к > Нк или к * к > Н* к, то в (22) или (23) необходимо соответственно принять к*к = Нк или к*к = Н*к (области вторичной пластичности в сечении х перекрывают области первичной пластичности).
После подстановки выражения для момента (20) в уравнение (1) получим уравнение движения в прогибах
w,tt = я(х,Ь) — Ь(х; w), Ь > 0, 0 < х < I, (24)
где
1 д2 р(х Ь)
1(х] 'Ш) = ~Щх) дх2 [А (х’ W,xx , W,Xx) 'Ш,хх —В (х^хх , 'Ш*хх , w*x*x)} , Я(х, Ь) = щх) (25)
Я(х) определена в (2) или (3).
Уравнение (24) с учетом (25) является параболическим уравнением в частных производных второго порядка по времени Ь и четвертого порядка по пространственной переменной х. Для однозначного интегрирования этого уравнения необходимо использовать начальные
w(x, 0) = 1щ(х), w,t (х, 0) = Уо(х), Ь = 0 (26)
и общеизвестные граничные условия [7], задаваемые на концах стержня х = 0, I. При этом
следует учитывать выражения для момента (20) и перерезывающей силы
Q М ,х (A(x, W,xx , W ,хх )‘^^,хх В(х, 'Ш,хх , w,xx , W ,хх )) ,х ■ (27)
В частности, при шарнирном опирании концов стержня имеют место граничные условия
1^(0, Ь) = 1^(1,1) = 0, IV,хх (0,Ь) = IV,хх (1,Ь) = 0, Ь > 0, (28)
а при защемлении концов
w(0,t) = w(l,t) = 0, w,x (0,t) = w,x (l,t) = 0, t > 0. (29)
(Стержень может иметь опоры не в концевых точках, а в промежуточных, в рамках настоящего исследования это не принципиально.)
Если в некоторых сечениях (x = xc) стержня заданы сосредоточенные перерезывающие силы Qc или моменты Mc, то дополнительно необходимо использовать условия сопряжения решения
w(xc - 0,t) = w(xc + 0,t), w,x (xc - 0,t) = w,x (xc + 0, t),
(30)
Q(xc - 0,t) + Q(xc + 0,t) = Qc(t), M(xc - 0,t) + M(xc + 0,t) = -Mc(t) (t > 0),
где для задания M и Q нужно учесть выражения (20), (27).
Так как используемая в настоящем исследовании диаграмма растяжения — сжатия фазовых материалов (5) имеет излом в точке £ = £sk, то в этой точке производная dak/d£ не определена. Поэтому производные по x от функций A и B в (25), (27) не определены в сечениях x = xc(t), в которых в крайних “волокнах” (|z| = H/2) k-й фазы композиции достигается впервые предел текучести, т. е. £sk = H lw,xx (xc,t)l /2. (Такие сечения xc являются границами между областями чистоупругого и упругопластического деформирования k-й фазы композиции.) В точках xc должны выполняться условия сопряжения (30)
(при этом Qc = 0, Мс = 0, если хс не совпадает с точками приложения сосредоточенных сил и моментов.)
По тем же причинам условия сопряжения (30) должны выполняться и в сечениях х = хс(Ь), разграничивающих области упругопластического деформирования и разгрузки, области разгрузки и вторичной пластичности, области вторичной пластичности и вторичной разгрузки материала к-й фазы.
2. Метод численного интегрирования сформулированной задачи
Для приближенного интегрирования уравнения (24) по времени Ь будем использовать обобщенные методы Рунге — Кутты [4]. С этой целью уравнение (24) перепишем в виде системы
w(x,t),t = ь(х,Ь), ь(х,Ь),1 = д(х,Ь) — Ь(х; 1^). (31)
Численно систему (31) проинтегрируем с помощью двустадийного обобщенного метода Лобатто ША — метода трапеций (наименование метода использовано в соответствии с принятой в [8] терминологией), имеющего второй порядок точности по т, где т — шаг по времени Ь. Согласно основной идее обобщения методов Рунге — Кутты [4] этот метод применительно к системе (31) реализуется так:
wn+l = ^ип + т (VI + V) /2, Уп+1 = Vй + т (дп — Ь(х; Ж:) + дп+1 — Ь(х; ^)) /2,
Ж: = 1^п + 0 ■ т (V: + У2), Ж2 = 1^п + т (V: + У2) /2, V: = vn + 0 ■ т (дп — Ь(х; Ж:) +
+ дп+1 — Ь(х; Ж2)) , V = vn + т (дп — Ь(х; Ж:) + дп+1 — Ь(х; Ж2)) /2, откуда после несложных преобразований получим
1ип+1 = Ж2(х), 1ип = Ж:(х), ьп+1 = V2(x), vn = Vl(x), vn+1 =
= 2(ги]п+1 — 1ип)/т — vn, vn+1 = ьп + т (дп — Ь(х; 1^п) + дп+1 — Ь(х; wn+1)) /2; (32)
т2Б(х; 1ип+1) + Ш(х)1ип+1 = т2 (рп + рп+1 — Б(х; wn)) + 4К(х) (туп + 1^п) , (33)
где
л / ч Л2 { . { Л2т^ d2wn / Л2тп
Б (х; wn) = К (х) Ь (х; wn) = —- [А[х, 2 ) —т—2-В [х,
(34)
Лх2 \ \ ’ Лх2 ) Лх2 \ ’ Лх2
дп(х) = д(х,Ьп), рп(х) = р(х,Ьп), 1^п(х) = 1^(х,Ьп), Vй(х) = v(x,tn)
Ьп+1 = Ьп + т, п = 0,1, 2■■■, Ьо = 0;
шаг по времени т > 0 может быть переменным (т = тп); Ж:, Ж2, V, V2 — вспомогательные функции.
Если решение на п-м слое известно (известны wn, vn), то уравнение (33) задает решение на следующем слое (п + 1). Определенным недостатком уравнения (33) является то, что для вычисления его правой части необходимо применить дифференциальный оператор четвертого порядка Б(х; ■) к известной функции 1ип (см. (33), (34)). Чтобы избежать этого дифференцирования, введем в рассмотрение функции
Рп(х) = т2Б(х; 1^п) + 4К(х)1Пп, п = 0,1, 2... (35)
Тогда разрешающее уравнение (33) примет вид
т 2Б(х; 1^п+1) + 4К(х)1ип+1 = Рп+:(х), где правая часть известна и определяется по рекуррентной формуле
Рп+1(х) = —Рп(х) + 4К(х)(пп + 2мГ‘) + т2(рп + рп+1),
(36)
(37)
вытекающей из сравнения (35) и правой части (33). В момент времени Ь0 = 0 в силу начальных условий (26) из (35) получим известную функцию
Р0(х) = т2Б (х; w0(x)) + 4К(x)w0(x) ^0(х) = w0(x), V0(х) = v0(x))■
В частности, при нулевых начальных условиях
Р0(х) = 0, V0 = v0(x) = 0, -иР = 1и0(х) = 0,
а из (37), (39) следует
Р1(х) = т2{р>0 + р1)
(38)
(39)
(40)
Таким образом, для определения прогиба на (п + 1)-м слое необходимо проинтегрировать уравнение (36) с известной правой частью (37) — (40) при граничных условиях, которые полностью соответствуют граничным условиям исходной начально-краевой задачи. (В частности, в граничных условиях (28), (29) и условиях сопряжения (30) с учетом (20),
(27) необходимо заменить w на wn+1.)
Уравнение (36) формально совпадает с уравнением статического упругопластического изгиба стержня на упругом основании. Для интегрирования граничной задачи, соответствующей уравнению (36), можно использовать известные методы статики.
Для удобства численного интегрирования указанной граничной задачи (например, (36),
(28) или (36), (29) и т.п.) обыкновенное дифференциальное уравнение четвертого порядка (36) (см. (34)) перепишем в виде системы двух обыкновенных дифференциальных уравнений второго порядка
А х,
d2wm
—т
Лх2
Л2М т Лх2
d2wm ( d2wm\ , ч
— ВЬчхг) + М т М = 0'
+ 4К (х) wm = Рт(х), т = п + 1,
(41)
где Мт(х) = М(х,Ьт) — подлежащий определению изгибающий момент на т-м слое.
Первое уравнение системы (41) нелинейное, для его линеаризации можно использовать следующий итерационный процесс, качественно аналогичный методу переменных параметров упругости, широко используемый при решении упругопластических задач статики
[6]. Пусть w(m), Мт) — известные г-е приближения разыскиваемых функций wm, Мт, а следующие приближения этих функций определим из системы линейных уравнений
А х,
d2w
{г)
Лх2
Лх2
Л2Мт)
—т2^411 + 4К (х) wm
Лх2
Лх2
Рт(х), т = п +1, 5 = г + 1.
т
В качестве начального приближения для wm можно выбрать функцию
wm)(x) = wn(x) + ^п(х), т = п +1, (43)
получающуюся по формуле Тейлора в предположении, что на предыдущем п-м слое решение задачи уже известно.
Заменив в (42) вторые производные конечными разностями на равномерной сетке по х с шагом к, получим для трехточечного шаблона конечно-разностные аналоги уравнений (42) со вторым порядком аппроксимации по к. Для решения соответствующей СЛАУ можно использовать матричную прогонку [9], которая устойчива.
Как показывают многочисленные расчеты, проведенные авторами, итерационный процесс (42), (43) всегда сходится, поэтому предложенная схема численного интегрирования граничной задачи для системы (41) имеет точность порядка к2.
Замечание 1. Граничную задачу для уравнения (36) можно приближенно интегрировать методом Бубнова — Галеркина с привлечением итерационного процесса, аналогичного приведенному выше. При этом в качестве базисных функций можно использовать известные балочные функции [10]. При таком подходе удается избежать дифференцирования коэффициентов А и В, а значит, не нужно ставить дополнительные условия сопряжения решения (30) (при Qc = 0, Мс = 0) в сечениях х = хс, разграничивающих зоны чистоупругого и упругопластического деформирования к-й фазы композиции, зоны разгрузки и пластической загрузки и т.п.
В случае упругопластического деформирования фаз композиции разрешающее уравнение динамики стержня (24) (с учетом (25)) является квазилинейным уравнением параболического типа относительно прогиба w. Известно [11], что для квазилинейных дифференциальных уравнений общая теория устойчивости и сходимости конечно-разностных схем разработана недостаточно, поэтому основным критерием доверия той или иной конечноразностной схеме служат приближенные решения для тестовых (модельных) задач, аналитические решения которых известны.
Авторам пока не удалось доказать устойчивость численной схемы (32), (36) в общем
случае, когда оператор О (см. (34)) нелинеен (упругопластический изгиб). Но в пользу
устойчивости этой схемы говорят физическая корректность (непротиворечивость) результатов многочисленных расчетов, проведенных авторами, и хорошее совпадение численных решений с известными аналитическими решениями в случае упругих колебаний однородного стержня (см. рис. 1 и комментарии к нему в разд. 3).
Замечание 2. При постановке задачи динамического изгиба стержня в настоящем исследовании была использована приближенная связь между деформацией и прогибом
(4). Вместо (4) можно использовать более точное выражение деформации е через кривизну к нейтрального слоя стержня [7]
е(х,г,Ь) = —гк(х,Ь), к(х,Ь) =-)ХХ ( , ) , ■ (44)
(1+ «<,1 (х,1))г/2 '
При этом во всех формулах настоящей работы w,xx необходимо заменить на к из (44), а все численные схемы остаются в силе без принципиального их усложнения. Единственное отличие будет заключаться лишь в том, что вместо коэффициента А в (20) и во всех последующих формулах, где он встречается, нужно использовать
А(х,к,к *)
Ак (х, к, к ) -
(1 + /ш2х (х,г))3/2
Кроме того, если интенсивности армирования Uk(x,y,z) задавать в виде ступенчатых функций по переменным y, z, то в рамках используемого здесь подхода можно получить численные схемы интегрирования задачи неупругой динамики слоистых стержней с нерегулярной структурой армирования.
3. Обсуждение результатов расчетов неупругой динамики стержней
В качестве примера исследуем динамический упругопластический изгиб изотропных и армированных стержней прямоугольного поперечного сечения. Ширина поперечных сечений предполагается постоянной (b(x) = b = const) и одинаковой для всех конструкций, высота сечений может быть постоянной (H(x) = H* = const) или переменной и определяется формулой
H(x) = nH* + п sin(nx/l)(1 — n)H*/2 (0 < x < l, n > 0). (45)
(При однородном армировании стержней (uk = const) и задании высоты сечения в виде (45) расход материала k-й фазы композиции будет такой же, как и в стержне постоянного поперечного сечения, так как объем стержней одинаков и равен H*lb. При 0 < n < 1 из (45) следует, что высота концевых сечений стержня меньше, чем центрального сечения; при n = 1 получаем стержень постоянного сечения.) Длина стержней l = 1 м, высота стержня постоянного сечения H* = 0, 03 м. Концы стержней могут быть шарнирно оперты или защемлены. Внешняя распределенная нагрузка является нагрузкой взрывного типа, не зависит от x и затухает по экспоненциальному закону
p(x,t) = p(t) = p0 exp(-at) (a =12 с-1, t > 0), (46)
где p0 — значение нагрузки в начальный момент времени t = 0. (Может быть задан другой закон изменения взрывной нагрузки по времени [1, 2], в рамках настоящего исследования это непринципиально.) Изотропные стержни изготовлены из алюминиевого сплава АДН; армированные стержни изготовлены из того же сплава АДН и усилены однородно и прямолинейно (вдоль x) борными волокнами. (Механические характеристики фазовых материалов приведены в табл. 1.)
Будем подбирать максимальные значения р0 в (46) так, чтобы за характерный период времени T > 0, в течение которого исследуется процесс колебаний, напряженное состояние в фазах композиции стержня не превышало некоторого предельного состояния. Через ps обозначим maxр0/b, при котором хоть в одной из фаз композиции напряжение (5) достигает предела текучести ask, но пластичность не развивается (при р0 < bps — упругие колебания, при p0 > bps — упругопластические колебания); через pss обозначим maxp0/b, при котором хоть в одной из фаз композиции достигается вторичная пластичность, но вторичная пластичность не развивается (при bps < p0 < bpss — первичная пластичность и
Т аблицаї
Механические характеристики фазовых материалов [12]
Материал Ek, ГПа ask, МПа ayk, МПа 5 к, %
Сплав АДН 71.0 100.0 150.0 6.0
Волокна бора 416.5 — 3150.0 0.2
разгрузка, при p0 > bpss — вторичная пластичность); через psss обозначим maxp0/b, при котором хоть в одной из фаз композиции достигается третичная пластичность, но не развивается (при bpss < p0 < bpsss — вторичная пластичность и разгрузка, при p0 > bpsss — третичная пластичность); через py обозначим значение p0/b, при котором в момент b = const достижения центральным сечением стержня (x = l/2) максимального прогиба напряжение (5) хоть в одной из фаз композиции достигает предела временного сопротивления ayk (начальное разрушение; при t < tmax — первичная пластичность, при t > tmax — первичная разгрузка в опасном сечении стержня).
В качестве характерного периода T выберем время, в течение которого нагрузка (46) уменьшается в 1000 раз (T = Тюоо = — ln(0.001)/a = 0.576 с). При t > Тюоо внешняя нагрузка ничтожно мала (p(t) < 0.001p0) и колебания можно считать установившимися, по крайней мере при упругом изгибе.
В расчетах характерный период времени будем разбивать на 1000 слоев (т = Т/1000), а по длине стержня введем 401 узел.
Рассмотрим сначала динамический изгиб изотропного, шарнирно опертого алюминиевого стержня постоянного поперечного сечения. Значения ps, pss, psss, py для этого случая приведены в табл. 2, где n — значение коэффициента в формуле (45), используемое при расчете, а ер — отношение пластической деформации к упругой деформации в опасном сечении стержня. Как видно из табл. 2, значение pss на 42.5 %, а значение psss на 83.3 % больше ps. Следовательно, при повторном нагружении стержня нагрузкой с начальной амплитудой p0 < bpss в нем (за счет возникшей при первом нагружении приспособляемости материала) не будут возникать пластические деформации. При значениях p0 > bpss в стержне при повторных нагружениях будут возникать пластические деформации, что может привести к быстрому разрушению стержня за счет малоцикловой усталости. Таким образом, для рассматриваемого стержня путем предварительного его нагружения с интенсивностью p0 = bpss можно предельную по критерию пластичности нагрузку увеличить на 42.5 %. (Для сравнения отметим, что в [6] показано: путем автофретирования (автокрепления) вращающихся дисков их несущую способность также можно увеличить примерно
Т аблица2
Экстремальные значения параметров для шарнирно опертых профильных стержней
Упругий изгиб до появления первичной пластичности Упругопластический изгиб до появления вторичной пластичности Упругопластический изгиб до появления третичной пластичности Предельное динамическое состояние
N Ps, кПа Єр n Pss, кПа Єр n psss} кПа Єр n Py, кПа Єр
Алюминиевый, шарнирно опертый стержень постоянного поперечного сечения
1.0 61.6 0.0 1.0 87.8 0.67 1.0 112.9 1.87 1.0 241.0 27.4
Алюминиевый, жестко защемленный стержень постоянного поперечного сечения
1.0 98.2 0.0 1.0 134.7 0.64 1.0 174.9 3.85 1.0 382.9 27.4
Алюминиевый, шарнирно опертый стержень переменного поперечного сечения
0.54 91.6 0.0 0.6 254.0 7.16 0.6 254.0 7.16 0.54 386.7 27.4
Армированный, шарнирно опертый стержень постоянного поперечного сечения
1.0 234.3 0.0 1.0 307.9 0.31 1.0 439.4 0.87 1.0 1197.0 4.11
Армированный, шарнирно опертый стержень переменного поперечного сечения
0.56 I 351.5 I 0.0 0.53 473.0 0.31 0.5 1 718.5 1.00 0.57 1795.7 1 4.11
Рис. 1.
Рис. 2.
в 1.5 раза без последующего возникновения пластичности при повторных нагружениях.)
На рис. 1, 2 изображены кривые, характеризующие колебания центрального сечения (х = 1/2) стержня при разных уровнях нагружения. По оси абсцисс отложено время t (чтобы не загромождать рисунок, колебания изображены не на всем характерном отрезке времени Tl000 = 0.576 с, а примерно на его половине), а по оси ординат — безразмерный прогиб W = H w/(2l2ey\), где єу\ = ayl/El > 0 — предельная упругая деформация борных волокон, ayl, El — предел временного сопротивления и модуль упругости упругохрупких борных волокон (см. табл. 1).
Сплошная кривая на рис. 1 характеризует численное решение при p0 = bps (упругий
изгиб), а пунктирная соответствует аналитическому решению [10]
w (х, t)
(nmx\
siny—)
a
exp (-at) — cos (Qmt) + sin (Qmt)
' 'm.
(47)
где
-l
Qm = (nm/l)2 sjA/R, Cm = 4p0 [(nm/l) (Ra2 + A (nm/l)4)]
A = EbHl/12 = const, R = p0bH = const
(см. (3), (21) при h0 = H = H*, u0 = 1, uk = 0, k = 1, 2,...,N), суммирование осуществляется по нечетным m = 1, 3, 5... Сравнение этих кривых показывает, что численное решение почти не отличается от аналитического (в последнем случае в (47) удерживалось 50 слагаемых), наблюдается лишь незначительное запаздывание приближенного решения от численного. Такое хорошее совпадение численного и аналитического решений позволяет доверительно относиться к результатам расчетов, полученных на основе разработанного в настоящем исследовании метода.
Кривая 1 на рис. 2 соответствует значению ps (упругий изгиб) и совпадает со сплошной линией на рис. 1, пунктирная линия 2 соответствует значению pss (упругопластический изгиб без появления вторичной пластичности), кривая 3 — значению psss (упругопластический изгиб без появления третичной пластичности). Поведение кривой 3 показывает, что после прохождения точки A, лежащей на этой линии, амплитуда колебаний изменяется скачком. Это связано с тем, что точке A предшествует появление и развитие вторичной пластичности, которой соответствует диссипация кинетической энергии стержня (кинетическая энергия переходит в пластическую составляющую энергии деформаций), вследствие чего амплитуда колебаний изменяется скачком. Кривая 3 левее точки A в окрестности локальных экстремумов ведет себя достаточно гладко, правее же точки A в окрестности локальных экстремумов наблюдается нарушение гладкости (в частности, в окрестности этих точек кривизна кривой (w,tt) меняет знак). Это нарушение гладкости кривой 3 объясняется тем, что после возникновения и развития вторичной пластичности (после прохождения точки A) жесткость стержня резко уменьшается и с приближением прогиба к точкам экстремумов в конструкции повторно достигается и развивается первичная и вторичная пластичность (третичная пластичность не достигается). Но развитие этой пластичности столь незначительно, что практически не сказывается на амплитуде колебаний, однако влияет на кривизну кривой 3 в окрестности экстремумов правее точки A.
Как видно из рис. 2, к концу рассматриваемого периода времени (t = 0, 3 с) кривая 1 колеблется около нулевого уровня, а кривые 2, 3 — около ненулевых уровней, определяемых остаточными прогибами стержня.
На рис. 3 кривая 1 соответствует кривой 3 на рис. 2, но изображена при других масштабах по обеим осям координат. Линия 2 на рис. 3 получена при значении p0 = bpy (см. табл. 2); правой точке этой кривой соответствует начало разрушения центрального сечения стержня. Кривая 3 на рис. 3 рассчитана при нагрузке p0/b = 200 кПа. После прохождения максимума этой кривой в стержне возникают вторичные пластические деформации, столь значительные, что, не совершив даже полупериода колебаний, в стержне возникает начальное разрушение (максимальное напряжение достигает предела временного сопротивления). Линия 4 на рис. 3 получена при p0/b = 176 кПа. В таком стержне развивается вторичная пластичность, что приводит к возникновению начального разрушения
W(U2, /)
0.6
0.5
0.4
0.3
0.2
0.1
0
О 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
Рис. 3.
менее чем за 1.5 периода колебаний. (В двух последних случаях имеет место малоцикловая усталость и начальному разрушению соответствуют правые точки кривых 3, 4.)
Так как до настоящего времени неупругий динамический изгиб стержня рассматривался в основном в рамках жесткопластической модели [1, 2], то естественно провести сравнительные расчеты для диаграмм растяжения — сжатия, “близких” к диаграммам идеального жесткопластического материала. С этой целью будем увеличивать модуль упругости материала по сравнению с реальным значением, предел временного сопротивления приравняем пределу текучести (ay0 = as0 = 100 МПа, см. табл. 1), а остаточное удлинение при разрыве увеличим в 1000 раз (чтобы допустить к рассмотрению большие деформации, которым в рамках жесткопластической модели соответствуют пластические шарниры). Кривая 4' на рис. 3 получена при нагрузке p0/b = 176 кПа и 10-кратном увеличении модуля упругости сплава АДН (E0 = 710 ГПа), а линия 4'' — при 100-кратном увеличении (E0 = 7100 ГПа). Правые точки кривых 4', 4'' соответствуют достижению максимального прогиба центральным сечением стержня. Правее этих точек начинаются низкоамплитудные, но высокочастотные колебания стержня, причем почти сразу достигается либо начальное разрушение, либо вторичная пластичность. (В используемом на рис.
3 масштабе эти колебания нецелесообразно изображать, так как они выглядят как вертикальные штрихи.) Сопоставление кривых 4, 4', 4'' позволяет проследить за тенденцией изменения остаточных прогибов от увеличения модуля упругости. Чем больше модуль упругости, тем больше величина остаточных прогибов и тем меньше время достижения максимального прогиба. При этом завышение остаточных прогибов и занижение времени “остановки” движения стержня могут быть значительными (в 2-3 раза, как на рис. 3, и более). В предельном случае при использовании жесткопластической схемы деформирования (E0 ^ то) остаточный прогиб и время остановки движения будут еще больше отличаться от соответствующих реальных величин.
На рис. 4, 5 сплошными линиями изображены эпюры безразмерных прогибов в моменты времени, когда прогиб центрального сечения стержня достигает максимума. Кривые 1 — 3 на рис. 4 рассчитаны при тех же входных данных, что и линии 1 — 3 на рис. 2 соответственно; линии 1', 2, 3' получены для жестко защемленного стержня (подробнее
Щха)
0.045 -|
0.04
0.035 -
0.03
0.025
0.02 -
0.015
0.01 -
0.005
о 4 о
Рис. 4.
0.6 -|
0.5 -0.4 0.3 -0.2 -
0.1 -
О ^
о
Рис. 5.
см. ниже). Кривые 1 — 4, 4', 4'' на рис. 5 получены при тех же входных данных, что и 1 — 4, 4', 4'' на рис. 3. Пунктирные кривые на рис. 4, 5 соответствуют первой собственной форме упругих колебаний стержня [10]. (Для кривых 1, 2, 1', 2 на рис. 4 и линии 1 на рис. 5 соответствующие штриховые линии визуально не отличаются от рассчитанных сплошных кривых.) Сопоставление сплошных и штриховых линий на рис. 4, 5 показывает, что при р0 = Ьр3 и р0 = Ьр33 (линии 1, 2, 1', 2' на рис. 4), т. е. при упругом изгибе и упругопластическом изгибе без появления вторичной пластичности, стержень колеблется почти неотличимо от первой собственной формы упругих колебаний и для таких случаев можно использовать более простую расчетную схему [14]. При более интенсивном нагружении стержня (кривые 3, 3' на рис. 4 и 2 — 4 на рис. 5) эпюра прогибов значительно отличается от первой собственной формы. Однако, как видно из рис. 4, 5, эпюру прогибов можно
удовлетворительно аппроксимировать суммой 2-3 первых форм собственных колебаний стержня. А значит, если в соответствии с замечанием 1 граничную задачу для уравнения (36) интегрировать методом Бубнова — Галеркина, то достаточно ограничиться несколькими первыми членами разложения решения по формам собственных упругих колебаний стержня.
Эпюры прогибов 4', 4'' на рис. 5 соответствуют правым точкам кривых 4', 4'' на рис. 4. Линии 4', 4'' на рис. 5 по форме близки к ломаным, что соответствует возникновению в центральном сечении стержня напряженного состояния, близкого к пластическому шарниру. Следовательно, наблюдается качественное сходство между расчетом по жесткопластической схеме, когда в шарнирно опертом стержне в центральном сечении вводят пластический шарнир, и расчетом по упругопластической схеме для идеального высокомодульного материала. Однако, в отличие от жесткопластической схемы, на начальных участках кривых 4', 4'' рис. 3 соответствующие эпюры прогибов не имеют “изломов” в сечении x = //2. С увеличением времени форма этих эпюр эволюционирует, проходя последовательно стадии, качественно схожие с кривыми 1, 2, 3 на рис. 4, затем с линиями 1, 4, 3, 2 на рис. 5, и лишь при развитии больших пластических деформаций в центральном сечении стержня возникает “излом”.
Замечание 3. Строго говоря, значения py, приведенные в табл. 2 для алюминиевых стержней, и кривые 2 — 4, 4', 4'' на рис. 3 и 5 рассчитаны не вполне корректно с точки зрения механики, так как прогибы при этом превышают высоту стержня, т. е. являются большими, и соответствующие расчеты необходимо проводить по теории гибких стержней, учитывая выражение для деформации (44). Однако указанные кривые дают определенное качественное представление об особенностях упругопластических колебаний стержней при интенсивных нагрузках. Прогибы же, изображенные на рис. 2, 4, являются малыми и соответствующие расчеты корректны с точки зрения механики.
Проследим за изменением значений ps, pss, psss, py при изменении закрепления стержня. В табл. 2 приведены значения указанных величин для жестко защемленного стержня постоянного поперечного сечения, выполненного из алюминиевого сплава АДН. Кривые 1', 2, 3' на рис. 4 характеризуют эпюры прогибов стержня в моменты времени, когда прогиб центрального сечения достигает максимума, при значениях ps, pss, psss соответственно. Сопоставление данных, приведенных в табл. 2, показывает, что для защемленного стержня значения ps, pss, psss, py в 1.5 — 1.6 раза больше, чем для шарнирно опертого. Это объясняется большей жесткостью защемленного стержня и возникновением в нем трех зон пластичности: в окрестностях концов стержня и в окрестности центрального сечения (при шарнирном опирании стержня возникает лишь одна пластическая зона в окрестности центрального сечения). Предельные напряженные состояния в защемленном стержне достигаются на его концах, а при шарнирном опирании — в центральном сечении.
Выше исследовались стержни постоянного поперечного сечения. Варьируя же форму сечения, можно управлять несущей способностью таких конструкций. Рассмотрим совокупность шарнирно опертых стержней прямоугольного поперечного сечения постоянной ширины b = const, высота которых определяется формулой (45). Выбор параметра n в (45) задает форму профиля стержня. Выберем из этой совокупности стержней такие, в которых максимума достигают значения ps(n), pss(n), psss(n), py (n). Эти экстремальные значения приведены в табл. 2. Сопоставление значений ps, pss, psss, py для шарнирно опертых профилированных стержней и стержня постоянного поперечного сечения показывает (см. табл. 2), что за счет рационального профилирования стержня значение ps можно увеличить в 1.49 раза, pss — в 2.89 раза, psss — в 2.25 раза, py — в 1.6 раза. Следовательно,
целесообразно осуществлять рациональное профилирование упруго и неупруго деформируемых стержней, что позволяет в разы увеличить их несущую способность. (Значения р««, р«««, приведенные в табл. 2 для профилированных стержней, совпадают, так как при нагрузке, большей чем р.. = р888, в стержне не достигается третичная пластичность, а возникает начальное разрушение за счет первичной и вторичной пластичности.)
Помимо профилирования стержня его несущую способность можно увеличить за счет армирования. В табл. 2 приведены значения р«, р88, р888, ру для шарнирно опертого армированного стержня постоянного поперечного сечения, изготовленного из сплава АДН и усиленного волокнами бора с интенсивностью армирования = 0.6. При этом значения р«, р««, р««« характеризуются достижением первичной, вторичной и третичной пластичности в связующем материале соответственно, а значение ру характеризуется началом разрушения упругохрупких волокон. Сопоставление указанных величин для армированного и изотропного шарнирно опертых стержней (см. табл. 2) показывает, что за счет армирования величина р« увеличилась в 3.8 раза, р88 — в 3.5 раза, р888 — в 3.89 раза, ру — в 4.97 раза.
Дополнительно несущую способность армированного стержня можно увеличить за счет рационального профилирования. Задавая высоту стержня формулой (45), определим максимальные значения р8(п), р88(п), р888(п), ру (п). Эти экстремальные значения для шарнирно опертого армированного (^ = 0.6) стержня приведены в последней строке табл. 2. Сравнение величин р«, р88, р888, ру в двух последних строках табл. 2 показывает, что за счет рационального профилирования армированного стержня его несущую способность можно увеличить в 1.5 ... 1.64 раза.
В заключение настоящего исследования остановимся на вопросе выбора характерного периода времени Т, в течение которого целесообразно рассматривать упругопластические колебания при нагрузках взрывного типа (46) и по истечении которого колебания можно считать установившимися. Выше в качестве характерного периода было выбрано время, за которое внешняя нагрузка (46) уменьшается по сравнению с начальным значением ро в 1000 раз (Т = Т1000 = 0.576 с). За это время шарнирно опертый алюминиевый стержень постоянного поперечного сечения успевает совершить примерно 40 колебаний (см. рис. 2), защемленный стержень — порядка 100 колебаний, профилированные и армированные стержни — еще большее количество колебаний. К концу этого периода колебания визуально не отличаются от установившихся. Если в качестве Т выбрать время 100-кратного уменьшения нагрузки, то из (46) получим Т = Т100 = 0.384 с. Расчеты, выполненные при Т = Т100, приводят к завышению на 3-5 % значений р88, р888 по сравнению с указанными в табл. 2. (На значения р«, ру выбор времени Т не оказывает влияния, необходимо лишь, чтобы время Т было больше времени начала полной первичной разгрузки и времени достижения центральным сечением максимального прогиба соответственно.) Если в качестве Т выбрать время 20-кратного уменьшения нагрузки (Т = Т20 = 0.250 с), то рассчитанные значения р««, р««« будут на десятки процентов превосходить те, что указаны в табл. 2. Следовательно, при нагрузках взрывного типа (46) в качестве Т вполне достаточно выбрать время Т1000.
Заключение
Таким образом, предложенный в настоящей работе численный метод позволяет эффективно решать задачи неупругой динамики стержней без излишней идеализации конструкции
(как это сделано в [2]). Некоторые результаты, полученные с применением этого метода, являются, насколько известно авторам, оригинальными (в частности, прохождение вторичной пластичности при упругопластических колебаниях; см. кривую 3 на рис. 2). Максимальные значения относительной пластической деформации , приведенные в табл.
2, позволяют заключить, что при значениях нагрузки pss, пластические деформации
сопоставимы с упругими (исключение составляет лишь изотропный, рационально профилированный стержень). Естественно предположить, что пластические деформации будут также малы в случае неупругой динамики пластин при уровнях нагрузки, равных pss, (когда в пластине не возникает вторичная или третичная пластичность соответственно). А значит, при таких уровнях нагружения для расчета упругопластического изгиба пластин можно использовать деформационную теорию пластичности [13]. (В пользу допустимости применения деформационной теории говорит и тот факт, что при нагрузках взрывного типа с интенсивностями pss, стержни колеблются по формам, почти не отличимым
от первой собственной формы упругих колебаний (см. рис. 4). Естественно предположить, что пластины при таких нагрузках также будут колебаться по формам, близким к первым собственным формам. Но при таком изгибе пластин применение деформационной теории вполне оправдано [13].) В рамках деформационной теории предложенный в настоящей работе численный метод можно модифицировать для расчета неупругой динамики изотропных и армированных пластин, исследование которой находится в настоящее время в самом начале [5].
Список литературы
[1] Комаров К.Л., НЕмировский Ю.В. Динамика жесткопластических элементов конструкций. Новосибирск: Наука, 1984. 236 с.
[2] Ананенко Л.А., Комаров К.Л. Динамика неупругих балок. Новосибирск: Наука, 1999. 151 с.
[3] Галин М.П. Поперечные колебания балок и плит за пределом упругости под действием импульсных и ударных нагрузок // Изв. АН СССР. ОТН. 1958. № 3. С. 42-50.
[4] Немировский Ю.В., Янковский А.П. Численное интегрирование двумерных краевых задач с большими градиентами решения // Вычисл. технологии. 2000. Т. 5, № 4. С. 82-96.
[5] Кошур В.Д., Немировский Ю.В. Континуальные и дискретные модели динамического деформирования элементов конструкций. Новосибирск: Наука, 1990. 200 с.
[6] МАлинин Н.Н. Прикладная теория пластичности и ползучести. М.: Машиностроение, 1968. 400 с.
[7] Феодосьев В.И. Сопротивление материалов. М.: Наука, 1986. 516 с.
[8] Деккер К., Вервер Я. Устойчивость методов Рунге — Кутты для жестких нелинейных дифференциальных уравнений. М.: Мир, 1988. 334 с.
[9] КАлиткин Н.Н. Численные методы М.: Наука, 1978. 512 с.
[10] РАБотнов Ю.Н. Механика деформируемого твердого тела. М.: Наука, 1979. 744 с.
[11] Самарский А.А. Теория разностных схем. М.: Наука, 1989. 616 с.
[12] КомпозиционныЕ материалы: Справочник. Киев: Наук. думка, 1985. 592 с.
[13] Ильюшин А.А. Пластичность. Основы общей математической теории. М.: Изд-во АН СССР, 1963. 272 с.
[14] Немировский Ю.В. Динамический изгиб упругопластических балок // Науч. вестн. НГТУ. 2001. № 1(10). С. 75-86.
Поступила в редакцию 25 ноября 2003 г., в переработанном виде —15 апреля 2004 г.
