Двухсеточное моделирование цилиндрических оболочек и панелей переменной толщины Текст научной статьи по специальности «Физика»

УДК 539.3 А.Д. Матвеев, А.Н. Гришанов

ДВУХСЕТОЧНОЕ МОДЕЛИРОВАНИЕ ЦИЛИНДРИЧЕСКИХ ОБОЛОЧЕК И ПАНЕЛЕЙ ПЕРЕМЕННОЙ ТОЛЩИНЫ*

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

Ключевые слова: композиты, упругость, цилиндрические оболочки, панели, метод конечных элементов, двухсеточные элементы.

A.D. Matveev, A.N. Grishanov

DOUBLE-GRID MODELING OF CYLINDRICAL SHELLS AND PANELS OF VARIABLE THICKNESS

Double-grid curvilinear elements that describe three-dimensional stress state of homogeneous and composite cylindrical shells and panels of variable thickness are offered, their attaching and loading complex nature is taken into account. Double-grid elements take into account the heterogeneous structure of panels and shells and give raise to discrete models of small dimension.

Key words: composites, elasticity, cylindrical shells, panels, finite element method, double-grid elements.

Введение. Как известно [1-3], при расчете однородных и композитных цилиндрических оболочек и панелей используют приближенные теории, в основе которых лежат гипотезы о характере распределения перемещений и напряжений. Это порождает неустранимую погрешность в решениях, обусловленную выбранной математической моделью. Существующие теории оболочек и панелей не учитывают сложный характер их закрепления и нагружения (например, закрепление и нагружение оболочки, панели частично по толщине), не всегда точно описывают поведение панелей и оболочек под действием локальных нагружений. При дискретизации цилиндрических панелей и оболочек с учетом их неоднородной структуры, сложного вида закрепления и нагружения возникает необходимость использовать мелкие разбиения. В этом случае дискретные модели имеют высокую размерность, что создает трудности при численной реализации метода конечных элементов (МКЭ) [4, 5].

В данной работе для расчета однородных и композитных цилиндрических панелей и оболочек по МКЭ с применением мелких разбиений предложены криволинейные двухсеточные конечные элементы (КЭ) 1-, 2- и 3-го порядков. Двухсеточные КЭ (ДвКЭ) формы прямоугольного параллелепипеда для анализа деформирования упругих тел неоднородной структуры рассмотрены в работах [6, 7]. В данном случае для построения ДвКЭ, имеющего неоднородную структуру, применяем две вложенные сетки: мелкую и крупную. Мелкая сетка порождена базовым разбиением ДвКЭ, которое учитывает его неоднородную структуру. На мелкой сетке определяется крупная сетка. Вначале неизвестные во всех внутренних узлах мелкой сетки ДвКЭ исключаются с помощью метода конденсации [5]. С помощью аппроксимирующих функций перемещений, построенных на крупной сетке, неизвестные узлов мелкой сетки, которые лежат на границе ДвКЭ, выражаются через узловые неизвестные крупной сетки.

Матрицы жесткости и узловые усилия предлагаемых криволинейных ДвКЭ определяются в локальных декартовых системах координат, а системы уравнений МКЭ для дискретных моделей оболочек и панелей - в глобальных декартовых системах координат. Связь между локальными и глобальными декартовыми системами координат осуществляется с помощью матриц вращений, которые определяются только для векторов узловых перемещений криволинейных ДвКЭ. При построении ДвКЭ используем однородные односе-точные КЭ. Здесь рассматриваем криволинейные однородные односеточные КЭ и композитные двухсеточные КЭ. На рисунке 1 показан криволинейный КЭ Уе 1-го порядка размерами Ъех х Ъеу х ке2, где ае - угол

* Работа выполнена при финансовой поддержке РФФИ (код проекта 11-01-00053).

