Шабров Н.Н. - Метод конечных элементов в расчётах деталей тепловых двигателей (1061803), страница 21
Текст из файла (страница 21)
Входными параметрами подпрограммы являются массив адресов диагональных элементов ХРР, массив БСт, в котором хранятся элементы профиля матрицы Ь, массив ЭО Х диагональных элементов матрицы Р и параметр ИЬК, Х имеющий то же назначение, Х что и в подпрограмме 8СРЕС. /$()"Х,— ПЙ Симе. Выходным параметром яв- Х ляется массив РО, в котором хранится решение Ь) ~Ч % системы уравнений (8.1). Х Х Х Х Х ХХ В этой подпрограмме при ХХ ХХ значении параметра цикла ХХХХ М8Т = 1 осуществляется Рис. 8.3. Схема перемыожеиия профиля римов ход а пр значе" матрицы на вектор нии )ПЯТ = 2 — обратный. Суммирование в (8.10) и (8.11) выполняется также по укороченным циклам.
При решении нестационарпых задач теплопроводности формирование новой правой части системы уравнений теплового баланса на каждом временнбм шаге требует выполнения операции перемножения некоторой симметричной матрицы на вектор начальных значений. Поскольку хранятся только элементы профиля матрицы, то обычный алгоритм перемножения должен быть модифицирован с учетом этого обстоятельства. Можно заметить, что любой элемент профиля дважды участвует в перемножении. Рассмотрим элемент'1~0 профиля матрицы К, представленный на рис.
8.3. Во-первых, этот элемент участвует в перемножении на компоненту Л1 как элемент, расположенный в строке Ь Во-вто- 5ф 131 рых, этот жс элемент участвует в перемножении на компоненту Х; как симметричный элемент, расположенный в строке Л. Подпрограмма МЮСгХ реализует перемножение симметричной матрицы, профиль которой хранится в массиве Яб, на вектор Х. Результат перемножения размещается в массиве г'. Параметру ХЫ соответствует порядок матрицы. 8.3. ОРГАНИЗАЦИЯ МЕТОДА ХОЛЕЦКОГО С УЧЕТОМ ДОСТУПА К ПЕРИФЕРИЙНОЙ ПАМЯТИ ЭВМ Как отмечалось ранее, профиль матрицы системы уравнений МКЭ может содержать до нескольких сотен тысяч чисел и не может быть размсщен целиком в оперативной памяти ЭВМ, по- 10 15 1В 17 лот ИЯ л'сХ 15 14 15 15 Е= 147 Рпс.
8.4. Пример сегментаиии профиля исходной ма- трицы этому он разделяется на сегменты, в которых содержатся полные строки. Размер сегментов устанавливается таким образом, чтобы любые два сегмента не перекрывали друг друга в оперативной памяти ЭВМ. На рис, 8.4 показана симметричная часть глобальной матрицы конечно-элементной модели, составленной из плоских четырех!Зя 5 7 17 15 74 5а 5 41 55 55 54 75 75 1254 55755 Ю125455755 Л7 1 7 Ю 5 5 7 5 5 Ю г 5 5 5 75 угольных элементов. Нумерация узлов в целях иллюстрации выбрана заведомо неудачной. Профиль матрицы содержит 147 элементов, что соответствует адресу последнего диагонального элемента в массиве МРР. Допустим для определенности, что максимальный размер одного сегмента ограничен 50-ю числами.
Тогда профиль матрицы должен быть разделен на четыре сегмента прн условии, что в каждом из сегментов содержатся полные строки. Работа алгоритма метода Холецкого с учетом посегмептного представления профиля матрицы сводится к последовательному а) 6 в,! Рис. 5.5. Организация метода Холсцного прн сегментном хранении профиля матрицы: 1, 11 — пассивные сегменты; 1!1 — активные сегменты выполнению Л', шагов, где Л", — общее число сегментов. К началу выполнения 1, шага, где 1', — номер текущего сегмента, элементы которого подлежат определению вместе с соответствующими элементами диагональной матрицы Р, уже определены элементы профиля матрицы Е и Р исключительно до 1, сегмента.
При этом г,-й сегмент называется ггктавным, а все предыдущие сегменты называются пассивныхи. Процедура декомпозиции при определении элементов 1-й строки матрипы Е в соответствии с (8.8) требует перемножения элементов двух строк с номерами 1 и 1, одна из которых принадлежит пассивному ссгменту. Диапазон изменения номера 7' в больших задачах обычно перекрывает несколько пассивных сегментов.
Следовательно, для вычисления всех вспомогательных величин К1; на 1-м шаге декомпозиции должен быть организован цикл по пассивным сегментам. Для отметки начала цикла по пассивным сегментам каждому активному сегменту дополнительно к двум целым числам ХР1т и М1 1т„введенным для указания номеров первой и последней строк, ставится в соответствие третье число МГ8, обозначающее помер первого пассивного сегмента.
Первый пассивный сегмент для текущего активного сегмента должен содержать строку, номер которой совпадает с минимальным из номеров столбцов первых ненулевых элементов активного сегмента. Рпс. 8.5 иллюстрирует схему алгоритма метода Холецкого на третьем шаге декомпозиции матрицы из примера на рис, 8.4. 133 Третий сегмент профиля матрицы содержит две строки. Номер первой строки этого сегмента ЫГК = 18, номер последней строки М[.К 19. Первый столбец третьего сегмента содержит ненулевой элемент, поэтому цикл по пассивным сегментам должен начинаться с первого сегмента профиля и соответственно МГ5 = 1. Заполнение активного сегмента вспомогательными величинами К;~ происходит в диапазоне номеров столбцов от МРС до И1.С, где ХгС вЂ” номер первой строки, а Х1 С вЂ” номер последней строки пассивного сегмента.
В конце цикла по пассивным сегментам, когда активный сегмент совпадает с пассивным, вспомогательные величины К,1 преобразуются в элементы 1;; матрицы 1. с параллельным вычислением диагональных элементов д; матрицы О. 8.4. МЕТОД СОПРЯЖЕННЫХ ГРАДИЕНТОВ Метод сопряженных градиентов является прямым методом решения системы линенных алгебраических уравнений [46[. Однако при решении упомянутых систем уравнений на ЭВМ этот метод ведет себя как итерационный.
Это связано с нарушением ортогональности некоторых векторов вследствие ошибок округления. Рассматривая метод сопряженных градиентов как итерационный метод для решения больших систем уравнений с редкой матрицей, можно обнаружить некоторые полезные свойства. Так, например, быстрая сходимость для хорошо обусловленных задач позволяет получать с достаточной точностью итерационное решение за сравнительно небольшое число итераций. Реализация алгоритма метода сопряженных градиентов без непосредственной сборки глобальной матрицы системы уравнений приводит к исключительно простой вычислительной процедуре.
Метод сопряженных градиентов предназначен для решения систем уравнений вида (8.1). Существует несколько вариантов этого метода, обзор которых приведен в [46). Остановимся на одном из них. Выберем в качестве начального приближения вектор 13о и получим соответствующий вектор невязки г, го = г — К1-[о. (8.12) Полагаем для начального значения вектора сопряженных направлений р, = г,.
Тогда, пользуясь рекуррентными соотношениями для~ ° 0,1,2, а; =(г;, г;Ир;, Кр;); (8.13) 1),„= Ц+а;р;; (8. 14) г~„— — г~ — а;Кр;; (8.15) Ь, = (г;„, г;,~)~(г;, г;); (8.16) (8.17) находим реКторы Ц,», г1,~, р;,~ и скалярьи а1 и Ь1. 134 В приведенном алгоритме наиболее трудоемкой операцией является многократное вычисление произведения матрицы К на вектор р;, для которого введем обозначение (1 = Кр.
(8.18) Остановимся подробнее на анализе произведения вида (8,18). Ранее было принято, что глобальная матрица К системы уравнений (8.1) формируется на основе соответствующих матриц К(е] макроэлементов. Справедливо соотношение К ~,' а(~) К(~) а(~), Е (8.19) где суммирование происходиг по всем макроэлементам, а(е)— матрица кинематических связей макроэлемента, которая устанавливает соответствие между глобальным и локальным номерами узловой точки сетки макроэлемента.
Подставим (8.19) в (8.18) и получим (е)'К(е) (е) Е Введем 'обозначения (е) (е) й (е) К(е) (е) . В (1( ) — а( ) =а у Тогда (8.20) запишется так (8.21 (8.22) (8.23) (8.24) Операция перемножения (8.21) эквивалентна логическим операциям пересылки соответствующих глобальных компонент вектора р в вектор меньшей размерности Х(е) на место, определенное локальным номером этих компонент. Поскольку симметричная часть матрицы макроэлемента К(е) полностью хранится в оперативной памяти ЭВИ1 в виде профиля, то перемножение в (8.22) осуществляется подпрограммой М1.ЗОХ.
Операция перемножения в (8.23) является обратной по отношению к (8.21) и эквивалентной логическим операциям пересылки компонент вектора у(е) в вектор большей размерности (1(е) на место, определенное глобальным номером этих компонент. Разумно совместить операции пересылки в (8,23) с параллельным накоплением результата в глобальном векторе (1, 5УВКОБТ1МЕ М(]50Х(!ЧРР,5О,Х,У,Х!.К) Р!МЕЬ]5!ОМ ЫРР(1),5О(1),Х(1),У(!) РО 1И 1=1,]Ч1.Р М= !+1 — МРР(1) 1Г (1.ОТ. 1) М==М+Ь]РР(1 — 1) РО 1И )=М,! Я е=)чРР(1) У(1)=У(1)+ВО(ЫРГ) е Х(1) 1Г (1,!ЧЕ.З) 4 т())=У(!) — '60(ЯРЕ) е Х(!) 1И СОЧТ1М1.!Е КЕТБКМ ЕМЭ Рис. 8.6 показывает схему вычислений в (8.21) — (8.23) для четвертого макроэлемента из примера на рис. 7.3.
г( (((е/ а(((' (а(4 Л'Л'Ы ЯРЕ Х тф (Я Рис. 8.6. Схема шага вычисления произведений глобальной матрицы на вектор в методе сопряженнын градиентов Очевидно, при такой схеме перемножения матрицы на вектор нет необходимости Формировать массивы Х(~)1. и ХМТ номеров узловых точек сетки макроэлемента и в качестве локальных номеров использовать пос ледовательность натуральных чисел в процессе обхода узлов сетки слева направо и сверху вниз (см. рис. 7.1). 00 1И 1=1,)ЧРЯ 00 1И !Ч--1,ЫВР МР(.=ИВЕ е (1 — 1)-(-1Ч ЫРЙ=ХЭЕ э ()Ч.")С(1) — 1)+М 1И ХЕ(г!Р()=Р(ХРО) СЛ1Л.