Расчет подкрепленных тонкостенных конструкций методом конечного элемента

_______УЧЕНЫЕ ЗАПИСКИ Ц А Г И__________

Том III Т9 72 № 1

УДК 629.7.015.3

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

Ю. И. Иванов

Изложена теория одного варианта метода конечного элемента для расчета подкрепленных тонкостенных конструкций широкого класса (плоские и искривленные панели, крылья и оперения различных схем, отсеки фюзеляжа). В качестве иллюстрации возможностей метода и программы для ЭЦВМ приводятся примеры расчета ряда конструкций. Дается сравнение с экспериментом.

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

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

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

На каждую поверхность обшивки наносится сетка, образованная двумя семействами линий аир, которые, как правило, совпадают с направлением стенок и подкрепляющих элементов. Точки пересечения линий а н р называются узловыми. Отрезок прямой, соединяющий соответствующие узловые точки на двух поверхностях, называется узлом. Поверхности задаются координатами узловых точек в правой системе координат хуг. В каждом узле вводится правая система трех взаимно ортогональных векторов а, Ь, п. Направления векторов (а, Ь, «),- называются главными направлениями в узле у. Вектор п1 всегда направлен вдоль узла у. Начало системы векторов (а, Ь, п)) находится в средней точке узла у'-

Поверхность, расположенная от начала в положительном направлении П), называется положительной, противоположная ей поверхность — отрицательной.

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

щие на узел у, раскладываются соответственно на компоненты А), В., /V;, Ма/, МьПоложительные направления перемещений и внешних усилий показаны на фиг. 1. К любому узлу по направлению любого компонента перемещения может быть приложена жесткая внешняя связь. В этом случае данный компонент перемещения равен нулю.

Элементами расчетной схемы конструкции являются пояса, работающие на растяжение — сжатие, панели обшивки и участки

6 5 /3

/—узел /; 2—стенки; /—панели обшивки;

гэсл у; 2—стенки; 3—отрицательная поверхность; панели обшивки; 5— положительная поверхность; б—пояса

Фиг. I. Отсек типовой конструкции.

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

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

— расчетная схема более естественным образом описывает действительную конструкцию, так как панели обшивки и стенки, как правило, имеют трапециевидную форму. Условия совместности деформаций для элементов выполняются более точно;

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

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

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

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

