Сегерлинд Л. Дж. - Применение метода конечных элементов (1051193), страница 15
Текст из файла (страница 15)
Например, коэффициент — ь/м заключенный в квадрат матрицы (7.7), находится на пересечении второй строки и четвер- Уравнения (7.6), очевидно, не идентичны уравнениям (6.19в). Чтобы полученная матрица соответствовала точной матрице жесткости третьего элемента, ее нужно переформировать и расширить. Алгоритм переформирования и расширения матрицы несложен. Строкам и столбцам сокращенной матрицы элемента приписываются номера глобальных степеней' свободы. Порядок расположения степеней свободы соответствует обходу элемента против часовой стрелки, начиная от ь-го узла.
Матрицы элементов в задаче о кручении имеют только одну степень свободы (искомую величину) в каждом узле, поэтому функции формы в (7.1) упорядочены так же, как и глобальные степени свободы. Используя указанную нумерацию для строк и столбцов матрицы (7.4), запишем Глава 7 108 того столбца глобальной матрицы жесткости. Коэффициент — /а, 1! заключенный в треугольник, находится на пересечении четвертой строки и пятого столбца. Расположение всех коэффициентов матрицы элемента в глобальной матрице жесткости показано на 1 2 3 4 Фиг. 7.1. Расширение и переформирование сокращенной матрицы жесткости элемента. а — сомращениан матрица; б — расмиреннан н аереформнронаннан матрица.
фиг. 7.1. Незаполненные прямоугольники для 1л1а11 соответствуют нулевым элементам. Сопоставление фиг. 7.1 и соотношения (6.19в) показывает, что мы имеем теперь точную матрицу жесткости данного элемента. Метод прямой жесткости построения глобальной матрицы жесткости является очень важным алгоритмом реализации метода конечных элементов на ЭВМ, потому что он значительно сокращает загрузку запоминающего устройства. В частности, он исключает необходимость запоминания больших матриц элементов, которые содержат всего несколько ненулевых коэффициентов.
Число строк и число столбцов сокращенной матрицы жесткости элемента равны числу степеней свободы элемента. 7.2. Система линейных уравнений При использовании метода конечных элементов получается система линейных уравнений, которая должна быть решена относительно неизвестных узловых параметров.
Решение этих уравнений является очень важным аспектом задачи в целом, потому что си- 109 реалиэаиил метода конечных элементов на ЭВМ стема уравнений обычно очень большая. Методы решения систем с малым или большим числом уравнений мало отличаются друг от друга. Реализация этих методов, однако, зависит от технических возможностей ЭВМ. Во второй главе, где рассматривался процесс дискретизации сплошной среды, было отмечено, что путем надлежащей нумерации узлов можно контролировать расположение коэффициентов в ~ О О О О Сь О О О с см о о С С СммО С С О С С С О С С С С С С С С С С О О С С С С с с о с с О~ С С О О~ С О О О С С С С С С С С С С С О С О О~ С С С С С С О С С О О О Фиг.
7.2, Общий вид системы уравнений, получаемой при использовании методе конечных элементов. глобальной матрице жесткости. Напомним, что при разумной нумерации узлов получается матрица ленточного типа вместо полной матрицы. Ленточная матрица характеризуется тем, что все ее ненулевые коэффициенты располагаются вблизи главной диагонали, а все коэффициенты за пределами некоторой полосы, ограниченной линиями, параллельными главной диагонали, равны нулю. Схема.- тически это проиллюстрировано на фиг.
7.2, где ширина полосы ленточной матрицы показана штриховыми линиями. Через С обозначены ненулевые члены. Вообще говоря, нулевые коэффициенты могут встречаться и внутри полосы. Два свойства результирующей системы уравнений делают ее идеальной: симметрия и положительная определенность матрицы. Наличие симметрии означает, что приблизительно половину ненулевых членов матрицы можно не запоминать. Положительная определенность означает, что коэффициент, стоящий на главной диагонали, всегда положителен и обычно много больше по вели- ыо Глава 7 чине, чем любой другой коэффициент соответствующей строки или столбца. В случае симметричной положительно определенной матрицы ленточного типа значительно сокращается объем вычислений, необходимых для получения решения системы уравнений. К тому же уменьшается вероятность больших ошибок округления.
Существование симметрии в матрице ленточного типа позволяет значительно сократить объем памяти, требуемой для хранения глобальной матрицы. Обычно при программировании предусматривается превращение матрицы, изображенной на фиг. 7.2, в прямоугольный массив, ширина которого совпадает с шириной полосы матрицы, а длина равна числу уравнений. Чтобы проиллюстрировать преимущество такого представления матрицы, допустим, что мы решаем задачу, которая включает 200 узловых неизвестных.
Обычно при этом получается глобальная матрица жесткости, для хранения которой требуется 200р,'200, т. е. 40000 единиц машин,ной памяти. Однако, если эта ленточная матрица имеет ширину полосы, равную 40, и хранится в виде прямоугольного массива, требуется уже только 8000 единиц машинной памяти для запоминания 40 столбцов по 200 элементов в каждом. Таким образом, загрузка машинной памяти сокращается на 20а/а по сравнению с загрузкой, требуемой при хранении квадратной матрицы.
Решение системы уравнений может быть проведено с помощью алгоритмов, которые обсуждаются во многих книгах, посвященных численному анализу. Следует подчеркнуть, что обращение матрицы — очень неэффективная процедура решения системы уравнений. Эта неэффективность может быть объяснена двумя причинами. Обращение матрицы эквивалентно решению системы У уравнений с У неизвестными.
Если при этом рассматривается ограниченное число столбцов правых частей (глобальный вектор нагрузки), то вычисление обратной матрицы мало оправдано. Кроме того, в результате обращения ленточной матрицы получается матрица неленточного типа. Процедура обращения матрицы неэффективна также еще и с точки зрения экономии машинной памяти.
7.2.1. Преобразование системы уравнений Результирующая система уравнений. имеет вид (К((% =(Г)' (7.8) она получается суммированием уравнений для всех элементов. Эта система должна быть преобразована, если некоторые составляющие (Ф) известны, что является скорее правилом, чем исключением. В большинстве задач теории поля некоторые граничные значения искомой величины заданы; во всех задачах теории упругости должны быть фиксированы некоторые перемещения с тем, чтобы Реализанин метода конезнззх элементов на ЭВМ исключить перемещение среды как жесткого тела. В задачах механики деформируемых сред матрица жесткости (К] будет сингулярной до тех пор, пока не заданы некоторые перемещения.
Цель этого раздела — обсуждение и иллюстрация процедуры преобразования [К) и (Р) таким образом, чтобы получить правильный ответ, не изменяя размеры (К) и (г), ибо это повлечет за собой трудности при программировании. Если фиксирована одна степень свободы узлового параметра (Ф), то преобразование системы уравнений представляет собой двухшаговую процедуру. Пусть, например, известно значение Фз, преобразование сводится тогда к следующему: 1.
Все коэффициенты пятой строки, за исключением диагональных, приравниваются нулю. Диагональный член остается неизменным. В форме равенства это выглядит как Кз;=0 при 1=1, ..., и и 1Ф5. Соответствующая компонента Рз вектора(г) заменяется на произведение КззФз' 2. Все остальные уравнения преобразуются вычитанием произведения Кз Ф из Р; н подстановкой Кзз=О, 1= 1,..., и, 1Ф5. Пример 44.
Требуется преобразовать следующую систему уравнений, ес- ли известно, что Ф1=150 и Фз=40: 55Ф, — 46Ф, + 4Ф, = 500, — 46Ф + 140Ф,— 46Фз =2000, 4Ф 46Фз+ 110Фз 46Фз+ 4Ф =1000 — 46Ф, + 142Ф,— 46Ф,=2000, 4Фз 46Фз+ 65Фз 900. На первом этапе приравняем нулю все коэффициенты в первой и пятой строках, за исключением диагональных членов, которые оставим неизменными. Компоненты Р1 и Рз в(Р) заменим соответ- ственно на КпФ1 и КззФз. В результате будем иметь 55Ф, = 8250, — 46Фз+ 140Фз — 46Фз =2000, 4Фз 46Фз+110Фз — 46Фз+ 4Фз=1000 — 46Фз+ 142Фз — 46Ф,=2000, 65Фз = 2600.
112 Глава 7 Второй этап состоит в исключении столбцов матрицы, коэффициенты которых умножаются на Ф~ и Фз. Это осуществляется переносом членов, содержащих Ф~ и Фм в правую часть системы. Например, величина Рз становится равной 2000+46Фь или 8900. Завершая второй шаг, получим 55Фз =8250, 140Фз — 46Фз — — 8900, — 46Ф, + 11ОФ вЂ” 46Фз =240, — 46Ф + 142Ф4 — — 3840, 65Ф, = 2600. Описанная выше процедура совершенно проста и легко поддается программированию.