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

4./2011 ВЕСТНИК _7/202J_МГСУ

О ПОСТРОЕНИИ МАТРИЦ ЖЕСТКОСТИ ТРЕХМЕРНОГО ДИСКРЕТНО-КОНТИНУАЛЬНОГО КОНЕЧНОГО ЭЛЕМЕНТА С ЧЕТЫРЕХУГОЛЬНЫМ ПОПЕРЕЧНЫМ СЕЧЕНИЕМ МЕТОДОМ БАЗИСНЫХ ВАРИАЦИЙ

CONSTRUCTION OF STIFFNESS MATRICES OF THREE-DIMENSIONAL DISCRETE-CONTINUAL FINITE ELEMENT WITH QUADRANGULAR CROSS-SECTION BY METHOD OF BASIC VARIATIONS

П.А. Акимов, М.Л. Мозгалева, B.H. Сидоров Pavel A. Akimov, Marina L. Mozgaleva, Vladimir N. Sidorov

ГОУ ВПО МГСУ

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

Correct algorithm of construction of stiffness matrices of three-dimensional discrete-continual finite element with quadrangular cross-section by method of basic variations are under consideration in the distinctive paper.

Введение.

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

1. Дискретно-континуальная аппроксимирующая модель конструкции.

Пусть x3 - переменная, соответствующая направлению регулярности физико-геометрических параметров рассматриваемой конструкции (основное направление), причем пусть вдоль x3 физико-геометрические параметры конструкции изменяются кусочно-постоянно (заметим, что по переменным x1, x2 соответствующие параметры могут изменяться произвольно); Qk, k = 1,2,..., nk -1 - подобласти постоянства физико-геометрических параметров конструкции с поперечным (по отношению к основному направлению) сечением S„'; <як, k = 1,2,..., nk -1 - соответствующие расширенные области, окаймляющие исходные Qk, k = 1,2,..., nk -1, причем Qk <^a>t.

Принимается дискретно-континуальная модель объекта следующего типа: по основному направлению (вдоль оси Ox3) конструкции задача остается континуальной, по поперечным (по отношению к основному направлению) направлениям (вдоль осей Oxx и

0х2) производится сеточная аппроксимация с использованием стандартной техники метода конечных элементов [3, 6].

Расширенные области сок, к = 1,2,..., пк -1 выбираются стандартными в виде параллелепипедов. Аппроксимация поперечного сечения стандартной области состоит в задании сетки, топологически эквивалентной прямоугольной таким образом, чтобы она как можно лучше соответствовала очертаниям поперечного сечения конструкции (рис. 1.1, 1.2). Понятие топологической эквивалентности в данном случае означает, что она может быть получена из прямоугольной сетки в результате некоторой невырожденной деформации ячеек последней без их «перекручивания». Выбор такого класса сеток, с одной стороны, дает возможность аппроксимировать большое количество разнообразных конструкций, а с другой - позволяет использовать простую регулярную нумерацию узлов (двухиндексную), что приводит в дальнейшем к удобным математическим формулам, эффективным вычислительным схемам и алгоритмам, а также существенно упрощает сбор исходной информации и вывод результатов. Переход к прямоугольной сетке с единичным шагом легко осуществляется локальной заменой координат внутри сеточной ячейки. При «выпрямлении» сетки общий вид соответствующих уравнений сохраняется - меняются лишь элементы матрицы коэффициентов исходной системы. Заметим, что во многих случаях эффективно непосредственное применение прямоугольных сеток, поскольку это значительно упрощает и ускоряет алгоритмы, а, следовательно, позволяет использовать большое количество узлов, что приводит в итоге к увеличению точности.

1 2

N.

Рис 1.1. Пример выбора сетки, аппроксимирующей «поперечное» по отношению к основному направлению сечение конструкции

дискретно-континуальный конечный элемент 0)А ;

