4.txt

文件1.txt,2.txt,3.txt和5.txt为用Fortran编写的有限元程序 4.txt为用c＋＋编写的钢筋混凝土异形柱的全过程非线性分析源程序

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                     杆系大挠度分析程序
!	 
! (取自《杆系与箱形梁桥结构分析及程序设计》，赵振铭 陈宝春编著，华南理工大学出版社)
!
!变量简介:
!NJ   节点数			NE   单元数
!NPJ  节点荷载数		NZC  支撑数
!E    弹性模量			NL   单元左端点数组
!NR   单元右端点数组	X    节点横坐标数组
!Y    节点纵坐标数组	EF   节点截面积数组
!EJ   单元截面抗弯惯矩	P    荷载列阵  (不平衡力列阵)
!P2   
!ZK   总刚矩阵			DK	 单元刚度矩阵
!DP   位移数组			STS  杆端力数组
!DK0  弹性刚度矩阵      DK1  大位移矩阵
!DK2  几何刚度矩阵		TX   轴力
!
!last modified by ZhangTaike. ^_^
!
		  PROGRAM PPRE
		  IMPLICIT INTEGER*2(I-N)
		  DIMENSION NL(50),NR(50),X(50),Y(50),EF(50),EJ(50)	   !节点只有50个
		  DIMENSION P(150),P1(150),PJ(30,2),NZ(20),P2(150)
		  DIMENSION ZK(150,150),DK(6,6),AL(6,6),AT(6,6)
		  DIMENSION STS(50,6),D1(6),D2(6),TX(150), &
		         &  DP(150),F(150),Q(150)
          COMMON/C1/E
		  COMMON/C2/I0,J0,SI,CO,EL
		  COMMON/C3/AL,AT

		  !打开数据输入文件，输出文件
		  OPEN(7,FILE='FIN.TXT',STATUS='OLD')
		  OPEN(8,FILE='FOUT.TXT',STATUS='UNKNOWN')
 		  OPEN(10,FILE='FOR.DBG',STATUS='UNKNOWN')	  !for DEBUG
													   

		  READ(7,*)KN  				 !线性分析非线性开关，1为考虑非线性。
		  IF(KN.NE.0)  READ(7,*)ROL	 !判敛条件
		  READ(7,*)NJ,NE,NPJ,NZC,E
		  WRITE(8,10)NJ,NE,NPJ,NZC,E
10        FORMAT(3X,'NJ=',I4,3X,'NE=',I4,3X,/,3X, &
              &  'NPJ=',I4,3X,'NZC=',I4,3X,' E=',F20.4)
          READ(7,*)(NL(I),NR(I),I=1,NE)
		  !write(8,'(" Element  Node1 Node2")')
		  WRITE(8,20)
20        FORMAT(/,3X,'NO=',10X,'NL=',10X,'NR=')
          WRITE(8,30)(I,NL(I),NR(I),I=1,NE)
30        FORMAT(/,3X,I4,10X,I4,10X,I4)

          READ(7,*)(X(I),Y(I),I=1,NJ)
		  WRITE(8,40)
40        FORMAT(/,3X,'NO=',10X,'X=',10X,'Y=')
          WRITE(8,50)(I,X(I),Y(I),I=1,NJ)
50        FORMAT(/,3X,I4,5X,F12.4,5X,F12.4)

          READ(7,*)(EF(I),EJ(I),I=1,NE)
		  WRITE(8,60)
60        FORMAT(/,3X,'NO=',10X,'EF=',10X,'EJ=')
          WRITE(8,70)(I,EF(I),EJ(I),I=1,NE)
70        FORMAT(/,3X,I4,F12.4,3X,F12.4)

          READ(7,*)((PJ(I,J),J=1,2),I=1,NPJ)
		  WRITE(8,80)
80        FORMAT(/,3X,'NO=',10X,'PJ(I,1)=',10X,'PJ(I,2)=')
          WRITE(8,90)(I,(PJ(I,J),J=1,2),I=1,NPJ)
90        FORMAT(/,3X,I4,5X,F12.4,5X,F12.4)
          READ(7,*)(NZ(I),I=1,NZC)	   !约束数据
		  WRITE(8,820)
820       FORMAT(/,5X,'NO=',10X,'NZ(I)=')
          WRITE(8,830)(I,NZ(I),I=1,NZC)
830       FORMAT(/,5X,I4,12X,I4)

		  !求解初始化
          NN=3*NJ		             !NN 总自由度
		  CALL CLEAR(TX,150,1)
		  CALL CLEAR(P,150,1)
		  CALL CLEAR(DP,150,1)
		  CALL CLEAR(ZK,150,150)
		  CALL CLEAR(STS,50,6)

		  DO 200 MJ=1,NPJ			  !
			  MI=PJ(MJ,2)			  !PJ(*,2)  节点荷载施加自由度
			  P(MI)=P(MI)+PJ(MJ,1)	  !形成荷载列阵
