Определение критических сил и форм потери устойчивости панели с трещиной методом конечного элемента Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Т о м XII 19 8 1

№ 1

УДК 624.073.518.61:502.831

ОПРЕДЕЛЕНИЕ КРИТИЧЕСКИХ СИЛ И ФОРМ ПОТЕРИ УСТОЙЧИВОСТИ ПАНЕЛИ С ТРЕЩИНОЙ МЕТОДОМ КОНЕЧНОГО ЭЛЕМЕНТА

В. П. Башкиров, В. Д. Чубань

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

Конечно-элементная модель пластины представлялась с помощью лагранжевых изопараметрических конечных элементов с общим числом неизвестных = ?000. В.сего было найдено 32 формы, соответствующие наименьшим по модулю критическим силам. Из них одна соответствовала растяжению, остальные — сжатию пластины.

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

Качественное объяснение этого явления довольно просто. При растяжении пластины вдоль оси х происходит раскрытие трещины; вследствие этого в зоне трещины появляются сжимающие

\у

^— 1 Т - \ 1 —

-—1 "Ь —

П С —► ^

—

——

•* •

-— т:1Х-

а

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

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

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

Следуя [1], представим тензор деформаций пластины в виде

е<7 = еи + \ + ш,т) {е]п + ®уя) ётп, (1)

где

ши— у ^«|)>

£тп— метрический тензор, уги} — ковариантная производная вектора,

ди) ,,/• . и! ~ ,1 Цг ’

х: — криволинейные координаты пластины.

Символы Кристоффеля 2-го рода ТТц можно представить': в виде

Г/гг = г5Г/,„

1 I деи , д§4 \ т, , .

где Гн £ = ~г~ I -2Н. _[---------^ — символы Кристофеля 1-го рода.

1 \ дх1 дх дх /

Пусть нагрузка на тело равна критической. Тогда для него будут возможны два бесконечно близких положения равновесия.

Будем обозначать компоненты вектора перемещений до потери устойчивости и<°>. Компоненты вектора перемещений пластины, потерявшей устойчивость, можно представить в виде

и. = и(°) + ),

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

Введем гипотезы:

— исходное состояние равновесия описывается уравнениями линейной теории упругости;

8—.Ученые записки ЦАГИ“ № 1 ИЗ

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

— изменениями размеров и формы тела в исходном состоянии равновесия можно полностью пренебречь, т. е. тело считается напряженным, но не деформированным;

— силы, действующие на тело, не зависят от его деформаций;

— закон Гука выполняется как для исходного, так и для конечного положения равновесия.

При этих предположениях можно доказать |1], что решение задачи о форме потери устойчивости и критической силе является экстремальной точкой функционала:

стпи = стп!‘ = с‘1'тп — упругие константы материала.

В этих выражениях индексом принимающим значения 0 и 1, отмечены величины, относящиеся к полям перемещений Я«Ч и м(1> соответственно.

Для численного решения задачи воспользуемся методом конечных элементов (МКЭ) в перемещениях [2].

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

Рассмотрим пространство с декартовой системой координат х.

где

(3)

атп (7) __ с,ппЧ е< 0.

2

Пусть V—область, занимаемая конечным элементом. Введем в ней криволинейную систему координат ?{;:?*, £2, £3}:

Х=*х{\), X £ Q, (4)

где 2 — область изменения 1-координат.

Система координат £ выбирается так, чтобы ее область изменения на элементе Q была достаточно простой (рис. 2). С помо-

А дхг*

щью матрицы Якоби преобразования (4) Д = —— метрический тен-зор gij определится из соотношения

__ rk /к

Sij — Jl J] •

Представим закон аппроксимации поля перемещений в элементе в виде

щ = Bus uas, i = 1, 2, 3, а = 1, 2, ... m, s = 1, 2, . . ., q, (5)

где ut — ковариантный компонент вектора перемещений точки среды в ^-системе координат; Bi0LS = B^s (X) — интерполянт; U*s — s — обобщенное перемещение узла а элемента; т — число узлов элемента; q — число обобщенных перемещений в каждом узле элемента.

Связь между тензором etj и обобщенными перемещениями элемента с учетом выражений (1) и (5) может быть представлена в виде

