Оценка наибольшего собственного числа глобальной матрицы жесткости регулярной сетки Текст научной статьи по специальности «Математика»

УДК 539.3:519.612+512.643

ОЦЕНКА НАИБОЛЬШЕГО СОБСТВЕННОГО ЧИСЛА ГЛОБАЛЬНОЙ МАТРИЦЫ ЖЕСТКОСТИ РЕГУЛЯРНОЙ СЕТКИ

П. М. Винник1, Т. В. Винник2

1 Балтийский государственный технический университет «ВОЕНМЕХ» им. Д. Ф. Устинова, Санкт-Петербург, Россия

2 Санкт-Петербургский государственный технологический институт (Технический университет), Санкт-Петербург, Россия

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

Ключевые слова: метод конечных элементов, регулярные сетки, матрица жесткости, собственные числа, теорема Гершгорина

Для цитирования: Винник П. М., Винник Т. В. Оценка наибольшего собственного числа глобальной матрицы жесткости регулярной сетки // Аэрокосмическая техника и технологии. 2024. Т. 2, № 1. С. 170-179. DOI 10.52467/2949-401X-2024-2-1-170-179. ЕРЫ EZDVJI

ESTIMATION OF THE GREATEST EIGENVALUE

OF THE GLOBAL STIFFNESS MATRIX FOR REGULAR MESH

P. M. Vinnik1, T. V. Vinnik2

1 Baltic State Technical University "VOENMEH", Saint Petersburg, Russia

2 Saint Petersburg State Technological Institute, Saint Petersburg, Russia

Abstract. The paper deals with condition of global stiffness matrices of regular mesh of finite elements. The estimation from above the greatest eigenvalue of such matrix offeres. The estimation is under construction on a local stiffness matrix of any finite element, hence, depends only from the size and the form of such element and does not depend on quantity of the final elements making a regular grid. At estimation construction the Gerschgorin's theorem and that fact are used that local matrixes of rigidity of final elements of regular grids differ from each other only shift of blocks. On a numerical example it is shown that the constructed estimation possesses split-hair accuracy and at a considerable quantity of the elements entering into a grid, it is possible to consider it almost coinciding with the greatest own number. The behaviour of an estimation is shown at change of quality of the form of final elements.

© Винник П. М., Винник Т. В., 2024

Keywords: a finite elements method, regular meshes, a stiffness matrix, eigenvalues, the Gerschgorin's theorem

For citation: Vinnik P. M., Vinnik T. V. Estimation of the greatest eigenvalue of the global stiffness matrix for regular mesh. Aerospace Engineering and Technology. 2024. Vol. 2, no. 1, pp. 170-179. DOI 10.52467/2949-401X-2024-2-1 -170-179. EDN EZDVJI (In Russian)

Введение

Одним из основных методов решения задач механики сплошных сред, в том числе для задач аэрокосмической области - прочности, газодинамики, является метод конечных элементов (МКЭ), в рамках которого задача сводится к решению системы линейных алгебраических уравнений (СЛАУ). Матрица СЛАУ предопределяется (до учета граничных условий) видом глобальной матрицы жесткости. Проблема собственных значений для матриц различного специального вида, хотя и решена в некоторых простых случаях [1, с. 178-180; 2], глобально продолжает оставаться актуальной. В ряде случаев вместо нахождения всех собственных чисел матрицы специального вида приходится ограничиться теми или другими оценками и локализацией спектра или его части [3-7].

Для построения СЛАУ МКЭ сначала вычисляются локальные матрицы всех конечных элементов (КЭ) е, затем каждая из локальных матриц вкладывается в матрицу [к(е)] большого порядка (для двумерного случая 2N х 2N, где N -количество узлов), и суммированием этих матриц большого порядка [8] находится вырожденная (нерегуляризованная) глобальная матрица жесткости:

[ КО ] = ^ [ к (е)],

ееЕ