200       CONTINUE

	  !for DEBUG
	  write(10,'("Load Vector:")')
      write(10,*)( P(j), j = 1,NN)

	  Iter = 0
810       CALL CLEAR(ZK,150,150)		!总刚清零 IF(R.GT.ROL.AND.KN.NE.0)GOTO 810
		  
			  Iter = Iter + 1
			  write(10,'("第 ",I3,"次跌代")')Iter
  			  WRITE(10,'("形成总刚阵求解")')

  			  CALL FZK(NJ,NE,NL,NR,X,Y,ZK,DP,NN,TX,EF,EJ)  ! 形成总刚
			  CALL FNZC(NZC,ZK,NN,P,NZ)	
			  				   ! 边界处理
			  DO 750 J=1,NN
				 P1(J)=P(J)			  !外荷载 / 不平衡力
750	          CONTINUE


			  WRITE(10,*)"求解荷载 "
			  DO 9751 I=1,NN
9751			 WRITE(10,'(E12.3)')P(I)
			  WRITE(10,*)

			  CALL GAS(NN,ZK,P)	      !解方程


			  DO 812 J=1,NN			 
				 F(J)=P(J)			  !当前位移，F()用于判敛
812	          CONTINUE

			  DO 710 J=1,NN			  
				 DP(J)=DP(J)+P(J)	  !累计位移
710	          CONTINUE

			  WRITE(10,*)"位移解 "
			  WRITE(10,*)"累计位移  当前位移"
			  DO 9755 I=1,NN
9755			 WRITE(10,'(2E12.3)')DP(I),P(I)
			  WRITE(10,*)

			  CALL CLEAR(TX,150,1)	  !轴力清零，注意位置
			  CALL CLEAR(STS,50,6)

			  DO 640 K=1,NE			  !对单元循环	 更新矩阵
				  CALL CHL(K,NJ,NE,NL,NR,X,Y)
				  CALL MLA					   !方向矩阵
				  CALL CLEAR(DK,6,6)
   				  WRITE(10,'("形成单刚阵求单元杆端力")')

				  IF(KN.EQ.0)THEN						  !线性
					  CALL CLEAR(Q,150,1)
					  CALL FDK(K,NE,DK,EF,EJ,TX(K),Q,NN)
	 			  ELSE
					  CALL FDK(K,NE,DK,EF,EJ,TX(K),DP,NN)  !非线性
				  END IF
				  DO 650 J=1,3
					  D1(J)=DP(I0+J)           !取单元节点位移
					  D1(J+3)=DP(J0+J)
650				  CONTINUE
				  CALL MATMU(AL,D1,D2,6,6,1)   !整体to局部
				  CALL MATMU(DK,D2,D1,6,6,1)   !局部

				  TX(K)=D1(4)				   !轴力

				  DO 950 J=1,6				   !
					 STS(K,J)=STS(K,J)+D1(J)   !杆端力
950		          CONTINUE
640	          CONTINUE


   			  WRITE(10,'("形成总刚阵计算杆端力列阵，计算不平衡力")')
			  CALL CLEAR(ZK,150,150)
			  CALL FZK(NJ,NE,NL,NR,X,Y,ZK,DP,NN,TX,EF,EJ)	   !重新组集总刚
			  CALL FNZC(NZC,ZK,NN,Q,NZ)

			  WRITE(10,*)"Globla Stiff Matrix:"
			  DO 101 I=1,6
				WRITE(10,103)(ZK(I,J),J=1,6)
103			  FORMAT(1X,6e12.3,';')
			  WRITE(10,*)
101           CONTINUE

			  CALL CLEAR(P2,150,1)
			  DO 751 I=1,NN
				  DO 752 MM=1,NN
					  P2(I)=P2(I)+ZK(I,MM)*P(MM)		  !形成杆端力列阵，当前位移or累计位移？
752		          CONTINUE
751		      CONTINUE

			  DO 754 I=1,NN
				  P(I)=P1(I)-P2(I)		!形成不平衡力  			
754	          CONTINUE

			  WRITE(10,*)"不平衡力计算  P ＝ P1 － P2"
			  WRITE(10,*)"P      P1     P2"
			  DO 9753 I=1,NN
9753			 WRITE(10,'(3E12.3)')P(I),P1(I),P2(I)
			  WRITE(10,*)


			  !判敛处理
			  CALL CLEAR(Q,150,1)
			  DO 756 I=1,NN
				  IF(DP(I).NE.0.0)Q(I)=ABS(F(I)/DP(I))