eij = Dijas Uas, (6)

где

Dijas = (Vi Вjas -f- уj Bias),

dB:„

_ D H*

Vic/M --- J‘ drm

j-lm _ dim

dxl

Связь между тензором (о,у и обобщенными перемещениями представим в виде

e>,. = Qy„ir, (7)

ии

где Qijas — — (Vi ^/«ts Vj Bias)-

Поле напряжений o<°> представим в виде:

oS3) = (8)

здесь X —параметр нагружения; — поле начальных напряжений, соответствующих Х=1.

С учетом (3), (6), (7) и (8) функционал (2) может быть записан в виде

F = -L (U* KU + W НИ), (9)

где К = [/С«?г] — матрица жесткости элемента,

К,,pr = Г Cm>“>Vm Bno.S V/

V

Н=[На5?г\ — матрица геометрической жесткости элемента, Н^г= ] Q^/„G//И,', Ятп?гс1У.

к. V

Здесь 0ЦтП =эV>)imglnt и = [£/“]— вектор обобщенных перемещений.- ;;-

Первая вариация функционала (9) по и должна обращаться в нуль:

(К+УЛ)и = 0. (10)

Таким образом, применение МКЭ позволило привести задачу поиска экстремальных точек функционала (9) к определению собственных чисел и векторов уравнения (10).

Представим уравнение (10) в виде

КУ + НУ\ = 0,

где

10 .1

У—матрица, столбцы которой соответствуют формам потери устойчивости, — собственные числа.

Пусть спектр собственных чисел выражения (10) записан в порядке возрастания собственных чисел по абсолютной величине

. V 0<|Х1|<|Х2|<. . . <|Х„|< . . . <|^|,

где N — общее число обобщенных перемещений конечно-элементной модели.

Этим собственным значениям соответствуют собственные вектора

^1> ^25 * • * ) ■ ’ * ) и N.

Отметим, что матрицы К и Н являются симметричными, матрица К — положительно определенной.

Как правило, число необходимых первых собственных значений много меньше общего числа собственных значений

й<А'. (11)

О' Метод, позволяющий наиболее эффективно находить собственные значения с учетом условия (11) (неполная проблема собственных значений), а также малой заполненности матриц К и И, основан на идее одновременной итерацпн в /г-мерном подпространстве первых п собственных форм [3—5].

Пусть заданы п линейно-независимых векторов (/?, и\, ■ ■ ■ ,

и„, которые являются начальными приближениями первых п собственных форм ии и2, ип. Предположим, что найдена г-я

итерация для векторов и{, По, ■ ■ ■, И1п, которые образуют матрицу ЫХп:

и1 = \и[, Щ, . . ., и‘\.

Алгоритм нахождения (г 4- 1)-итерацип для этих векторов может быть представлен в следующем виде.

1. Построим матрицу Рн'1 размерностью ДгХп

Г41 = ни1.

Мб

2. Найдем матрицу Х1+ї размерностью ЛОО из следующего соотношения:

что соответствует решению системы линейных уравнений для п правых частей.

3. Построим матрицы Л, + 1 и к размерностью п'Хп, являющиеся энергетическими эквивалентами матриц Н и К соответственно в «-мерном подпространстве ‘

после чего на главной диагонали матрицы А‘+1 будут стоять 1.

5. Решим полную проблему собственных значений и векторов для уравнения: .

где х1+1 ■—матрица размерностью составленная из столбцов

собственных векторов уравнения (12); Л'+1 — матрица собственных чисел уравнения (12) размерностью «Хи-

6. Ортогонализируем векторы А1!4-1, Х12+\ . . ., X]п+1 в норме К и соответственно Н

Можно доказать [4], что при I оо имеют место соотношения

Свойства (13) позволяют в ходе итераций следить за точностью получаемых форм и критических сил потери устойчивости.

В качестве объекта расчета рассмотрим пластину с размерами ребер &=■12, а = 24 и толщиной Л == 0,02, изготовленную из материала с модулем упругости Е = 9,827- 10е и коэффициентом Пуассона [а = 0,3. Длина трещины 1=3 (рис. 1).

Конечно-элементная модель (рис. 3) как на этапе определения а<°>, так и на этапе определения собственных чисел и векторов представлена 124 лагранжевыми изопараметрическими конечными элементами (рис. 2) с общим- числом переменных ^3000. В зоне трещины расчетная сетка более густая, чем у краев пластины.

КХІЛ = Ґ+\

где

1

0

0

(12)

и1+1 Хі+1хі+:

(13)

Для этого элемента в качестве обобщенных перемещений в узле выбираются линейные перемещения узла и2, а3 и углы поворота узла «р1, Тогда

£/* 1 = и1, гГ2 = и2, иа3 = и3, £7“4 = <?\ иа5--=ср2.

Аппроксимация поля перемещений внутри элемента строится на основе подхода Тимошенко—Рейсснера: исходная нормаль к срединной поверхности пластины в процессе деформации остается прямой, но не обязательно нормальной к ней, т. е. допускается наличие в пластине поперечных сдвигов.

С учетом этого интерполянт В для лагранжевого элемента с поликвадратичной аппроксимацией поля перемещений имеет вид:

Вы3 - /а и ;2); Вы = — Ц (&\ V);

Вы5 = ;3): 5=1,2, 3,

где /а1, Л2, /Г3 — единичные векторы ортогонального базиса в узле а; — толщина элемента в узле а; £«(?*, ;2) — интерполяцион-

ные полиномы, которые вычисляются как произведение полиномов Лагранжа по направлениям н ;2 соответственно.

Геометрическое положение элемента в пространстве определяется соотношением:

здесь /?а — радиус-вектор узла а элемента; /?(£\ ?2, ?3) —радиус-вектор точки элемента в х-системе координат, имеющей координаты (I1, £2, I3) в криволинейной системе координат.

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

Пластина нагружена равномерно распределенной единичной силой, приложенной к коротким сторонам.

Расчет проводился в следующем порядке.

1. Определение напряженно-деформированного состояния пластины при заданных нагрузках и закреплениях. На рис. 4 представлено деформированное состояние конструкции.

2. Формирование на основе полученных напряжений матрицы геометрической жесткости И.

3. Определение первых 32 собственных чисел и векторов выражения (10). Следует отметить, что первые двадцать четыре собственных числа получились отрицательными, что соответствует потере устойчивости от сжимающих усилий, причем первое значение критической нагрузки (Я, =— 1,338) хорошо совпадает с Эйле-ровским решением задачи о потере устойчивости шарнирно опертого сжатого стержня

Двадцать пятая форма потери устойчивости, соответствующая растяжению, представлена на рис. 5. Для большей наглядности сетка недеформированной конструкции вычерчена более крупно. Критическая сила при данной форме потери устойчивости Р->ь — = 657,3.

На рис. 6 приведены эпюры изгибающего напряжения ауу вдоль линий АВ и ВС соответственно. Другие компоненты напряжений

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

Время центрального процессора ЭВМ для определения поля напряжений о<°) и поиска критических сил и форм потери устойчивости. составило соответственно 50 и 130 мин.

В заключение отметим, что все расчеты и чертежи на графопостроителе выполнялись с помощью комплекса программ „Сис-тема-4“, разработанного для ЭВМ БЭСМ-6 А. Б. Кудряшовым, Т. В. Снисаренко, Ю. А. Шевченко и авторами.

ЛИТЕРАТУРА

1. Новожилов В. В. Основы нелинейной теории упругости.

М.—Л., Гостехиздат, 1948.

2. Зенкевич О. Метод конечных элементов в технике. М., .Мир*, 1975.

3. Jennings A., Orr D. R. L. Application of the simultaneous

iteration method to undamped vibration problems. „Int. J. for Num. Meth. in Eng.“, vol. 3, 1971.

4. С orr B. R., Jennings A. Implementation of simultaneous ite-

ration for vibration analysis. „Computers and Siructures”, vol. 3, 1973.

5. Bathe K.-J., Wilson E. L. Solution methods for eigenvalue problems in structural mechanics. „Int. J. for Num. Meth. in Eng.“, vol. 6, 1973.

Рукопись поступила 26\VI 1979 г.