CC BY
1Математическое моделирование и его применение в научных исследованиях
УДК 519.673
DOI:10.25729/ESI.2024.35.3.012
Модификация алгоритма минимальной степени для повышения эффективности вычислений при итерационном решении контактной задачи теории упругости методом конечных элементов Дудаев Михаил Алексеевич
Иркутский государственный университет путей сообщения,
Россия, Иркутск, [email protected]
Аннотация. В статье предложена модификация классического алгоритма минимальной степени, предназначенного для упорядочения структуры глобальной матрицы коэффициентов системы линейных алгебраических уравнений метода конечных элементов, с целью повышения эффективности вычислительного алгоритма решения контактной задачи теории упругости при моделировании сопряжений деталей специальными контактными элементами сопряжения конструкций, применение которых обусловливает итерационное решение задачи с постепенным уточнением их состояния. Предложенная модификация реализована совместно с прямым методом Холецкого решения систем линейных алгебраических уравнений, вычислительная эффективность которого проявляется в наибольшей степени при решении систем с симметричной, положительно определенной матрицей коэффициентов, преимущественно имеющей в методе конечных элементов, разреженную структуру. Важной проблемой в методе Холецкого является заполнение «треугольного множителя» по сравнению с исходной матрицей в процессе разложения; уменьшение заполнения реализуется обычно с использованием известных эвристических алгоритмов, среди которых существенное распространение получил алгоритм минимальной степени. Модификация этого алгоритма, предложенная в работе, симметричной перестановкой строк и столбцов упорядочивает структуру исходной глобальной матрицы таким образом, что позволяет в процессе итерационного решения контактных задач теории упругости, проводить частичное разложение Холецкого на повторных итерациях с экономией времени расчета по сравнению с полным разложением. В работе показано, что структура матрицы, упорядоченная модифицированным алгоритмом, в целом, является менее предпочтительной, нежели при использовании классического алгоритма минимальной степени, так как претерпевает большее заполнение, однако, эффективность применения вычислительного алгоритма достигается с проведением повторных итераций при постепенном уточнении состояния контактных конечных элементов сопряжения конструкций. С использованием тестовой конечноэлементной модели сборного ротора, а также модели роторной системы реального авиационного газотурбинного двигателя, определены и показаны портреты исходных матриц и их треугольных множителей, проведены оценки количества арифметических операций и условного времени их выполнения при использовании классического и модифицированного алгоритма минимальной степени.
Ключевые слова: метод конечных элементов, СЛАУ, контактная задача, метод Холецкого, разложение, факторизация, алгоритм минимальной степени, скалярные произведения, заполнение
Цитирование: Дудаев М.А. Модификация алгоритма минимальной степени для повышения эффективности вычислений при итерационном решении контактной задачи теории упругости методом конечных элементов / М.А. Дудаев // Информационные и математические технологии в науке и управлении, 2024. - № 3(35). - С. 133-147. - Б01: 10.25729/ЕБ1.2024.35.3.012.
Введение. При решении задач механики твердого деформируемого тела широко используется метод конечных элементов (МКЭ) [1]. Его применение в задачах статики и динамики во многом связано с построением симметричных, положительно определенных и разреженных матриц [2] (жесткости, масс и демпфирования), являющихся матрицами коэффициентов в разрешающих системах линейных алгебраических уравнений (СЛАУ). Хранение таких матриц в памяти ЭВМ обычно реализуется с использованием компактных схем хранения [3].
Численное решение СЛАУ в МКЭ в задачах статики и динамики обычно выполняется с использованием прямых методов [4], выполняющих его с доступной машинной точностью за конечное и предсказуемое количество арифметических операций. Наибольшее распространение при решении СЛАУ с положительно определенными симметричными разреженными матрицами получил метод Холецкого [4-6], обеспечивающий решение за минимальное количество арифметических операций. В этом случае матрица коэффициентов [А] представляется в виде:
[ 4 = [ ¿][ ¿Г, (1)
где [Ь] - нижнетреугольная матрица, называемая треугольным «множителем» Холецкого, а процедура ее вычисления - разложением (факторизацией).
В работе [5] исследованы несколько алгоритмов вычисления треугольного множителя, однако, при использовании компактных схем хранения матриц [А] и [Ь] преимущество имеет алгоритм в форме скалярных произведений. В таком случае элементы I треугольного множителя Холецкого [Ь] определяются выражениями:
1-1
(
\
lü=^\aü-E '1 = П lji = J a J, "E lk 'ljk J = 1 + J •П
(2)
\ к=1 Ч( V к=1 У
где г и] - номер строки и столбца элемента соответственно; п - количество строк матриц [А] и [Ь]; а - элемент матрицы [А].
При выполнении первичного присваивания [Ь] = [А] порядок обращения к элементам матрицы [Ь] в процессе ее разложения в форме скалярных произведений показан на рис. 1.
Рис. 1. Порядок обращения к элементам матрицы [Ь] в процессе ее разложения
Недостатком прямых методов решения СЛАУ является «заполнение» матрицы [7], при котором в процессе разложения в матрице [Ь] появляются дополнительные ненулевые значения. Эта особенность при использовании компактных схем хранения в памяти ЭВМ не позволяет хранить треугольный множитель [Ь] в массивах, выделенных под исходную матрицу [А], а также увеличивает время решения СЛАУ из-за дополнительных арифметических операций.
Заполнение треугольного множителя зависит от структуры исходной матрицы [А], которую можно изменить симметричной перестановкой строк и столбцов [4, 5], что в математической интерпретации означает умножение исходной матрицы [А] слева и справа на матрицу перестановок [Р]:
[В] = [[А] [Р] (3)
В матрице перестановок [Р] в каждой строке присутствует всего один ненулевой элемент, равный единице; умножение слева на матрицу [Р]т обеспечивает перестановку строк исходной матрицы [А], а умножение справа на матрицу [Р] - перестановку столбцов.
Для уменьшения заполнения в треугольном множителе обычно применяются эвристические алгоритмы «минимальной степени» (АМС), «вложенных сечений», обратный алгоритм «Катхилла-Макки» и др., а также их комбинации [4, 5].
Классический алгоритм минимальной степени. При анализе динамического поведения многокаскадных сборных роторных систем турбомашин [8] с применением специально разработанного комплекса программ [9], реализованного на языке программирования FORTRAN 90 [10], используется АМС, показавший в этой задаче более эффективный результат уменьшения заполнения треугольного множителя [L], и, как следствие, позволяющий существенно экономить память ЭВМ [4], а также получивший наибольшее распространение [5].
Принцип работы алгоритма поясняется на примере конечноэлементной конструкции, показанной на рис. 2а, с использованием элементов теории графов [11, 12]. Соответствующий исходной структуре матрицы граф смежности показан на рис. 2б, где черным цветом указаны номера вершин графа (узлов КЭ), а красным - степени вершин.
в)
д)
Рис. 2. Исключение вершин графа смежности в разной очередности при его ассоциации с
узлами конечноэлементной модели Порядок исключения вершин графа без перестановки строк и столбцов показан на рис. 2в. После исключения вершины Xi добавляются ребра («клики») так, чтобы все вершины из множества Adj(Xi) оказались попарно смежными. Случай совпадения добавленного ребра с уже существующим свидетельствует о проведении математических операций с числом в процессе разложения, а случай появления новых ребер - о заполнении. Структура (портрет) матрицы жесткости при использовании начальных номеров узлов в качестве номеров строк и столбцов показана на рис. 2г, где ненулевые внедиагональные элементы исходной матрицы обозначены символом « х », а поддиагональные заполнения, получаемые в ходе разложения -« ® ». В соответствии с рис. 2г, без перестановки строк после разложения в матрице появляется 5 новых элементов под главной диагональю.
На рис. 2д показано исключение вершин в соответствии с АМС. После предварительного расчета степеней всех вершин производится исключение вершин 3 и 6, имеющих минимальную степень, равную двум. После их исключения добавляются ребра 2-5 и 1-4, которые совпадают с уже существующими ребрами. Степени оставшихся вершин пересчитываются и
в приведенном примере оказываются равными. В следующую очередь происходит исключение вершины 1 и добавляется ребро 2-4 (клика), поскольку каждая из этих вершин являлась смежной с исключенной, что свидетельствует о заполнении элемента /42. Процесс повторяется до исключения всех вершин. Структура матрицы, совмещенная с множителем Холецкого после упорядочения АМС показана на рис. 2, е.
Особенности решения контактной задачи теории упругости прямыми методами в МКЭ. При решении контактной задачи анализа динамического поведения многокаскадных сборных роторных систем турбомашин, математическая модель которой представлена в работах [13-15], используется специальный контактный конечный элемент сопряжения конструкций - КЭСК (рис. 3).
Рис. 3. Контактный конечный элемент сопряжения конструкций Данный КЭ может находиться в одном из двух состояний: «открытом» и «закрытом». Оценка состояний КЭСК проводится на основании невязки перемещений Ди его узлов г и у.
Аи = щ - щ + Щ , (4)
где иг и щ - перемещения узлов г и ] КЭСК соответственно в направлении оси х локальной системы координат; ио - начальный зазор/натяг КЭСК.
При условии Ди > 0 КЭСК находится в открытом состоянии, а при Ди < 0 - в закрытом. В зависимости от состояния КЭСК изменяется его матрица жесткости и вектор контактных сил:
[ к ]; =
[к]; -[к];
.-[кг [к];
[ к ]- =
[к]- -[к]--[к]- [к]-
к; о о
; [ к ]; = о к; о ; {f};= о
о о к
" к- о о "
] = о к- о ; {F}- = к--¿и -
о о
(5)
l о о -1 о о
(б)
где [ К и [ К ]- - матрицы жесткости КЭСК в открытом и закрытом состоянии соответственно; к* и к- - коэффициенты продольной жесткости КЭСК в открытом и закрытом состоянии соответственно [15]; к+у , к+г и к- , к- - коэффициенты поперечной жесткости КЭСК
в открытом и закрытом состоянии соответственно; {Р} и {Р} - вектор сил КЭСК в открытом и закрытом состоянии соответственно.
Смена состояний КЭСК в процессе статического или динамического анализа обусловливает итерационное решение контактной задачи с постепенным уточнением на каждой итерации состояния всех КЭСК. Блок-схема работы алгоритма решения контактной задачи МКЭ в статической постановке представлена на рис. 4.
Рис. 4. Блок-схема алгоритма решения контактной задачи теории упругости методом конечных элементов в статической постановки
Основную трудоемкость вычислительного алгоритма при итерационном решении контактной задачи составляет решение СЛАУ, а именно вещественное разложение Холецкого, проводимое на каждой итерации ввиду изменений в глобальной матрице жесткости из-за смены состояний КЭСК.
Для демонстрации указанной особенности рассмотрен пример, показанный на рис. 5 в виде системы из пяти стержневых конечных элементов (КЭ) с номерами 1, 2, 4-6 (номера КЭ обведены окружностями) с одинаковой продольной жесткостью к и одним КЭСК (№ 3) с осевой жесткостью кк, величина которой изменяется в зависимости от его состояния.
Рис. 5. Конечноэлементная модель стержневой системы с одним КЭСК
Матрица жесткости представленной системы с учетом закрепления узлов 1 и 7 в горизонтальном направлении (номера узлов показаны на рис. 5) имеет вид [1] (строки № 1 и 7 матрицы вычеркнуты):
2k -k О О О
-k к ^ kK -kK О О
О -kK к ^ kK -k О
О О -k 2к -k
О О О -k 2к
[ * ] =
Треугольный множитель Холецкого [L] для матрицы [K] имеет вид
О
(7)
[ L] =
42k -к >/2к
О О
О О
2kТ2
' к ^ 2ктт
\к ( к + ЗкК )
к ^ 2ктт
к (к + 2кК ) к (к + 4кК )
О
О
к ^ Зктт
О
к ^ Зктт
к (к + ЗкК ) к (к + 5кК )
к + 4k,.
к + 4k,.
(S)
В представленном примере видно, что изменение величины контактной жесткости кк повлечет за собой изменение величин элементов треугольного множителя, начиная с l22 и расположенных ниже и правее него под главной диагональю, поскольку во все эти элементы входит величина kK. Таким образом, даже в примере с одним КЭСК изменение его состояния приведет к необходимости практически полного повторного разложения Холецкого. При наличии множества КЭСК их состояния уточняются в процессе итерационного решения, поэтому разложения Холецкого будут вычисляться на каждой итерации, что существенно сказывается на времени расчета даже в статической постановке задачи. В задачах динамики [1] решение контактной задачи проводится на каждом шаге интегрирования по времени [13-15], что может неприемлемо увеличивать время расчета КЭ модели, особенно, на этапах ее отладки, когда решение задачи требуется проводить многократно.
Модификация алгоритма минимальной степени. В процессе отладки итерационного решения контактной задачи теории упругости МКЭ при анализе динамического поведения многокаскадных сборных роторных систем турбомашин был принят во внимание тот факт, что изменения коэффициентов жесткостей на каждой итерации происходят не во всех строках и столбцах глобальной матрицы жесткости, а только тех, в которые входят коэффициенты жесткости КЭСК. В связи с тем, что для упорядочения структуры глобальной матрицы жесткости уже был программно реализован АМС, возникла идея его модификации таким образом, чтобы при проведении упорядочения строки и столбцы, связанные с узлами КЭСК, симметрично переставлялись в конец матрицы (вниз и вправо). Для этого к реальным степеням узлов (соответствующим вершинам графа смежности), связанных с КЭСК, достаточно прибавить большое число, например, соответствующее размерности n матрицы. Это приводит к тому, что степени таких вершин графа смежности (узлов КЭСК) становятся больше, чем у остальных вершин. В таком случае продолжение работы АМС приводит к «принудительной» симметричной перестановке строк и столбцов матрицы, связанных с узлами КЭСК в ее конец.
Работа модифицированного таким образом АМС рассматривается на уже представленном примере (рис. 5). Соответствующий КЭ модели граф смежности показан на рис. 6а.
закреплено
исходный граф
а)
1 7 7 2 1
I этап 11 этап III этап
(2И?К?НТН® (¿-(¿KD
6 7 1 6 6 5
б)
6
х
3
х
Рис. 6. Исключение вершин графа смежности с использованием модифицированного
алгоритма минимальной степени Узлы 1 и 7 в КЭ модели закреплены, поэтому соответствующие им строки вычеркнуты из матрицы, а соответствующие узлам вершины - исключены из графа смежности до начала работы АМС. На первой итерации АМС определяются естественные степени вершин графа, показанные на рисунке над вершинами красным цветом. После этого модифицированным АМС определяются вершины, ассоциированные с узлами КЭСК, и к степеням этих вершин прибавляется размерность матрицы (равная пяти) так, чтобы полученные величины степеней были самыми большими в графе. Модифицированные степени показаны под графом смежности малиновым цветом. Затем, в соответствии с работой АМС, из графа исключаются вершины 2 и 6, степень которых минимальна и равна единице. На второй итерации модифицированного АМС пересчитываются реальные степени оставшихся вершин; к степеням вершин, ассоциированных с узлами КЭСК, прибавляется размерность матрицы; после этого исключается вершина 5, степень которой минимальна и равна единице. Аналогично проводится третья итерация, в ходе которой исключаются последние вершины 3 и 4. Портрет матрицы после упорядочения модифицированным АМС показан на рис. 6б.
При разложении Холецкого представленная в примере матрица не претерпевает заполнений как в исходном, так и в упорядоченном варианте, что обусловлено линейной структурой графа смежности. Однако в преобразованном виде матрица имеет вид:
[ * ] =
2k 0 0 -k 0
0 2k -k 0 0
0 -k 2k 0 -k
-k 0 0
0 0
-k
а ее треугольный множитель Холецкого [L]:
[ L] =
V2k 0 0
0 42k 0
0 -1 JF
-k 0 0
42k
0 0 -f
k ^ k^ k^ k ^ k^
0 0
(9)
0 0
2kГ2
k (k + 5kK )
'k + 2kK V3 (k + 2kK )
(10)
В таком представлении видно, что жесткости контактных элементов встречаются в матрице и треугольном множителе Холецкого, начиная только с элемента 44 (в соответствии с номерами строк и столбцов матрицы) и далее (снизу и справа от этой позиции).
Это обстоятельство, а также порядок обращения к элементам матрицы при ее разложении в форме скалярных произведений (рис. 1) свидетельствует о том, что при изменении состояний КЭСК повторное разложение следует проводить не для всей матрицы [K], а только для ее части, связанной с узлами КЭСК. В приведенном примере с использованием классического АМС при повторном разложении пришлось бы вычислить величины семи элементов треугольного множителя (8), а с использованием модифицированного АМС - только для трех (10).
Представленный подход был реализован в комплексе проблемно-ориентированных программ [9], разработанном для анализа динамического поведения многокаскадных сборных роторов турбомашин.
Численный эксперимент и анализ эффективности вычислительного алгоритма на конечноэлементных моделях сборных роторных систем. Для тестирования эффективности модифицированного АМС была разработана КЭ модель сборного ротора, эскиз которой показан на рис. 7а. Моделирование основных деталей ротора реализовано с применением КЭ типа «изопараметрический гексаэдр», контактных сопряжений - с применением КЭСК, моделирование крепежа деталей - с применением КЭСК, балочных КЭ Тимошенко и КЭ типа «невесомый стержень», а моделирование упругих опор - с применением упруго-демпферных КЭ [8]. КЭ модель представлена на рис. 7б.
Рис. 7. Эскиз и конечноэлементная модель тестового сборного ротора
Всего представленная КЭ модель включает в себя: 1888 узлов (48 из которых закреплены по трем поступательным степеням свободы) и 1200 КЭ: 864 восьмиузловых КЭ типа «изопараметрический гексаэдр», 48 упруго-демпферных КЭ, 248 КЭСК, 24 балочных КЭ Тимошенко, и 16 КЭ типа «невесомый стержень». Доля узлов КЭСК составляет 27% от общего числа незакрепленных узлов модели. С учетом «вычеркивания» закрепленных степеней свободы матрица жесткости имеет размерность n = 5520.
На рис. 8а показан портрет исходной матрицы жесткости (главная диагональ и поддиа-гональные элементы), упорядоченный с использованием классического АМС, содержащий 147.912 элементов, а на рис. 8б - портрет треугольного множителя Холецкого, содержащий 933.585 элементов. Построение портретов матриц реализовано с применением языка программирования Python [16] и библиотеки «Pillow» для работы с изображениями; каждый пиксель черного цвета обозначает наличие вещественного ненулевого числа в матрице. Всего представленная матрица 5520*5520 на главной диагонали и под ней содержит 15.237.960 элементов.
При разложении Холецкого было подсчитано количество арифметических операций: 5.520 извлечений квадратного корня, 130.686.523 умножений, 130.686.523 сложений и вычитаний, 922.546 делений - всего 262.301.112 операции.
б)
.."■■у ' \v,w ¡к . -
Рис. 8. Портреты матриц с использованием классического и модифицированного АМС С использованием дополнительных тестовых примеров было установлено время простых арифметических операций, приведенное в секундах для 1010 операций в табл. 1.
Таблица 1. Время выполнения 1010 арифметических операций (FORTRAN), сек
Операция Сложение Вычитание Умножение Деление Извлечение квадратного корня
Опыт 1 3,50 3,54 3,50 9,06 12,74
Опыт 2 3,45 3,48 3,45 9,03 12,81
Опыт 3 3,47 3,51 3,65 9,04 12,75
Опыт 4 3,47 3,50 3,44 9,05 12,75
Опыт 5 3,44 3,50 3,46 9,04 12,75
Среднее значение 3,49 9,04 12,76
Таким образом, условное время выполнения одной итерации разложения при использовании классического АМС (без учета других операций, например, присваивания, поиска элемента массива и т.п.) составляет 92,078 миллисекунды.
На рис. 8в показан портрет исходной матрицы жесткости, упорядоченный с использованием модифицированного АМС, содержащий 147912 элементов, а на рис. 8г - портрет треугольного множителя Холецкого, содержащий 1577346 элементов.
Количество арифметических операций при полном разложении Холецкого составляет: 5520 извлечений квадратного корня, 430524817 умножений, 430524817 сложений и вычитаний, 1566307 делений - всего 862621461 операции.
Условное время выполнения первой итерации (полного) разложения при использовании модифицированного АМС без учета других операций составляет 301,99 миллисекунды. Таким образом, время, затраченное при разложении Холецкого на первой итерации с упорядочением
структуры матрицы модифицированным АМС, в 3,28 раз превышает время, затраченное при разложении с упорядочением классического АМС.
Однако количество арифметических итераций и время повторных разложений оказываются несоизмеримо меньшими времени первой итерации. На рис. 8в и рис. 8г соответственно дополнительно выделены красным цветом портреты повторно раскладываемой части исходной матрицы и треугольного множителя Холецкого.
При этом в исходной матрице на главной диагонали и под ней содержится 21840 ненулевых элементов, а в треугольном множителе - 832416 ненулевых элементов. Количество арифметических операций при повторном (неполном) разложении Холецкого составляет: 1488 извлечений квадратного корня, 1022 умножений, 1022 сложений и вычитаний, 829441 делений - всего 832973 операции. Условное время выполнения одной итерации неполного разложения составляет 0,753 миллисекунды.
Таким образом, эффективность применения модифицированного АМС по сравнению с классическим АМС проявляется при решении задачи за четыре и более итераций (рис. 9).
„ 2000
с* 1500 т е
S 1000
а р
я
е р
В
500 0
1 1 1 1 1 1 Классический AMC Mодифицированный AMC
0 2 4 6 8 10 12 14 Количество итераций
1б 18 20
Рис. 9. Зависимость времени разложения от количества итераций при использовании разных методов упорядочения структуры исходной матрицы При анализе динамического поведения сборных роторов турбомашин сходимость решения контактной задачи, оцениваемая по относительному количеству КЭСК с неизменным состоянием по сравнению с предыдущей итерацией, достигается за 6-10 итераций [14]. Кроме того, в некоторых случаях возникают ситуации неустойчивого поведения КЭСК, при котором несколько таких конечных элементов изменяют свое состояние от итерации к итерации, не позволяя завершить решение. В таких случаях количество итераций решения контактной задачи ограничено (двадцатью) и может быть изменено на этапе настройки решения. Таким образом, использование модифицированного АМС для упорядочения дает ощутимую экономию времени расчета (в представленном примере - в 1,81 раз при 6 итерациях, в 2,98 раза при 10 итерациях и в 5,82 раза при 20 итерациях).
Также с использованием тестовой КЭ модели было исследовано изменение времени выполнения полного и частичного разложения Холецкого при изменении доли общего числа КЭСК, но с сохранением общего количества узлов и размерности матриц. Для этого независимо от физической стороны задачи было проведено разложение для КЭ модели с 32, 56.. .392 КЭСК с шагом 24 элемента (соответствующему количеству узлов в окружном направлении ротора) и оценено количество итераций, на которых одно полное и несколько неполных разложений при применении модифицированного АМС дают выигрыш во времени по сравнению с несколькими разложениями при применении классического АМС (рис. 10). Оценка количества итераций определялась на основе неравенства:
Т + (N-1) , < Nr , (11)
где T и t - соответственно время одного полного и неполного разложения Холецкого с применением модифицированного АМС; N - количество итераций; т - время одного полного разложения Холецкого с применением классического АМС.
Рис. 10. Количество итераций решения контактной задачи, при котором упорядочение модифицированным АМС выгоднее упорядочения классическим АМС в зависимости от
относительного количества КЭСК в модели Исходя из полученных результатов можно сделать вывод, что эффективность упорядочения структуры матрицы жесткости модифицированным АМС может проявляться даже при большом относительном количестве КЭСК в модели. Однако на практике могут встречаться случаи, в которых модифицированный АМС приводит к большим затратам времени, вызванным вмешательством в работу основного АМС и неоптимальным (с точки зрения данного алгоритма) упорядочением.
Представленный в работе модифицированный АМС был применен для упорядочения структуры матрицы жесткости при анализе динамического поведения многокаскадной сборной роторной системы реального авиационного ГТД (рис. 11).
Рис. 11. Эскиз и КЭ модель, совмещенная с твердотельной моделью с вырезом одной четверти многокаскадной сборной роторной системы реального авиационного ГТД Всего представленная КЭ модель включает в себя: 46596 узлов (104 из которых закреплено по трем поступательным степеням свободы), 30488 КЭ: 22568 восьмиузловых КЭ типа «изопараметрический гексаэдр», 104 упруго-демпферных КЭ, 3900 КЭСК, 1200 балочных КЭ Тимошенко, 2340 КЭ типа «невесомый стержень» и 376 КЭ типа «сосредоточенная масса». Доля узлов КЭСК составляет 25,6% от общего числа незакрепленных узлов модели. С учетом закрепленных степеней свободы матрица жесткости имеет размерность n = 139476.
На рис. 1 2а показан портрет исходной матрицы жесткости (главная диагональ и поддиа-гональные элементы), упорядоченный с использованием классического АМС, содержащий 3791880 элементов, а на рис. 12б - портрет треугольного множителя Холецкого, содержащий 31239756 элементов.
Рис. 12. Портреты матриц с использованием классического и модифицированного АМС для
КЭ модели роторной системы реального ГТД На рис. 12в показан портрет исходной матрицы жесткости, упорядоченный с использованием модифицированного АМС, содержащий 3791880 элементов, а на рис. 12г - портрет треугольного множителя Холецкого, содержащий 52541100 элементов.
Условное время выполнения полного разложения на первой итерации при использовании модифицированного АМС составляет 23,6151 с., время неполного разложения на последующих итерациях - 0,025014 с., а время полного разложения на каждой итерации при упорядочении классическим АМС составляет 6,4247 с. Таким образом, время, затраченное на разложение Холецкого в первой итерации при использовании упорядочения модифицированным АМС, в 3,68 раз, превышает время, затраченное при упорядочении классическим АМС. Эффективность применения модифицированного АМС, по сравнению с классическим АМС, проявляется при решении задачи за четыре и более итераций (рис. 13). 160
с
Í 120 чсра 80
е р
m
40
- ...... Классический AMC Ыодифицированный AMC
0 2 4 6 8 10 12 14 Количество итераций
1б 18 20
Рис. 13. Зависимость времени разложения от количества итераций при использовании разных методов упорядочения структуры исходной матрицы для КЭ модели роторной
системы реального ГТД
0
Заключение. При сравнении результатов, полученных в работе, видно, что применение классического АМС упорядочивает структуру матрицы коэффициентов СЛАУ таким образом, что существенно позволяет уменьшить ее заполнение при разложении Холецкого, и, как следствие, уменьшить время разложения, что повышает эффективность вычислительного метода. Однако решение контактных задач итерационным способом можно проводить более эффективно (с меньшими вычислительными затратами) за счет проведения неполного разложения Холецкого на повторных итерациях при условии, что количество узлов, связанных с КЭСК, существенно меньше общего количества узлов КЭ модели; упорядочение структуры матрицы коэффициентов СЛАУ в таком случае выполняется модифицированным АМС. Таким образом, при решении контактной задачи за несколько итераций применение модифицированного АМС в совокупности с частичным разложением может существенно повышать эффективность вычислительного алгоритма. Поскольку при этом структура матрицы упорядочивается не самым эффективным с точки зрения классического АМС способом, может возникать чрезмерное заполнение. Поэтому в программной реализации решения [9] проводятся два символических разложения Холецкого (с применением классического и модифицированного АМС) и оценивается меньшее время решения задачи. Однако даже для КЭ модели роторной системы реального ГТД с большим количеством элементов множителя Холецкого и арифметических операций для его получения более высокая эффективность достигается при применении модифицированного АМС с учетом решения контактной задачи за 6-10 итераций [14].
Список источников
1. Bathe K.J. Finite element procedures. Upper Saddle River, New Jersey, Prentice hall, 1996, 1038 p.
2. Гантмахер Ф.Р. Теория матриц / Ф.Р. Гантмахер. - Изд. 5-е. - Москва : Физматлит, 2004. - 559 с.
3. Расчеты машиностроительных конструкций методом конечных элементов: Справочник / В.И. Мяченков, В.П. Мальцев, В.П. Майборода и др.; под общ. ред. М.И. Мяченкова. - М.: Машиностроение, 1989. - 520 с.
4. Фиалко С.Ю. Прямые методы решения систем линейных уравнений в современных МКЭ-комплексах /С.Ю. Фиалко. - М.: Издательство СКАД СОФТ, Издательство Ассоциации строительных вузов (АСВ), 2009. - 160 с.
5. Джордж А. Численное решение больших разреженных систем уравнений: Пер. с англ. / А. Джордж, Дж. Лю. - М.: Мир, 1984. - 333с., ил.
6. Исследование структурных свойств алгоритма разложения Холецкого: от давно известных фактов до новых выводов / А.В. Фролов, В.В. Воеводин, И.Н. Коньшин, А.М. Теплов // Вестник Уфимского государственного авиационного технического университета, 2015. - Т. 19. - № 4(70). - С. 149-162.
7. Пирова А.Ю. Параллельные алгоритмы и программные средства для уменьшения заполнения фактора симметричных разреженных матриц: диссертация на соискание ученой степени кандидата технических наук / Пирова Анна Юрьевна, 2022. - 143 с.
8. Дудаев М.А. Математическое моделирование динамического поведения двухроторных систем турбома-шин с межвальным подшипником / М.А. Дудаев, А.А. Пыхалов // Труды МАИ, 2024. - № 2.
9. Свидетельство о государственной регистрации программы для ЭВМ № 2019617798. Конечноэлементный решатель задачи роторной динамики одно- и двухвальных турбомашин с контактным взаимодействием деталей и межвальными связями: заявка № 2019616638 от 11.06.2019 РФ / М.А. Дудаев (РФ) - Зарегистрировано в Реестре программ для ЭВМ 20.06.2019 г. (РФ).
10. Бартеньев О.В. Современный Фортран. - 3-е изд., доп. и перераб. / О.В. Бартеньев. - М.: ДИАЛОГ-МИФИ, 2000. - 449 с.
11. Белоусов А.И. Дискретная математика / А.И. Белоусов, С.Б. Ткачев. - 7-е издание, исправленное. - Москва: Издательство МГТУ им. Н. Э. Баумана, 2021. - 704 с. - ISBN 5-7038-1769-2
12. Герб А.Р. Применение теории графов в алгебраических многосеточных методах для решения разреженных СЛАУ / А.Р. Герб, Г.А. Омарова // Проблемы информатики, 2022. - № 3(56). - С. 77-89.
13. Дудаев М.А. Контактная задача в анализе динамического поведения сборных роторов турбомашин / М.А. Дудаев, А.А. Пыхалов // Вестник НГТУ, 2015. - № 3. - С. 113-129.
14. Милов А.Е. Контактная задача динамики сборных роторов турбомашин : диссертация на соискание ученой степени кандидата технических наук / Милов Александр Евгеньевич, 2007. - 174 с.
15. Пыхалов А.А. Контактная задача статического и динамического анализа сборных роторов турбомашин: диссертация на соискание ученой степени доктора технических наук / А.А. Пыхалов Анатолий Александрович, 2006. - 405 с.
16. Бизли Д. Python: Подробный справочник / Д. Бизли. - 4-е изд.. - Санкт-Петербург: Символ, 2010. - 858 с. Дудаев Михаил Алексеевич. Старший преподаватель Иркутского государственного университета путей сообщения. Основные направления исследований автора: механика деформируемого твердого тела, динамика, метод конечных элементов. AuthorID: 688571, SPIN: 7575-9125, [email protected], 664074, г. Иркутск, ул. Чернышевского, 15.
UDC 519.673
DOI:10.25729/ESI.2024.35.3.012
The minimum degree algorithm modification for improving the computational efficiency in the iterative solution of the contact problem of elasticity theory using the finite element method Mikhail A. Dudaev
Irkutsk state transport university, Russia, Irkutsk, [email protected]
Abstract. The article proposes a modification to the classical minimum degree algorithm, designed to optimize the structure of the global coefficient matrix in the system of linear algebraic equations used in the finite element method. The aim is to improve the efficiency of the computational algorithm when solving contact problems in elasticity theory, specifically when modeling the connection of parts using special contact elements that determine the iterative solution with a gradual refinement of the state. The proposed modification is implemented in conjunction with the direct Cholesky method for solving systems of linear algebraic equations. The computational efficiency of this method is demonstrated to the greatest extent when applied to systems with a symmetric, positive-definite coefficient matrix, which often has a sparse structure when using the finite element method. A significant problem in the Cholesky decomposition method is the filling of the "triangular multiplier" relative to the original matrix during the decomposition process. The reduction of this filling is typically achieved using well-known heuristic algorithms, among which the minimum degree algorithm has become widely used. The modification of the algorithm proposed in this work, by symmetric permutation of rows and columns of the initial global matrix, reorders the structure of this matrix in a way that, during the iterative solution of contact problems of elasticity theory, allows for partial Cholesky factorization to be performed in repeated iterations, thereby saving calculation time compared to performing a full factorization. The paper demonstrates that overall the matrix structure ordered by the modified algorithm is less preferable than that produced by the classical minimum-degree algorithm, since it is undergoing more filling. However, the efficiency of the computational method is achieved through repeated iterations that gradually specified the state of the contact finite elements of the structures coupling. Using a test finite element model of a rotor, as well as a model of the real aviation gas turbine engine, portraits of the initial matrices and their triangular multipliers are determined and presented, estimates of the number of arithmetic operations and the estimated time required for their execution are calculated using both a classical and a modified minimum degree algorithm.
Keywords: finite element, finite element method, contact problem, Cholesky method, minimum degree algorithm, scalar products, filling, symmetric permutation
References
1. Bathe K.J. Finite element procedures. Upper Saddle River, New Jersey, Prentice hall, 1996, 1038 p.
2. Gantmaher F.R. Teorija matric [Matrix theory]. Izd. 5-e., Moskva, Fizmatlit [Ed. 5th., Moscow, Fizmatlit], 2004, 559 p.
3. Mjachenkov V.I., Mal'cev V.P., Majboroda V. P. et al. Raschety mashinostroitel'nyh konstrukcij metodom konech-nyh jelementov: spravochnik [Calculations of machine-building structures by the finite element method: Handbook]. M., Mashinostroenie [Mechanical engineering], 1989, 520 p.
4. Fialko, S. Ju. Prjamye metody reshenija sistem linejnyh uravnenij v sovremennyh MKJe-kompleksah [Direct methods for solving systems of linear equations in modern FEM complexes]. M., Izdatel'stvo SKAD SOFT, Iz-datel'stvo Associacii stroitel'nyh vuzov (ASV) [M., Publishing house SKAD SOFT, Publishing house of the Association of construction universities (ACS)], 2009, 160 p.
5. Dzhordzh A. Chislennoe reshenie bol'shih razrezhennyh sistem uravnenij: Per. s angl. [Numerical solution of large sparse systems of equations: transl. from English]. M., Mir, 1984, 333p.
6. Frolov A.V., Voevodin V.V., Kon'shin I.N., Teplov A.M. Issledovanie strukturnyh svojstv algoritma razlozhenija Holeckogo: ot davno izvestnyh faktov do novyh vyvodov [Study of structural properties of cholesky decomposition: from well-known facts to new conclusions] Vestnik Ufimskogo gosudarstvennogo aviacionnogo tehnich-eskogo universiteta [Bulletin of the Ufa state aviation technical university], 2015, vol. 19, no. 4(70), pp. 149-162.
7. Pirova A.Ju. Parallel'nye algoritmy i programmnye sredstva dlja umen'shenija zapolnenija faktora simmetrichnyh razrezhennyh matric: dissertacija na soiskanie uchenoj stepeni kandidata tehnicheskih nauk [Parallel algorithms and software tools for reducing the filling factor of symmetric sparse matrices : dissertation for the degree of Candidate of technical sciences], 2022. - 143 p.
8. Dudaev M.A., Pyhalov A.A. Matematicheskoe modelirovanie dinamicheskogo povedenija dvuhrotornyh sistem turbomashin s mezhval'nym podshipnikom [The mathematical modeling of the dynamic behavior of twin-rotor turbomachine systems with an inter shaft bearing]. Trudy MAI, 2024, no. 2.
9. Svidetel'stvo o gosudarstvennoi registratsii programmy dlia EVM №2019617798. Konechnoelementnyi reshatel' zadachi rotornoi dinamiki odno- i dvukhval'nykh turbomashin s kontaktnym vzaimodeistviem detalei i mezh-val'nymi sviaziami: zaiavka №2019616638 ot 11.06.2019 RF [Certificate of state registration of the computer program No. 2019617798. Finite element solver of the problem of rotor dynamics of single- and two-shaft turbomachines with contact interaction of parts and inter-shaft connections: application No. 2019616638 dated 06/11/2019 of the Russian Federation]. Dudaev M.A. Zaregistrirovano v reestre programm dlia EVM 20.06.2019 g. (RF).
10. Barten'ev O. V. Sovremennyj Fortran [Modern Fortran]. M., DIALOG-MIFI, 2000, 449 p.
11. Belousov A.I., Tkachev S.B. Diskretnaja matematika [Discrete mathematics]. 7-e izdanie, ispravlennoe. Moskva, Izdatel'stvo MGTU im. N. Je. Baumana [7th edition, revised. Moscow, Publishing house of Bauman Moscow state technical university], 2021, 704 p., ISBN 5-7038-1769-2.
12. Gerb A.R., Omarova G.A. Primenenie teorii grafov v algebraicheskih mnogosetochnyh metodah dlja reshenija razrezhennyh SLAU [Application of graph theory in algebraic multigrid methods for solving sparse SLAEs]. Prob-lemy informatiki [Problems of informatics], 2022, no. 3(56), pp. 77-89.
13. Dudaev M.A., Pyhalov A.A Kontaktnaja zadacha v analize dinamicheskogo povedenija sbornyh rotorov turbomashin [The contact problem in the analysis of the dynamic behavior of modular turbomachine rotors]. Vestnik NGTU, 2015, no. 3. pp. 113-129.
14. Milov A.E. Kontaktnaja zadacha dinamiki sbornyh rotorov turbomashin : dissertacija na soiskanie uchenoj stepeni kandidata tehnicheskih nauk [The contact problem of the dynamics of prefabricated rotors of turbomachines : dissertation for the degree of candidate of technical sciences]. Milov Aleksandr Evgen'evich, 2007, 174 p.
15. Pyhalov A.A. Kontaktnaja zadacha staticheskogo i dinamicheskogo analiza sbornyh rotorov turbomashin : dissertacija na soiskanie uchenoj stepeni doktora tehnicheskih nauk [The contact problem of static and dynamic analysis of combined rotors of turbomachines : dissertation for the degree of doctor of technical sciences]. Pyhalov Anatolij Aleksandrovich, 2006, 405 p.
16. Bizli D. Python: Podrobnyj spravochnik [Python: A detailed reference guide]. Sankt-Peterburg, Simvol [Symbol], 2010, 858 p.
Dudaev Mikhail Alekseevich. Assistant professor Irkutsk state transport university, The main areas of research of the author: mechanics of a deformable solid body, dynamics, finite element method. AuthorID: 688571, SPIN: 75759125, [email protected], 664074, Irkutsk, Chernyshevsky str., 15.
Статья поступила в редакцию 13.05.2024; одобрена после рецензирования 01.10.2024; принята к публикации 10.10.2024.
The article was submitted 05/13/2024; approved after reviewing 10/01/2024; accepted for publication 10/10/2024.
