4.txt
来自「文件1.txt,2.txt,3.txt和5.txt为用Fortran编写的有限元」· 文本 代码 · 共 569 行 · 第 1/2 页
TXT
569 行
PZ(1)=P(I0+1)
PZ(2)=P(I0+2)
PZ(3)=P(I0+3)
PZ(4)=P(J0+1)
PZ(5)=P(J0+2)
PZ(6)=P(J0+3)
CALL MATMU(AL,PZ,PB,6,6,1)
CALL FDK1(DK1,PB(2),PB(5),PB(3),PB(6),EL,EF(K))
CALL FDK2(DK2,EL,TN)
DO 210 I=1,6
DO 210 J=1,6
DK(I,J)=DK0(I,J)+DK1(I,J)+DK2(I,J)
210 CONTINUE
WRITE(10,102)K
102 FORMAT(1X,'Element ',I2,' Stiff Matrix(大位移)')
DO 101 I=1,6
WRITE(10,103)(DK1(I,J),J=1,6)
103 FORMAT(1X,6e12.3,';')
write(10,*)
101 CONTINUE
WRITE(10,402)K
402 FORMAT(1X,'Element ',I2,' Stiff Matrix(弹性+大位移)')
DO 401 I=1,6
WRITE(10,403)(DK0(I,J)+DK1(I,J),J=1,6)
403 FORMAT(1X,6e12.3,';')
write(10,*)
401 CONTINUE
WRITE(10,201)K
201 FORMAT(1X,'Element ',I2,' Stiff Matrix(几何)')
WRITE(10,204)TN
204 FORMAT(1X,'轴力=',e12.3)
DO 202 I=1,6
WRITE(10,203)(DK2(I,J),J=1,6)
203 FORMAT(1X,6e12.3,';')
write(10,*)
202 CONTINUE
WRITE(10,301)K
301 FORMAT(1X,'Element ',I2,' Stiff Matrix(切线刚阵)')
DO 302 I=1,6
WRITE(10,303)(DK(I,J),J=1,6)
303 FORMAT(1X,6e12.3,';')
write(10,*)
302 CONTINUE
RETURN
END
!-------------------------子程序FDK0-------------------------
!计算弹性刚度矩阵Ke
!
!
!
SUBROUTINE FDK0(DK0,A,TL,TI)
IMPLICIT INTEGER*2(I-N)
REAL DK0(6,6),A,TL,TI
COMMON/C1/E
CALL CLEAR(DK0,6,6)
DK0(1,1)=E*A/TL
DK0(1,4)=-DK0(1,1)
DK0(2,2)=12*E*TI/TL/TL/TL
DK0(2,3)=6*E*TI/TL/TL
DK0(2,5)=-DK0(2,2)
DK0(2,6)=DK0(2,3)
DK0(3,3)=4*E*TI/TL
DK0(3,5)=-DK0(2,3)
DK0(3,6)=2*E*TI/TL
DK0(4,4)=E*A/TL
DK0(5,5)=DK0(2,2)
DK0(5,6)=-DK0(2,3)
DK0(6,6)=DK0(3,3)
DO 100 I=1,5
DO 100 J=I+1,6
DK0(J,I)=DK0(I,J)
100 CONTINUE
! DO 101 I=1,6
! DO 101 J=1,6
! WRITE(10,103)I,J,DK0(I,J)
!103 FORMAT(1X,'DK0(',I3,',',I3,')=',E12.3)
!101 CONTINUE
RETURN
END
!-------------------------子程序FDK1-------------------------
!计算大位移刚度矩阵
!
!
SUBROUTINE FDK1(DK1,VI,VJ,THI,THJ,TL,TA)
IMPLICIT INTEGER*2(I-N)
REAL DK1(6,6),DK11(6,6),DK12(6,6),A,B,C,D, &
& AA,BB,CC,AB,AC,BC,TL,TA,VI,VJ,THI,THJ
COMMON/C1/E
D=VJ-VI
A=6*D/5/TL-(THI+THJ)/10
B=-D/10-2*TL*THI/15-TL*THJ/30
C=D/10+TL*THI/30-2*TL*THJ/15
AA=72.0*D*D/35/TL/TL/TL+ &
& 3*(THI*THI+THJ*THJ)/35/TL-18*D*(THI+THJ)/35/TL/TL
BB=3*D*D/35/TL+2*TL*THI*THI/35+ &
& TL*THJ*THJ/210+D*(THI-THJ)/70-THI*THJ/70
CC=3*D*D/35/TL+TL*THI*THI/210- &
& 2*TL*THJ*THJ/35-D*(THI-THJ)/70-TL*THI*THJ/70
AB=-9*D*D/35/TL/TL+6*THI*D/35/TL+(THI*THI-THJ*THJ)/140
AC=9*D*D/35/TL/TL-4*THJ*D/35/TL+ &
& (THI*THI+2*THI*THJ-THJ*THJ)/140
BC=D*(THI+THJ)/70+(THI*THI+THJ*THJ)/140-TL*THI*THJ/105
CALL CLEAR(DK1,6,6)
CALL CLEAR(DK11,6,6)
CALL CLEAR(DK12,6,6)
DK11(1,2)=A
DK11(1,3)=-B
DK11(1,5)=-A
DK11(1,6)=C
DK11(2,4)=-A
DK11(3,4)=B
DK11(4,5)=A
DK11(4,6)=-C
DO 110 I=1,5
DO 110 J=I+1,6
DK11(J,I)=DK11(I,J)
110 CONTINUE
DO 120 I=1,6
DO 120 J=1,6
DK11(I,J)=DK11(I,J)*E*TA/TL
120 CONTINUE
DK12(2,2)=AA
DK12(2,3)=-AB
DK12(2,5)=-AA
DK12(2,6)=AC
DK12(3,3)=BB
DK12(3,5)=AB
DK12(3,6)=-BC
DK12(5,5)=AA
DK12(5,6)=-AC
DK12(6,6)=CC
DO 130 I=1,5
DO 130 J=I+1,6
DK12(J,I)=DK12(I,J)
130 CONTINUE
DO 140 I=1,6
DO 140 J=1,6
DK12(I,J)=E*TA*DK12(I,J)
140 CONTINUE
DO 150 I=1,6
DO 150 J=1,6
DK1(I,J)=DK11(I,J)+DK12(I,J)
150 CONTINUE
RETURN
END
!-------------------------子程序FDK2-------------------------
!计算几何刚度矩阵
!
!
SUBROUTINE FDK2(DK2,TL,TN)
IMPLICIT INTEGER*2(I-N)
REAL DK2(6,6),TL,TN
CALL CLEAR(DK2,6,6)
DK2(2,2)=6*TN/5/TL
DK2(2,3)=TN/10
DK2(2,5)=-DK2(2,2)
DK2(2,6)=DK2(2,3)
DK2(3,3)=2*TN*TL/15
DK2(3,5)=-DK2(2,3)
DK2(3,6)=-TN*TL/30
DK2(5,5)=DK2(2,2)
DK2(5,6)=-DK2(2,3)
DK2(6,6)=DK2(3,3)
DO 160 I=1,5
DO 160 J=I+1,6
DK2(J,I)=DK2(I,J)
160 CONTINUE
RETURN
END
!-------------------------子程序GAS--------------------------
! 高斯法求解方程组
! N 方程组规模；A 方程组；B 输入为荷载，输出为位移
!
SUBROUTINE GAS(N,A,B)
IMPLICIT INTEGER*2(I-N)
DIMENSION A(150,150),B(150)
DO 310 K=1,N
T=A(K,K)
DO 320 J=K,N
WRITE(*,*)A(K,J),T
320 A(K,J)=A(K,J)/T
B(K)=B(K)/T
DO 330 I=K+1,N
T=A(I,K)
DO 340 J=K+1,N
340 A(I,J)=A(I,J)-A(K,J)*T
330 B(I)=B(I)-B(K)*T
310 CONTINUE
DO 350 I=N-1,1,-1
T=0.0
DO 360 J=I+1,N
360 T=T+A(I,J)*B(J)
B(I)=B(I)-T
350 CONTINUE
!RETURN
END
!-------------------------子程序CHL--------------------------
!单元角度正余弦
!
SUBROUTINE CHL(K,NJ,NE,NL,NR,X,Y)
IMPLICIT INTEGER*2(I-N)
DIMENSION NL(NE),NR(NE),X(NJ),Y(NJ)
COMMON/C2/I0,J0,SI,CO,EL
I=NL(K)
J=NR(K)
I0=3*(I-1)
J0=3*(J-1)
CO=X(J)-X(I)
SI=Y(J)-Y(I)
EL=SQRT(CO*CO+SI*SI)
CO=CO/EL
SI=SI/EL
RETURN
END
!-------------------------子程序MLA--------------------------
!计算方向矩阵 AL
!转置矩阵 AT
!
!
SUBROUTINE MLA
IMPLICIT INTEGER*2(I-N)
DIMENSION AL(6,6),AT(6,6)
COMMON/C2/I0,J0,SI,CO,EL
COMMON/C3/AL,AT
CALL CLEAR(AL,6,6)
AL(1,1)=CO
AL(1,2)=SI
AL(2,1)=-SI
AL(2,2)=CO
AL(3,3)=1.0
DO 400 I=1,3
DO 400 J=1,3
AL(I+3,J+3)=AL(I,J)
400 CONTINUE
DO 410 I=1,6
DO 410 J=1,6
AT(I,J)=AL(J,I)
410 CONTINUE
RETURN
END
!-------------------------子程序MATMU------------------------
SUBROUTINE MATMU(A,B,C,M,N,L)
!Multi Array
!A(M,N)*B(N,L) = C(M,L)!
DIMENSION A(M,N),B(N,L),C(M,L)
DO 420 I=1,M
DO 420 J=1,L
C(I,J)=0.0
DO 420 K=1,N
C(I,J)=C(I,J)+A(I,K)*B(K,J)
420 CONTINUE
RETURN
END
!-------------------------子程序CLEAR------------------------
SUBROUTINE CLEAR(A,M,N)
!
!Clear Array(M,N)!
DIMENSION A(M,N)
DO 430 I=1,M
DO 430 J=1,N
A(I,J)=0.0
430 CONTINUE
RETURN
END
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?