Рис. 1.2. Общий вид ДККЭ с четырехугольным поперечным сечением 2. Некоторые предварительные обозначения.

Введем обозначения: (г, у) - двухиндексный номер рассматриваемого дискретно-континуального конечного элемента (ДККЭ) с четырехугольным поперечным сечением; со.. - область, занимаемая ДККЭ; З1'-1 - область, занимаемая поперечным (по от-

ношению к основному направлению) сечением рассматриваемого ДККЭ (постоянна вдоль х3); 13 - длина ДККЭ по основному направлению (х3 е [0, /3]); х\к, к = 1,..., пк

- координаты сечений, в которых может происходить «скачкообразное» (разрывы первого рода) изменение физико-геометрических параметров конструкции). Рассматриваемый ДККЭ а>.. (рис. 1.2) определяется в виде:

„к-1

= и®к« , где а>к .. = ^Х2,Х3): С^х2) е ^ ; Х3 е [ х1к, хз,к+1 ]}. (2.1)

3. Локальная система координат дискретно-континуального конечного элемента. Восполнение неизвестных на элементе.

Рассмотрим произвольный (¡, у)-й ДККЭ. Пусть (х1(',7), х2',7)) - координаты (¡,у) -го узла сетки дискретно-континуальных конечных элементов, аппроксимирующих конструкцию. Тогда рассматриваемый ДККЭ имеет узлы (узловые линии) со следую-

щими номерами: (>, х2'-.о, (х;-+1^>, х2,+и)), (х(

„(¡,7+1) „(¡,7+1)

) и (

(¡+1,7+1) „(¡+1,7+1)

В поперечном (по отношению к основному направлению) сечении ДККЭ вводится локальная система координат 0^ и 012 (рис. 3.1), при этом е [0,1]; 12 е [0,1].

?Ь,7+1)

7(1+1,7+1)

Рис. 3.1. Переход к локальной системе координат на ДККЭ

Имеет место соответствие глобальных и локальных координат узлов элемента: (х('7>, х27>) ^ (0,0); (хГи', х™) ^ (1,0);

(х^, х2',7+1)) ^ (0,1); (хГ1,7+1), х2'+1,7+1)) ^ (1,1). (3.1) Пусть 1 = [ ^ t2Y и х = [ х1 х2]г - векторы координат произвольной точки дискретно-континуального конечного элемента в локальной и исходной глобальной системах координат. Формула преобразования координат на элементе (рис. 3.1):

х(11,12) = х(',л + 11А1 х + 12Д2х +1112А12х , (3.2)

где А,х(''7) = !('+и) -Р7; А2х(и) = I7 -'; Д.,х?7 = !('+и+1) -!('+и) - А2х . (3.3)

^ 1 п п ' 2 п п ' 12 п п 2 V /

Здесь х(', К'1"1) - векторы координат (', 7) -го узла элемента в глобальной и локальной системах координат соответственно.

В результате замены переменных вычисление производных по х1 и х2 производится согласно формулам дифференцирования сложной функции:

дф дф д 1к

дх ы д 1 дх

5 = 1,2,

(3.4)

где ф - некоторая функция.

к=1

2

2

Величины д гк /дхг из (3.4) образуют якобиан системы. Матрица Якоби:

л

пх

8 г1/ 8 х1 8 г1/ 8 х2 8 и / 8 х 8 и / 8 х.

В(г„ г 2) =

8 г„

а™^, г2) \г„ г2) «2!;у)(г1, г2) ^\гх, г2)

Элементы матрицы функционального определителя 0 = Пх / Пг находятся как

Д^, г 2) '(-ц г 2) ^2Г)(г1, г 2) ^) (-1, - 2) где ^'(-цг2) = Ддхр" + -з_9Д12х*;1р = 1,2; д = 1,2;

р = 1,2;

= Л;„.(г1,г2), где «к, = ^ ,(3.5)

ох

(3.6)

(3.7)

1 р

р = 1,2; А,х{,,]^ = х1 - х{;,]^ * 2 р р

(;',"_ „(¡+1,1+1)

Л,, х(У > = х

- х

р р ^) -л х(;1>

р = 1, 2 .

12 2 Формула вычисления определителя матрицы В(г1, г2):

/;, 1 (г1, г2) = ёе1[Л;,. (г, г2)] = /3^\гх, г2)р%\гх, г2)-р^\гх, г2)р% \ги г2). Матрица Якоби Л(г1, г2) имеет вид:

Л;, 1 (г^ г2) = в:\ (г,, г 2) =

(;,.) 1,1 (;,1),

(г„ г 2) «м лгц г 2) (г1 , г2 ) ^2,2 (г1, г2 )

(3.8)

(3.9)

(3.10)

где аЦ) =а17^х1, х2) =

Р^Л-и г 2)

Л. (г„ г г)

(;,л,

®1, 2 ' ~ ®1, 22 (х1, х2) _

ЛУЧ^, г 2) Л, (г„ г 2)

В(;,')(г г)

сс 21 — сс 21 (х1, х2) — — ; сх,

в(и]и г )

= Ч х1, х2) . (3.11)

¿и (г„ гг) 2,2 ' 1 JiJ (г„ г2)

В качестве основных неизвестных в узлах принимаются составляющие перемещений Щ км2км3к) и их производные у1< к\ v2k\ уЗ к) попеременной х3, т.е. для (р, д) -го

узла это и\

где

(к,р,д) ,,(к,р,д) ,,(к,р,д.) ,,(к,р,д) ,,(к,р,д) ,,(к,р,д)

V ', V и

23

(к,р,д)

(к,р,д) (к,р,д)

(хз) =

V = V

Формула восполнения перемещений на элементе:

1

(к,р,д) 2

(к,р,д)

3

и соответственно векторы неизвестных

(к,р,д)

-(к,р,д) _ Т7(к,р,д)

(хз) =

1

(к,р,д) 2

(к,р,д)

3

где \и(k,;,1у = и(к,;+1,1) - и

и {к)(гх, г2) = и(к,;,Л + г1А1и <k,;J' + г2 А2и <k,;J' + г1г2А12и

11 а2й(к,;-Л=м -Й(к,;'1

(к )

Формула восполнения производных от перемещений по х3 на элементе:

V(к )(г1, г2) = V) + г1А1у(к,;,]') + г2 А 2У(к,;^') + г1г2А12У

Т(к

(3.12)

(3.13)

(3.14)

(3.15)

где Ау(к,;,]} = V(к) - V{к,;,]}; Д,V= V(к,;- V

'(к,;,1).

А^(к,;'Л = V(k,;+lJ'+1' - V(к,;+1^) - Д^(к,;'Л . (3.16)

4. Вычисление частных производных от перемещений, деформаций и напряжений на элементе.

Частные производные от перемещений по г1 и г2 определяются по формулам:

7(к)

5г„

(г,г2) = Д й(к,;,л + -з_р^к,;,1\ Р = 1,2;

(4.1)

(x<,x2'e0k

2

А и {к,;," = й(к" - Й- А,Й{к,;,"

12

2

Формулы определения частных производных от перемещений по х1, х2 и х3:

7<к)

дх„

дх,

(11,0 = £

7<к)

Э1„

(11,12)<;Л(11,12), р = 1,2;

( х1,х2)е®к .¡.

7(к)

(11,12) = V

(к)

1( Х1,Х2)ЕЮк '

( Х1,Х2)ЕЮк '

Формула для компонент тензора деформаций на элементе:

471, о =

/к)

дх

(11,0 +

( х1,х2)ейк .¡.

,( к)

&

(11, 12)

( х1,х2)еак .¡. 7

р = 1,2,3; д = 1,2,3.

(4.2)

(4.3)

(4.4)

(4.5)

... ....................(4.6)

где 8 - символ Кронекера [4, 5]; Хк 7 =вклА; ~Цк=вк вк. - характеристическая функция соответствующего дискретно-континуального конечного элемента.

5. Вычисление квадратичной части функционала энергии на элементе.

Квадратичная часть функционала энергии на (¡, 7) -м элементе на фрагменте конструкции сок определяется по формуле:

1 1

Формула для объемной деформации:

^ >(!„ 12)=]г ?ркг )(11, и).

Р=1

Формула для компонент тензора напряжений на элементе:

Лк,',7)/

12) = 8р,Ли£<к^)(1„ 12) + )(1„ 12), р = 1,2,3; д = 1,2,3

М] М2 3 3