где Е — множество всех КЭ. Матрица [КО] является вырожденной, имеющей ровно три нулевых собственных числа и (2N - 3) ненулевых (рассматривается двумерный случай, т. е. каждый узел имеет две координаты).

Рассмотрим треугольные КЭ. В [9] показано, что процедура регуляризации матрицы [КО] - перехода от [КО] к невырожденной матрице путем задания трех перемещений - является неканонической, и обоснована целесообразность применения в качестве меры обусловленности СЛАУ отношения наибольшего собственного числа матрицы [КО] к наименьшему положительному. В [2] установлено, что указанное отношение применимо с такой целью только асимптотически.

Для регулярных сеток КЭ (составленных из КЭ, имеющих одни и те же форму и размеры) вид матрицы [КО] предопределен видом матриц [к(е)]. Каждая матрица [К (е)] имеет только шесть ненулевых строк и столбцов. Составим

подматрицу из элементов матрицы [ К (е)], взятых в том же порядке, расположенных на пересечении ненулевых строк и столбцов:

[ К о(е)] =

41

21

к

V 31

12 22 ^32 лТ

13

23

^33 у

где к.. - матрица размера 2x2, причем [К0(е)]Т = [Ко(е)].

кп + к12 + к13 = 3,

к 21 + к 22 + к 23 = 3 ;

к31 + к32 + к33 = 3 :

(1)

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

Рис. 1. Фрагмент регулярной КЭ-сетки

Если вычислить сумму матриц [К(е1)], [К(е2)], [К(е3)], [К(е4)], [К(е5)], [К(еб)]

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

к11 + к 33 к12 к 32 к13 + к31 3 3 3 1

к 21 к33 + к 22 3 к 23 + к 32 к 31 3 3

к 23 3 к 22 + к11 к12 + к 21 3 к13 3

к13 + к 31 к 23 + к 32 к12 + к 21 2(кп + к22 + к33 ) к12 + к 21 к 23 + к32 к13 + к31

3 к13 3 к12 + к 21 к + к 22 11 3 к 23

3 3 к 31 к 23 + к32 3 к 33 + к 22 к 21