Р — аР, (О

то в силу ортогональности матрицы а обратное соотношение будет

Р = а'Ё.

Здесь и в дальнейшем штрих означает операцию транспонирования. Если некоторые координаты вектора Г тождественно равн^1 нулю, то их можно исключить, понизив размерность вектора Р. При этом вычеркиваются соответствующие строчки в матрице а и соответствующие столбцы в матрице а'. Матрицы а и а' становятся прямоугольными.

Матрица жесткости К» для элемента ц устанавливает соотношение между силами в узлах и перемещениями этих узлов

Данное соотношение в системе 7 (пусть это будет местная система координат, связанная с элементом) запишется в виде

В системе 7 (система главных направлений) соотношение (3) с учетом (1) и (2) запишется как

где

(4>

(5>

Матрица кесть матрица жесткости в системе 7. Если матрица к1Х симметричная, то из (5) следует, что симметричной будет и матрица к^.

Пояс. Для участка пояса (фиг. 2) соотношение (3) относительно оси £ запишется в виде [2]

I

1 - 1 - 1 I

Р*

Я

В системе главных направлений вводятся векторы

п -\Ъ\

Рп~ \ р, I '

Р-1Р< \ Р(

Л-

На основании фиг. 2 можно получить

И, !

в% V,

м ТО» 1

А/ Л V '?<*” 1

)

(6>

Плюс соответствует элементу положительной поверхности, а минус — элементу отрицательной поверхности. Матрица жесткости пояса в системе главных направлений определяется по формуле (5) с учетом (6) и (9).

При известных перемещениях усилие в средней точке пояса по его длине определится по формуле

ЕЕ

Р = ~{~ \ 1 )а„р„.

Трапециевидная панель обшивки. Для трапециевидной панели обшивки выбирается местная система координат Ь] (фиг. 3); ось т\ направлена вдоль параллельных сторон. Усилия и перемещения в панели в системе Ь] представляются в виде

где для V = к, /, т

Р, *=>

Л

Рф

х,

у.

I - \ І р>

— I Рп I Р т

р= \ -1/,

Панель нагружена по граням потоками касательных сил, для которых имеют место следующие соотношения [3|:

Потенциальная энергия трапециевидной панели вычисляется по формуле [3]

<725

\У

205

С,

где

5-у (/,-№*,

40

С— I + — — (<ё"’ /1 + tg /Л х> + tg‘ /,).

З Е

Потоки касательных сил (11) статически эквивалентным образом приводятся в узлах к усилиям X и У. В результате получим

где £ имеет выражение

и — {—1 ——1 І-, 1 I, 1 - цу,

“ —1~ + tg Хг> ^>2 — Хї-

Работа сил (13) на перемещениях (10) равна

А Яор0=у Ьд^р0.

Па равенства работы внешних сил (15) и потенциальной энергии (12) следует

<38Ь .

9=21С Рг (

После подстановки (16) в (13) получим

Р» — , Ро —: к„ Рц,

где

4SC

4SC

есть матрица жесткости трапециевидной панели в местной системе координат.

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

0 л

#V=

pJ Pi

р* i Pt\ p.i

Pi Pi

рт Pm

а п =

О

0

Здесь подвекторы в Р„ и р„ определены выражениями (8) для v = = /, к, I, т\ матрицы с» имеют вид

Н, Я,

cos

//- АЛ

с. =

_c-°s

cos

(’i- «О ± -у cos (t), bj

Относительно знаков справедливо замечание, сделанное выше. Матрица жесткости панели в главных направлениях определяется по формуле (5) с учетом (17) и (18).

Трапециевидная стенка. Стенка трапециевидной формы рассматривается как случай панели обшивки, у которой грань /'/и совпадает с узлом я, а грань — с узлом/ конструкции. При этом

Л = я„ /, = //,.

Из допущения о нерастяжимости конструкции по направлениям узлов следуют соотношения

'и, = г>т= г>3. ^кшшю1 = Ъй.

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

п. -

и.,. — — и,=

- - и,-

и,= - ик =

В результате векторы усилий и перемещений для стенки в местной системе координат можно записать в виде

у ' S V*

Ж

S ; Рс = 1 S

У, V,

м., Ь і

Переход от (10) к (22) и обратный переход осуществляются по формулам

яс = с'я0> р0 = ОсРс, где с учетом

а,. —

0 - HJ2 0 0

1 0 0 0

0 0 0 — Н,! 2

0 0 1 0

0 0 0 и,і 2

0 0 1 0

0 Н 2 S 0 0

1 0 0 0

Матрица жесткости стенки относительно (22) с учетом (17) и (23) принимает вид

кс —Оск0ас — Lc.L t,

4 SC

где

Lc =* Ос L.

Здесь L получается из (14) при условии (19).

С учетом (16), (23) и (25) средний поток в стенке определяется по формуле

В системе главных направлений векторы усилий и перемещений узлов s и t имеют вид (7), (8). Матрица ас перехода от (22) к (7) сохраняет выражение (9), где для

cos(i), аУ) cos(y), b.,) cos(t), «,) 0 0

ООО —cos

Матрица жесткости стенки в главных направлениях определяется по формуле (5) с учетом (9), (24), (27).

3. Система уравнений равновесия и ее решение. Система уравнений равновесия узлов конструкции в матричной форме может быть записана в виде

Кг = R,

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

Л, --

определенных по главным направлениям в соответствии с формулами предыдущего раздела. Задача автоматического формирования матрицы К для любой конструкции, относящейся к классу рассматриваемых, при произвольных кинематических граничных условиях является одной из основных при разработке алгоритма расчета. Такой алгоритм разработан и реализован в виде программы для ЭЦВМ типа БЭСМ.

Если столбцы г и /? расчленить на группы и в каждую группу отнести компоненты, относящиеся к узлам, расположенным на одной линии сетки, то в большом числе случаев система (28) приводится к квазитрехленточному виду [4]. Система такого типа достаточно эффективно решается методом последовательного группового исключения неизвестных [5]. Стандартная программа решения такой системы для ЭЦВМ типа БЭСМ приводится в работе 16].

После определения перемещений усилия в элементах конструкции определяются по формулам разд. 2.

4. Примеры расчета. В качестве иллюстрации возможностей метода приводятся некоторые результаты расчетов различных конструкций, выполненных с помощью универсальной программы, составленной для ЭЦВМ типа БЭСМ.

На фиг. 4 показано распределение интенсивности напряжений

3, = + «5 + Зх»,

в стенке силовой нервюры с вырезом. Рядом с кривыми на фиг. 4 указаны значения интенсивности напряжений, отнесенные к наибольшей величине интенсивности в стенке. Нервюра нагружена в левом сечении парой сосредоточенных сил, приложенных к поясам. Эти силы уравновешиваются по контуру нервюры потоками касательных сил, приходящих со стороны обшивки и стенок лонжеронов. Выбрана расчетная сетка со следующими параметрами: относительный шаг но высоте нервюры 0,05; относительный шаг но длине нервюры переменный и изменяется от 0,025 до 0,05; число узлов на половине нервюры 278. Число неизвестных составляет 538. Для расчета одного варианта нагружения требуется 26 мин машинного времени.

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

сечений даны распределения напряжений при кососимметричном нагружении, отнесенных к напряжениям при симметричном нагружении, •/. = огг- кс/згг. с- Там же показана принятая расчетная сетка, содержащая на одной консоли 89 узлов. Число неизвестных составляет 425. При первом решении задачи на один случай нагружения требуется 22 мин машинного времени. Каждый последующий расчет на другую нагрузку выполняется за 8 мин.

На фиг. 6 показаны некоторые результаты расчета консольной круговой оболочки, нагруженной сосредоточенной радиальной силой. Данные эксперимента и расчета этой оболочки методом Хоффа приводятся в работе [7]. Сплошные линии соответствуют методу конечного элемента, пунктирные — методу Хоффа, штрих-пунктирные — элементарному решению, основанному на допущении, что шпангоуты абсолютно жесткие. Эпюры изгибающих моментов в нагруженном шпангоуте показаны кривыми /; экспериментальные значения обозначены кружочками; решения по методу конечных элементов и методу Хоффа для изгибающих моментов в этом шпангоуте практически совпали. Кривыми 2 показаны эпюры потоков касательных сил в обшивке в сечении 2, а кривыми 3—эпюры потоков касательных сил в сечении 3. Экспериментальные значения напряжений для сечения 2 обозначены крестиками,а для сечения 3— треугольниками. Решение Хоффа получено в предположении,, что обшивка не работает на нормальные напряжения в окружном направлении. В решении методом конечного элемента в жесткость шпангоутов включалась часть площади прилегающих полос обшивки.

Фиг. б

Этим объясняется расхождение решений, полученных двумя методами. При решении методом конечного элемента выбрана расчетная сетка со следующими параметрами: шаг в направлении образующих равен расстоянию между шпангоутами, шаг в окружном направлении соответствует Д'-р = 15°, число узлов на половине оболочки — 52. Число неизвестных составляет 157. Время решения задачи 9 мин.

ЛИТЕРАТУРА

1. Аргирис Дж. Современные достижении п методах расчета конструкций с применением матриц, пер. с англ. под ред. А. Ф. Смирнова. М., Стройиздат, 1968.

2. Современные методы расчета сложных статически неопределимых систем. Сб. переводов под ред. Филина А. П. Л., Судпром-гиз, 1961.

3. И в а и о в Ю. И. Расчет крыльев малого удлинения методом перемещений. Сб. .Прочность летательных аппаратов". Труды ЦАГИ, вып. 1069, 1967.

4. Д л у г а ч М. И., Ш и и к а р ь А. И. Решение на электронных вычислительных машинах симметричных систем линейных алгебраических уравнений строительной механики и теории упругости.

ДАН УССР, № 4, 196!.

6. И в а и о в Ю. И. Программа решения систем линейных алгебраических уравнений с квазитрехленточной симметричной матрицей. Труды ЦАГИ, вып. 1331. 1971.

7. .Вопросы прочности цилиндрических оболочек". Сб. переводов под ред. Даревского В. М. М., Оборонгиз, 1960.

Рукопись поступила 1П\/ 1971 г.