756	          CONTINUE
			  R=0.0
			  DO 758 I=1,NN
				  IF(R.LT.Q(I))R=Q(I)
758	          CONTINUE

          !IF(R.GT.ROL.AND.KN.NE.0.AND.Iter.LT.3)GOTO 810		  !跌代循环判敛
          IF(R.GT.ROL.AND.KN.NE.0)GOTO 810		  !跌代循环判敛

		  write(8,'("跌代总次数 ",I3)')Iter

		  !位移输出
		  WRITE(8,620)
620       FORMAT(/,20X,'THE DISPLACEMENTS',/, &
              &  3X,'NO=',15X,'U=',15X,'V=',15X,'ACF=')
          DO 930 J=1,NJ
		  WRITE(8,630)J,DP(3*J-2),DP(3*J-1),DP(3*J)
630       FORMAT(3X,I4,5X,E12.3,5X,E12.3,5X,E12.3)
930       CONTINUE

		  !内力输出
          WRITE(8,969)
969       FORMAT(/,10X,'FORCE COMPONENTS OF ELEMENTS',/)
          WRITE(8,970)
970       FORMAT(1X,'NO=',7X,'NI=',7X,'QI=', &
              &  7X,'MI=',7X,'NJ=',7X,'QJ=',7X,'MJ=')
          DO 1000 M=1,NE
		  WRITE(8,980)M,STS(M,1),STS(M,2),STS(M,3), &
		      &       STS(M,4),STS(M,5),STS(M,6)
980       FORMAT(1X,I4,1X,E10.3,1X,E10.3,1X,E10.3, &
              &  1X,E10.3,1X,E10.3,1X,E10.3)
1000      CONTINUE
          END
!-------------------------子程序FNZC-------------------------
!边界条件
          SUBROUTINE FNZC(NZC,ZK,NN,P,NZ)
		  IMPLICIT INTEGER*2(I-N)
		  DIMENSION ZK(150,150),P(150),NZ(20)
		  DO 600 I=1,NZC
		  IZ=NZ(I)
		  DO 610 J=1,NN
		  ZK(IZ,J)=0.0
		  ZK(J,IZ)=0.0
610       CONTINUE
          ZK(IZ,IZ)=1.0
		  P(IZ)=0.0
600       CONTINUE
          RETURN
		  END
!-------------------------子程序FZK--------------------------
!总刚阵组集
!
!
          SUBROUTINE FZK(NJ,NE,NL,NR,X,Y,ZK,P,NN,TX,EF,EJ)
		  IMPLICIT INTEGER*2(I-N)
		  DIMENSION ZK(150,150),NL(50),NR(50),X(50),Y(50),P(150)
		  DIMENSION EF(50),EJ(50),TX(50)
		  DIMENSION DK(6,6),AL(6,6),AT(6,6),AN(6,6)
		  COMMON/C1/E
		  COMMON/C2/I0,J0,SI,CO,EL
		  COMMON/C3/AL,AT
		  DO 500 K=1,NE
			  CALL CHL(K,NJ,NE,NL,NR,X,Y)
			  CALL MLA
			  CALL FDK(K,NE,DK,EF,EJ,TX(K),P,NN)
			  CALL MATMU(DK,AL,AN,6,6,6)
			  CALL MATMU(AT,AN,DK,6,6,6)
			  DO 510 I=1,3
				  DO 510 J=1,3
					  ZK(I0+I,I0+J)=ZK(I0+I,I0+J)+DK(I,J)
					  ZK(I0+I,J0+J)=ZK(I0+I,J0+J)+DK(I,J+3)
					  ZK(J0+I,I0+J)=ZK(J0+I,I0+J)+DK(I+3,J)
					  ZK(J0+I,J0+J)=ZK(J0+I,J0+J)+DK(I+3,J+3)
510			  CONTINUE

500       CONTINUE
          RETURN
		  END
!-------------------------子程序FDK--------------------------
!形成单元切线刚度矩阵  DK
!K 单元索引；NE 单元总数；DK 切线刚阵；EF 单元截面积；EJ 单元抗弯惯矩
!TN 初始轴力；P ？；NN 结构总自由度
!
!
!
          SUBROUTINE FDK(K,NE,DK,EF,EJ,TN,P,NN)
		  		  IMPLICIT INTEGER*2(I-N)
		  DIMENSION DK(6,6),DK0(6,6),DK1(6,6),DK2(6,6)
		  DIMENSION EF(50),EJ(50),P(150),PZ(6),PB(6),AL(6,6),AT(6,6)
		  COMMON/C1/E
		  COMMON/C2/I0,J0,SI,CO,EL
		  COMMON/C3/AL,AT
		  CALL CLEAR(DK,6,6)
		  CALL FDK0(DK0,EF(K),EL,EJ(K))