{ 3 3 3 к13 + к 31 к32 к12 к11 + к33 у

Заметим, что четвертая блочная строка матрицы к содержит все ненулевые элементы глобальной матрицы [КО], соответствующие перемещениям узла 5 на рис. 1. Поскольку все внутренние узлы рассматриваемой регулярной сетки входят в состав таких же шести КЭ, что и узел 5, то для каждого из внутренних узлов все ненулевые элементы глобальной матрицы [КО] содержатся в точно такой же блочной строке.

Рассмотрим теперь узлы, находящиеся на контуре области. Возможны три варианта. Во-первых, это угловые узлы (3 и 7 на рис. 1), входящие в состав единственного КЭ. Во-вторых, это угловые узлы (1 и 9 на рис. 1), входящие в состав ровно двух КЭ. В-третьих, это узлы, расположенные не в углах (2, 4, 6 и 8 на рис. 1), и входящие в состав ровно трех КЭ.

Для первого из этих случаев элемент блочной строки на главной диагонали матрицы [КО] - это к22, остальные ненулевые блоки - к21 и к23. Для второго случая элемент блочной строки на главной диагонали матрицы [КО] - это (к11 + к 33), остальные ненулевые блоки - к12, к32 и (к13 + к31). Для третьего случая элемент блочной строки на главной диагонали матрицы [Ш] - это (к11 + к22 + к33), остальные ненулевые блоки - к12, к21, (к13 + к31) и (к23 + к32) либо к23, к32, (к13 + к31)

и (к 21 + к12 ) .

Уточним, что номера выбираемых блоков для узлов, находящихся на контуре, зависят от локальной нумерации узлов в КЭ. Выбранные блоки соответствуют порядку узлов 2-5-1 на рис. 1.

Обусловленность СЛАУ с невырожденной матрицей оценивается по числу обусловленности (ЧО), которое для положительно-определенной матрицы можно определить, как отношение наибольшего собственного числа к наименьшему. Матрица [КО] является вырожденной. В [9] показано, что процесс перехода от [КО] к невырожденной матрице путем задания трех перемещений является неканоническим, т. е. его результат зависит от того, какие именно перемещения задаются. В [9] обосновано, что для устранения этой неканоничности можно расширить понятие ЧО на неотрицательно-определенные матрицы (т. е. матрицы, часть собственных чисел которых положительна, а часть - нулевая), полагая ЧО равным отношению наибольшего собственного числа к наименьшему положительному. Показано, что такое ЧО будет наименьшим из возможных получаемых при задании трех произвольных перемещений.

Для оценки этого ЧО необходимо оценить наибольшее собственное число матрицы [KG] сверху и наименьшее положительное - снизу. Простейшим инструментом оценки собственных чисел является теорема Гершгорина [10, 11]. В силу равенств (1) все круги Гершгорина содержат нуль, поэтому оценить снизу наименьшее положительное собственное число матрицы [КО] с их помощью напрямую невозможно. Однако оценить сверху наибольшее собственное число можно.

Во-первых, все собственные числа матрицы [КО] вещественны, т. е. лежат на отрезках, являющихся пересечением оси ОХ с кругами Гершгорина. Во-вторых, центрами кругов Гершгорина являются диагональные элементы матрицы [ КО], а радиусами - суммы модулей недиагональных элементов одной строки. Поэтому наибольшее собственное число не превосходит наибольшей координаты правых концов отрезков. Различных блочных строк всего пять (одна для внутренних узлов, две для угловых узлов и две разных (из-за перестановки индексов) для неугловых узлов), поэтому различных кругов Гершгорина всего 10, независимо от количества КЭ в регулярной сетке.

Вычислим 10 чисел - правые концы указанных отрезков, т. е. суммы координат центров и радиусов кругов Гершгорина:

А1 = 2( ^П + А22 + А33 )п + |(кП + А22 + А33 | + ^^ + А23 )п | + 1^32 + А23 ^2 | + + |(к12 + А21 )п| +|(к12 + А21 + |(А13 + А31 )п| + |(к13 + А31 )12 |) ;

А2 = 2( (к11 + А22 + А33 )22 + |(А11 + А22 + А33 )2^ + 1^32 + А23 )2^ + ^^ + А23 )2^ + + |(к12 + А21 )2^ + |(к12 + А21 )2^ + |(А13 + А31 )2^ + |(А13 + А31 )22 |) ;

