Алгоритмы формирования матриц коэффициентов и векторов правых частей граничных условий при использовании дискретно-континуального метода конечных элементов для расчета трехмерных конструкций

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

ALGORITHMS OF AUTOMATIC CREATION OF MATRICES OF COEFFICIENTS AND RIGHT-SIDE VECTORS OF BOUNDARY CONDITION WITHIN DISCRETE-CONTINUAL FINITE ELEMENT METHOD FOR THREE-DIMENSIONAL PROBLEMS OF STRUCTURAL ANALYSIS

П.А. Акимов, М.Л. Мозгалева Pavel A. Akimov, Marina L. Mozgaleva

ФГБОУ ВПО «МГСУ»

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

Effective computer-oriented algorithms of automatic creation of matrices of coefficients and right-side vectors of boundary conditions within discrete-continual finite element method for three-dimensional problems of structural analysis are under consideration in the distinctive paper. These algorithms are based on so-called method of extended domain, proposed by Alexander B. Zolotov.

1. Постановки задач.

Пусть х3 - координата, соответствующая основному направлению объекта; х34к -

координата рассматриваемой граничной точки, в которой задаются граничные условия; N - число узлов сетки дискретно-континуальных конечных элементов, аппроксимирующих поперечное (по отношению к основному направлению) сечение объекта по направлению осей йх1, Ох2 соответственно; Пк(х3) - искомая вектор-функция неизвестных, структура которой определяется следующими формулами:

ик = ик (Хз) = [U )T (V )T г

где

U = u(x2) = [(U^y (urv ... (ur-y ...

... (и(кд'2))т (Un(kA2))T ... (Un(k'N"2))T ... (U'nkXN2')T (U'nk'2'N2')T Vt = V(Х2) = [(v^V (vnk,2,V ... (v^-y ...

• •• (V<k12))T (V(kA2))T ... (v<kNi'2))T ... (vnk1N2))T (v<kAN2^T

)T

(u

(v

(k,N,,N2))T ]T

(k,N,, N2КT -|T.

)T ]T ;

n

и{к^> = м(

п п

(*з) =

.

(к,-,]) 2

. ) 3

—км) _ у

пп

(*з) =

1

,(к ) 2

,(к,/, ]) 3

м^,'-]), и?,;-1), и

нии и

значения неизвестных компонент перемеще-и их производных у^), —2к), у(к) попеременной х3 в (-',])-м узле дискретно-континуальной модели.

Пусть ст^х1, х2, х(), р = 1,2,3; д = 1,2,3 - компоненты тензора напряжений на

подобласти а>к с поперечным (по отношению к основному направлению) сечением ¡лк (х1, х2) и Лк (х1, х2) - параметры Ламе на подобласти а>к; [стР4]0,]) - значение

компоненты тензора напряжений а '

(-', ]) -м узле дискретно-

континуальной модели. Рассмотрим граничные условия типа «свободный край» и «идеальный контакт» [2].

Свободный край. Без ограничения общности положим, что граничные условия типа «свободный край» задаются в граничной точке х^ = 0 (т.е. к = 1). Тогда имеем соответствующие граничные условия следующего вида:

ви

где В1+, В~ - заданные матрицы коэффициентов граничных условий, квадратные 6NN2 -го порядка; , - заданные 6N1И2 -мерные векторы правых частей граничных условий;

Условия типа «свободный край» по поперечному к основному направлению сечению х( = х3ь1 задаются, как известно, следующими формулами:

<

Дискретно-континуальный аналог (1.6) имеет вид: [^ ) (хЬ,1 + 0) = 0, [о-™ Г,]) (хЬ,1 + 0) = 0, [<( ) (хЬ,1 + 0) = 0,

- = 1,2,..., N1 -1, ] = 1,2,..., N2 -1. (1.7) Идеальный контакт. Граничные условия типа «идеальный контакт», задаваемые в граничной точке х34к, представимы в виде:

в;йк_1

где В~, В+к - заданные матрицы коэффициентов граничных условий, квадратные 6N1N2 -го порядка; ~§+к - заданные 6N1N2 -мерные векторы правых частей граничных условий.

Условия типа «идеальный контакт» по поперечному к основному направлению сечению х3ьк задаются, как известно, следующими формулами:

,(к-1)/

,(к-1)/

^1<,Г) (х1, х2 , х(к - 0) = (х1, х2 , хС,к + 0) :

г(к-1) /

г( к-1'Ч

Дискретно-континуальные аналоги (1.9), (1.10) имеют вид:

2

2,(

ик) (хЬ,к - 0) = и^) (хЬ,к + 0), ик-1,(') (х3Ьк - 0) = и™) (х^ + о),

и?™^ -0) = '

КГ Г,я (х1 - 0) = [<]"-л (хЬ,к + 0), г>(хЬ,к - 0) = [< Г">(хЬ,к + 0),

[ст3

Требуется представить методику формирования матриц коэффициентов и векторов правых частей граничных условий в соответствии с их формами представления (1.5) и (1.8) [5].

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

Рассмотрим произвольный (', у) -й дискретно-континуальный конечный элемент (рис. 2.1). В его поперечном (по отношению к основному направлению) сечении вводится локальная система координат О11 и О12, при этом ^ е [0,1]; 12 е [0,1].

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

Имеет место соответствие глобальных и локальных координат узлов элемента: (х?Л х' ^ (0,0); (х;-+1,л, х2' ^ (1,0);

Пусть 1 = [ ^ t2Y и х = [ х1 х2 ]г - векторы координат произвольной точки дискретно-континуального конечного элемента в локальной и исходной глобальной системах координат.

Формула преобразования координат на элементе (рис. 2.1):

х

где Д.х' = !

^ 1 п п ' 2 п п ' 12 п п 2\/

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

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

и к)

где А1м(к,'' = и(к,'- и(к,"; А2М(к,'' = и (к,','+1) - и ;

Д12м

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

V(к '(г.,г2) = V (к,;,л + г1А1У(к,;,]' + г2 А 2У<к,;,]) + г.г2 Д12у

(к,;,]) _ - (к ,¡+1,.,') _ - (к,;,]) . д - (к,/, ]) _ - (к,/, ]+1) - (к,;,]) .

( к ,;,Л

где А1у (к,',]) = V(к,''+и) - V , ; А2У

= V4' - у

Д12Vк,-' = V

Функциональная матрица В(г., г2) имеет вид:

^

В(г., г2) =

где

А. х= х(,+и) - х, А2х= х' - х(-,Л

1 р р р ' 2 р р р

ДГ'

о,])

)('1, г2) = Д9х^> + г,А12х\ р = 1,2; д = 1,2;

р,д V 1 > 2 / д р 3-д 12 р

(;',]') _ Л+1,]+1) _ V(¡+1,]') _ л л

д12хр',]) = Xр+и+1) - хр,+1,]) - А2х-,'>, р = 1,2 . (2.10) Формула вычисления определителя матрицы В(г1, г2) :

,]'

Матрица Якоби Л(г., г2), определяемая формулой

л ¡, ]

имеет вид

"а™

Л ;,, ('1, О =

2 1,2 1 2 <)('1, О «2,2 )('1, '2 )_

Л?/ Чг., '2)

вы ]

(;,]) _ О',]')/ ' 2,2 \г1'г2/ , е,]') _ е,]') / \ _ /~1,

где 1 — сс.1 (х. , х2) — ; сс.2 — сс.2 (х. , х2) —

3;,/

/?<;,у)

(¡V) _ <-,Л( ) _ ' 2,1 У'1>'2/ ; <;,у) _ <;,j) ( ) _ И\,\ 2 '

&2 1 — 1 (X., х2 ) — ; 2 — £^2 2 (X., х2 ) — "

3;,]

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

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

7<к)

дг„

(г.,г2) = Д и<к,;,]) + гъ_рАпм(к,;,]), р = 1,2;

<х1 ,х2 )ек.;.

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

7<к)

5х„

(гl, г2) =

< х1 ,х2 )ек.;.

7<к)

5г„

5 _

дх,

-и

< к)

< х1 ,х2)ек.;. ] (г., г2) = V(к"

|< х1,х2)ейк.;.

-!< х1 ,х2 )ек.;'. ]

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

(г- -1 [— —1 л

^ )(г1, г 2) = -

р,д V 1' 2/ ^

5 <к) -м

5х р

(г., г2) +

< х1,х2)еОк.;'. ]

-и

дхр д

< к)

< х1,х2)ек.;.

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

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

,

(1г, Ч) = к?(1г, 1г).

?=1

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

о-'1,12) = 12) + 12), р = 1,2,3; ; = 1,2,3,

