УДК 539.3
ОБ ОДНОЙ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ДЕФОРМИРОВАНИЯ УПРУГИХ МНОГОСЛОЙНЫХ АРОЧНЫХ КОНСТРУКЦИЙ ПРИ БОЛЬШИХ ПЕРЕМЕЩЕ НИЯХ И УГЛАХ ПОВОРОТА
В.Г. Дмитриев, О.В. Егорова, С.И. Жаворонок, Л.Н. Рабинский
Представлена техника численного решения задач о нелинейном деформировании многослойных гибких арок, описываемых обобщенной моделью Кирхгоффа. Вариационно-разностным методом получены уравнения дискретного аналога модели арки. Изложен алгоритм модифицированного варианта метода квазидинамического установления, получены оценки параметров итерационного процесса, обеспечивающие его устойчивость и сходимость численного решения.
Ключевые слова: панели многослойные, арки упругие, перемещения большие, метод вариационно-разностный, квазидинамическое установление, устойчивость, сходимость.
Рассматриваются упругие многослойные арочные конструкции, а также панели в условиях плоской деформации, под действием статической системы краевых и/или поверхностных нагрузок общего или локального характера. Для описания процессов деформирования упругих арок (панелей) при больших перемещениях точек координатной поверхности и неограниченных углах поворота нормали к ней в качестве неизвестных принимаются изменения декартовых координат х,у (рис.1) [1,2].
Для многослойного пакета в целом принимаются условия жесткого контакта слоев без взаимного отрыва и проскальзывания, при этом в общем случае рассматриваются слои- переменной толщины Ит; т - индекс
Рис. 1. Кинематика криволинейной панели
слоя: 1<m<M, M - число слоев (рис.2). В качестве координатной поверхности принимается срединная поверхность одного из слоев, поверхность контакта слоев или линия, проходящая через центры тяжести поперечных сечений.
До деформации элемент недеформированной образующей длиной ds0 имеет начальные координаты х0,у0, кривизну и начальный угол 60 между осью х и нормалью к образующей (рис. 1). Полагая координаты х0,у0 заданными функциями линейной координаты вдоль образующей х0=х0^0), у0=у0^0), для недеформированного состояния можно записать
где R10 - радиус кривизны в направлении образующей (рис. 1). В результате деформации элемент будет иметь длину ds, координаты х,у, угол поворота 6 и кривизну ^. Рассматривая координаты х,у как функции линейной координаты s вдоль деформированной образующей х=х^), у=у^), для деформированного состояния имеем соотношения, аналогичные (1)
где 6 - угол поворота в соответствии с гипотезой Кирхгоффа-Лява о "жесткой" нормали. Компоненты деформации вдоль образующей E11, а также изменение кривизны линии приведения определяются как
Б„ = -1; = ^ • (1 + Бп) -^0. (3)
Компоненты деформации для волокна оболочки на расстоянии z от координатной поверхности с учетом z•k10<<1 можно выразить в виде
30
Рис. 2. Структура слоев криволинейной панели
(1)
(2)
Е?1 = Бц + z • Ки. (4)
Силовые факторы в многослойном элементе - продольная сила Т11 и изгибающий момент М11 - выражаются через компоненты деформации координатной поверхности
Тп = БцЕц + ЛПКП; Мц = ЛцЕц + БПКП. (5)
Жесткостные коэффициенты Л11,Б11 для случая панельной конструкции определяются через упругие характеристики слоев и их толщины [2,3]
М 7ш Е(т) М 7ш Е(т) • z
Б11 = Т | 1 п (Ш)п (ш) Л11 = Т | Т-Пшошг
ш=12Ш-1 1 - П12 П 21 ш=1гШ-1 1 - п12 V 21
М 7ш Е(ш) _2
в ^ ш Е •z
11 = Т I 1 П(ш)п(ш) ^ (6)
ш=12ш-1 1- П12 П 21
где
Е(ш) у (ш) = Е(ш) V(ш),
4 >п2^ = Е2 (7)
и где Е((ш),Е2ш), п(ш), V 2ш)- модули Юнга и коэффициенты Пуассона. Полагая физико-механические характеристики материала ш-го слоя неизменными в пределах слоя и отсчитывая координату z от внутренней поверхности панели, формулы (6) можно привести к виду
М Е(ш) (7 - 7 ) 1 М Е(ш) (72 - 72 )
Б — 1 V ш ^ш-1А Л = 1 ^ ш ^ш-1А
Б11 ^ 1 п (ш)п (ш) ; Л11 1 п (ш) п (ш) ;
ш=1 1- п12 П 21 2 ш=1 1- п12 п 21
1 М Е(ш)(73 - 73 )
В11 =1Т Е1 (7ш ^ 7ш-1). (8)
11 3ш=1 1- п(ш)п 2ш)
Для случая однослойной конструкции, принимая в качестве координатной срединную поверхность, имеем
Т11 = Б11Е11; М11 = БПКП, (9)
где
Е Ь Е Ь3
Б11 = Г-^^; °11 = Т2ТТЕ1Ь-^. (10)
1 П12П 21 12(1 П12П 21)
Соотношения (6)-(10) относятся к случаю деформирования упругих панелей при их цилиндрическом изгибе. Физические соотношения для арочных конструкций могут быть получены из (6)-(10) путем соответствующих упрощений. Уравнения равновесия деформированного элемента арки (панели) записываются в виде [1,2]
^ + ВДц + д! = 0; ^ - к1Т11 + = 0; ds ds
01, = (11)
ds
где q1,q2 - компоненты поверхностной нагрузки. На краях рассматриваются различные варианты граничных условий [2].
При численном решении системы нелинейных дифференциальных уравнений (1)-(11) дискретизация по пространственным переменным осуществляется методом конечных разностей [4-8]. В области непрерывного изменения аргумента б0 (в - для деформированного состояния) вводятся две сетки: основная с целочисленными индексами ¿±1 и вспомогательная с индексами ¿±1/2. Сеточные функции обобщенных перемещений х(в),у(в) соотносятся с узлами основной сетки ¿±1, а сеточные функции компонент деформаций и напряжений - с узлами вспомогательной ¿±1/2. Шаг сетки для начального (недеформированного) состояния А0=Ав0 определяется как
1 о =—, (12)
0 N -1
где Ьо - начальная длина образующей конструкции (область изменения аргумента б0: 0< б0< Ьо), N - число точек дискретизации: 1< К N.
Для деформированного состояния шаг сетки Аб является переменной величиной, поэтому для упрощения формы записи дискретизирован-ных уравнений в исходных соотношениях (2)-(12) целесообразно перейти к дифференцированию по координате б0 вдоль недеформированной образующей оболочки, используя первую формулу из (3). Тогда для аппроксимации производных могут быть использованы разностные операторы второго порядка точности 0( 120) на сетке с постоянным шагом 10. Конечно-разностные аналоги уравнений равновесия (11), дискретизированные относительно узловой точки основной сетки в операторной форме могут быть представлены как
[Ь 1о(ик)] + = 0, (13)
где [Ьх (ик)]; - соответствующие конечно-разностные операторы для вектора ик сеточных функций обобщенных перемещений: и1=х1 (к=1), и2=у1 (к=2); - сеточные функции обобщенных компонент поверхностной и краевой нагрузки.
Для численного решения системы сеточных уравнений (13) используется квазидинамическая форма метода установления с переходом к эволюционной задаче путем замены уравнений (13) на уравнения, совпадающие по форме с уравнениями движения оболочки в вязкой среде [4,5]. Поскольку уравнения равновесия (11) и их конечно-разностные аналоги (13) записаны в проекциях для компонент тангенциального и и нормального w перемещения точек, то нестационарные уравнения метода установления запишутся как
[Ь 10 (ик)]1 + ^к )1 = (тки к)1 +(е ки к)1, (14) где ек - параметры удельной вязкости среды, р - плотность, тк=рИ (к=1,2), и1=и1 (к=1), U2=w1 (к=2) - сеточные функции обобщенных перемещений. Аппроксимация уравнений (14) на временной сетке с шагом А1=сопб1 с ис-
32
2
пользованием разностных операторов второго порядка точности О(Д t ) позволяет получить в явном виде выражения для скоростей [ú k ](n+1/2) на
(n+1/2)
временном слое t и сеточные функции обобщенных перемещений
Г -|(n+1) ,(n+1)
[úk ]i на временном слое t
[ú ](n+1/2) = [2mk -ekDt] [ú ](n-1/2) + 2Dt •[Li0(Uk) + qk |n), 1 [2mk +ekDt\ 'LúkJl [2mk +ekDt\ '
[úk ](n+1) =[úk f +Dt •[ú k Г/2). (15)
Таким образом, разностная аппроксимация нестационарных уравнений (14) приводит к итерационному процессу (15) нахождения решения исходной стационарной задачи (13). Сеточные функции перемещений xl,yl узловых точек в системе координат x,y определяются через сеточные функции компонент перемещений úl,wl по формулам перехода как
Xi = (Хо)1 + Dx¡; yi = (Уо)1 +Ду^ (16)
где
Дх1 = úl • sin0l - wl • cos0l; Ду1 = úl • cos0l + wl • sin0l. (17)
Параметры итерационного процесса - удельные вязкости среды ek(i) и шаг по времени Dt - определяются из условия ускорения сходимости и устойчивости разностной схемы [4]. Для случая стационарного итерационного процесса оценочные формулы с учетом структуры уравнений (14) запишутся в виде [2, 5]
e k = 2öe,(k)
m k m',(k>m 2'(t>; Dtk = 2at(k)
mi,(k) + m 2,(k) У
mk
mi,(k) 2,(k)
(18)
где mi,(k) и m2,(k) - наименьшие и наибольшие собственные числа для соответствующих разностных операторов в уравнениях (14); ae,(k) и at,(k) - близкие к единице поправочные коэффициенты. Шаг по времени Dt для всей разностной схемы определяется из условия вида:
Dt = min Dtk.
k k
Для нелинейных задач точное определение границ спектров разностных операторов связано со значительными математическими трудностями, поэтому m1;(k) и m2,(k) оцениваются в рамках линейных соотношений при соответствующих упрощениях в исходных уравнениях. Оценочные формулы для m1,(k) и m2,(k) для случая однослойной ортотропной конструкции, описываемой физическими соотношениями (9),(10), можно представить следующим образом:
mi,(i) = 4-f •sm2 ^Г0; mi,(2) = 1^-rfsm4 --0- + k120Bn; (19)
70 2 —0 7 0 2 —0
наименьшие собственные числа;
m 2,(2) = 16 cos4 p10 + ki0Bn, (20)
- наибольшие собственные числа. Здесь к10 - характерное значения кривизны.
При расчете реальных конструкций отношение т 2(к) / т1(к) обычно
велико, поэтому метод установления, построенный в форме оптимального линейного итерационного процесса (15), позволяет сократить число итераций примерно в д/т2 / т1 раз по сравнению с вариантом метода установления, когда в правых частях уравнений (14) учитывается только второе слагаемое [4,5].
Развитие метода применительно к нелинейно деформируемым тонким композиционным оболочкам представляется необходимым для решения задач об их нестационарном взаимодействии с потоками жидких или газообразных сред [9-13], где ограничение линейной моделью нередко приводит к недопустимым погрешностям.
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты РФФИ №№ 16-01-00751, 17-0801461).
1. Григолюк Э.И., Шалашилин В.И. Проблемы нелинейного деформирования: Метод продолжения решения по параметру в нелинейных задачах механики твердого деформируемого тела. М.: Наука, 1988. 232 с.
2. Нелинейное деформирование многослойных композитных оболочек вращения при больших перемещениях и углах поворота нормали / Дмитриев В.Г., Бирюков В.И., Егорова О.В., Жаворонок С.И., Рабинский Л.Н. // Известия вузов. Авиационная техника. 2017. №2. С. 8-15.
3. Васильев В.В. Механика конструкций из композиционных материалов. М.: Машиностроение, 1988. 272 с.
4. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Физматлит. Лаборатория Базовых Знаний, 2001. 632 с.
5. Dmitriev V.G., Egorova O.V., Rabinsky Lev N. Solution of nonlinear initial boundary-value problems of the mechanics of multiply connected composite material shells on the basis of conservative difference schemes // Composites: Mechanics, Computations, Applications: An Int. J. 2015. V. 6, Issue 4. P.
6. Трушин С. И., Князев А. А., Жаворонок С. И. Решение нелинейной задачи устойчивости многослойной оболочки вращения из композиционного материала с низкой сдвиговой жесткостью // Механика композиционных материалов и конструкций. 2001. Т. 7. № 3. С. 363-373.
Список литературы
265-277.
7. Trushin S. I., Zhavoronok S. I. Nonlinear analysis of multilayered composite shells using finite difference energy method // Space Structures 5: Proc. of the Fifth Int. Conf. on Space Struct., Univ. of Surrey, Guildford, UK. 2002. Thomas Telford Ltd. 2002. P. 1527-1533.
S. Дмитриев В. Г., Жаворонок С. И., Коровин Е. К., Москвитин Г. В. Оптимальные вычислительные технологии в математическом моделировании нелинейных задач механики деформируемого твердого тела // Инженерная физика. 200S. № б. С. 2-5.
9. Плоская задача дифракции акустической волны давления на тонкой ортотропной панели, помещенной в жесткий экран / Горшков А. Г., Жаворонок С. И., Медведский А. Л., Рабинский Л. Н. // Известия РАН. МТТ. 2004. № 1. С. 209-220.
10. Жаворонок С. И., Рабинский Л. Н. Осесимметричная задача нестационарного взаимодействия акустической волны давления с упругой оболочкой вращения // Механика композиционных материалов и конструкций. 200б. Т. 12, № 4. С. 541-554.
11. Нестационарная задача дифракции цилиндрической акустической волны давления на тонкой оболочке в форме эллиптического цилиндра / Горшков А. Г., Жаворонок С. И., Медведский А. Л., Рабинский Л. Н. // Int. J. Comput. Civil Struct. Engng. 2007. V. 3. No. 2. P. S2-93.
12. Егорова О. В., Жаворонок С. И., Рабинский Л. Н. Взаимодействие оболочки средней толщины с акустической волной // Вестник Московского авиационного института. 2010. Т. 17. № 2. С. 127-135.
13. Жаворонок С.И., Куприков М.Ю., Медведский А. Л., Рабинский Л.Н. Численно-аналитические методы решения задач дифракции акустических волн на абсолютно твердых телах и оболочках. М.: ФИЗМАТЛИТ, 2010. 223 с.
Дмитриев Владимир Георгиевич, д-р техн. наук, проф., vgd21D5@mail.ru, Россия, Москва, Московский авиационный институт (Национальный исследовательский университет),
Егорова Ольга Владимировна, канд. физ.-мат. наук, доц., janus olgaamail.ru, Россия, Москва, Московский авиационный институт (Национальный исследовательский университет),
Жаворонок Сергей Игоревич, канд. физ.-мат. наук, доц., Zhavoro-nok a iam. ras. ru, Россия, Москва, Институт прикладной механики Российской академии наук,
Рабинский Лев Наумович, д-р физ.-мат. наук, проф., декан факультета, ra-binskiy@mail.ru, Россия, Москва, Московский авиационный институт (Национальный исследовательский университет)
A MA THEMA TICAL MODEL OF THE DEFORMING OF ELASTIC LAMINA TED ARCH STRUCTURES CONSIDERING FINITE TRANSLATIONS AND ROTATION ANGLES
V.G. Dmitriev, O.V. Egorova, S.I. Zhavoronok, L.N. Rabinskiy
A new numerical simulation technique for nonlinear deforming flexible laminated arches is presented. The equations of the discretized model corresponding to the generalized Kirchhoff theory are constructed on the background of the variation-difference approach. The modified pseudoviscosity method for solution of the nonlinear problem is briefly described, and the estimation of optimum parameters of the numerical algorithm securing stability and convergence are obtained.
Key words: laminated panels, elastic arches, finite displacements, variation difference approach, pseudoviscosity method, numerical stability, convergence.
Dmitriev Vladimir Georgievich, doctor of technical sciences, professor, vgd2105@mail.ru, Russia, Moscow, Moscow Aviation Institute (National Research University),
Egorova Olga Vladimirovna, candidate of physical and mathematical sciences, do-cent, janus olgaamail.ru, Russia, Moscow, Moscow Aviation Institute (National Research University),
Zhavoronok Sergey Igorevich, candidate of physical and mathematical sciences, do-cent, Zhavoronok@iam. ras. ru, Russia, Moscow, Institute of Applied Mechanics of Russian Academy of Sciences,
Rabinskiy Lev Naumovich, doctor of physical and mathematical sciences, professor, dean of faculty, rabinskiy@mail. ru, Russia, Moscow, Moscow Aviation Institute (National Research University)