reconfea2005.f90

康清良老师编制的非线性有限元程序（2）

   write(11,"(a)",advance='yes')
   write(11,*)
   read(10,"(a)")inputstr
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") '粘结单元材料类型向量：'
   write(11,"(a)") "============================================================================================"  
   
   FType=1
   if(NFmat>1) then       
      read(10,*) FType	  
   end if
   do i=1,NFelem     
      write(11,"(i2)",advance='no') FType(i)      
   end do
   write(11,"(a)",advance='yes')
   write(11,*)
   read(10,"(a)") inputstr
   read(10,*) G_Pcoord
   
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "节点坐标信息："
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)")"节点号           X坐标           Y坐标"
   do i=1,Npoin
      write(11,"(i5,5x,3e16.3)") i,(G_pcoord(j,i),j=1,ndime)
   end do
   write(11,*)
   read(10,"(a)") inputstr
   read(10,*) G_ECJtno

   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "混凝土单元的节点信息："
   write(11,"(a)") "============================================================================================" 
   
   if (type3==1) then
      write(11,"(a)")"单元号         节点1     节点2     节点3"
      do i=1,NCelem
         write(11,"(i5,5x,3i10.5)") i,(G_ECJtno(j,i),j=1,3)
      end do
   else 
      write(11,"(a)")"单元号         节点1     节点2     节点3     节点4"
      do i=1,NCelem
         write(11,"(i5,5x,3i10.5)") i,(G_ECJtno(j,i),j=1,4)
      end do
   end if
   write(11,*)

   read(10,"(a)") inputstr
   read(10,*) G_ESJtno
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "钢筋单元的节点信息："
   write(11,"(a)") "============================================================================================" 
   write(11,"(a)")"单元号         节点1     节点2"
   do i=1,NSelem
      write(11,"(i5,5x,2i10.5)") i,(G_ESJtno(j,i),j=1,2)
   end do
   write(11,*)

   read(10,"(a)") inputstr
   read(10,*) G_EFJtno
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "粘结单元的节点信息："
   write(11,"(a)") "============================================================================================" 
   write(11,"(a)")"单元号         节点1     节点2"
   do i=1,NFelem
      write(11,"(i5,5x,2i10.5)") i,(G_EFJtno(j,i),j=1,2)
   end do
   write(11,*)


   G_pdof=1
   read(10,"(a)") inputstr
   if(Nvfix>0) read(10,*) (k,G_pdof(:,k),i=1,nvfix)
   read(10,"(a)") inputstr
   if(load_nodes>0) then
      allocate(ptload(load_nodes))  
	  read(10,*) ptload

   end if
   call form_G_pdof(G_pdof)
   Ntotv=maxval(G_pdof)
   allocate(Kcol(Ntotv),diag(Ntotv))
   allocate(G_load(0:Ntotv),All_dis(0:Ntotv),Ub_G_load(0:Ntotv))   
   do i=1,NCelem    
      ECJtno=G_ECJtno(:,i); 
      call EJointNo_to_Edof (ECJtno,G_Pdof,Ecdof)  
      G_Ecdof(:,i)=Ecdof  ! 形成单元定位向量
   end do 

   do i=1,NSelem    
      ESJtno=G_ESJtno(:,i); 
      call EJointNo_to_Edof ( ESJtno ,G_Pdof , ESdof)  
      G_ESdof(:,i)=ESdof  ! 形成单元定位向量
   end do

   do i=1,NFelem    
      EFJtno=G_EFJtno(:,i); 
      call EJointNo_to_Edof ( EFJtno ,G_Pdof , EFdof)  
      G_EFdof(:,i)=EFdof  ! 形成单元定位向量
   end do

   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "节点定位向量："
   write(11,"(a)") "============================================================================================" 
   if (ndofn==2) then
      write(11,"(a)")"节点号          自由度1     自由度2"
      do i=1,Npoin
         write(11,"(i5,5x,2i12.5)") i,(G_pdof(j,i),j=1,ndofn)
      end do
   else
      write(11,"(a)")"节点号          自由度1     自由度2     自由度3"
      do i=1,Npoin
         write(11,"(i5,5x,3i12.5)") i,(G_pdof(j,i),j=1,ndofn)
      end do
   end if
   write(11,*)
   
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "混凝土单元定位向量："
   write(11,"(a)") "============================================================================================"    
   write(11,"(a)")"节点号          自由度："
   selectcase(nevab)
   case(6)
      do i=1,NCelem
         write(11,"(i5,5x,6i10.5)") i,(G_Ecdof(j,i),j=1,nevab)
      end do
   case(8)
      do i=1,NCelem
         write(11,"(i5,5x,8i10.5)") i,(G_Ecdof(j,i),j=1,nevab)
      end do
   case(12)
      do i=1,NCelem
         write(11,"(i5,5x,12i10.5)") i,(G_Ecdof(j,i),j=1,nevab)
      end do
   end select
   write(11,*)

   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "钢筋单元定位向量："
   write(11,"(a)") "============================================================================================"    
   write(11,"(a)")"节点号          自由度："
   do i=1,NSelem
      write(11,"(i5,5x,4i10.5)") i,(G_ESdof(j,i),j=1,4)
   end do
  
   write(11,*)

   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "粘结单元定位向量："
   write(11,"(a)") "============================================================================================"    
   write(11,"(a)")"节点号          自由度："
   do i=1,NFelem
      write(11,"(i5,5x,4i10.5)") i,(G_EFdof(j,i),j=1,4)
   end do  
   write(11,*)
   call setmatband(Kcol,Ntotv,NCelem,NSelem,NFelem,G_Ecdof,G_ESdof,G_EFdof)
   write(11,"(a)") "============================================================================================"  
   write(11,"(a,i5,a)") "本工程的方程个数为:", ntotv,"个"
   write(11,"(a)") "============================================================================================"    
   
   !到此为止没有错
   
   
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   
      
   
   
   
   
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   Ub_G_load=0.0
   G_load=0.0
   all_dis=0.0
   cstate=11
   SState=1  !钢筋单元是弹性
   FState=1  !粘结单元是弹性
   E=Eo     !混凝土的E都是初始弹性模量
   pu=poisson !混凝土的泊松比都是初始泊松比
   FYB=0.0   !粘结单元的初始应变都是0.0
   CYL=0.0
   sic=-0.0033
   siu=0.0
   CYB=0.0
   CYL=0.0  
   do i=1,NLstep
      destroynum=0
      write(*,"(a,i,a)") "现在是算第",i,"步."
	  write(11,"(a,i,a)") "现在是算第",i,"步."
	  numb=1
      write(*,"(a,i,a,i,a)") "现在是算第",i,"步的第",numb,"小步."
	  write(11,"(a,i,a,i,a)") "现在是算第",i,"步的第",numb,"小步."
	  do j=1,Ntotv  
	     kcol(j)%row=0.0
      end do
!集装总刚开始      
      select case(type3)	  
	  case(1)
	     do j=1,NCelem
            ECJtno=G_ECJtno(:,j)
            Eccoord=G_Pcoord(:,ECJtno)						
	        Ecstif=0.0
			thera=CYL(6,1,j)
            call shape_T3_B(Eccoord,Bmatx,area)
			call D_matrix_2d_plane(type2,Dmatx,e(1,1,j),e(2,1,j),pu(1,j))
            call shape_YL_translate_matrix(thera,T)
            Dmatx=matmul(matmul(T,dmatx),transpose(T))		
	        smatx=matmul(dmatx,bmatx)
	        Ecstif=con_t*area*matmul(transpose(bmatx),smatx)			             									
			Ecdof=G_Ecdof(:,j) 					  
            call shape_Astif(kcol,Ecstif,Ecdof,Nevab)
		 end do
      endselect

	  do j=1,NSelem
	     ESJtno=G_ESJtno(:,j)
         ESCoord=G_Pcoord(:,ESJtno)
		 ESStif=0.0
		 call shape_ste_stif(ESCoord,ESStif)		 		 
		 if (SState(j)==1) then
		    Es=steprop(1,SType(j))
		 else
		    Es=10
		 end if	     
		 ESStif=Es*steprop(2,SType(j))*ESStif		 
	     ESdof=G_ESdof(:,j) 						  
         call shape_Astif(kcol,ESStif,ESdof,4)
	  end do

	  do j=1,NFelem
	     EFJtno=G_EFJtno(:,j)
         
		 EFStif=0.0
		 tmp1=abs(FYB(1,j))	 
		 tmp2=pi*feltprop(1,FType(j))*feltprop(2,FType(j))		 
		 if (FState(j)==1) then
		    if (curvetype==1) then
			   EsL=9.81e2-11.5e4*tmp1+2.51e6* tmp1**2
			else
			   Esl=(5.3e2-5.04e4*tmp1+17.6e5*tmp1**2-21.9e6*tmp1**3)/sqrt(Fcu*0.67/40.7)
			end if
			KH=Esl*tmp2
			Kv=feltprop(4,FType(j))*feltprop(2,FType(j)) !kv＝Ec×Lm
		else
		   KH=10.0
		   kv=feltprop(4,FType(j))*feltprop(2,FType(j)) !kv＝Ec×Lm
	    end if        
		c=cos(feltprop(3,FType(j))*pi/180)
		s=sin(feltprop(3,FType(j))*pi/180)
        call shape_flet_stif(c,s,kh,kv,EFStif)	
						    
        EFdof=G_EFdof(:,j) 						  
        call shape_Astif(kcol,EFStif,EFdof,4) 
	  end do
!集装总刚结束
      do j=1,ntotv
	     if (kcol(j)%row(j)<1e-4) then
		    write(*,"(a,i5,a)") "试件在第",i,"级荷载下破坏。由刚度矩阵为0引起" 
			write(11,"(a,i5,a)") "试件在第",i,"级荷载下破坏。由刚度矩阵为0引起"
		    stop
		 endif
	  end do	  
	  G_load=Ub_G_load	  
      Ub_G_load=0.0
	  if (selfweight>0 .and. i==1) then
	     selectcase(type3)
	     case(1)
		    do j=1,NCelem
			   ECJtno=G_ECJtno(:,j)
               Eccoord=G_Pcoord(:,ECJtno)
			   Area=Eccoord(1,1)*Eccoord(2,2)+Eccoord(1,2)*Eccoord(2,3)+Eccoord(1,3)*Eccoord(2,1)
		       Area=Area-Eccoord(1,1)*Eccoord(2,3)-Eccoord(1,2)*Eccoord(2,1)-Eccoord(1,3)*Eccoord(2,2)
			   Area=abs(0.5*Area)
			   Ecdof=G_Ecdof(:,j)
			   G_load(Ecdof(2))=g_load(Ecdof(2))-Area*rzh*Con_t/3.0
               G_load(Ecdof(4))=g_load(Ecdof(4))-Area*rzh*Con_t/3.0
			   G_load(Ecdof(6))=g_load(Ecdof(6))-Area*rzh*Con_t/3.0
			end do
		 end select
	  end if   
	  if (load_nodes>0) then
	     do j=1,load_nodes
	        G_load(G_pdof(1,ptload(j)%Jno))=G_load(G_pdof(1,ptload(j)%Jno))+ptload(j)%xload/NLstep
			G_load(G_pdof(2,ptload(j)%Jno))=G_load(G_pdof(2,ptload(j)%Jno))+ptload(j)%yload/NLstep
         end do
	  end if	     
	  call Varbandsolv(Kcol,G_load,Ntotv)	  
	  G_load(0)=0.0
	  All_dis=All_dis+G_load
	  write(11,*)"==================="
	  write(11,*) all_dis(124)	 
	     do j=1,ntotv
		    if (abs(All_dis(j))>1000) then
			   write(*,"(a,i5,a)") "试件在第",i,"级荷载下破坏，位移引起" 
			   write(11,"(a,i5,a)") "试件在第",i,"级荷载下破坏，位移引起"
			   stop
			end if
		 end do         
		 do j=1,NSelem   !钢筋单元的分析
		    ESJtno=G_ESJtno(:,j)
		    ESdof=G_ESdof(:,j)
			ESdis(:)=G_load(ESdof)
            ESCoord=G_Pcoord(:,ESJtno)
		    if (SState(j)==1) then  !弹性阶段
		       Es=steprop(1,SType(j))
		    else
		       Es=10             !水平段了
		    end if
            c=ESCoord(1,2)-ESCoord(1,1)
			s=ESCoord(2,2)-ESCoord(2,1)
			L=sqrt(c*c+s*s)
			c=c/L
			s=s/l
			SYL(j)=SYL(j)+Es*(-c*ESdis(1)-s*ESdis(2)+c*ESdis(3)+s*ESdis(4))/L !存第二点轴应力,拉为正
			     
		    if (abs(SYL(j))<steprop(3,SType(j))) then			   
			   SState(j)=1
			else
			   if (sstate(j)==1) then
			      destroynum=destroynum+1
			   end if
			   SState(j)=2
			   
               if (SYL(j)>0) then
			      tmp1=(SYL(j)-steprop(3,SType(j)))*steprop(2,SType(j))