Главная страница сайта
Российские промышленные издания (узловые агрегаты)
1 ...
21 22 23 [
24 ]
25 26 27 Блок-схема программы GDATA
Блок-схема подпрограммы LOAD
Считывание и печать тнтрольныя! данных
Считывание и печать харантеристик мат^иала
Считывание ноординат узлов
Считывание информации о свяои между элементами и типе алвмента
Считывание граничных условий
Нет
Нужно ли печатать | исходные даш1ые ? |
| | |
Печать ноординат узлов |
Печать информации о сВяви между злементами
Печать граничных условий
Возврат в основную программу
С С С
С С С
|
| | |
Задание нулевого вектора нагрузки |
| | |
Печать заголовка |
| | |
| |
Считывание инд)ормаиии о нагрузке с пердзонарты и печать |
| | |
Засылка значений в.вектор нагрузки |
| | |
Меньше ли номер узла махсималькр возможного ? |
| | нет |
Возврат в основную программу |
Считывание информации об элементе
READ(5,3) (N,(NOP(N,M).M=. 1,3),IMAT(N).N=. 1,NE)
Считывание граничных условий
READ(5,4) (NBC(r).NFrX(I).r=. 1,NB) ГР(П.1ЧЕ.О) GO TO 500
С С С
Печать введенных данных WRrTE(6,102)
WRITE(6.2) N.(CORD(N,M),M=. 1,2).N= 1,NP) WRITE(6,103)
WRITE(6,3) (N,(NOP(N,M).M=- l,3),IMAT(N),N=i l.NE) WRrTE(6,104)
WRITE(6,4) (NBC(r),NFIX(I),I=i l.NB)
40 42
52 53
500 CONTINUE
1 FORMAT(915) - 2 FORMAT(I10,2FI0.3)
3 FORMAT(615)
4 FORMAT(215)
7 FORMAT(12A6)
8 FORMAT(I10,2F10.2) 100 F0RMAT(1HI,12A6)
102 FORMAT(20H0 NODAL POINTS)
103 FORMAT(20H0 ELEAIENTS)
104 FORMAT(21H0 BOUNDARY CONDITIONS) 108 FORMAT(lH0,20H MATERIAL PROPERTIES)
RETURN END
Программа 20-3
С С С
SUBROUTINE LOAD
C0MM0N/C0NTR/T1TLE(12),NP,NE,NB,NDF,NCN.NLD,NMAT, NSZF,LI,NT4
COMMON CORD(100,2),NOP(200,4),IMAT(200),ORT(25,2).NBC(25),
NFIX(25) 1,R1(200).SK(200,40) 2,R(3)
Задание нулевого массива нагрузки
DO 160 J= l.NSZF 160 R1(J) = 0.
WRITE(6,100) TITLE.LI WR1TE(6,109)
I Считывание с перфокарты, печать и засылка в память ; информации о нагрузке
165 CONTINUE READ(5,9) 1 NQ,(R(K),K=1,NDF) WRITE(6,9)
=1,NDF) : 1,NDF 1C = (NQ- 1).NDF + K 170 Rl(lC) = R(K) + Rl(tC)
1 NQ,(R(K),K DO 170 K =
С С С С
Если помер узла ие равен максимальному NP, то возврат и продолжение считывания
1F(KQ.LT.NP) GO ТО 165 9 FORMAT(110,3F10.2) 20 FORMAT(10X,4I5)
100 F0RMAT(1H1,12A6,5X, lOHLOAD CASE, 15) )09 FORMAT(lH0,6H LOADS)
RETURN
58 71 72
75 77
80 81 82 83
10 11 12 13
20 21 22 23 24 25 26 27
9 28 39 40 41 42 43 44
20.4. Формирование матрицы жесткости
Основная задача подпрограммы формирования матрицы жесткости - составление матрицы жесткости элемента для исследуемой задачи. Такая подпрограмма обычно использует всю необходимую информацию из общих массивов памяти и в конце работы либо возвращает матрицу жесткости элемента в вызывающую программу, либо отсылает ее в периферийную память. Структура подпрограммы в большой степени зависит от математического описания жесткости элемента и, в частности, от того, необходимо ли для элемента численное интегрирование или возможно точное интегрирование.
Для простых элементов основными операциями являются: а) описание элемента в локальных координатах; . б) составление матрицы В, связывающей деформации с перемещениями (или ее эквивалента);
в) составление матрицы 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)
Всюду обозначает номер элемента Параметры, определяющие связи элемента; позже используются как счетчики цикла Локальные координаты треугольника Матрица, связывающая деформации
с перемещениями Матрица, связывающая напряжения с деформациями; в дальнейшем используется при построении матрицы жесткости элемента; Матрица напряжений для обратного хода
С С С
С С С
Программа 20-4
SUBROUTINE STIFT2(N) 1
C0MM0N/C0NTR/TITLE(I2),NP,NE,NB,NDF,NCN,NLD,NMAT, NSZF,LI,NT4
COMMON CORD(100,2),NOP(200,4),rMAT(200),ORT(25,2),NBC(25),
NFIX(25) 1,R1(200),SK(200,40) 2, ESTIFM(12,12),A(3,6),B(3,9)
Определение связей элемента
I = N0P(N,1) 11
J = N0P(N,2) 12
k; = N0P(N,3) 13
L = IMAT(N) 14
Образование системы локальных координат
AJ = C0RD(J,1)-C0RD(I,1) 15
AK = C0RD(K,1)-C0RD{1,1) 16
BJ = C0RD(J,2) - C0RD(I,2) 17
BK = C0RD(K,2) - C0RD(I,2) 18
AREA = (AJ.BK-AK.BJ)/2 19
IF(AREALE.O.) GO TO 220 20
Формирование матрицы связи деформаций с перемещениями
A(1,1)=.BJ-BK з
А(1,2) = 0. 24
А(1,3) = ВК 25
А(1,4) = 0. 26
A(1.5) = -BJ 27
А1,6) = 0. 28
С С С
A(2,l) = 0. А(2,2) = АК - AJ А(2,3) = 0. А(2,4) = - АК А(2,5) = 0. А(2,6) = AJ A(3,l) = AK-AJ А(3,2) = BJ - ВК А(3;3) = -АК А(3,4) = ВК; А(3,5) = AJ А(3,6) --BJ
Формирование матрицы связи напряжений с деформациями COMM = ORT(L,I)/((l. + ORT(L,2)).(l.-ORT(L,2).2).AREA) ESTIFM(1,1) = С0ММ.(1. - 0RT(L,2)) ESTIFM(1,2)= C0MM.0RT(L,2) ESTIFM(1,3) = 0. ESTIFM(2,1) = ESTIFM(1,2) ESTIFM{2,2) = ESTIFM(1,1) ESTIFM(2,3) = 0. ESTIFM(3,l) = a ESTIFM(3,2) = 0.
ESTIFM(3,3) = ORT(L.l)/(2.(l-+ORT(L,2).AREA)
В - матрица напряжений для обратного хода; хранится иа магнитной ленте
DO 205 I = 1,3 DO 205 J= 1,6 B(1,J) = 0. DO 205K= 1,3 205 B(I,Jl = B{I,J)-(-ESTIFM(I,K)/2..A(K,J) WRITE(NT4)N,((B(1,J),J = I,6),I = 1,3)
ESTIFM -матрица жесткости
DO 210 1 = 1,6 DO 210 J = 1,6 ESTIFM(I,J) = 0.
DO 210 K= 1,3 210 ESTIFM(I,J) = EST1FM(I,J) -t- B(K,I)/2..A(K,J) RETURN
Выход из программы при ошибке в задании связей^ 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. Составление ансамбля и решение уравнений
В любой программе, реализующей метод конечных элементов, ключевой является подпрограмма решения систем уравнений. Выбор метода решения зависит от числа уравнений задачи и от типа используемой вычислительной машины.
В простейшем случае уравнения формируются полиостью. При таком подходе требуется Л' ячеек памяти, где Л^ -количество уравнений, и его можно использовать только для небольших задач и при применении больших ЭВМ. При 100 уравнениях, например, требуется 10 000 ячеек памяти.
Для решения систем линейных уравнений хорошо разработаны два метода: а) прямой метод, позволяющий получить точное (в пределах ошибки округления) решение; б) итерационный метод, использующий сходящийся к точному решению процесс последовательных приближений.
В качестве прямого метода часто используется метод исключения Гаусса, а в качестве итерационного - метод Гаусса - Зей-деля.
20.5.1. Метод исключения Гаусса
Пусть система Л/ уравнений представлена а блочном виде
(20.1)
где АГц -матрица размерности 1 X 1,
K\i - матрица размерности 1 X (Л^ - 1). К21 - матрица размерности (Л/ - 1) X U К22 - матрица размерности (Л/ - 1) X (ЛГ - 1), 6 -вектор неизвестных величин, F - вектор известных величин. Метод исключения Гаусса позволяет уменьшить размерность матрицы К и получить матричное уравнение размерности - 1 в виде
Гд = Г, (20.2)
где
К* = - К2\К\1 К\2, F = р2 - KiiKwFi
(20.3) (20.4)
Разбивая точно так же на блоки матрицу К*, этот прием можно првторить. Основной операцией прн этом будет вычисление тройного произведения KiiKTiKn. Так как матрица Км имеет размерность 1 X 1, количество операций пропорционально (Л^- 1). Когда размерность матрицы К будет сведена к 1 X 1, последнюю неизвестную 8 можно будет найти непосредственно.
Используя обратный ход'), остальные неизвестные можно
) Процесс нсключення Гаусса обычно называют прямым ходом, а решение системы с треугольной матрицей- обратным ходом. -Ярнл. ред.
определить из уравнения типа
Si = KiiFi - КчКА.
(20.5)
Описанная процедура представляет собой простейший прямой метод, при использовании которого все элементы матрицы хранятся в оперативной памяти и участвуют в вычислениях. Используя свойство симметрии матриц жесткости, можно сократить число операций почти вдвое, производя все действия только с частью матрицы, расположенной выше диагонали. При этом для решения необходимо выполнить около /еЛ операций.
Типичная матрица жесткости ансамбля элементов содержит много нулей, в частности, на некотором расстоянии от диагонали находятся только нулевые члены. Такая матрица называется ленточной, а расстояние от диагонального члена элемента до последнего ненулевого элемента этой же строки называется половиной ширины ленты. Ленточный характер этой матрицы можно продемонстрировать, записывая ее в блочном виде, как показано ниже, где нулевая подматрица заменяет часть подматрицы Кц, содержащую нули:
Кп Ki2 О
Ки О
К22 К23
К23 Кзг
(20.6)
В силу симметрии матрицы будут иметь размерности
К (1X1), КАХт), KiiimXm):-
К23 (m X (ЛГ - 1 - /п)), KaMN - 1 - m) X (iV - 1 - т)).
При исключении лодматрицы Ки изменяется только К22, так как нулевая подматрица не изменяет Кгз и Кзз- Количество операций исключения пропорционально т', а общее количество операций - величине
или приблизительно ЧгМт, где /Птах - максимальная величина половины ширины ленты.
На практике половина ширины ленты обычно меньше /ю размерности матрицы, и описанный способ использования свойства ленточности матрицы позволяет уменьшить число арифметических операций почти до 3% от числа выполняемых операций прн использовании метода, не учитывающего ленточного характера матрицы.
Другим преимуществом является экономное использование памяти машины, так как матрица может быть помещена в прямоугольный массив размерности X max, как показано иа фиг. 20.1. Однако при этом общее число уравнений, которое может быть решено, ограничивается размерами прямоугольного массива.
55
Kll К12
К22 К23
Кзз Кз4
/Cg5 О
О 2i О
О
о
Фиг. 20.1. Хранение матрицы при решении уравнений методом ленточных
матриц.
Ленточная матрица, задаваемая формулой (20.6), особенно удобна при решении больших систем уравнений, так как в процессе исключения операции производятся только над К22 и поэтому только эта матрица нужна в оперативной памяти. Однако истинная эффективность достигается, если составление уравнений для ансамбля и исключение производятся параллельно, т. е. если уравнение исключается сразу после его составления. Это легко осуществить, если в качестве первого узла каждого элемента всегда использовать узел с наименьшим номером, а элементы располагать так, чтобы номера их первых узлов располагались последовательно. При таком способе процесс составления уравнений ансамбля будет автоматически прекращаться, как только номер первого узла станет больше номера рассматриваемого узла. После исключения Кц оставшаяся матрица сдвигается так, что первый элемеит измененной матрицы К22 занимает позицию (1,1), а остальные - соответствующие им места. Исключенные уравнения могут временно пересылаться в промежуточное запоминающее устройство, а затем переписываться в виде блока в накопитель до следующего вызова. Для небольших задач объем промежуточной памяти может оказаться достаточным для того, чтобы вместить в себя все исключенные уравнения, при этом делать пересылку нет никакой необходимости.
Новая матрица может быть вновь выражена через подматрицы Км, Ki2 и К22, после чего процесс исключения повторяется.
На каждом этапе исключения требуется (m+l)X(m + 2) .ячеек памяти или при наличии симметрии - половина этого
числа. С помощью этого метода можно решить практически неограниченное число уравнений при условии конечности максимальной ширины ленты. В большинстве случаев при решении за-дачиЧ1етодом конечных элементов точность ие является проблемой. При необходимости можно решать уравнения с двойной точностью, однако при этом увеличиваются используемый объем памяти и время решения. В некоторых случаях точность можно повысить, если подставить полученные значения в первоначальную систему уравнений и вычислить невязки, а затем, используя эти невязки в качестве новых правых частей, снова решить систему. Сумма этих двух решений даст уточненное решение. Однако чаще эти невязки используются для оценки ошибок округления.
20.5.2. Экономное распределение памяти для решения ленточных систем
Чтобы уменьшить используемый объем памяти вычислительной машины, в программе решения больших систем FESS (в настоящей главе ие приводится) для храиерия в памяти верхней
X-элементы, не хранящиеся вптттц
Фнг. 20.2. Порядок размещения в памяти уравнений в программе FESS.
треугольной части матрицы жесткости, необходимой при исключении самого верхнего уравнения, используется одномерный массив. Пример расположения в памяти элементов матрицы показан иа фиг. 20.2. В этом примере элементы матрицы располагаются по столбцам, так как это удобнее, чем построчное расположение. В программе FESS матрицы, необходимые для обратного хода [уравнение (20.5)], хранятся в верхней части одномерного массива. Если матрица жесткости и эти матрицы перекрывают друг друга, содержимое остатка переписывается иа ленту и процесс повторяется. Такой способ достаточно эффективен, так как позволяет использовать длинные записи иа ленту. В некоторых случаях удается решить уравнения без использования виешией памяти.
20.5.3. Итерационный метод Гаусса - Зейделя
В общем случае п-е уравнение системы N уравнений может быть записано в виде
а-1 N
1=\ i=n+l
Из этого уравнения можно найти
(20.7)
(20.8)
Если процесс итераций таков, что в правой части используются последние приближения б,-, то для те-й итерации имеем
:i==k-AF,-tknibT- т jn,6r-.
I 1 = 1 i=n-l )
(20.9)
В этом заключается итерационный процесс Гаусса - Зейделя.
Часто для уточнения решения используется прием, состоящий в умножении разности между двумя итерациями для б на некоторый коэффициент и представлении уточненной величины б в виде
б7 = бГ+р(бГ-бГ).
(20.10)
где б™* - величина, вычисленная в соответствии с (20.9), а р - коэффициент верхней релаксации, значение которого обычно лежит между 1 и 2. Установлено, что во многих практических случаях самым подходящим является значение, близкое к 1,8.
Итерационный метод Гаусса - Зейделя легко программируется. Матрица жесткости хранится в компактной форме без нулевых членов вместе с матрицей-указателем номеров столбцов, в которых находятся ее элементы. Каждое уравнение итерируется в соответствии с (20.9), и найденное значение уточняется в соответствии с (20.10). Процесс повторяется столько раз, сколько необходимо для получения приемлемого решения, причем сходимость обычно оценивается путем вычисления разности между двумя последовательными приближениями.
Итерационные методы удобны для решения нелинейных задач, так как при их использовании обычно требуется решение ряда сходных задач, поэтому для вектора решения всегда есть удачное начальное приближение. Кроме того, если на промежуточных этапах ограничиться невысокой точностью решения, то число арифметических операций уменьшается.
Недостатком использования итерационных методов является необходимость повторения основного цикла для всех уравнений.
Поэтому при использовании для хранения уравнений внешней памяти процесс решения занимает много времени. Однако компактное расположение матрицы жесткости дает возможность решать большое число уравнений с использованием только оперативной памяти, например 1000 дополнительных уравнений на 32К слов памяти. Главный же недостаток итерационных методов состоит в том, что при решении, как правило, может быть рассмотрен только один вектор нагрузки, так как перемещения занимают почти всю память.
В приложении 20А приведены две подпрограммы: FORMK (формирование матрицы жесткости и соответствующей матрицы-указателя) и SOLVE (решение системы уравнений итерационным методом Гаусса - Зейделя). Это упрощенные подпрограммы, в которых заданы число циклов и пределы сходимости. Кроме того, в них используется взаимно однозначное соответствие между матрицей жесткости и матрицей-указателем. Эти подпрограммы совместимы с остальными подпрограммами, приведенными в этой главе. Заметим, что подпрограмма FORMK вычисляет члены, лежащие как выше, так и ниже диагонали компактной матрицы жесткости, и этим она отличается от другой подпрограммы FORMK (программа 20-5), которая используется при применении прямого метода решения.
20.5.4. Некоторые другие прямые методы
Видоизменением метода исключения, используемого в случае ленточных матриц, является так называемый метод редкозапол-ненных матриц, часто применяемый в экономических расчетах. В соответствии с этим методом подматрица записывается в компактной форме вместе с матрицей-указателем, назначение которой состоит в указании номера строки и столбца полной матрицы. При решении задач методом конечных элементов с одинаковым числом степеней свободы в каждом узле матрица-указатель может быть записана в компактной форме для того, чтобы соответствовать узловым подматрицам.
В методе редкозаполненных матриц исключение производится точно так же,-как описано выше, но при этом в качестве индекса используется матрица-указатель. Преимущество, получаемое за счет избавления от некоторых операций, необходимо сравнить с дополнительными затратами времени на отыскание элемента в полной матрице с помощью матрицы-указателя. Требуемый объем памяти зависит от числа ненулевых подматриц на каждом шаге процесса исключения и не зависит от ширины ленты. Следует отметить, что при исключении на местах, где раньше стояли нули, появляются ненулевые подматрицы. Для максимально эффективного использования памяти и времени
уравнения следует упорядочить так, чтобы в матрице в каждый момент содержалось минимальное число элементов. Определение оптимального порядка является сложной задачей динамического программирования, и для ее решения могут использоваться другие приближенные методы [5]. С помощью аналогичного упорядочения уравнений можно минимизировать число операций в ленточных методах. Однако при этом ширина ленты может сильно увеличиться, что сведет на нет все преимущества метода.
В приложении 20Б приведены подпрограмма FORMK составления матрицы жесткости (верхней треугольной части в форме прямоугольника) вместе с соответствующей матрицей-указателем и подпрограмма SOLVE решения систем уравнений методом редкозаполненных матриц.
Более существенными отличиями обладает фронтальный метод решения. Как и в методе редкозаполненных матриц, матрица жесткости записывается в компактной форме вместе с матрицей-указателем. Однако здесь исключение производится по элементам независимо от нумерации узлов. Своим названием метод обязан тому, что в процессе решения создается некоторый фронт, проходящий через систему узлов. Метод применяется главным образом при наличии дополнительных узлов на сторонах элементов, так как он может привести к значительной экономии памяти по сравнению с ленточным методом.-
Преимущество использующей метод фронтального исключения системы FINESSE особенно заметно при использовании элементов высоких порядков и в том случае, когда основное внимание уделяется характеристикам элементов. .Поскольку, например, предпочтительна поэлементная выдача, число расчетных параметров нагружения, которое может быть помещено в оперативной памяти вычислительной машины, значительно увеличивается. В каждый момент времени в процессе исключения в оперативной памяти машины находятся только сведения о векторе нагрузки.
Фронтальный метод решения наиболее эффективен при решении больших задач с применением трехмерных элементов.
20.5.5. Некоторые специальные приемы
Для улучшения более крупных программ могут использоваться некоторые специальные приемы:
а) Выделение подконструкций. Вместо того чтобы решать полную систему уравнений, иногда бывает удобно строить -новую матрицу меньшей размерности. Этот прием позволяет построить матрицу жесткости для совокупности элементов, в которой все точки, кроме внешних, исключены.
б) Включение дополнительной жесткости. В дополнение к выделению подконструкций используется введение в полную матрицу жесткости специальных дополнительных жесткостей. Это позволяет непосредственно учитывать, например, граничные условия, описывающие упругое опирание.
Эти два приема особенно полезны при решении нелинейных задач, в которых большая часть конструкции остается линейно-упругой. В таких случаях линейно-упругая часть конструкции заменяется некоторой подконструкцией в виде линии или поверхности, примыкающей к исследуемой области нелинейности, а затем решается нелинейная задача при упругом, граничном условии. Получающаяся нелинейная задача значительно меньше по объему.
20.5.6. Учет граничных условий
При применении метода конечных элементов граничные условия в напряжениях автоматически учитываются вектором нагрузки; если граница свободная, элементы вектора нагрузки, соответствующие расположенным на границе узлам, равны нулю.
При заданных перемещениях (например, осадка фундамента и др.) матрица жесткости и вектор нагрузки должны быть модифицированы. В общем случае для упрощения- индексации размерность матрицы не должна изменяться, т. е. никакие строки или столбцы не должны исключаться. Существуют два способа, с помощью которых, не нарушая указанных требований, можно удовлетворить граничным условиям, не учтенным при формировании элемента. Положим, что имеется некоторая система N уравнений
12 22
п2 пл nN
(20.П)
и, скажем, щ = а.
В соответствии с первым способом столбец нагрузки видоизменяется так, что F{ = Fi - kna {i = 2, N) и = a. Тогда соответствующая строка и столбец становятся нулевыми, а диагональный член - единичным. В частном, но часто встречающемся случае, когда а = О (т.е. опора неподвижна), необходимо изме-
НИТЬ матрицу описанным выше способом, оставляя матрицу нагрузки неизменной, кроме члена Fi= 0.
Второй способ состоит в умножении соответствующего диагонального элемента матрицы на некоторое большое число, скажем 10*, перед модификацией соответствующего коэффициента нагрузки. В рассматриваемом случае мы бы получили
fe = fell-108, f = fe -108-a,
kij - kij (за исключением случая, когда г = 1, / = 1),
В полученном решении i будет почти равно а. Этот способ пригоден для любых методов решения.
Во многих случаях возникает необходимость определять силы реакции в граничных точках, например реакции мостовых опор или контрфорсов плотин.
Если требуется определить реакции, то при обратном ходе необходимо преобразовать соответствующие уравнения так, чтобы вместо определения в точке опоры неизвестного перемещения по известной нагрузке отыскивалась неизвестная нагрузка по известному перемещению. В соответствии с первым из уравнений (20.11) получаем следующее выражение для силы реакции:
20.5.7. Пример подпрограммы
Ниже приведены блок-схемы и тексты двух подпрограмм. Подпрограмма FORMK используется для построения ленточной прямоугольной матрицы жесткости и учета граничных условий первым способом, описанным в предыдущем разделе (при а=0). Подпрограмма SOLVE применяется для решения систем алгебраических уравнений методом ленточных матриц ) Блок-схемы приведены на стр. 488 и 489.
Обозначения переменных в подпрограмме FORMK
NBAND Максимально возможная в про-
грамме ширина ленты
NROWB, NCOLB, NCOL Переменные, определяющие поло-
жение элемента матрицы жесткости
) В приложении 20В приведены другие подпрограммы, несовместные с си-, стемой FFSS.
NR, NX
ESTIFM* (12,12)
Переменные, используемые для
записи граничного условия Матрица жесткости элемента
Обозначения переменных в подпрограмме SOLVE NBAND Максимально возможная в про
С
R1 * (200)
грамме ширина ленты Счетчик числа уравнений для исключения и обратного хода Рабочая переменная для процесса
исключения Вектор правых частей; в конце работы программы на его место помещается решение
Программа 20-5 SUBROUTINE FORMK
С формирование верхнего треугольника матрицы жесткости
С
С
C0MM0N/C0NTR/TITLE(12),NP,NE,NB,NDF,NCN,NLD,NMAT, NSZF,LI,NT4
COMMON CORD(100,2),NOP(200,4),IMAT(200),ORT(25,2),NBC(25),
NFIX(25) 1,R1(200),SK(200,40) 2,ESTIFM(12,12)
С
С Ввод максимальной ширины ленты и информации о количестве
С решаемых уравнений
С
NBAND = 40
С
С Задание нулевой матрицы жесткости С
DO 300 N = l.NSZF DO 300 М = 1,NBAND 300 SK(N,M) = 0.
С Цикл по элементам С
DO 400 N = 1,NE CALL STIFT2(N)
С
С Возврат к ESTIFM, как к матрице жесткости С
с с с с
с с с
Блок-схема подпрограммы FORMK
Задание нулевой
Начало цикла по алементам
Составление оюсткости | матрицы злемента |
| | |
Запоминание матрицы , жвсткост алемента в прямоугольной матрице |
Конец цикла по злементам
Учет ераничных условий
возврат в основную гроврамму
Засылка ESTIFM в массив SK
Обход по строкам
DO 350 JJ= l.NCN
NROWB = (NOPfN.JJ) - I ).NDF
DO 350 J= !,NDF
NROWB = NROWB + 1
I = (JJ- l).NDF +J
Затем no столбцам
DO 330 KK = I.NCN
NCOLB = (NOP(N,KK)- n.NDF
DO 320 К = I.KDF
L = (KK- I).NDF+K
NCOL = NCOLB + К + I - NROWB
Блок-схема подпрограммы SOLVE
Начало цикла по уравнениям (Ы)
Модиринация членов, находящихся внутри ленты уравнений (20
Модификация вектора наврузки
Конец цикла ло уравнениям
Начало цикла обралтоао хода по уравнениям
т
Решение уравнения (20.2) при обратном ходе
Конец цикла по уравнениям
I
Возврат в основную программу
С С
С
Если элемент ниже ленты, то пропуск засылки
IF(NCOL) 320,320,310 310 SK(NROWB,NCOL) = SK(NROWB,NCOLJ + ESTiFM(I,L) 320 CONTINUE 330 CONTINUE З.ЧО CONTINUE 400 CONTINUE
1 ...
21 22 23 [
24 ]
25 26 27