ф

(к,¡,7) _

ИХ! / (к,',л(С

2 М1М2 ,=1 ,=1 р=1 ,=1

(5.1)

где

/(к^(1„ 12) = 12); (5.2)

р=1 ,=1

М1 и М2 - количество точек интегрирования вдоль оси 011 и оси 012 соответственно (как правило, М1 = М2); ^' и - координаты точек интегрирования,

--0г-',2,...,-М^ 4}1 -и,...,М-; (5-3)

компоненты 12) и ,?ркд',7)(11,12) тензоров напряжений и деформаций определя-

ются по формулам (4.6) и (4.4) соответственно.

6. Связь локальной и глобальной индексации.

Вектор неизвестных на рассматриваемом (¡, 7) -м ДККЭ имеет вид:

и

(Ш) = и (7 х3) =

(к ,¡,7) '(к ,¡,7)

где

= ~<к,'^)(х3) =

V

—(к ,¡,7)' ~{к ,¡,7)

-(к ,¡,7)

Г(к,¡,7) _ к,¡,7)

(х3) =

Т(к ,¡,7) 1

Т(к ,¡,7)

Т(к ,¡,7)

(6.1) (6.2)

2

ВЕСТНИК 4/2011

Т(к,'',1) _ 77<к,''1

(Х3) =

(к,'',1) р

(к,'+1,1') Р

№,'',1+1) Р

(к,''+1,1+1) Р

7(к,'',1) _ 77(к,'',1)

(=

Р

у(к,''+1,1) Р

Р

у(к,''+1,1+1) Р

Р = 1,2,3

(6.3)

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

¿/(к->(Хз) « У(Хз), где У(Хз) = [ у у ... у24]т . (6.4)

Пусть ' - глобальный индекс элемента вектора (6.4). Поставим ему в соответствие локальные индексы.

Номер ]п соответствующего узла (условно вдоль оси Ох2):

]„ = ]■

' -1 _

12

где здесь и далее запись типа [п] обозначает целую часть числа п . Номер гп соответствующего узла (условно вдоль оси Ох1):

"' -12(]п - ') -1]

(6.5)

(6.6)

Индикатор типа неизвестной (= 1 - неизвестная компонента перемеще ния; = 2 - неизвестная производная от компоненты перемещения по х3):

" 1, -12('п - ') - 6('п -'') -1]

М -1 +

3

(6.7)

(6.8)

Номер р неизвестной компоненты

Р =-12(1 - Л - 6('п " О - 3(1п^ "1). Очевидным следствием формулы (6.8) является следующая:

= 12(]п -1) + 6('п -') + 3(1па -1) + Р . (6.9)

7. Алгоритм формирования матрицы жесткости дискретно-континуального конечного элемента.

Формирование матрицы жесткости ДККЭ производится методом базисных вариаций [4, 5]. Формула для определения элементов матрицы имеет вид:

(К(;,]) )ч, к = Ф(км > (ё + ён) - Ф(км' (ё) - Ф(к,',]) (ё1е) + Ф( км) (е0),

(7.1)

(7.2)

= 1,2,..., 24; = 1,2.....24,

где ёр, е0 - 24-х мерные векторы, элементы которых определяются по формуле

(вр\ =Зра; (во\ = 0, д = 1, 2,..., 24 ;

с>Рд - символ Кронекера [4, 5].

Таким образом, можно сформулировать описанный ниже алгоритм формирования матрицы жесткости ДККЭ.

1. На основании вышеприведенных формул (3.9), (3.11), (4.1)-(4.6), (6.5)-(6.8), (7.2) вычисляется значение Ф(к,',1}(ё0).

2. Последовательно перебираются ' = 1,2, ...,24. Для каждого фиксированного зна-

чения I выполняются следующие действия:

Р

Р

Р

Р

6

2.1. Вычисляется значение Ф<k,',-')(e¡ ) по формулам (3.9), (3.11), (4.1)-(4.6), (6.5)-(6.8), (7.2);

2.2. Последовательно перебираются у = 1,2, ...,24. Для каждого фиксированного

значения выполняются следующие действия:

2.2.1. Вычисляется значение Ф(к,',-"(ё. + ё. ) по формулам (3.9), (3.11), (4.1)-(4.6), (6.5)-(6.8), (7.2);

2.2.2. Вычисляется значение Ф(к,',л(ёу ) по формулам (3.9), (3.11), (4.1)-(4.6), (6.5)-(6.8), (7.2);

2.2.3. Вычисляется значение элемента матрицы жесткости ДККЭ (К'к',л). . по формуле (7.1).

Структура получаемой в результате матрицы жесткости К',у) ДККЭ следующая:

Гк1

КЦ,Л =

к 2

к;, К, к 4

-1,1 К! КЦ кЦ К Ц кЦ кЦ

"1,1 КЦ кЦ кц кЦ кЦ кЦ кЦ

"2,1 КЦ кц кЦ к:ц, кц кЦ кЦ

"2,1 к* к-: кЦ К"'I кЦ кЦ кЦ

"3,1 к* кЦ к£ кЦ кЦ кЦ

"3,1 К'1 кц кЦ КЦ кЦ кЦ кЦ

"4, 1 КЦ к£ кЦ КЦ кЦ к!'1 кЦ

"4,1 кЦ к£ кЦ кЦ кЦ кк кЦ

(7.3)

Здесь КЦ, КЦ, КЦ , К£™, /, т = 1,2,3, 4 - матрицы третьего порядка.

На основе матрицы жесткости ДККЭ (7.3) формируются соответствующие поэле-

ментные матрицы К

КЦ =

КЦ =

К

К

К

КЦ , кЦ' , КЦ > И кц

-1,1 кЦ КЦ КЦ

"2,1 кЦ кЦ КЦ

"3,1 КЦ кЦ к '1

"4,1 кЦ кЦ кц

1,1 кЦ кЦ кк1,,4,

2,1 кц; кЦ ккк2,,, 4,,

3,1 к ,, 1, кЦ К '1 ккк ,,, 4,,

КЦ кЦ ккк4,,, 4,,

12-го порядка

К^ =

ЧУ =

К1 К К К

К К

к

К

К Кк Кк К

1,3

2,3 3,3 4,3

К К К К

К1,

кЦ

К'3 км

КЦ кЦ кЦ кЦ

; (7.4)

(7.5)

Замечания. Исследования проводились в рамках следующих работ:

1. Грант НШ-8684.2010.8 Президента Российской Федерации для государственной поддержки ведущих научных школ Российской Федерации «Многоуровневые численные, аналитические и экспериментальные методы исследования прочности зданий и сооружений с учетом конструктивных и физических особенностей» на 2010-2011 гг.

2. Грант 2.3.9 Российской академии архитектуры и строительных наук «Разработка и исследование дискретно-континуальных методов для расчета строительных конструкций с кусочно-постоянными физико-геометрическими параметрами по одному из направлений» на 2011-2013 гг.

2

1

2

1

2

2

2

к

ВЕСТНИК 4/2011

3. НИР «Разработка теории и алгоритмов построения корректных аналитических решателей многоточечных краевых задач применительно к расчетам строительных конструкций», выполняемой по аналитической ведомственной целевой программе «Развитие научного потенциала высшей школы (2009-2011 годы)» (проект 2.1.2/12148).

Литература

1. Акимов П.А. Корректный метод точного аналитического решения многоточечных краевых задач расчета конструкций для систем обыкновенных дифференциальных уравнений первого порядка с кусочно-постоянными коэффициентами. // Вестник МГСУ, №1, 2011, с. 11-16.

2. Акимов П.А. Корректный метод точного аналитического решения многоточечных краевых задач расчета конструкций для систем обыкновенных дифференциальных уравнений второго порядка с кусочно-постоянными коэффициентами. // Вестник МГСУ, №4, т. 1, 2010, с. 24-28.

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

4. Золотов А.Б., Акимов П.А., Сидоров B.H., Мозгалева М.Л. Дискретно-континуальные методы расчета сооружений. - M.: Издательство «Архитектура - С», 2010. - 336 с.

5. Золотов А.Б., Акимов П.А., Сидоров B.H., Мозгалева М.Л. Дискретно-континуальный метод конечных элементов. Приложения в строительстве. - M.: Издательство ACB, 2010. - 336 с.

6. Секулович M. Метод конечных элементов. - M.: Стройиздат, 1993. - 664 с.

References

1. Akimov P.A. Correct Analytical Solution of Multipoint Boundary Problems of Structural Analysis for Set of First-order Differential Equations with Piecewise-constant Coefficients. // Bulletin MSUSE, #1, 2011, pp. 11-16 (in Russian).

2. Akimov P.A. Correct Analytical Solution of Multipoint Boundary Problems of Structural Analysis for Set of Second-order Differential Equations with Piecewise-constant Coefficients. // Bulletin MSUSE, #4, vol. 1, 2010, pp. 24-28 (in Russian).

3. Zienkiewicz O.C. The Finite Element Method in Engineering Science, London, Mc.Graw-Hill, 1971, 521 pages.

4. Zolotov A.B., Akimov P.A., Sidorov V.N., Mozgaleva M.L. Discrete-continual methods of structural analysis. Moscow, "Arkhitectura - S", 2010, 336 pages (in Russian).

5. Zolotov A.B., Akimov P.A., Sidorov V.N., Mozgaleva M.L. Discrete-continual finite element method. Applications in Construction. Moscow, "ASV", 2010, 336 pages (in Russian).

6. Sekulowicz M. The Finite Element Method. Moscow, "Stroyizdat", 1993, 664 pages.

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

Key words: multipoint boundary problem, set of first-order differential equations, piecewise-constant coefficients, correct analytical solution, discrete-continual methods, structural analysis.

Авторы:

1. Акимов Павел Алексеевич, доктор технических наук, профессор, член-корреспондент РААСН (ГОУ ВПО МГСУ); 129337, Россия, г. Москва, Ярославское шоссе, дом 26; тел./факс: +7(499) 183-59-94; e-mail: pavel.akimov@gmail.com.

2. Мозгалева Марина Леонидовна, кандидат технических наук, доцент, профессор (ГОУ ВПО МГСУ); 129337, Россия, г. Москва, Ярославское шоссе, дом 26; тел./факс: +7(499) 183-5994; e-mail: marina.mozgaleva@gmail.com.

3. Сидоров Владимир Николаевич, доктор технических наук, профессор, советник РААСН, заведующий кафедрой (ГОУ ВПО МГСУ); 129337, Россия, г. Москва, Ярославское шоссе, дом 26; тел./факс: +7(499) 183-59-94; e-mail: sidorov.vladimir@gmail.com.

Рецензент: Белостоцкий A.M., профессор, д.т.н., генеральный директор ЗАО «Научно-исследовательский центр СтаДиО»