reconfea2005.f90

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

                  SYL(j)=steprop(3,SType(j))
			   else
			      tmp1=(SYL(j)+steprop(3,SType(j)))*steprop(2,SType(j))
				  SYL(j)=-steprop(3,SType(j))
			   end if
	           ESF(1)=-tmp1*c
			   ESF(2)=-tmp1*s
			   ESF(3)=tmp1*c
			   ESF(4)=tmp1*s
			   Ub_G_load(ESdof)=Ub_G_load(ESdof)+ESF	
			   continue     
			end if               
 		 end do
		 do j=1,NFelem !粘结单元
		    EFJtno=G_EFJtno(:,j)
		    EFdof=G_EFdof(:,j)
			EFdis(:)=G_load(EFdof)
  		    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)
            EFYB(1)=-c*EFdis(1)-s*EFdis(2)+c*EFdis(3)+s*EFdis(4)  !增量
			EFYB(2)=s*EFdis(1)-c*EFdis(2)-s*EFdis(3)+c*EFdis(4)              			
			FYB(1,j)=FYB(1,j)+EFYB(1)  !粘结单元的应变H或称位移  总的
			FYB(2,j)=FYB(2,j)+EFYB(2)   !粘结单元的应变v
            !此fyb是某个荷载步的应变，不是总应变,这样才是总应变           
            if (Fstate(j)==1) then
			   FForec(1,j)=FForec(1,j)+KH*EFYB(1)  !!粘结单元的力H
			   FForec(2,j)=FForec(2,j)+Kv*EFYB(2) !!粘结单元的力V		    
			   FYL(j)=FForec(1,j)/tmp2  !H方向的应力
             end if
			!以下是粘结单元的分析
			if (curvetype==1) then
			   s0=0.011   
			   tmax=4.96
			else
			   s0=0.03
			   tmax=4.64*sqrt(Fcu*0.67/40.7)
			end if            
			tmp1=abs(FYB(1,j))  ! 刷新应变
			if (tmp1<s0) then
			   FState(j)=1
			   if(curvetype==1) then
			      R_FYL(j)=9.81e2*tmp1-5.74e4*tmp1**2+0.837e6*tmp1**3
			   else
			      R_FYL(j)=(5.3e2*tmp1-2.52e4*tmp1**2+5.87e5*tmp1**3-5.47e6*tmp1**4)*sqrt(Fcu*0.67/40.7)
			   end if
			   if (FYB(1,j)<0) R_FYL(j)=-R_FYL(j)   
			   
			   tmp3=FForec(1,j)-R_FYL(j)*tmp2               
			   
			   FForec(1,j)=R_FYL(j)*tmp2
               FYL(j)=R_FYL(j)
			else
			   if (fstate(j)==1) then

			      destroynum=destroynum+1
			   end if
			   FState(j)=2			  
			   if (FYB(1,j)<0) Tmax=-tmax
               tmp3=FForec(1,j)-tmax*tmp2               
			   FForec(1,j)=tmax*tmp2
               FYL(j)=tmax
		    end if
			EFF(1)=-tmp3*c
			EFF(2)=-tmp3*s
			EFF(3)=tmp3*c
		  	EFF(4)=tmp3*s
			Ub_G_load(EFdof)=Ub_G_load(EFdof)+EFF  	
					            
		 end do
         
         do j=1,NCelem !混凝土单元单元分析
		    ECJtno=G_ECJtno(:,j)
		    Ecdof=G_Ecdof(:,j)
			ECdis(:)=G_load(Ecdof)
            Eccoord=G_Pcoord(:,ECJtno)
			call shape_T3_B(Eccoord,Bmatx,area)
			ECYB=matmul(Bmatx,ECdis)   !单元当前步的应变
			CYB(:,1,j)=CYB(:,1,j)+ECYB  !所有单元当前步的应变

			tmp1=0.5*sqrt((CYB(1,1,j)-CYB(2,1,j))**2+CYB(3,1,j)**2)			
			
			C_main_YB(1,1,j)=(CYB(1,1,j)+CYB(2,1,j))*0.5+tmp1
            
			C_main_YB(2,1,j)=(CYB(1,1,j)+CYB(2,1,j))*0.5-tmp1  !主应变
            
			thera=CYL(6,1,j)        
			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))								
			ECYL=matmul(Dmatx,ECYB)			
			CYL(1:3,1,j)=CYL(1:3,1,j)+ECYL  !累加的应力			
			tmp1=CYL(1,1,j)+CYL(2,1,j)			
			tmp2=sqrt((CYL(1,1,j)-CYL(2,1,j))**2+4*CYL(3,1,j)**2)			
			CYL(4,1,j)=(tmp1+tmp2)*0.5
			CYL(5,1,j)=(tmp1-tmp2)*0.5			
			Tx=CState(1,j)/10
			Ty=CState(1,j)-10*Tx
			oldstate=cstate(1,j)							
			if(Tx>=2 .or. Ty>=2) then
			   angle=cyl(6,1,j)
			   goto 1000  !非弹性状态11
            end if
            tmp3=CYL(2,1,j)-CYL(5,1,j)
			if (abs(tmp3)<1e-4) then			   
			   angle=0.0
			else			  
			  angle=90-57.296*atan(CYL(3,1,j)/tmp3)
			end if            !求角度
1000        tmp1= ECYL(1)+ECYL(2)
            tmp2=sqrt((ECYL(1)-ECYL(2))**2+4*ECYL(3)**2)			
			tmp3=(tmp1+tmp2)*0.5
			tmp4=(tmp1-tmp2)*0.5 
            if ( abs(cyl(6,1,j)-angle)<45 .or. abs(cyl(6,1,j)-angle+180)<45 .or. abs(cyl(6,1,j)-angle-180)<45) then
			   siu(1,1,j)=siu(1,1,j)+tmp3/E(1,1,j) 
			   siu(2,1,j)=siu(2,1,j)+tmp4/E(2,1,j)  !等效单向应变的累加
            else
			   tmp5=siu(1,1,j)
			   siu(1,1,j)=siu(2,1,j)+tmp3/e(1,1,j)
			   siu(2,1,j)=tmp5+tmp4/e(2,1,j)   
			end if
            cyl(6,1,j)=angle			
			if (CState(1,j)==11) then
			   call kupf(Fcu,CYL(:,:,j),Fz(:,:,j))  !算破坏应力
			end if						
	    fc=-Fcu*0.667		   
            if (CYL(5,1,j)<0.0 .and. CYL(4,1,j)<0.0) then  !双压
			   if (cstate(1,j)/=11) then
			      write(111,*) cstate(1,j)
			   end if
			  ! if (cstate(1,j)==11) then
			      if (abs(CYL(5,1,j))<abs(fz(2,1,j))) then			      
			         call saenz(1,Fz(:,:,j),Fc,sic(:,:,J),siu(:,:,J),E(:,:,j),R_CYL(:,:,J),Eo,cyl(:,1,j)) !算E和真实应力
				     CState(1,j)=11
				     pu(1,j)=0.2
				     sl(1)=R_CYL(1,1,j)
				     sl(2)=R_CYL(2,1,j)
				  else
				     do k=1,2
				        if (abs(siu(k,1,j))<abs(4*sic(k,1,j)) .and. abs(siu(k,1,j))>=abs(sic(k,1,j))) then
					       sl(k)=fz(k,1,j)-(0.8*fz(k,1,j)*(siu(k,1,j)-sic(k,1,j)))/(3.0*sic(k,1,j))      			            					
						   CState(1,j)=22												
					    else if (abs(siu(k,1,j))>abs(4*sic(k,1,j))) then
                           sl(k)=0.0					
                           CState(1,j)=33						
					    end if			
				     enddo
				     E(:,1,j)=10
				     pu(1,j)=0.2
                  end if
               !elseif(cstate(1,j)==22) then
			    !  continue
			   !elseif(cstate(1,j)==33) then
			    !  continue
               !elseif(cstate(1,j)==12) then
			    !  continue
               !elseif(cstate(1,j)==13) then
			    !  continue
               !elseif(cstate(1,j)==21) then
			    !  continue
               !elseif(cstate(1,j)==23) then
			    !  continue
               !elseif(cstate(1,j)==31) then
			    !  continue
			   !elseif(cstate(1,j)==32) then
			    !  continue
			   !end if

			   
			elseif (CYL(4,1,j)>=0.0 .and. CYL(5,1,j)<0.0) then  !拉压
			   if (abs(CYL(4,1,j))<abs(fZ(1,1,j))) then
		          call saenz(2,Fz(:,:,j),Fc,sic(:,:,J),siu(:,:,J),E(:,:,j),R_CYL(:,:,J),Eo,cyl(:,1,j))
                  R_CYL(1,1,j)=CYL(4,1,j)
				  sl(1)=R_CYL(1,1,j)
				  sl(2)=R_CYL(2,1,j)
				  CState(1,j)=11
				  pu(1,j)=0.2
			   else 
			      Eiu=abs(siu(1,1,j))
			      Etu=0.00015
				  if (eiu<etu .and. eiu>=abs(Fz(1,1,j)/Eo)) then
				     sl(1)=fz(1,1,j)							           
					 tx=2					 
                  else if(eiu>=etu) then
				     sl(1)=0.0				
					 tx=3				     					 
				  end if  			
			      E(1,1,j)=10
				  pu(1,j)=0.2
				  !以下是受压
				  Eiu=abs(siu(2,1,j))
				  Ecu=0.0033
				  Ec0=0.002
				  if (eiu<Ec0) then
				     Es=abs(Fc/Ec0)
					 q=eiu/Ec0
					 s=1+(Eo/es-2.0)*q+q**2
					 E(2,1,j)=eo*(1.0-q**2)/s**2
					 sl(2)=-eo*eiu/s
					 ty=1
				  elseif (eiu>=ec0 .and. eiu<ecu)then
				     e(2,1,j)=10
					 sl(2)=fc-0.8*fc*(eiu-ec0)/(ecu-ec0)
					 ty=2					 
                  else
				     sl(2)=0.0
					 e(2,1,j)=10
					 ty=3								    
				  end if					 			   
			   end if			   
               pu(1,j)=0.2
			   CState(1,j)=10*tx+ty			
			elseif (CYL(4,1,j)>=0.0 .and. CYL(5,1,j)>=0.0) then  !双拉
			   if (abs(CYL(4,1,j))<abs(fz(1,1,j))) then
			      sl(1)=CYL(4,1,j)
				  sl(2)=CYL(5,1,j)
				  CState(1,j)=11
			   else
			      eiu=abs(siu(1,1,j))
				  etu=0.00015
				  if (eiu>=abs(fz(1,1,j)/Eo) .and. eiu<etu) then
				     sl(1)=fz(1,1,j)
					 tx=2					 
				   elseif (eiu>=etu) then
				      sl(1)=0.0
					  tx=3					 
				   end if				 
				   e(1,1,j)=0.0
				   pu(1,j)=0.2
				   !小主应力方向
				   eiu=abs(siu(2,1,j))
				   if (abs(CYL(5,1,j))<abs(fz(2,1,j)) .or. eiu<abs(fz(2,1,j)/Eo)) then
				      sl(2)=CYL(5,1,j)
					  ty=1
					  e(2,1,j)=eo  !这句不写我觉得更好
					  pu(1,j)=0.2
				   elseif (eiu>=abs(fz(2,1,j)/eo) .and. eiu<etu) then
				      sL(2)=fz(2,1,j)
					  ty=2					 
					  E(2,1,j)=10
					  pu(1,j)=0.2
				   else
				      sl(2)=0
					  ty=3					  
					  E(2,1,j)=10
					  pu(1,j)=0.2
                   end if
			      CState(1,j)=10*tx+ty			
               end if
			end if	
			if (oldstate/=cstate(1,j)) then
			   destroynum=destroynum+1
			end if			
			thera=cyl(6,1,j)
			call shape_YL_translate_matrix(thera,T)														   		  
		    yl4(1)=CYL(4,1,j)-sl(1)
		    yl4(2)=CYL(5,1,j)-sl(2)
		    yl4(3)=0.0
		    yl1=matmul(T,yl4)
	   	    call shape_T3_B(Eccoord,Bmatx,area)     
		    yl3=area*Con_t*matmul(transpose(bmatx),yl1) 
            Ub_G_load(Ecdof)=Ub_G_load(Ecdof)+yl3		    
			yl4(1)=sl(1)
		    yl4(2)=sl(2)
		    yl4(3)=0.0            
			yl1=matmul(T,yl4)
		    CYL(1:3,1,j)=yl1
		    CYL(4,1,j)=sl(1)
		    CYL(5,1,j)=sl(2)  !这里为什么不进一步算角度 我求	
			IF (CState(1,j)==11)THEN
			   tmp3=CYL(2,1,j)-CYL(5,1,j)
			   if (abs(tmp3)<1e-4) then			   
			      angle=0.0
			   else			  
			      angle=90-57.296*atan(CYL(3,1,j)/tmp3)
			   end if
			   CYL(6,1,J)=ANGLE
		    END IF			
		 end do		      
200   Ub_G_load(0)=0.0      
      if (i==8) then
	  continue
	  end if
      wuca=max(abs(maxval(Ub_G_load)),abs(minval(Ub_G_load)))	  
	  !if (destroynum>0  .or. (wuca>error .and. numb<100)) then
	  if (wuca>error .and. numb<100) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         destroyNum=0
	     numb=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