Главная страница сайта  Российские промышленные издания (узловые агрегаты) 

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 [ 83 ] 84 85 86 87 88

57 X(I) = XAUX(1)/EIG DO 67 1 = 1,NV

IF(ABS((X(I) - XUX(I)/X(I)) - TEST)67,67,82 67 CONTINUE

Достаточно

GO TO 50

Повторение

82 ITS = ITS - 1 IF(ITS)21,21,25

21 WRITE(5,22)

22 FORMAT(26H ITERATION COUNT EXCEEDED) GO TO 50

25 DO 26 I = 1,NV

26 XUX(I) = X(I)

42 FORMAT (4E16.8)

GO TO 14 50 0M = SQRT(1./EIG)

WRITE(5,13) II.OM

WRITE(5,42) (X(IU = 1,4) 13 FORMAT (15,E16.8)

Формирование видоизменений матрицы

CALL TPRD(X,W,XUX,NV,1,NVJ CALL MPRD(XUXAXAUX,1,NV,1) AA = EIG/XAUX(1) DO 68 I = 1,NV 68 XAUX(I) = X(I).AA

CALL MPRD (XAUX,XUX,EAUX,NV,1,NV) DO 110 1 = 1,NV DO 110J=1,NV 110 EGG(I,J) = EGG(I,J) - EAUX(I,J) 1 CONTINUE RETURN END

Программа 20-10

subroutine mprd(d,B,db,l,m,n) dimension d(4,41,B(4,4),db(4,4) DB(l,n) = d(l,m).BiM,n)

x)0 110 J=1,n do 110 1= 1,l DB(i,J) = 0 do 110 k = 1,m

110 db;i,J) = db(i,J) + d(i,k).b(k,J)

return end

Программа 20-11

SUBROUTINE TPRD (D,B,DB,M,L,N) DIMENSION D(4,4),B(4,4),DB(4,4)

; DB(L,N) = (TpaHcnoHHpOBaHHafl D(M,L)).B(M,N)

DO 110 J = 1,N DO 110 I=-1,L DB(I,J)=0 DO 110 K = 1,M 110 DB(I,J) = DB(I.J) + D(K,I)B(K,J) RETURN END

20.10. Заключительные замечания

В этой главе рассмотрен весь процесс реализации метода конечных элементов и приведены образцы программ. Изложенный материал не содержит каких-либо утонченных приемов, которые могли бы оказаться непонятными для начинающих. Напротив, авторы привели достаточно простые, но весьма эффективные рабочие программы, отметив в то же вреявозможности использования других методов и учета особенностей исследуемой системы.

Программа решения задач о собственных значениях, приведенная в этой главе, очень проста и не использует свойство симметрии матриц жесткости и масс (или геометрической жесткости). Читателям, интересующимся применением более совершенных методов и возможностями экономии памяти машины, следует обратиться к работе Андерсона [8], посвященной задачам о колебаниях и устойчивости.

ЛИТЕРАТУРА

1. Clough R. W., The Finite Element in Plane Stress Analysis, Proc. 2nd ASCE Conf. on Electronic Computation, Pittsburgh, Pa., Sept. 1960.

2. Irons B. M., Economical Computer Techniques for Numerically Integrated Finite Elements, Ini. S. Num. Meth. Eng., 1, 201-203 (1969).

3. Martin H. C, Introduction to Matrix Methods of Structural Analysis, McGraw-Hill, 1966.

4. Clough R. W., Johnson C. p., A Finite Element Approximation for the Analysis of Thin Shells, Int. J. Solids Struct. 4, 43-60 (1968).

5. King I. P., An Automatic Recording Scheme for Simultaneous Equations Derived from Network Systems (будет опубликовано).

6. Irons В. М., A Frontal Solution Program for Finite Element Analysis, Ini. J. Num. Meth. Eng., 2, 5-32 (1970).

7. Zienkiewicz O. C, King I. P., Discussion on «The Analysis of a Four-Span Bridge Using an Electrical Analogue Computers, Proc. Inst. Civ. Eng., 37, 819-820 (1967).

8. Anderson R. G., A Finite Element Eigenvalue Solution System, Ph. D. Thesis, Univ. of Wales, Swansea, 1968. .. --

ПРИЛОЖЕНИЕ 20А

Приведенные программы являются подпрограммами формирования и решения уравнений итерационным методом Гаусса - Зейделя.

с с с

с с с

с с с

с с с

Программа 20-12

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,RH200),S1(200,20),ISP(200,20) 2,ESTIFM(12,12)

Задание максимального числа членов

NMAX = 20 NOFF = 20

Задание нулевых массивов

DO 300 N = l.NSZF

DO 280 М = l.NMAX 280 SI(N,M) = 0

DO 290 M = 2,N0FF 290 ISP(N,M) = 0 300 ISP(N,1) = N

Обход no всем элементам

DO 400 N = l.NH CALL STIFT2(N)

Возврат к ESTIFM, как к матрице жесткости

Засылка ESTIFM в массив SI н члена-указателя в ISP

С С С

Обход сначала по строкам 1 = 0

DO 350 JJ= l.NCN

NROWB = (NOP(N,JJ) - I).NDF

DO 350 J=I,NDF

305 310

с с с

с с с

NROWB = NROWB -f- 1

I = I-f 1

Затем обход по столбцам 11=0

DO 330 КК = l.NCN

NCOLB = (NOP(N,KK) - 1)»NDF

DO 330 К = l.NDF

NCOLB = NCOLB -b 1

II = II-f 1

Поиск в ISP номера сто-тбца DO 310 M = 1,N0FF

IF(ISP(NROWB,M) - NCOLB) 305,320,305 IF(ISP(NROWB,M)) 315,315,310 CONTINUE

Поиск свободного места для хранения NCOLB ISP(NROWB,M) = NCOLB

Запись ESTIFM

SI(NROWB,M) = ESTIFM\I,II) -j- SI(NROWB,M)

Конец цикла по столбцам CONTINUE

Конец цикла по строкам CONTINUE

Коаец цикла по элементам CONTINUE

Учет граничных условий

DO 500N= 1,NB NX= 10..(NDF- 1) I = NBC(N)

NROWB = (I- 1).NDF

Проверка каждой степени свободы

DO490 М = I.NDF NROWB = NROWB -f 1 ICON = NFIX(N)/NX IF(ICON) 450,450,420

С С С

Запоминаине большого номера на диагонали

420 SI(NROWB,1) = SI(NROWB,1)>10..20 NFIX(N) = NFIX(N) - NX.ICON

450 NX = NX/10

490 CONTINUE

500 CONTINUE RETURN END

Программа 20-13

SUBROUTINE SOLVE

С Подпрограмма предиазначеиа для решения уравнений ятерацнон-С иым методом

C0MM0N/C0NTR/TITLE(12),NP,NE,NB,NDF,NCN,NLD,NMAT,

COMMONa3RD(100,2),NOP(200,4),IMAT(200),ORT(25,2),NBC(25), NFIX(25)

1,R(200),A(200.20),ITEM(200,20),DIS(200) NT = 20

Задание коэффициента релаксации и точности

T0LER = .lE-3 RELAX =1.8

Прн отрицательном числе уравнений пропуск задания начальных данных

IF(NEQ.LT.0)G0 ТО 310

DO 300 N = l.NEQ DO 250 М= 1,NT IF(ITEM(N,M).NE.O) GO TO 250

Массив 1TEM(N,1) содержит счетчик расстояния от диагонали

с с с с

с с с с

с с с

ITEM(N,1) = M- 1

GO ТО 260 250 CONTINUE 260 CONTINUE 300 CONTINUE 310 NEQ=5=IABS(NEQ)

Задание максимального количества циклов

NCYC = NEQ/2 IF(NCYC.LT.25)NCYC = 25

С С С С

С С С С

I Задание нулевого массива неизвестных

DO 320 N = l.NEQ IF( А( N, 1 ).NE.O.)A( N,1) = 1 ./A(N, 1) 320 DIS(N) = 0.

I Начало цикла по циклам

DO 500 NC = l.NCYC SUM = 0. SUMD = 0.

I Начало цикла по уравнениям

DO 450N=l,NEQ FX = R(N) NUM = ITE.M(N,1) DO 330 M = 2,NUM L = ITEM(N,M) 330 FX = FX-A(N,M).DIS(L)

FX - общее отклонение от RHS DX - измененное значение

DX = A(N,1).FX-D1S(N) DIS(N) = DlS(N) + RLAX.DX

SU.Vl и SUMD - параметры, характеризующие сходимость

С С С

с с с

SUM = SUM + ABS(DX) SUMD = SUMD + ABS(DIS(N)) 450 CONTINUE

Выход аз цикла арн достижении сходимости ND = NC

IF(SUM.LT.SUMD.TOLER) GO ТО 550 500 CONTINUE

Пересылка окончательных результатов в массив R

550 DO 600 N = l.NEQ 600 R(N) = DIS(N)

Печать последнего значения суммы а т. д.

WRITE(6,10),ND,SUM,SUMD FORMAT(20H0 LAST CYCLE NO. = ,110 1, / 20H (SN-SN-1)/SN =,EI0.3

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 [ 83 ] 84 85 86 87 88