Основная задача подпрограммы формирования матрицы жесткости - составление матрицы жесткости элемента для исследуемой задачи. Такая подпрограмма обычно использует всю необходимую информацию из общих массивов памяти и в конце работы либо возвращает матрицу жесткости элемента в вызывающую программу, либо отсылает ее в периферийную память. Структура подпрограммы в большой степени зависит от математического описания жесткости элемента и, в частности, от того, необходимо ли для элемента численное интегрирование или возможно точное интегрирование.
Для простых элементов основными операциями являются: а) описание элемента в локальных координатах; . б) составление матрицы В, связывающей деформации с перемещениями (или ее эквивалента);
в) составление матрицы D, связывающей напряжения и де формации;
г) получение матричного произведения BDB;
д) интегрирование по пдощади элемента произведения ма триц (в случае плоского напряженного состояния эта операция сводится к простому умножению иа площадь треугольника);
е) выполнение при необходимости обратного преобразования полученной матрицы к глобальным координатам.
Выбор системы локальных координат обычно зависит от иС пользуемого способа получения матрицы жесткости. В простейших случаях (например, треугольный элемент при плоской деформации) система локальных координат может либо просто совпадать с глобальной, либо может быть получена путем пере носа начала координат в один из узлов или в центр тяжести элемента. На этой стадии, если необходимо, могут применяться /.-координаты или криволинейные координаты, описанные в гл. 7 и 8.
Подпрограмма формирования матрицы жесткости может ис пользоваться для построения матрицы напряжений, после умножения которой иа соответствующие узловые перемещения полу чаются напряжеиия элемента. Эта матрица обычно строится попутно при формировании матрицы жесткости, и получение ее не требует большого машинного времени. При этом возможны два варианта. Первый - составлять матрицу напряжений одновременно с матрицей жесткости й хранить ее до дальнейшего использования в накопителе (т. е. сначала вычислять произведение DB, а затем BDB), а второй - отдельно вычислять матрицу DB непосредственно перед использованием. Выбор того или другого способа зависит от скорости выборки данных из
накопителя и от времени, необходимого для повторного проведения некоторых вычислений.
Для сложных элементов, например изопараметрического типа, указанный порядок вычислений, как правило, неэкономичен. В этом случае можно использовать специальные приемы, описанные, в частности, Айронсом [2].
В зависимости от характера задачи основная программа может изменяться. Некоторые возможности описаны ниже.
а) При использовании элементов высоких порядков с дополнительными узлами на сторонах или групп малых элементов одной и той же конфигурации исходные данные могут задаваться только в некоторых определенных узлах, а для полного описания геометрии элемента могут использоваться специальные подпрограммы. Иногда элемент может содержать узел, не связанный с другими элементами (например, внутреннюю степень свободы). В таких случаях соответствующую степень свободы можно исключить, а в ансамбле использовать сокращенную матрицу жесткости.
б) Координатные оси можно выбирать по направлению перемещений (например, координаты, связанные с узлом, которые описаны в работе [3]), что влечет за собой необходимость преобразования матриц от узла к узлу. Применение таких координат приводит к тому, что матрицы жесткости элементов записываются в разных системах координат. Указанный подход удобен в следующих случаях:
1) при наличии узлов, в которых под некоторым углом к глобальным осям направлена одна (или несколько) удерживающая связь, что сильно затрудняет непосредственный учет граничных условий;
2) при учете симметрии или антисимметрии с целью уменьшения общего числа уравнений; например, при исследовании /а вместо пластины в случае двойной симметрии или при использовании неосесимметричных элементов в осесимметричных задачах [3].
3) при исследовании оболочек двойной кривизны, когда ось z в каждом узле направляется по внешней нормали к поверхности оболочки, так что для каждого узла требуется всего лишь пять степеней свободы [4].
Пример программы. Приведенная программа строит матрицу жесткости размерностью 6X6 для элемента, используемого при расчете плоской деформации, и записывает матрицу напряжений на магнитную ленту для дальнейшего использования при вычислении напряжений элемента. Блок-схема программы приведена на стр. 477.
Обозначения переменных в подпрограмме STIFT2 (N)
I, J, К -
AJ, BJ, АК, ВК А(3,6) *
ESTIFM (12,12)*
В (3,9)
Всюду обозначает номер элемента Параметры, определяющие связи элемента; позже используются как счетчики цикла Локальные координаты треугольника Матрица, связывающая деформации
с перемещениями Матрица, связывающая напряжения с деформациями; в дальнейшем используется при построении матрицы жесткости элемента; Матрица напряжений для обратного хода
Выход из программы при ошибке в задании связей 220 WRITE(6,I00)N
100 FORMAT(33HIZERO OR NEGATIVE AREA ELEMENT N014/21 HOEXECUTION ITERMINATED) STOP END
С С С С
с с с
29 30 31
32 33 34 35 36 37 38 39 40
42 43 44
46 47 48
50 51 52 S3 54
55 56 57 58 59 63
66 67 68
Блок-схема подпрограммы STIFT2
Начало
Размещение связей уолов
Вычисление размеров злементов
Проверка правильности нумции
Правильно
Непрцвально
Составление
матрицы связи
напряжений с деформациями
Вычисление мат
чццы напряокеиий
Запись матрицы напряженай на маенитную ленту NT4
Вычисление жесткости
матрицы злемента
Выдача сигнала ошибни
Останов
Возврат в основную программу
20.5. Составление ансамбля и решение уравнений
В любой программе, реализующей метод конечных элементов, ключевой является подпрограмма решения систем уравнений. Выбор метода решения зависит от числа уравнений задачи и от типа используемой вычислительной машины.
