N3, 2000

Использование эффективных алгоритмов решения больших систем линейных алгебраических уравнений в задачах геотехники.

К.Г.Шашкин

При решении геотехнических задач методом конечных элементов довольно часто встречается ситуация, когда расчетную схему невозможно свести к плоской или осесимметричной постановке и пренебрежение пространственным характером работы основания и конструкций здания приводит к грубым ошибкам. Однако пространственная постановка задачи для получения удовлетворительной точности в большинстве случаев требует решения больших систем линейных алгебраических уравнений (СЛАУ), состоящих из сотен тысяч уравнений. Причем задачи, моделирующие работу основания, много сложнее задач аналогичной размерности при решении стержневых или пластинчатых систем, которые используются для расчета конструкций. Действительно, например, при дискретизации расчетной области основания восьмиузловыми гексаэдрами каждый узел расчетной схемы оказывается связан с 26 окружающими узлами, что приводит к возникновению в матрице жесткости системы 26 ненулевых блоков в каждой строке. Для пластинчатых систем число таких блоков оказывается равным 8, а для стержневых систем оно зависит от числа стержней сходящихся в одну точку, но обычно не превышает 6. Кроме того, решение задач деформирования основания в большинстве случаев невозможно без учета нелинейного характера работы грунта. Решение таких задач приводит к необходимости многократного решения СЛАУ с меняющимся вектором сил или матрицей жесткости. При разработке программы "FEM models" была поставлена задача создать возможность применения трехмерных расчетных схем с учетом нелинейности в практических расчетах при существующем уровне развития вычислительной техники. Это означает, что такие расчеты должны выполняться, во-первых, на персональном компьютере среднего (достижимого для проектной фирмы) уровня, во-вторых, такие расчеты должны производится в ограниченное время. Иными словами, при современных требованиях к вариантному проектированию должна быть обеспечена возможность выполнить расчеты по нескольким вариантам расчетной схемы в течение одного рабочего дня. В противном случае расчеты с учетом пространственной работы основания и конструкций останутся единичными, применимыми только для анализа работы уникальных сооружений. Для достижения поставленной цели при создании модуля решения СЛАУ программы "FEM models" был произведен практический анализ большого числа алгоритмов решения СЛАУ. Следует отметить, что в связи с развитием вычислительной техники несколько изменился критерий оценки применимости многих алгоритмов. Действительно, в настоящее время объем оперативной памяти компьютера является скорее вопросом финансовым, чем техническим, поэтому определяющим является время решения задачи. Вопрос об объеме занимаемой памяти встает чаще всего в этом же аспекте, поскольку для обработки большого количества памяти требуется много времени. Применение различных методов решения СЛАУ тесно связано с принципом хранения матрицы системы в оперативной памяти. Традиционно для задач МКЭ применяется хранение матрицы в виде ленты, ширина которой зависит от максимальной разности номеров узлов схемы. Однако, как показывает практика, при решении пространственных задач такой принцип приводит к достаточно большой ширине полосы, большая часть которой заполнена нулевыми элементами. Другим способом хранения матрицы в оперативной памяти является хранение списков ненулевых элементов строк матрицы. В практических задачах такой принцип обеспечивает экономию памяти, однако он не столь эффективен, поскольку требует хранения большого количества индексной информации (например, о номерах столбцов элементов списка). В программе "FEM models" принят компромиссный принцип хранения матрицы системы. Матрица хранится по блокам, соответствующим узлам конечно-элементной сетки в виде списка со связями по строкам и по столбцам матрицы (рис. 1). Хранение в виде такого списка удобно для алгоритмов разложения матрицы на множители.

Рис. 1. Схема хранения матрицы.