раствора КЭ Уе, 01 х1 у1- локальная декартовая система координат; г101 у1 - плоскость симметрии; ей - ось цилиндрической оболочки; Я{, Яе2 - радиусы нижней и верхней поверхностей элемента Уе; Не - толщина; Ьу - длина КЭ Уе; Иех =аеЯ{, узлы отмечены точками. Прямоугольники размерами Ьг2 х Ьеу есть боковые грани, криволинейные прямоугольники со сторонами Ье2, аеЯе, аеЯе2 - торцевые грани КЭ. Форма КЭ Уе есть прямая призма высотой Ьу, основанием которой является криволинейный прямоугольник. ДвКЭ Уа (рис. 2) имеет такую же геометрическую форму, как и КЭ Уе (рис. 1). Так как при мелком разбиении углы раствора криволинейных КЭ малы (угол ае мал, рис. 1), то их формы мало отличается от форм прямоугольных параллелепипедов.

В связи с этим для аппроксимации перемещений криволинейных КЭ используем известные полиномы 1-, 2- и 3-го порядков (которые применяются при построении характеристик КЭ формы прямоугольного параллелепипеда [4, 5]) и уравнения трехмерной задачи теории упругости, записанные в локальных декартовых системах координат данных КЭ (рис. 1, 2).

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

1. Односеточные криволинейные элементы. Рассмотрим односеточный однородный криволинейный КЭ Уе 1-го порядка со сторонами ЬХ х Ьу х Ьг2 (см. рис. 1). Перемещения, деформации и напряжения элемента Уе удовлетворяют соотношения Коши и закон Гука [4], которые отвечают трехмерной задаче теории упругости [8], т. е. во всей области КЭ Уе реализуется трехмерное напряженное состояние. Для аппроксимирующих функций перемещений ые, \е, we элемента Уе используем полином 1-го порядка, записанный в локальной декартовой системе координат 01 х1 у1(рис. 1)

ие, уе, we = а1 + ае х1 + а3 у1 + а4 21 + а5 х1 у1 + а6 21 х1 + а7 21 у1 + а8 х1 у 1. (1)

Обозначим: и1 - перемещение узла '' КЭ Уе в направлении оси 01 х1; х'1, у'1, г1 - координаты '' -го узла КЭ Уе, ' = 1,...,8, номера узлов отмечены на рисунке 1. Коэффициенты а1 в (1) находим из условия ие(х1,у1,21) = и1, ' = 1,...,8, которое представим [С]а = и, где а = {а1,...,а8}Т, и = [и\,...,и^}Т,

[С] - матрица размерности 8 х 8, Т - транспонирование. Откуда следует а = [С]-1 и. Используя найденный вектор а в (1), определяем функцию ие. Аналогично строим аппроксимирующие функции перемещений Уе, we КЭ Уе, которые запишем

8 8 8 ие = £ М1и], Уе = £ М^, Л^е = £ М^], (2)

'=1 '=1 '=!

где и1,у1,л1 - перемещения узла '; М = М (х1,у1, г1) - функция формы ' -го узла КЭ Уе; ' = 1,...,8 . Вектор функций перемещений ие = {ие, Уе, }т КЭ Уе представим

и е = [Ме Й , (3)

где бе - вектор неизвестных КЭ Уе, имеющий в локальной декартовой системе координат 01 х1 у1 структуру

бе = {и\,...,и\,у1,...,у1,л1,...л1]т , (4)

[Ме ] - матрица функций формы КЭ Уе .

Полную потенциальную энергию КЭ Уе запишем в матричной форме

Же = 2 У £Т в ейУ - \Уе иТ FedУ - ^ иТ Ц , (5)

Т Т

где £е = {ех! , £у, , , Уху , Уу^ , Ух^ } , ве = ^х^ , °гх , ахт , ^ , } - векторы функций деформаций и напряжений, найденных в локальной системе координат 01 х1 у1 г1 КЭ Уе; Fe, це -векторы объемных и поверхностных сил КЭ Уе; Уе, - область и поверхность КЭ Уе. Соотношения Коши и закон Гука для КЭ Уе имеют вид

£е = [Ве ]бе , ве = №е ]ве , (6)

где [Ве ], [Бе ] - матрицы деформаций и модулей упругости КЭ Уе. Подставляя (3), (6) в (5), получаем

Же (бе ) = 1У (бе )Т [Ве ]Т [Бе ][Ве ]б'ейУ & )Т [Ме ]Т FedУ ^ )Т [Ме ]Т Ц^ .

Из выполнения дЖе (бе)/ дбе = 0 следует

[Kl] 6i = Г] ,

где [Kl] = f [Be]T[De][Be]dV, Pl = f [Ne]TFedV - f [Ne]TqedS, [K¡] - матрица жесткости

JVe V e Se

и P] - вектор узловых сил КЭ Ve, которые определятся в локальной декартовой системе координат

° x y1 z1 КЭ Ve.

Процедуры построения однородных односеточных криволинейных КЭ 2-, 3-го порядков, имеющих такую же форму, как и КЭ Ve (рис. 1), аналогичны вышеописанной.

2. Композитные криволинейные двухсеточные элементы. Рассмотрим криволинейные трехмерные композитные ДвКЭ 1-, 2-, 3-го порядков. Изложим процедуру построения ДвКЭ 3-го порядка 2-го типа, который обозначим через Va (рис. 2), где Oxyz - локальная декартовая система координат ДвКЭ Va;

zOy - плоскость симметрии; ha - толщина; hya - длина ДвКЭ Va. Считаем, что между компонентами неоднородной структуры ДвКЭ связи идеальны. Функции перемещений, напряжений и деформаций компонентов ДвКЭ удовлетворяют закон Гука и соотношения Коши, которые отвечают трехмерной задаче теории упругости [8], т. е. во всей области ДвКЭ Va реализуется трехмерное напряженное состояние. Область

ДвКЭ Va представляем базовым разбиением Ra, которое состоит из однородных односеточных криволинейных элементов Ve 1-го порядка (рис. 1). Базовое разбиение Ra учитывает неоднородную структуру ДвКЭ Va и порождает мелкую трехмерную узловую сетку ha. На мелкой сетке ha определяем крупную сетку Ha . Узлы крупной сетки Ha на рисунке 2 отмечены точками, 32 узла. На сетке Ha строим аппроксимирующие функции ua, va, wa для перемещений u , v, w ДвКЭ Va, используя при этом полином 3-го порядка (записанный в локальной декартовой системе координат Oxyz, рис. 2)

ua, va, wa = a1 + a 2 x + a3 y + a 4 z + a5 xy + a6 zx + a7 zy + a8 xyz +

2 2 2 2 2 2 2 2 2 + a9 x + a10 y + a11 z + a12 x y + a13 x z + a14 y x + a15 y z + a16 z x + a17 z y +

3, 3, 3, 3, 3, 3, 3, 3, 3,

+ a18 x + a19 y + a20 z + a21 x y + a22 x z + a23 y x + a24y z + a25 z x + a26 z y +

2 2 2 3 3 3

+ a27 x yz + a28y xz + a29 z xy + a30 x yz + a31 y xz + a32 z xy .

По процедуре, аналогичной процедуре п. 1, на сетке Ha для ДвКЭ Va определяем аппроксимирующие функции перемещений ua, va, wa вида

32 32 32

ua =1 Niui, va =Z Nivi, wa Niw,, (7)

i=1 i=1 i=1

где ui, vi, wi - перемещения i -го узла крупной сетки Ha ; Ni = Ni (x, y, z) - функция формы i-го узла сетки Ha, i = 1,...,32; вектор 5a неизвестных крупной сетки Ha в локальной декартовой системе координат Oxyz имеет структуру

5 a = W^.. u32, v32, w32}T . (8)

Пусть ось O1 y1 КЭ Ve (рис. 1) (который расположен в базовой модели ДвКЭ Va) параллельна оси Oy (рис. 2) и между осями O1 x1 и Ox угол равен ве, вектор Pe, q e неизвестных КЭ Ve в декартовой системе координат Oxyz имеет структуру

6 е = к<, Vе!<, < }Т, (9)

где ие, уе,»е - перемещения узла / КЭ V, / = 1,...,8 .

Между векторами 6^, 6е имеем связь 6^ = [Те ]6е [3], где [Те ] - матрица вращений [4] (размерности 24 х 24) в силу (4), (9) имеет вид

[Те ] =

[ Е,] [Е0] [ Е2] [ Е0] [ Е ] [Е0] - [Е2] [Е0] [Е,]

(10)

[Е0],...,[Е2] - матрицы размерности 8х8; [Е0] - нулевая и [Е] - единичная матрицы, [Е,] = со*в [Е], [Ег] = мпв [Е].

Используя равенство 6^ = [Те ]6 е, получим соотношения [4]

[Ке ] = [Те ]Т [К!][Те ] , Ре = [Те ]Т Р\ , (11)

где [Ке ] - матрица жесткости и Ре - вектор узловых сил КЭ Уе, определяемые в декартовой системе координат Охуг (рис. 2).

На базовом разбиении Яа ДвКЭ Уа с помощью метода конденсации [5] строим суперэлемент О,. Используя (11), систему уравнений для разбиения Яа представим в виде

[ЛП]61 + [Лп]6 , = В,, [А21]61 + [А22]6 , = В 2, (12)

где 61, 6 , - векторы перемещений соответственно внутренних и граничных узлов сетки На базового разбиения Яа суперэлемента О,; [Ли], [Л22] ([Л12], [Л21]) - квадратные (прямоугольные) матрицы, В;, В 2 - векторы узловых сил.

Из (12) для суперэлемента О, получаем систему уравнений [К,, ] 6, = Р,, где [К,, ] - матрица жесткости, Р,, 6, - векторы узловых сил и неизвестных суперэлемента О,,

[К,] = [Л22] - [Л^Лц]-1^], Р, = В

2 - [Л21][Л11] 1В1. Полную потенциальную энергию П,

суперэлемента О, запишем в виде

П, = 2 6Т [К, ]6, - 6ТР,. (13)

Используя (7), между векторами 6,, 6а установим связь

6 , = [Ла ]6 а , (14)

где [Ла ] - прямоугольная матрица; 6а - вектор узловых неизвестных ДвКЭ V,.

Подставляя (14) в (13), из условия дП, / 56а = 0 получаем [Ка ] qа = Га, где

[Ка ]=[л: ]Т [К, ][Ла ], г = [Ла ]Т Р,, (15)

[Ка ] - матрица жесткости и Га - вектор узловых сил ДвКЭ Уа (2-го типа).

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

3. Результаты численных экспериментов. Рассмотрим в глобальной декартовой системе координат О0х0у0г0 (рис. 3) модельную задачу о деформировании композитной цилиндрической панели У0

(прямоугольной в плане) ступенчато переменной толщины, для которой: а0 = п/2 - угол раствора; О2х0у0 - плоскость симметрии панели.

Рис. 3. Панель У0

При у0 = 0 панель У0 закреплена. Внутренний радиус панели равен 24. Длина панели Ь = 45 . При у < Ь / 2 толщина панели равна 2И, внешний радиус равен 29. При Ь / 2 < у < Ь толщина равна И, где И = 2,5, внешний радиус равен 26,5. Базовая дискретная модель панели У0 состоит из однородных изотропных КЭ Уе 1-го порядка со сторонами Иех х Иеу х Иег (рис. 1), которая учитывает структуру и порождает криволинейную сетку ¡,к размерами 82х9\х\9. Для КЭ Уе имеем ае = а0/8\, Иех = аеЯ^, Иу = Ь0/8\ = 0,5, И1 = 2И/\8 = 0,278. Дискретная модель панели состоит из ДвКЭ Уа 3-го порядка с размерами 9Иех х 9Иеу х 9И%, Иах = 9Иех, Иау = 9Иеу, И"., = 9И1 (рис. 2). ДвКЭ Уа армирован непрерывными волоконами (с криволинейным поперечным сечением со сторонами Иех, ИI), направленными вдоль оси Оу, на рисунке 4 сечения волокон закрашены.

Рис. 4. Разбиение в сечении Охг ДвКЭ

Расстояния между волокнами по дуговой координате 2Иех , по толщине ДвКЭ Уа - 2Ъг2. Модуль Юнга связующего материала равен 1, волокон - 10, коэффициент Пуассона равен 0,3. На внешней поверхности панели У0 в узлах г,к действуют силы # = 0,1, где к = 10, г = 55 + 3(а -1), ] = 91,

а = 1,...,10.

Пусть ось О0у0 (рис. 3) параллельна оси Оу (рис. 2), а между осями О0х0 и Ох угол равен ва,

при этом вектор 5° узловых неизвестных ДвКЭ Уа в глобальной декартовой системе координат О0х0у0г0 имеет структуру