В1 = (к22 )11 + \(А22 + \(А21 ^ | + |(А23 ; В2 = (к22 )22 + |(А22 )2^ + ^^21 )2^ + |(к23 )2^ ;

С1 = ( (к11 + к33 )11 + |(к11 + к33 + |(к32 + |(к32 + + |(к12 + |(к12 + |(А13 + А31 + |(А13 + А31 )12 |)

С2 = ( (к11 + А33 )22 + |(А11 + А33 )2^ + ^^ )2^ + ^^32 )2^ + + |(к12 )2^ + |(к12 )2^ + |(А13 + А31 )2^ + |(А13 + А31 )22 ^

А = ^11 + А22 + А33 )и + |(А11 + А22 + А33 ^2 | + |(к32 | + |(к32 | + 1(^23 | + |(к23 ^2 | + + |(А12 + А21 ^ +|(А12 + А21 ^ + |(А13 + А31 ^ + |(А13 + А31 ^ ;

А2 = (А11 + А22 + А33 )22 + + А22 + А33 )2^ + ^^ )2^ + 1^32 )2^ + ^^23 )2^ + ^^23 )2^ +

+ |(к12 + А21 )2^ + |(к12 + А21 )2^ + |(А13 + А31 )2^ + |(А13 + А31 )2^ ;

А3 = ^11 + А22 + А33 + |(А11 + А22 + А33 )12 | + 1(^2 | + |(к12 ^2 | + |(к21 | + |(к21 ^2 | + + |(к23 + А32 )п I + |(к23 + А32 ^2 I + |(к13 + А31 I + |(к13 + А31 ^2 I ;

D4 = (к11 + к22 + к33 ^21 +|(к11 + к 22 + к33 ^2 | + |(к12 1 + |(к12 ^2 | + |(к21 | + |(к21 ^2 | + + |(к23 + к32 )2^ + |(к23 + к32 )2^ + |(к13 + к31 )2^ + |(к13 + к31 )2^ '

Таким образом, если Л1 >Л2 >... >Л2ы, то

Ж = тах(^1, А2,В1,В2,С1,С2,,D2,Б3,D4).

Для демонстрации качества построенной оценки Ж рассмотрим численный пример. Взята прямоугольная область размера 65x13 со сторонами, параллельными осям координат, разделенная прямыми линиями, параллельными сторонам прямоугольника и равноотстоящими друг от друга, на 65x13 одинаковых квадратов, затем все квадраты однотипной диагональю были разделены на два треугольных КЭ. Таким образом, всего построено 1690 одинаковых треугольных КЭ, причем каждый треугольник являлся прямоугольным равнобедренным с катетами, длина которых равна 1.

Локальная матрица жесткости КЭ имеет вид [12, с. 220]:

[К0(е)] = [Вт ][D][B]Ft,

где [В] — матрица градиентов перемещений; [D] — матрица упругих констант материала; F — площадь КЭ; / — толщина КЭ,

(/

[ В] = ^

2F

Ъ 0 Ъ, 0 Ък 01

0 с 0 су 0 ск

с Ъ с Ъу ск Ъ)

[ D] =

Е (1 -V)

(1 + v)(1 - 2v)

1 V

V

(1 -V) 0

(1 -V) 1

0

0 0

1 - 2у 2(1 -V),

где Ь, Ъу, Ък, сг, е}, ск понимаются традиционно для МКЭ; матрица [О] взята для

плоского деформированного состояния; V - коэффициент Пуассона; Е - модуль Юнга. При расчетах примем v = 1 / 3. Модуль Юнга будет присутствовать множителем у всех коэффициентов А1, А2, В1, В2, С1, С2, D1, D2, D3, D4 и собственных чисел матриц [Ш] и [К0(е)]. Тогда получаем:

Л1 = 11,88797407Е;

А1 = А2 = 12Е, В1 = В2 = С1 = С2 = 3Е, D1 = D2 = D3 = D4 = 6Е;

Ж = тах(А1, А2, В1, В2, С1, С2, D1, D2, D3, D4 ) = 12Е.

Относительная погрешность оценки составляет всего около 0,9 %. Необходимо отметить, что суть в росте наибольшего СЧ матрицы с увеличением количества КЭ: оценка Ж не меняется, так как вычисляется по мат-

рице [к0(*)] одного КЭ, а наибольшее СЧ увеличивается. Укажем, что, хотя в рассматриваемом примере матрица [ ш ] имеет порядок 1848 (размер 1848*1848), при реальных реализациях МКЭ порядок [ко] часто достигает сотен тысяч. При таких порядках матриц [Ш] следует ожидать дальнейшего уменьшения относительной погрешности. Это значит, что можно принять наибольшее СЧ равным этой оценке W.

Из практики применения теоремы Гершгорина [12, с. 31-32] следует, что строчные и столбцовые круги Гершгорина дают различные по точности оценки СЧ, причем для симметричных матриц оценки совпадают. Для рассматриваемых матриц жесткости регулярных сеток КЭ некоторым аналогом такого различия может являться случай, когда значения, например, констант А1, А2, не совпадают друг с другом. В отличие от строчных и столбцовых кругов Гершгорина это - плохой случай, так как предложенная оценка W - максимум из нескольких чисел.

Рассмотрим, когда такой случай будет иметь место, и что тогда получается. В численном примере предполагалось, что катеты всех КЭ параллельны координатным осям и равны друг другу по длине. Предположим теперь, что один катет имеет длину 1Д/т, а другой - 4т, и, что первый катет образует с осью ОХ угол, равный аг^(к). Тогда случай т = 1, к = 0 соответствует разобранному численному примеру. Поле значений оценки W при разных значениях величин т и к приведено на рис. 2.

Рис. 2. Значения оценки Ш/Е при различных значениях отношения т длины второго катета КЭ к длине первого и тангенса к угла поворота первого катета относительно оси ОХ

Вследствие того, что константы А1,А2 (как и В1,В2,С1,С2,Ц,Б2,Д,Д ) по-разному реагируют на изменение значений величин т и к (рис. 3), отсутствует симметричность поля W относительно прямой к = 0. Существенный рост оценки W при отходе величины т от значения т = 1 объясняется ухудшением качества формы таких КЭ.

Рис. 3. Значения величин A1, A2 при различных значениях отношения m длины второго катета КЭ к длине первого (а) и тангенса к угла поворота первого катета относительно оси ОХ (б)

1 - AJE; 2 - A2/E

Заключение

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

Показано, что оценка имеет очень малую относительную погрешность даже для достаточно малого количества КЭ в регулярной сетке.

Конфликт интересов / Conflict of interests

Авторы заявляют об отсутствии конфликта интересов / The authors declare no conflict of interests.

Библиографический список

1. Даугавет И. К. Теория приближенных методов. Линейные уравнения: уч. пособие. СПб: БХВ-Петербург, 2006. 288 с.

2. Винник П. М., Винник Т. В., Хакимов А. А. Доказательство существования предела точности МКЭ при применении одинаковых одномерных линейных конечных элементов // Вычислительные технологии. 2022. Т. 27, № 1. С. 39-51. DOI: 10.25743/ICT.2022.27.1.004

3. Батагаева Т. А., Казаков А. Л. Аналитическое и численное исследование числа обусловленности квазитеплицевой трехдиагональной матрицы матрицы с неограниченной размерностью // Вестник ИрГТУ. 2015. Т. 99, № 4. С. 122-130.

4. Alvarez-Nodarse R., Petronilho J., Quintero N. R. Spectral properties of certain tridiagonal matrices // Linear Algebra and its Applications. 2012. Vol. 436, no. 3, pp. 682-698. DOI: 10.1016/j.laa.2011.07.040

5. Noschese S., Pasquini L., Reichel L. Tridiagonal Toeplitz Matrices: Properties and Novel Applications // Numerical Linear Algebra with Applications. 2013. Vol. 20, Iss. 2, pp. 302-326. DOI: 10.1002/nla.1811

6. Veerman J. J. P., Hammondy D. K., Baldiviesoz P. E. Spectra of Tridiagonal Matrices. URL: https://arxiv.org/pdf/1801.04977.pdf (дата обращения 02.11.2023).

7. Stukopin V., Grudsky S., Voronin I., Barrera M. Asymptotics of the eigenvalues of seven-diagonal Toeplitz matrices of a special form. URL: https://arxiv.org/pdf/2111.07196v1.pdf (дата обращения 02.11.2023)

8. Зенкевич О. Метод конечных элементов в технике: пер. с англ. М.: Мир, 1975. 542 с.

9. Иванов К. М., Винник П. М., Иванов В. Н. Численное моделирование разделительных процессов обработки давлением // Вестник Самарского государственного аэрокосмического университета. 2012. Т. 33, № 2. С. 192-198. EDN: PWOBVL

10. Уилкинсон Дж. X. Алгебраическая проблема собственных значений: пер. с англ. М.: Наука, 1970. 564 с.

11. Икрамов Х. Д. Несимметричная проблема собственных значений. Численные методы. М.: Наука, 1991. 240 с.

12. Сегерлинд Л. Применение метода конечных элементов: пер. с англ. М.: Мир, 1979. 392 с.

Дата поступления: 04.03.2024 Решение о публикации: 11.03.2024

Контактная информация:

ВИННИК Петр Михайлович - д-р техн. наук, заведующий кафедрой (Балтийский государственный технический университет «ВОЕНМЕХ» им. Д. Ф. Устинова, Россия, 190005, Санкт-Петербург, 1-я Красноармейская ул., д. 1), [email protected]

ВИННИК Татьяна Викторовна - канд. техн. наук, доцент (Санкт-Петербургский Государственный Технологический Институт (Технический университет), Россия, 190013, Санкт-Петербург, Московский пр., д. 26), [email protected]

References

1. Daugavet I. K. Teoriya priblizhennyh metodov. Linejnye uravneniya: uch. posobie [Theory of approximate methods. Linear Equations: Study Manual] St. Petersburg: BHV-Peterburg, 2006, 288 p. (In Russian)

2. Vinnik P. M., Vinnik T. V., Khakimov A. A. Finite element method: proof of the existence for the accuracy limit when applying one-dimensional linear finite elements. Computational Technologies. 2022. Vol. 27, no 1, pp. 39-51. DOL10.25743/ICT.2022.27.1.004 (In Russian)

3. Batagaeva T. A., Kazakov A. L. Analytical and numerical study of condition number of quasi-toeplitz tridiagonal matrix with unlimited dimension. Proceedings of Irkutsk State Technical University.. 2015. Vol. 99, no. 4, pp. 122-130. (In Russian)

4. Alvarez-Nodarse R., Petronilho J., Quintero N.R. Spectral properties of certain tridiagonal matrices. Linear Algebra and its Applications. 2012. Vol. 436, no. 3, pp. 682-698. DOI: 10.1016/j.laa.2011.07.040

5. Noschese S., Pasquini L., Reichel L. Tridiagonal Toeplitz Matrices: Properties and Novel Applications. Numerical Linear Algebra with Applications. 2013. Vol. 20, Iss. 2, pp. 302-326. DOI: 10.1002/nla.1811

6. Veerman J. J. P., Hammondy D. K., Baldiviesoz P. E. Spectra of Tridiagonal Matrices. URL: https://arxiv.org/pdf/1801.04977.pdf (accessed 02.11.2023).

7. Stukopin V., Grudsky S., Voronin I., Barrera M. Asymptotics of the eigenvalues of seven-diagonal Toeplitz matrices of a special form. URL: https://arxiv.org/pdf/2111.07196v1.pdf (accessed 02.11.2023)

8. Zienkiewiez O. C. The finite element method in engineering science. London, 1971. 521 p.

9. Ivanov K. M., Vinnik P. M., Ivanov V. N. Numerical simulation of separation processes in mechanical working. Vestnik of Samara University. Aerospace and Mechanical Engineering. 2012. Vol. 33, no. 2, pp. 192-198. EDN: PWOBVL (In Russian)

10. Wilkinson J. H. The algebraic eigenvalue problem. Oxford: Clarendon Press, 1965, 662 p.

11. Ikramov H. D. Nesimmetrichnaya problema sobstvennyh znachenij. Chislennye metody [Unbalanced eigenvalue problem. Numerical methods]. Moscow: Nauka, 1991, 240 p. (In Russian)

12. Segerlind L. J. Applied finite element analysis. New York: Wiley, 1976, 422 p.

Date of receipt: March 4, 2024 Publication decision: March 11, 2024

Contact information:

VINNIK Petr Mikhailovich - Doctor of Engineering Sciences, Associate Professor, Head of Department (Baltic State Technical University "VOENMEH", Russia, 190005, Saint Petersburg, 1st Krasnoarmeyskaya ul., 1), [email protected]

VINNIK Tatyana Viktorovna - Candidate of Engineering Sciences, Associate Professor (Saint Petersburg State Technological Institute, Russia, 190013, Saint Petersburg, Moskovsky pr., 26), vinnik.tv92@gmail .com