где ; - символ Кронекера(если р = ;, то 8р; = 1; если р Ф ;, то с>р; = 0).

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

Приведем ниже формулы для определения напряжений и деформаций в узлах дискретно-континуальной модели (ниже всюду р = 1,2,3 и ; = 1,2,3):

- внутренние узлы

+ &.

к: г у' = хш (уж^Н+<^ (0,1)0,

'р,; (1,0Ж,Ы,У  ~р,;

+ я

к,'-1,'-1 ^ С р,;

>(1,0)0^,, + <,л <^ (0,1)0, -

(0,0)0, ,,');

к ,',у-1

+ £

(0,0)0к',');

- «левая вертикальная граница»

[СГ = (^^к;1"-1)(0,1)0к,1,'._1 + <;1")(0,0)0кД,'); &' = (С-1)(0,1)^к' +^;Рк;1"')(0,0)0к,1,'.);

- «правая вертикальная граница»

[^(к)-|(Ы1,'-) _„ (_<к N -1,у-1) (1 1)Д | ^(к N -1,Л(1 0)Д )•

[^р,; ] ~Хк,Ы1,' (^р,; (1,1)^к,Ы1 -1,'-1 + ^р,; (1,0)^к,Ы1 -1,' );

- «нижняя горизонтальная граница»

[О™ =^к,',1(^1Г,1)(1,0)0к,'-1,1 +^^";,1)(0,0)0к,',1);

р,^ ^к ,',1 ^ р,;

(1,0)0к,

р,; к,',1 <кл,!)(0,0)0к,'Д);

- «верхняя горизонтальная граница»

[<;г2)у (^рг,Ы2 -1)(1,1)^к,'-[*£г 2)=*к,', у 2 "чуж

к ,'-1,Ы2 -1 + ^р,;

.

рл

<к,' ,N2 -1)

(0Д)0к,',Ы2 -1);

(о,l)0k,¡,ы2 -1);

- «угловые узлы границы» в случае ' = 1; у = 1 имеем:

<0, 0)0к,1,1; [<;г =^к,1,1^1к;1,1)(0, 0)0.

в случае ' = 1; у = N2 имеем:

<к )П<1, N2)

_

(0,1)0кД„ _1; [е

р,; ] %кД^2 ^р,;

[а{к>]

в случае ' = ; у = 1 имеем:

[<]<^ = Х^^К?1 "1Д)(l,о)0k,Ыl-1,1; [<;]<

в случае ' = Ы1; у = N2 имеем:

„<к,1,Ы2 -1)

(0,1)0к,

,< к, N1 -1,1) /

<к)]<N1,N2) _

I- р,; -1

< к, N1 -1, N2 -1)

(1Д)0к,

-1;

< к)]< N1, N2) _ р,;

41,0)0^ -1

^ -1,N2 -1)

к, N1 -1, N2 -1 '

12

12

12

А/ЭПИ ВЕСТНИК

5. Алгоритм формирования матриц и векторов граничных условий типа «свободный край».

Очевидно, что при задании условий типа (1.5) следует положить

0; Гк = 0.

Алгоритм формирования матрицы В* при задании условий (1.7) описан ниже.

1. Последовательно задаются ^ = 1,2,..., NN. Для каждого фиксированного значения Iг выполняются действия, перечисленные ниже.

1.1. Последовательно перебираются р = 1,2, где р - индекс, определяющий тип неизвестной (при р = 1 рассматривается компонента вектора перемещений й1, а при р = 2 рассматривается компонента вектора производных от перемещений по переменной х3, т.е. вектора V,). Для каждого фиксированного значения р выполняются действия, перечисленные ниже.

1.1.1. Последовательно перебираются I = 1,2,..., N . Для каждого фиксированного значения I выполняются действия, перечисленные ниже.

1.1.1.1. Последовательно перебираются у = 1,2,..., М2. Для каждого фиксированного значения у выполняются действия, перечисленные ниже.

1.1.1.1.1. Последовательно перебираются д = 1,2,3 . Для каждого фиксированного значения д выполняются действия, перечисленные ниже.

1.1.1.1.1.1. Вычисляется глобальный индекс

у = д + 3

1.1.1.1.1.2. В качестве вектора неизвестных Ц1(х^ + 0) (см. формулу (1.5)) задается

х^ + 0) = ё]г, где

8} д - символ Кронекера.

1.1.1.1.1.3. В соответствии со структурой (1.1)-(1.4) вектора неизвестных по соответствующим формулам из числа представленных выше (4.1), (4.3), (4.5), (4.7), (4.9), (4.11)-(4.14) определяются приведенные к узлу компоненты [о-1<13)](',у), [о^]0,у) и

[о^]0,7' тензора напряжений.

1.1.1.1.1.4. Соответствующие элементы матрицы В1+ определяются по формулам

6. Алгоритм формирования матриц и векторов граничных условий типа «идеальный контакт».

Очевидно, что при задании условий типа (1.8) следует положить

£;= 0; 0.

Алгоритм формирования матриц Вк и В~ при задании условий (1.11)-(1.12) описан ниже.

1. Элементы матриц Вк и Вк определяются по формуле

В)р,д = Зр,д, В)р,д = -8ГЛ, р = 1,2,...,3^М2, д = 1,2,...,3^М2.

Заметим, что вычисление элементов матрицы Вк и Вк по формулам (6.2) соответствует заданию граничных условий (1.11).

2. Последовательно задаются ге = 1,2,..., NN. Для каждого фиксированного значения выполняются действия, перечисленные ниже.

2.1. Последовательно перебираются р = 1,2, где р - индекс, определяющий тип неизвестной (при р = 1 рассматривается компонента вектора перемещений й1, а при р = 2 рассматривается компонента вектора производных от перемещений по переменной х3, т.е. вектора у1). Для каждого фиксированного значения р выполняются действия, перечисленные ниже.

2.1.1. Последовательно перебираются , = 1,2,..., ^ . Для каждого фиксированного значения , выполняются действия, перечисленные ниже.

2.1.1.1. Последовательно перебираются у = 1,2,...,^. Для каждого фиксированного значения у выполняются действия, перечисленные ниже.

2.1.1.1.1. Последовательно перебираются д = 1,2,3 . Для каждого фиксированного значения д выполняются действия, перечисленные ниже.

2.1.1.1.1.1. Вычисляется глобальный индекс у

' = д + 3

2.1.1.1.1.2. В качестве векторов неизвестных Цк (х3ьк + 0) и Цк_1(х3ьк - 0) задаются

ик

2.1.1.1.1.3. В соответствии со структурой (1.1)-(1.4) вектора неизвестных по соответствующим формулам из числа представленных выше (4.1), (4.3), (4.5), (4.7), (4.9), (4.11)-(4.14) определяются приведенные к узлу компоненты у), \&{2l~\{,,у),

[^Р> и [о-1<кГ1)]<,,у\ \ [^3<")\<,,у' тензора напряжений.

2.1.1.1.1.4. Соответствующие элементы матриц В^ и Вк определяются формулами (В+) = [а<к) ]<,,у'); (В+) = [а<к(В+) = [а<к)]<í,'К

Заключение.

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

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

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

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

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

Б/2011 ВЕСТНИК

4. НИР «Разработка, исследование и развитие математических методов и программно-алгоритмических средств, реализующих безопасную технологию возведения уникальных зданий и сооружений сложной архитектурной формы», выполняемая в рамках Тематического плана ГОУ ВПО МГСУ (проект 5.1.11).

Литература

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

2. Варданян Г.С., Андреев В.И., Атаров Н.М., Горшков А.А. Сопротивление материалов с основами теории упругости и пластичности. - М.: Издательство АСВ, 1995. - 572 с.

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

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

5. Zolotov A.B., Akimov P.A., Sidorov V.N. Correct Discrete-Continual Finite Element Method for Three-Dimensional Problems of Structural Analysis. // Journal of Beijing University of Civil Engineering and Architecture. Vol. 25, No. 2, Jun. 2009.

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. Vardanyan G.S., Andreev V.I., Atarov N.M., Gorshkov A.A. Strength of Materials with Foundations of Theory of Elasticity and Plasticity. Moscow, "ASV", 1995, 572 pages (in Russian).

3. 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).

4. 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).

5. Zolotov A.B., Akimov P.A., Sidorov V.N. Correct Discrete-Continual Finite Element Method for Three-Dimensional Problems of Structural Analysis. // Journal of Beijing University of Civil Engineering and Architecture. Vol. 25, No. 2, Jun. 2009.

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

Keywords: discrete-continual finite element method, static analysis, plate analysis, piecewise constant physical and geometrical parameters, method of extended domain.

Авторы:

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.

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