5а = [п1..,п°32,у°,...,у°32,, (16)

где и0, V0, - перемещения узла г ДвКЭ Уа.

С помощью матрицы вращений [Та ] размерности 96 х 96 между векторами 5а, 5°а установим связь 5а = [Та ]5а. Матрица [Та ] в силу (8), (16) имеет вид

[Ta ] =

[ Ea ] [ Ea ] [ Ea ]

[ E^a ] [ Ea ] [ E^a ]

- [ Ea ] [ Ea ] [ Ea ]

(17)

где [E0a ],..., [E% ] - матрицы размерности 32 x 32; [Ea ] - нулевая и [Ea ] - единичная матрицы, [ Ea ] = cos в a [ Ea ] , [ E a ] = sin ва [E a ] .

Используя равенство 5 a = [Ta ]5 °a, построим соотношения [4]

[K0 ] = [Ta ]T [Ka ][Ta ] , F = [Ta ] Fa , (18)

где [K0] - матрица жесткости; Fa - вектор узловых сил ДвКЭ Va, которые определяются в глобальной декартовой системе координат O0х0y0z0 панели V0.

Для ДвКЭ Va находим угол ва и по формуле (17) вычисляем матрицу [Ta], a = I,...,N, N - общее число ДвКЭ Va, N = 135 . Используя (21), определяем матрицы жесткости и векторы узловых сил ДвКЭ Va, a = 1,...,N, с помощью которых строим систему уравнений МКЭ для дискретной модели панели V0.

Анализ результатов расчетов показывает, что максимальные перемещения двухсеточной wh = 55,605 и базовой w0 = 58,890 моделей отличаются на 5,58 %. Максимальные эквивалентные напряжения, возникающие в окрестности изменения толщины панели (при 22 < L < 23), двухсеточной ah = 2,693 и базовой а0 = 2,649 моделей отличаются на 1,66 %. Расчеты показывают, что для однородной изотропной панели V0 (рис. 3) погрешность для максимальных перемещений (по оси O0z0) равна

1,40 %, для максимальных эквивалентных напряжений - 0,55 %.

Число неизвестных базовой модели равно 321030, двухсеточной - 4500. Ширина ленты системы уравнений МКЭ базовой модели равна 9474, двухсеточной - 1206. Двухсеточная модель занимает в 560 раз меньше памяти ЭВМ, чем базовая. Время реализации МКЭ для двухсеточной модели в 16 раз меньше, чем для базовой модели.

Литература

1. Гольденвейзер А.Л. Теория упругих тонких оболочек. - М.: Наука, 1976.

2. Огибалов П.М., Колтунов М.А. Оболочки и пластины. - М.: Изд-во МГУ, 1969.

3. Андреев А.Н., Немировский Ю.В. Многослойные анизотропные оболочки и пластины. Изгиб, устойчивость и колебания. - Новосибирск: Наука, 2001.

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

5. Норри Д., Ж. де Фриз. Введение в метод конечных элементов. - М.: Мир, 1981.

6. Матвеев А.Д. Некоторые подходы проектирования упругих многосеточных конечных элементов / Ин-т вычислительного моделирования СО РАН. - Красноярск, 2000. - Деп. в ВИНИТИ № 2990-В00.

7. Матвеев А.Д. Многосеточное моделирование композитов нерегулярной структуры с малым коэффициентом наполнения // ПМТФ. - 2004. - № 3.

8. Самуль В.И. Основы теории упругости и пластичности. - М.: Высш. шк., 1982.

УДК 519.163 С.А. Тарасов, Ю.С. Тарасов

МЕТОД АНАЛИЗА ЛОГИСТИЧЕСКИХ МОДЕЛЕЙ

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

Ключевые слова: алгоритмы, анализ моделей, логистика, графы.

S.A. Tarasov, Yu.S. Tarasov

METHOD OF THE LOGISTIC MODEL ANALYSIS

The analysis method for the logistic models described in terms of the graph theory is offered. The model analysis algorithms, developed in accordance with the proposed method, lead to the considerable reduction of the set task computing complexity.

Key words: algorithms, model analysis, logistics, graphs.

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

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

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