📄 先处理法程序pff.f90.txt
字号:
!---------------------------------------------------------------------------------------
! STATIC ANALYSIS FOR PLANE FRAME
! METHOD OF FORMER TREATMENT
! 平面杆系结构静力分析程序
! 先处理法源程序
!
! 主要功能和特点：
! 可调数组
! 输入单元结点编号和结点位移编码；
! 先处理法形成总体刚度阵和结构综合结点荷载
! 存储总刚上三角元素；
! GAUSS消元法解线性代数方程组；
! 计算结点位移和单元杆端力
!
!------------------------------------------------------------------------------------
! Define the interface of assumed-shape arrays
! 定义可调数组接口模块
MODULE KECALL
INTERFACE
SUBROUTINE ESM(M,NE,E,A,ZI,BL,SI,CO,EK)
REAL::E,BL
INTEGER::M,NE
REAL::A(NE),ZI(NE),EK(:,:)
END SUBROUTINE
SUBROUTINE CTM(SI,CO,T)
REAL::SI,CO
REAL::T(:,:)
END SUBROUTINE
SUBROUTINE TTKT(EK,T)
REAL::EK(:,:),T(:,:)
END SUBROUTINE
END INTERFACE
END MODULE
! Main program reads the control parameters and organizes
! the whole calculation by calling subroutines.
! 主程序
! 输入控制参数，实现子程序调用和组集
PROGRAM PFF
implicit none
INTEGER::iAl_Err=0 ! 错误信息控制变量。若iAl_Err\=0，则输出出错信息。
REAL::E
INTEGER::NE,NJ,NR,NP,NF,N
REAL,ALLOCATABLE::X(:),Y(:),A(:),ZI(:),PJ(:,:),PF(:,:),TK(:,:),P(:) ! 定义可调数组
INTEGER,ALLOCATABLE::IJ(:,:),JN(:,:)
CHARACTER*12 INDAT,OUTDAT
WRITE(*,*) 'Please input primary data file name!'
READ(*,'(A12)') INDAT
WRITE(*,*) 'Please input calculation result file name!'
READ(*,'(A12)') OUTDAT
OPEN(1,FILE=INDAT,STATUS='OLD')
OPEN(2,FILE=OUTDAT,STATUS='NEW')
READ(1,*) NE,NJ,N,NP,NF,E ! 输入结构控制参数
! NE: 单元总数，整型
! NJ: 结点总数，整型
! N: 结点位移总数，整型
! NP: 结点荷载总数，整型
! NF: 非结点荷载总数，整型
! E: 弹性模量，实型
ALLOCATE(X(NJ),Y(NJ),A(NE),ZI(NE),PJ(NP,3),PF(NF,4), &
TK(N,N),P(N),IJ(NE,2),JN(NJ,3),STAT=iAl_Err)
! 为可调数组分配地址
IF (iAl_Err/=0) THEN ! 输出出错信息
WRITE(*,*) 'Inadequate heap RAM!'
WRITE(*,*) 'Program terminated!'
CLOSE(1)
CLOSE(2)
STOP
ENDIF
! 在文件中输出结构控制参数
WRITE(2,10) NE,NJ,N,NP,NF,E
10 FORMAT(3X, 'PLANE FRAME STRUCTURE ANALYSIS', /, 5X, 'NE=', I2, &
8X, 'NJ=', I2, 8X, 'N=', I2, /, 5X, 'NP=', I2, 8X, 'NF=', I2, &
8X, 'E=', E12.4)
CALL INPUT(NE,NJ,N,NP,NF,X,Y,IJ,A,ZI,JN,PJ,PF) ! 输入、输出原始数据
CALL TSM(NE,NJ,N,E,X,Y,IJ,A,ZI,JN,TK) ! 形成结构总刚
CALL JLP(NE,NJ,N,NP,NF,X,Y,IJ,JN,PJ,PF,P) ! 计算结构综合结点荷载列阵
CALL GAUSS(TK,P,N) ! GAUSS消元法计算结点位移
CALL MVN(NE,NJ,N,NF,E,X,Y,IJ,A,ZI,JN,PF,P) ! 计算并输出结点位移和单元杆端力
CLOSE(1)
CLOSE(2)
END
! Read and print primary data.
! 在数据文件中输入输出原始数据
SUBROUTINE INPUT(NE,NJ,N,NP,NF,X,Y,IJ,A,ZI,JN,PJ,PF)
implicit none
INTEGER::NE,NJ,N,NP,NF ! 带入哑元，整型变量
REAL::X(NJ),Y(NJ),A(NE),ZI(NE),PJ(NP,3),PF(NF,4) ! 带出哑元，实型数组
INTEGER::IJ(NE,2),JN(NJ,3) ! 带出哑元，整型数组
INTEGER::I,J ! 中间变量，整型
READ(1,*) (X(I),Y(I),I=1,NJ) ! 输入结点坐标
READ(1,*) (IJ(I,1),IJ(I,2),A(I),ZI(I),I=1,NE)
! 输入单元始、末端结点编码及单元面积和惯性矩
READ(1,*) ((JN(I,J),J=1,3),I=1,NJ)
! 输入结点位移分量编码，按总体编码的顺序输入，一行对应一个结点
IF(NP.GT.0) READ(1,*) ((PJ(I,J),J=1,3),I=1,NP) !输入结点荷载，
! 其中PJ(I，1)、PJ(I，2)、PJ(I，3)分别为第I个结点荷载作用的结点码、
! 作用方向代码及荷载数值。
! 作用方向代码沿整体坐标X方向为“1”，整体坐标Y方向为“2”，力偶为“3”。
! 荷载数值输代数值，即结点力与整体坐标正向一致为正，结点力偶以顺时针为正。
IF(NF.GT.0) READ(1,*) ((PF(I,J),J=1,4),I=1,NF) !输入非结点荷载
! 其中PF(I，1) : 第I个非结点荷载作用的单元号；
! PF(I，2): 荷载类型号
! PF(I，3): 荷载参数q, 以与单元局部坐标正向一致为正，力偶以顺时针为正
! PF(I，4): 位置参数c(以上参见表2-1)
! 输出结构原始参数，以便检查数据
WRITE(2,10) (I,X(I),Y(I),I=1,NJ)
WRITE(2,20) (I,IJ(I,1),IJ(I,2),A(I),ZI(I),I=1,NE)
WRITE(2,30) (I,(JN(I,J),J=1,3),I=1,NJ)
IF(NP.GT.0) WRITE(2,40) ((PJ(I,J),J=1,3),I=1,NP)
IF(NF.GT.0) WRITE(2,50) ((PF(I,J),J=1,4),I=1,NF)
10 FORMAT(/2X,'COORDINATES OF JOINT'/6X,'JOINT',12X,'X',&
12X,'Y'/(10X,I4,5X,2F12.4))
20 FORMAT(/2X,'INFORMATION OF ELEMENTS'/6X,'ELEMENT',&
4X,'JOINT-I',4X,'JOINT-J',10X,'A',14X,'ZI'/(2X,3I16,8X,2F12.6))
30 FORMAT(/2X,'INFORMATION OF JOINT DEGREES OF FREEDOM'/&
6X,'JOINT',6X,'u',7X,'v',6X,'ceta'/(6X,4I8))
40 FORMAT(/2X,'JOINT LOAD'/6X,'JOINT'8X,'XYM',12X,'LOAD'/ &
(6X,F5.0,10X,F5.0,10X,F12.4))
50 FORMAT(/2X,'NON-JOINT LOAD'/6X,'ELEMENT',8X,'TYPE',8X,&
'LOAD',12X,'C'/(6X,F6.0,18X,F6.0,6X,2F12.4))
END
! Calculate length,sine and cosine of member.
! 计算单元常数
SUBROUTINE LSC(M,NE,NJ,X,Y,IJ,BL,SI,CO)
implicit none
INTEGER::M,NE,NJ ! 带入哑元，整型变量
INTEGER::IJ(NE,2) ! 带入哑元，整型数组
REAL::X(NJ),Y(NJ) ! 带入哑元，实型数组
REAL::BL,SI,CO,DX,DY ! 代出哑元，实型变量
INTEGER::I,J ! 中间变量，整型
I=IJ(M,1) ! 单元始端结点编号
J=IJ(M,2) ! 单元末端结点编号
DX=X(J)-X(I)
DY=Y(J)-Y(I)
BL=SQRT(DX*DX+DY*DY) ! 单元长度
SI=DY/BL
CO=DX/BL
END
! Set up element location vector.
! 确定单元定位向量
SUBROUTINE ELV(M,NE,NJ,IJ,JN,LV)
implicit none
INTEGER::M,NE,NJ ! 带入哑元，整型变量
INTEGER::IJ(NE,2),JN(NJ,3) ! 带入哑元，整型数组
INTEGER::LV(6) ! 代出哑元，整型数组
INTEGER::I,J,K ! 中间变量，整型
I=IJ(M,1)
J=IJ(M,2)
DO K=1,3
LV(K)=JN(I,K)
LV(K+3)=JN(J,K)
ENDDO
END
! Calculate element stiffness matrix referred to global coordinate system.
! 生成整体坐标下的单刚阵
SUBROUTINE ESM(M,NE,E,A,ZI,BL,SI,CO,EK)
implicit none
INTEGER::M,NE ! 带入哑元，整型变量
REAL::E,A(NE),ZI(NE),BL,SI,CO ! 带入哑元，实型数组和实型变量
REAL::EK(:,:) ! 代出哑元，实型数组
REAL::C1,C2,C3,C4,S1,S2,S3,S4,S5,S6 ! 中间变量，实型
INTEGER::I,J ! 中间变量，整型
C1=E*A(M)/BL
C2=2.0*E*ZI(M)/BL
C3=3.0*C2/BL
C4=2.0*C3/BL
S1=C1*CO*CO+C4*SI*SI
S2=(C1-C4)*SI*CO
S3=C3*SI
S4=C1*SI*SI+C4*CO*CO
S5=C3*CO
S6=C2
EK=0.0
EK(1,1)=S1
EK(1,2)=S2
EK(1,3)=-S3
EK(1,4)=-S1
EK(1,5)=-S2
EK(1,6)=-S3
EK(2,2)=S4
EK(2,3)=S5
EK(2,4)=-S2
EK(2,5)=-S4
EK(2,6)=S5
EK(3,3)=2.0*S6
EK(3,4)=S3
EK(3,5)=-S5
EK(3,6)=S6
EK(4,4)=S1
EK(4,5)=S2
EK(4,6)=S3
EK(5,5)=S4
EK(5,6)=-S5
EK(6,6)=2.0*S6
FORALL(I=1:5,J=2:6,(J>I.AND.EK(I,J)/=0.0)) EK(J,I)=EK(I,J)
END
! Assemble total stiffness matrix.
! 组集总刚
! 按单元定位向量对号入座顺序存储单刚元素在总刚TK中
SUBROUTINE TSM(NE,NJ,N,E,X,Y,IJ,A,ZI,JN,TK)
USE KECALL ! 该处应用了可调数组，需调用模块KECALL
implicit none
INTEGER::NE,NJ,N,IJ(NE,2),JN(NJ,3) ! 带入哑元
REAL::E,X(NJ),Y(NJ),A(NE),ZI(NE) ! 带入哑元
REAL::TK(N,N) !代出哑元
REAL::BL,SI,CO,EK(6,6),T(6,6) ! 中间变量
INTEGER::I,J,K,L,M,LV(6) ! 中间变量
TK=0.0 !总刚赋初值
DO M=1,NE
CALL LSC(M,NE,NJ,X,Y,IJ,BL,SI,CO) ! 计算单元长度，方位角正弦和余弦
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -