📄 后处理法程序pfl.f90.txt
字号:
INTEGER::LV(6)
TK=0.0 ! 总刚阵置初值
DO M=1,NE ! 对单元循环
CALL LSC(M,NE,NJ,X,Y,IJ,BL,SI,CO) ! 计算单元长度，方位角正弦和余弦
CALL ESM(M,NE,E,A,ZI,BL,EK) ! 生成局部坐标下的单刚阵
CALL CTM(SI,CO,T) ! 计算单元坐标变换阵
CALL TTKT(EK,T) ! 计算整体坐标下的单刚阵
! 形成单元结点位移编号
DO K=1,3
LV(K)=3*(IJ(M,1)-1)+K
LV(3+K)=3*(IJ(M,2)-1)+K
ENDDO
! 单刚元素对号入座组集总刚
DO L=1,6
I=LV(L)
DO K=1,6
J=LV(K)
TK(I,J)=TK(I,J)+EK(L,K)
ENDDO
ENDDO
ENDDO
END
! Calculate element fixed-end forces.
! 计算单元固端力
SUBROUTINE EFF(L,PF,NF,BL,FO)
implicit none
REAL::BL
INTEGER::L,NF
REAL::PF(NF,4),FO(6)
REAL::Q,C,B,C1,C2,C3
INTEGER::NO
NO=INT(PF(L,2))
Q=PF(L,3)
C=PF(L,4)
B=BL-C
C1=C/BL
C2=C1*C1
C3=C1*C2
FO=0.0
SELECT CASE(NO) ! 按不同的荷载工况计算单元固端力
CASE(1) ! 跨中作用均布力
FO(2)=-Q*C*(1.0-C2+C3/2.0)
FO(3)=-Q*C*C*(0.5-2.0*C1/3.0+0.25*C2)
FO(5)=-Q*C*C2*(1.0-0.5*C1)
FO(6)=Q*C*C*C1*(1.0/3.0-0.25*C1)
CASE(2) ! 跨中作用集中力
FO(2)=-Q*B*B*(1.0+2.0*C1)/BL/BL
FO(3)=-Q*C*B*B/BL/BL
FO(5)=-Q*C2*(1.0+2.0*B/BL)
FO(6)=Q*C2*B
CASE(3) ! 跨中作用集中力矩
FO(2)=6.0*Q*C1*B/BL/BL
FO(3)=Q*B*(2.0-3.0*B/BL)/BL
FO(5)=-6.0*Q*C1*B/BL/BL
FO(6)=Q*C1*(2.0-3.0*C1)
CASE(4) ! 跨中作用线性分布力
FO(2)=-Q*C*(0.5-0.75*C2+0.4*C3)
FO(3)=-Q*C*C*(1.0/3.0-0.5*C1+0.2*C2)
FO(5)=-Q*C*C2*(0.75-0.4*C1)
FO(6)=Q*C*C*C1*(0.25-0.2*C1)
CASE(5) ! 跨中作用轴向均布力
FO(1)=-Q*C*(1.0-0.5*C1)
FO(4)=-0.5*Q*C*C1
CASE(6) ! 跨中作用轴向集中力
FO(1)=-Q*B/BL
FO(4)=-Q*C1
END SELECT
END
! Form total joint load vector.
! 形成结构综合结点荷载列阵.
SUBROUTINE JLP(NE,NJ,NP,NF,X,Y,IJ,PJ,PF,P,N)
implicit none
INTEGER::NE,NJ,NP,NF,N
REAL::X(NJ),Y(NJ),PJ(NP,3),PF(NF,4),P(N)
INTEGER::IJ(NE,2)
REAL::BL,SI,CO
INTEGER::I,J,L,M
REAL::FO(6),PE(6)
P=0.0
! 生成直接结点荷载列阵
IF(NP.GT.0) THEN
DO I=1,NP
J=INT(PJ(I,1))
L=3*(J-1)+INT(PJ(I,2))
P(L)=PJ(I,3)
ENDDO
ENDIF
! 生成综合结点荷载列阵
IF(NF.GT.0) THEN
DO L=1,NF
! 计算整体坐标下的等效结点荷载
M=INT(PF(L,1))
CALL LSC(M,NE,NJ,X,Y,IJ,BL,SI,CO)
CALL EFF(L,PF,NF,BL,FO)
PE(1)=-FO(1)*CO+FO(2)*SI
PE(2)=-FO(1)*SI-FO(2)*CO
PE(3)=-FO(3)
PE(4)=-FO(4)*CO+FO(5)*SI
PE(5)=-FO(4)*SI-FO(5)*CO
PE(6)=-FO(6)
! 生成综合结点荷载列阵
I=IJ(M,1)
J=IJ(M,2)
P(3*I-2)=P(3*I-2)+PE(1)
P(3*I-1)=P(3*I-1)+PE(2)
P(3*I)=P(3*I)+PE(3)
P(3*J-2)=P(3*J-2)+PE(4)
P(3*J-1)=P(3*J-1)+PE(5)
P(3*J)=P(3*J)+PE(6)
ENDDO
ENDIF
END
! Introduce support condition.
! 引入位移约束条件——主1副0法.
SUBROUTINE ISC(NR,JR,TK,P,N)
implicit none
INTEGER::NR,N
REAL::TK(N,N),P(N)
INTEGER::JR(NR,4)
INTEGER::I,J,L,K
DO I=1,NR
J=JR(I,1)
DO K=1,3
IF(JR(I,K+1).NE.0) THEN
L=3*(J-1)+K
TK(L,:)=0.0
TK(:,L)=0.0
TK(L,L)=1.0
P(L)=0.0
ENDIF
ENDDO
ENDDO
END
! Solution of simultaneous equations by the GAUSS elimination method.
! GAUSS消元法解线性代数方程组，求解结点位移.
SUBROUTINE GAUSS(A,B,N)
implicit none
INTEGER::N
REAL::A(N,N),B(N)
REAL::A1
INTEGER::I,J,K
DO K=1,N-1
DO I=K+1,N
A1=A(K,I)/A(K,K)
DO J=K+1,N
A(I,J)=A(I,J)-A1*A(K,J) ! 消元
ENDDO
B(I)=B(I)-A1*B(K)
ENDDO
ENDDO
B(N)=B(N)/A(N,N)
DO I=N-1,1,-1
DO J=I+1,N
B(I)=B(I)-A(I,J)*B(J) ! 回代
ENDDO
B(I)=B(I)/A(I,I)
ENDDO
END
! Print joint displacements.Calculate and print member-end forces of elements.
! 计算、输出单元结点位移和单元杆端力.
SUBROUTINE MVN(NE,NJ,NF,E,X,Y,IJ,A,ZI,PF,P,N)
USE KECALL
implicit none
REAL::E
INTEGER::NE,NJ,NF,N
REAL::X(NJ),Y(NJ),A(NE),ZI(NE),PF(NF,4),P(N)
INTEGER::IJ(NE,2)
REAL::BL,SI,CO
INTEGER::I,J,K,L,M
REAL::FO(6),F(6),D(6),TD(6),T(6,6),EK(6,6)
WRITE(2,10)
10 FORMAT(//2x,'JOINT DISPLACEMENTS'/5x,'JOINT',12X,&
&'u',14X,'v',11X,'ceta') ! 输出数据的文字标识
DO I=1,NJ
WRITE(2,15) I,P(3*I-2),P(3*I-1),P(3*I) ! 输出整体坐标下的结点位移
15 FORMAT(2X,I6,4X,3E15.6)
ENDDO
WRITE(2,25)
25 FORMAT(/2X,'MEMBER-END FORCES OF ELEMENTS'/4X,&
&'ELEMENT',13X,'N',17X,'V',17X,'M')
DO M=1,NE ! 计算并输出局部坐标下的单元杆端力
CALL LSC(M,NE,NJ,X,Y,IJ,BL,SI,CO)
CALL ESM(M,NE,E,A,ZI,BL,EK)
CALL CTM(SI,CO,T)
I=IJ(M,1)
J=IJ(M,2)
DO K=1,3
D(K)=P(3*(I-1)+K)
D(K+3)=P(3*(J-1)+K)
ENDDO
TD=0.0
TD=MATMUL(T,D)
F=0.0
F=MATMUL(EK,TD)
IF(NF.GT.0) THEN
DO L=1,NF
I=INT(PF(L,1))
IF(M.EQ.I) THEN
CALL EFF(L,PF,NF,BL,FO)
DO J=1,6
F(J)=F(J)+FO(J)
ENDDO
ENDIF
ENDDO
ENDIF
WRITE(2,80) M,(F(I),I=1,6)
80 FORMAT(2x,I8,4X,'N1=',F12.4,3X,'V1=',F12.4,3X,'M1=',&
& F12.4/14X,'N2=',F12.4,3X,'V2=',F12.4,3X,'M2=',F12.4)
ENDDO
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -