reconfea2005.f90

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

!============================================================================
Program RcConFea2005
!     Editor by Yiweiwen
!============================================================================
   use CoreCalculation
   implicit none
   integer(kind=ikind)::Nelem,Ntotv,Npoin,Nvfix,NSelem,NCelem,NFelem,Load_nodes,NLstep,istep
   integer(kind=skind)::Nnode,Ndofn,Nevab,Ngaus,type1,type2,type3,NSmat,NFmat,curvetype,LOCATECHAR
   integer(kind=skind),parameter::Ndime=2,Nstre=3

   logical if_file_exist,selfweight
   character(len=80):: filename,inputstr,refinestr
   real(kind=rkind) time_begin,time_end
   character(len=12)::Real_clock(3)
   integer(kind=ikind):: date_time(8)
   real(kind=rkind)::Fcu,Eo,Con_t,rzh,poisson !立方强度 初始弹性模量 单元厚度 容重 泊松比
   
!乱来变量   
   real(kind=rkind)::fc,Eiu,Etu,Ecu,Ec0,q 
   integer(kind=ikind)::i,j,k
   
   real(kind=rkind)::YL4(3),YL3(6),YL1(3)
!乱来变量   
    
   type(typ_Kcol),allocatable::Kcol(:)  !刚度
   type(Tptload),allocatable::ptload(:) !节点荷载

   integer(kind=ikind)::ESJtno(2),ESdof(4)  !钢筋的节点好，自由度号
   integer(kind=ikind)::EFJtno(2),EFdof(4)   !粘结单元的
   real(kind=rkind)::ESStif(4,4),ESCoord(2,2),ESdis(4),ESF(4) !钢筋的单元刚度，坐标，位移，力
   real(kind=rkind)::EFStif(4,4),EFCoord(2,2),EFdis(4),EFF(4) !粘结单元的   
   real(kind=rkind)::ECYB(3),ECYL(3)     !混凝土单元的应变和应力   
   real(kind=rkind)::T(3,3),Dmatx(3,3) !应力转换矩阵和弹性矩阵

   real(kind=rkind),allocatable::G_Pcoord(:,:) !节点荷载矩阵
   real(kind=rkind),allocatable::G_load(:)     !荷载向量
   real(kind=rkind),allocatable::All_dis(:)    !总位移
   real(kind=rkind),allocatable::Ub_G_load(:)  !不平衡荷载向量
   integer(kind=ikind),allocatable::G_Pdof(:,:) !节点定位向量
  
   real(kind=rkind),allocatable::ECstif(:,:),ECcoord(:,:),ECdis(:) !混凝土的单元刚度，坐标，位移，力
   
   integer(kind=ikind),allocatable::Ecdof(:),ECJtno(:) !混凝土单元的节点好，自由度号

   integer(kind=ikind),allocatable::G_Ecdof(:,:),G_ECJtno(:,:) !混凝土单元整体的节点好，自由度号
   integer(kind=ikind),allocatable::G_ESdof(:,:),G_ESJtno(:,:)  !钢筋的
   integer(kind=ikind),allocatable::G_EFdof(:,:),G_EFJtno(:,:) !粘结单元的

   integer(kind=skind),allocatable::SState(:),FState(:),CState(:,:) !钢筋，粘结混凝土单元的状态
   
   real(kind=rkind),allocatable::E(:,:,:),pu(:,:) !混凝土单元的E1，E2，和possion ratio
   real(kind=rkind),allocatable::diag(:),Bmatx(:,:),smatx(:,:) !一个是来分解方程，一个是应变矩阵

   integer(kind=skind),allocatable::SType(:),FType(:) !钢筋和粘结各单元的材料组号向量
   real(kind=rkind),allocatable::steprop(:,:),feltprop(:,:) !钢筋和粘结的材料向量

   real(kind=rkind),allocatable::SYL(:) !钢筋的应力

   real(kind=rkind),allocatable::FYB(:,:),FYL(:),FForec(:,:),R_FYL(:)
   !粘结单元的应变，用前一步算出的应力，力，用当前步位移算出的应力
   real(kind=rkind),allocatable::CYB(:,:,:),C_main_YB(:,:,:),CYL(:,:,:),R_CYL(:,:,:)
   !混凝土单元的应变，主应变，用前一步算出的应力，用当前步位移算出的应力

   real(kind=rkind),allocatable::Fz(:,:,:),sic(:,:,:),siu(:,:,:)
   !存两个方向的破坏应力，与破坏应力对于的应变，等效单向应变
 
 
!过程变量   
   real(kind=rkind)::area,thera,Es,Esl!混凝土单元面积，混凝土第一主应力的角度，钢筋的弹性模量,粘结单元的滑移模量
   real(kind=rkind)::s,c,L !一般为角度的sin和cos,和杆的长度
   real(kind=rkind)::KH,KV !粘结单元的水平和垂直滑移模量
   real(kind=rkind)::tmp1,tmp2,tmp3,tmp4,tmp5  !乱来的变量
   real(kind=rkind)::s0,tmax  !粘结单元的H方向最大应力对应的位移，和最大应力
   integer(kind=skind)::Tx,Ty
   real(kind=rkind)::error,wuca  !容许误差和计算过程中的误差
   real(kind=rkind)::sl(3)  !真正的应力
   integer(kind=ikind)::numb !在每一个荷载步迭代的数目
   real(kind=rkind)::EFYB(2)  !在某一步中的粘结单元的应变增量
   real(kind=rkind)::angle
   integer(kind=ikind)::destroyNum,oldstate
   !integer(kind=ikind),allocatable::destroyNelem(:)
   !iteger(kind=ikind)::Dnelem

   
   

   
!program begin

   write(*,*)  
   write(*,"(a)") '!==========欢迎使用ReConFea2005程序,本程序由华南理工大学易伟文编写==========!'
   write(*,"(a)") '!==========     RcConFea2005是一个钢筋混凝土的非线性有限元程序   ===========!'
   write(*,*)
   write(*,"(a)") '请输入数据文件的路径：'
   filename="ConFea.txt"
   write(*,"(a)") 'ConFea.txt'
   read(*,"(a)") filename
   if (filename=='') filename="ConFea.txt"
   inquire(file=filename,Exist=if_file_exist)
   if (.not.if_file_exist) then
      write(*,*) "Error,file does not exist!"
	  stop
   end if
   open (unit=10 , file = filename , status = 'old' ,    action ='read')
   write(*,"(a)") '请输入输出文件的路径：'
   filename="ConFeaout.txt"
   write(*,"(a)") 'ConFeaout.txt'
   read(*,"(a)") filename
   if (filename=='') filename="ConFeaout.txt"   
   open (unit=11 , file = filename , status = 'replace', action='readwrite') 
   
   open (unit=111 , file = "state.txt",status = 'replace', action='readwrite') 

   call date_and_time(real_clock(1),real_clock(2),real_clock(3),Date_time)
   call cpu_time(time_begin)
   
   READ(10,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) filename

   write(11,"(a)") "============================================================================================"
   write(11,"(a)") "这个是ReConFea2005的结果输出文件"
   write(11,"(a)") "============================================================================================"   
   write(11,"(A,14x,i5,A,i2.2,a,i2.2)") "Today is:",DAte_time(1),'-',date_time(2),'-',date_time(3)
   write(11,"(3(a,i2),a,i3.3)") 'The beginning time is:  ',Date_time(5),':',date_time(6),':',Date_time(7),'.',date_time(8)
   write(11,*)
   write(11,"(a)") "============================================================================================"  
   write(11,"(2a)") "本工程名称：",filename
   write(11,"(a)") "============================================================================================"  

   READ(10,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) type1
   if (type1==1) then
      write(11,"(a)") "下面进行的是弹性计算。"
   else
      write(11,"(a)") "下面进行的是弹塑性计算。"
   end if

   READ(10,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) type2
   if (type2==1) then
      write(11,"(a)") "此工程为平面应力问题。"
   else
      write(11,"(a)") "此工程为平面应变问题。"
   end if

   READ(10,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) type3
   selectcase(type3)
   case(1)
      write(11,"(a)") "本工程的混凝土单元使用T3单元，即常应变三角形单元。"
   case(2)
      write(11,"(a)") "本工程的混凝土单元使用Q4单元，即双线性四边形单元。"
   case(3)
      write(11,"(a)") "本工程的混凝土单元使用GC－Q6单元，即用广义协调条件修正的Willon广义协调四边形单元。"
   case(4)
      write(11,"(a)") "本工程的混凝土单元使用GR12单元，即具有旋转自由度的广义协调四边形单元。"
   case default
      write(11,"(a)") "本工程的混凝土单元使用GR12M单元，即在有旋转自由度的基础上引入泡沫位移的广义协调四边形单元。"
   end select  
   call confirm_parameter(type3,Nnode,Ngaus,Nevab,Ndofn)


   write(11,*)
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "本工程总体信息："
   write(11,"(a)") "============================================================================================"  

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) NCelem

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) NSelem

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) NFelem
   Nelem=NCelem+NFelem+NSelem

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) Npoin

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) Nvfix

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) NLstep

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) selfweight

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) error

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) Curvetype

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) NSmat

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) NFmat

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) load_nodes

   allocate(G_Pcoord(2,Npoin),G_Pdof(Ndofn,Npoin))
   allocate(Ecstif(nevab,nevab),Eccoord(Ndime,Nnode),ECdis(nevab))
   allocate(Ecdof(nevab),ECJtno(Nnode))
   allocate(G_ESdof(4,NSelem),G_ESJtno(2,NSelem))
   allocate(G_EFdof(4,NFelem),G_EFJtno(2,NFelem))
   allocate(G_ECJtno(nnode,NCelem),G_Ecdof(nevab,NCelem))
   allocate(SState(NSelem),FState(NFelem))
   allocate(bmatx(3,nevab),smatx(3,nevab))
   allocate(SType(NSelem),FType(NFelem))
   allocate(steprop(3,NSmat),feltprop(4,NFmat))
   allocate(SYL(NSelem))
   allocate(FYB(2,NFelem),FYL(NFelem),R_FYL(NFelem),FForec(2,NFelem))

      
   if (type3==1) then       
	  allocate(CState(1,NCelem))	 
	  allocate(E(2,1,NCelem),pu(1,NCelem))
	  allocate(CYB(3,1,NCelem),C_main_YB(2,1,NCelem),CYL(6,1,NCelem),R_CYL(2,1,NCelem))
	  allocate(Fz(2,1,NCelem),sic(2,1,NCelem),siu(2,1,NCelem))
   else      
	  allocate(CState(4,NCelem))	  
	  allocate(E(2,4,NCelem),pu(4,NCelem))
	  allocate(CYB(3,4,NCelem),C_main_YB(2,4,NCelem),CYL(6,4,NCelem),R_CYL(2,4,NCelem))
	   allocate(Fz(2,4,NCelem),sic(2,4,NCelem),siu(2,4,NCelem))
  end if
         
   read(10,"(a)") inputstr
   read(10,"(a)") inputstr
   read(10,*)  Fcu,Eo,Con_t,rzh,poisson
   
   
   write(11,*)
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "以下是本工程混凝土材料的特征参数：" 
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "材料号      抗压强度        初始弹性模量    单元厚度        容重            泊松比"
   write(11,"(i2,4x,5e16.4)") 1,Fcu,Eo,Con_t,rzh,poisson

   read(10,"(a)") inputstr
   read(10,"(a)") inputstr
   read(10,*) steprop

   write(11,*)
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "以下是本工程钢筋材料的特征参数：" 
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "材料号      弹性模量        钢筋面积        设计强度"   
   
   do i=1,NSmat
      write(11,"(i2,4x,3e16.4)") i,steprop(:,i)
   end do

   read(10,"(a)") inputstr
   read(10,"(a)") inputstr
   read(10,*) feltprop

   write(11,*)
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "以下是本工程粘结单元的特征参数：" 
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "材料号      钢筋直径和      控制长度        和x轴的角度     抗拉弹性模量"   
   
   do i=1,NFmat
      write(11,"(i2,4x,4e16.4))") i,feltprop(:,i)
   end do
   write(11,*)

   read(10,"(a)")inputstr
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") '钢筋单元材料类型向量：'
   write(11,"(a)") "============================================================================================"  
   
   SType=1
   if(NSmat>1) then       
      read(10,*) SType	  
   end if
   do i=1,NSelem     
      write(11,"(i2)",advance='no') SType(i)      
   end do