Хранение матрицы по блокам также удобно тем, что блоки матрицы могут иметь различную размерность. Это избавляет от недостатка, присущего большинству конечно-элементных программ, а именно бессмысленного резервирования памяти по максимальному числу степеней свободы в узле. Например, при сопряжении стержневой системы (6 степеней свободы в узле) с пространственными элементами (3 степени свободы в узле) большинство программ создают матрицу с 6 степенями свободы во всех узлах, что приводит к необходимости фиктивного закрепления узлов, не примыкающих к стержням, от поворота вокруг трех осей. В программе "FEM models" принят следующий принцип определения размеров блоков матрицы: каждый конечный элемент выдает информацию о степенях свободы в своих узлах, по опросу всех элементов формируются списки различных степеней свободы в узлах, затем по размеру этих списков определяются размеры блоков. При грамотном написании программная реализация данного принципа работает практически так же эффективно, как и традиционное составление матрицы системы с фиксированным размером блоков. Следует отметить, что блочная структура хранения матрицы позволяет (с некоторыми корректировками) использовать обычные (не блочные) алгоритмы решения СЛАУ. Как известно, алгоритмы решения СЛАУ делятся на прямые и итерационные. К первым относятся различные схемы разложения матриц на множители. Наиболее часто употребляемыми являются алгоритмы Гаусса (для несимметричных матриц) и Холецкого (для симметричных матриц). Эти алгоритмы позволяют представить исходную матрицу системы в виде произведения двух треугольных матриц. Имея такое разложение легко получить решение с любой правой частью. Размер используемой памяти и время решения с использованием этих алгоритмов сильно зависит от сортировки строк и столбцов матрицы, т.е. от нумерации узлов в схеме. Существует несколько достаточно удачных алгоритмов квазиоптимальной нумерации узлов схемы. Наиболее эффективно эта задача решается в терминах теории графов. Целью перенумерации является уменьшение числа ненулевых элементов (так называемых новых ненулевых элементов), возникающих в матрице при применении алгоритмов Гаусса или Холецкого разложения матриц на треугольные множители. Часть алгоритмов сводит ненулевые элементы матрицы ближе к диагонали, уменьшая ширину полосы ленточной матрицы. Однако при разложении матрицы на треугольные множители происходит заполнение матрицы новыми ненулевыми элементами в пределах ширины полосы, что приводит к значительным затратам памяти и машинного времени. В книге А. Джорджа и Дж. Лю [1] приводится сравнение различных методов сортировки матриц. Наиболее эффективным и быстродействующим является так называемый алгоритм минимальной степени. В сочетании с техникой построения фактор-графов этот алгоритм позволяет получить перестановку, дающую при разложении Холецкого близкое к минимальному число новых ненулевых элементов. Как показывает практика, число новых ненулевых элементов при использовании этого метода оказывается в 2 и более раз меньшим, чем при использовании методов приведения к ленточному виду. Положительными чертами прямых методов является их безусловная сходимость и возможность одновременного решения задач с разными правыми частями. Однако при увеличении размерности системы время счета с использованием прямых методов быстро возрастает (см. рис. 2). При использовании различных алгоритмов сортировки матрицы и вычислительных машин различного быстродействия время решения, естественно, будет различаться, однако характер зависимости (резкое возрастание при увеличении порядка системы) при этом сохраняется. В результате для решения пространственных геотехнических задач большой размерности прямые методы практически неприменимы, так как требуют слишком большого времени для решения. Возможность решать одновременно несколько систем с разными правыми частями при рассмотрении нелинейных задач также является сомнительным преимуществом, поскольку при этом исключается применение эффективных алгоритмов решения с использованием переменной матрицы системы. В связи с указанными недостатками прямых методов решения СЛАУ был проведен анализ существующих итерационных методов. Было рассмотрено более 10 различных линейных и нелинейных итерационных методов (в том числе методы Якоби, Ричардсона, Гаусса-Зейделя, метод наискорейшего спуска, метод сопряженных градиентов и др. [2]). Как оказала программная реализация этих методов, самой высокой скоростью сходимости для большинства практических задач обладает метод сопряженных градиентов. Однако в целом, ни один из методов в исходном виде, несмотря на выполнение теоретических предпосылок для сходимости, не обеспечивал пригодной для практических расчетов скорости сходимости. Это связано с тем, что обычно матрица системы МКЭ имеет достаточно широкий спектр, что и обусловливает низкую скорость сходимости в большинстве итерационных методов. Поэтому многие авторы рекомендуют перед применением итерационных методов производить предварительное преобразование системы к системе где матрица обладает свойствами, улучшающими сходимость. Такое преобразование называется предобусловливанием. В частности, в [3] предлагается такое преобразование, которое позволяет снизить количество итераций при решении методом сопряженных градиентов задач небольшой размерности. Однако при повышении размерности задачи скорость сходимости снова становится низкой. Очевидно, что скорость сходимости большинства итерационных методов будет высокой, если матрица близка к единичной матрице. Для получения матрицы, близкой к единичной мы должны вычислить матрицу, близкую к . Очень простым и эффективным способом получения приближения к обратной матрице является так называемое неполное разложение Гаусса или Холецкого, которое отличается от обычного разложения тем, что не предусматривает резервирования новых ненулевых элементов. В результате получается разложение на треугольные множители некоторой матрицы , при этом достаточно близка к . Действительно, неизменной характеристикой матрицы системы для практически любых задач МКЭ является то, что диагональный член матрицы всегда больше по абсолютному значению, чем остальные члены матрицы. При появлении нового ненулевого элемента в разложении Гаусса или Холецкого его величина вычисляется следующим образом . Поэтому новые ненулевые элементы, заполняя строки матрицы, постепенно убывают по абсолютной величине. С этим связана возможность пренебречь возникновением новых ненулевых элементов, получив при этом хорошее приближение к . Получение матрицы в явном виде для использования итерационных методов не нужно. В большинстве алгоритмов достаточно решения системы , что легко выполняется, поскольку представлена в виде треугольного разложения. Метод преобразования системы с помощью неполного разложения Холецкого для симметричных систем при использовании метода сопряженных градиентов изложен в [4]. Однако, как показывает произведенный анализ, при использовании неполного разложения Гаусса такой метод применим и к несимметричным системам. Кроме того, при применении неполного разложения улучшается сходимость не только метода сопряженных градиентов, но и многих других итерационных методов. При решении практических задач МКЭ часто встречается ситуация, когда члены матрицы системы сильно отличаются по абсолютной величине (например, при введении в схему элементов большой жесткости). В этом случае неполное разложение дает матрицу не достаточно близкую к и итерационный процесс затягивается. В этом случае представляется целесообразным применять гибкую схему неполного разложения Гаусса или Холецкого. При применении такой схемы заполнение новых ненулевых элементов происходит в случае, если новый элемент по абсолютному значению превосходит средний отбрасываемый ненулевой элемент. Применение такого принципа позволяет снизить негативные качества системы с большим разбросом величин членов матрицы и в несколько раз уменьшить количество необходимых итераций для получения требуемой точности. Сравнение прямых и итерационных методов решения СЛАУ приведено на рис. 2.

Рис. 2. Сравнение методов решения СЛАУ.

При решении пространственных задач начиная с сетки 10х10х10 элементов итерационные методы значительно выигрывают как по времени решения, так и по объему используемой памяти по сравнению с прямыми. В целом применение рациональных способов хранения матрицы системы и эффективных итерационных алгоритмов решения СЛАУ позволило сделать реальным применение расчетов сложных пространственных конечно-элементных задач, например, с порядком системы до 200000 уравнений (что отвечает конечноэлементной сетке 40х40х40), в повседневной геотехнической практике.        

  Литература:  

1. А. Джордж, Дж. Лю. Численное решение больших разреженных систем уравнений. М. "Мир", 1984.

2. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М., "Наука", 1984

3. Расчеты машиностроительных конструкций методом конечных элементов. Справочник. Под общей редакцией В.И.Мяченкова. М.,"Машиностроение", 1989.

4. Дж. Голуб, Ч. Ван Лоун. Матричные вычисления. Перевод с английского. М. "Мир", 1999.