📄 后处理法程序pfl.f90.txt
字号:
!----------------------------------------------------------------------------------------
! STATIC ANALYSIS FOR PLANE FRAME
! (METHOD OF LATTER TREATMENT)
! 平面框架结构静力分析程序
! 后处理法
!
! 主要功能：
! 输入单元结点编号,自动生成结点位移编号；
! 总刚元素按上三角阵存储；
! 主1副0法引入位移约束条件
! GAUSS消元法解线性代数方程组；
!
!----------------------------------------------------------------------------------------
! Define the interface of assumed-shape arrays 定义可调数组接口模块
MODULE KECALL
INTERFACE
SUBROUTINE ESM(M,NE,E,A,ZI,BL,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 PFL
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(:,:),JR(:,:)
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,NR,NP,NF,E ! 输入结构控制参数
! NE: 单元总数，整型
! NJ: 结点总数，整型
! NR: 结点约束总数，整型
! NP: 结点荷载总数，整型
! NF: 非结点荷载总数，整型
! E: 弹性模量，实型
N=NJ*3 ! 结点位移总数＝3×结点总数，整型
ALLOCATE(X(NJ),Y(NJ),A(NE),ZI(NE),PJ(NP,3),PF(NF,4),&
& TK(N,N),P(N),IJ(NE,2),JR(NR,4),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,NR,NP,NF,E
10 FORMAT(3X,'PLANE FRAME STRUCTURE ANALYSIS'/5X,'NE=',&
& I2,8X,'NJ=',I2,8X,'NR=',I2,/5X,'NP=',I2,8X,'NF=',I2,8X,'E=',E12.4)
CALL INPUT(NE,NJ,NR,NP,NF,X,Y,IJ,A,ZI,JR,PJ,PF) ! 输入、输出原始数据.
CALL TSM(NE,NJ,E,X,Y,IJ,A,ZI,TK,N) ! 形成结构总刚
CALL JLP(NE,NJ,NP,NF,X,Y,IJ,PJ,PF,P,N) ! 计算结构综合结点荷载列阵
CALL ISC(NR,JR,TK,P,N) ! 引入支座约束条件
CALL GAUSS(TK,P,N) ! GAUSS消元法计算结点位移
CALL MVN(NE,NJ,NF,E,X,Y,IJ,A,ZI,PF,P,N) ! 计算并输出结点位移和单元杆端力
CLOSE(1)
CLOSE(2)
END
! Read and print primary data.
! 由数据文件输入、输出原始数据.
SUBROUTINE INPUT(NE,NJ,NR,NP,NF,X,Y,IJ,A,ZI,JR,PJ,PF)
implicit none
INTEGER::NE,NJ,NR,NP,NF
REAL::X(NJ),Y(NJ),A(NE),ZI(NE),PJ(NP,3),PF(NF,4)
INTEGER::IJ(NE,2),JR(NR,4)
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,*) ((JR(I,J),J=1,4),I=1,NR) ! 输入结点约束信息，
! JR(I,1): 支座结点码；
! JR(I,2): 支座结点沿整体坐标x轴方向约束信息，有＝1，无＝0；
! JR(I,1): 支座结点沿整体坐标y轴方向约束信息，有＝1，无＝0；
! JR(I,1): 支座结点沿转角方向约束信息，有＝1，无＝0；
IF(NP>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>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) ((JR(I,J),J=1,4),I=1,NR)
IF(NP>0) WRITE(2,40) ((PJ(I,J),J=1,3),I=1,NP)
IF(NF>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'/(6X,I4,5X,2F12.4))
20 FORMAT(/2X,'INFORMATION OF ELEMENTS'/6X,'ELEMENT',&
&4X,'JOINT-I',4X,'JOINT-J',10X,'A',12X,'ZI'/(2X,3I10,6X,2F12.6))
30 FORMAT(/2X,'INFORMATION OF RESTRICTION'/6X,'RES.-JOINT',&
&7X,'XR',8X,'YR',7X,'CETA'/(4X,4I10))
40 FORMAT(/2X,'JOINT LOAD'/6X,'JOINT'8X,'XYM',12X,'LOAD'/ &
&(6X,F5.0,6X,F5.0,6X,F12.4))
50 FORMAT(/2X,'NON-JOINT LOAD'/6X,'ELEMENT',8X,'TYPE',8X,&
&'LOAD',12X,'C'/(6X,F6.0,6X,F6.0,4X,2F12.4))
END
! Calculate length,sine and cosine of member.
! 计算单元长度、方位角的sin和cos值.
SUBROUTINE LSC(M,NE,NJ,X,Y,IJ,BL,SI,CO)
implicit none
REAL::BL,SI,CO
INTEGER::M,NE,NJ
REAL::X(NJ),Y(NJ)
INTEGER::IJ(NE,2)
REAL::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
! Calculate element stiffness matrix referred to element
! 生成局部坐标下的单元刚度矩阵
SUBROUTINE ESM(M,NE,E,A,ZI,BL,EK)
implicit none
REAL::E,BL
INTEGER::M,NE
REAL::A(NE),ZI(NE),EK(:,:)
REAL::EA,EI
INTEGER::I,J
EK=0.0
EA=E*A(M)
EI=E*ZI(M)
EK(1,1)=EA/BL
EK(1,4)=-EK(1,1)
EK(2,2)=12.0*EI/BL/BL/BL
EK(2,3)=6.0*EI/BL/BL
EK(2,5)=-EK(2,2)
EK(2,6)=EK(2,3)
EK(3,3)=4.0*EI/BL
EK(3,5)=-EK(2,3)
EK(3,6)=2.0*EI/BL
EK(4,4)=EK(1,1)
EK(5,5)=EK(2,2)
EK(5,6)=-EK(2,3)
EK(6,6)=EK(3,3)
FORALL(I=1:5,J=2:6,((J>I).AND.(EK(I,J)/=0.0))) EK(J,I)=EK(I,J) ! 给下三角元素赋值
END
! Form coordinate transformation matrix.
! 形成单元坐标变换阵.
SUBROUTINE CTM(SI,CO,T)
implicit none
REAL::SI,CO
REAL::T(:,:)
INTEGER::I,J
T=0.0
T(1,1)=CO
T(1,2)=SI
T(2,1)=-SI
T(2,2)=CO
T(3,3)=1.0
FORALL(I=1:3,J=1:3,T(I,J)/=0.0) T(I+3,J+3)=T(I,J)
END
! Calculate element stiffness matrix referred to global coordinate system.
! 生成整体坐标下的单元刚度阵
SUBROUTINE TTKT(EK,T)
implicit none
REAL::EK(:,:),T(:,:)
REAL,ALLOCATABLE::TE(:,:)
INTEGER::iAl_Err=0
ALLOCATE(TE(UBOUND(EK,1),UBOUND(EK,2)),STAT=iAl_Err)
IF (iAl_Err/=0) THEN
WRITE(*,*) 'Inadequate heap RAM!'
WRITE(*,*) 'Program terminated!'
CLOSE(1)
CLOSE(2)
STOP
ENDIF
TE=0.0
TE=MATMUL(TRANSPOSE(T),EK) ! 调用FORTRAN库函数进行矩阵相乘
EK=0.0
EK=MATMUL(TE,T)
END
! Assemble total stiffness matrix.
! 组集总刚
SUBROUTINE TSM(NE,NJ,E,X,Y,IJ,A,ZI,TK,N)
USE KECALL ! 该处应用了可调数组，需调用模块KECALL
implicit none
REAL::E
INTEGER::NE,NJ,N
REAL::X(NJ),Y(NJ),A(NE),ZI(NE),TK(N,N)
INTEGER::IJ(NE,2)
REAL::BL,SI,CO
INTEGER::I,J,K,L,M
REAL::EK(6,6),T(6,6)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -