module.f90

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

!*************************************************************************
module NumKind
!*************************************************************************
      
   implicit none
   integer(kind=1),parameter::ikind=kind(1),rkind=kind(0.D0),skind=1
   real(rkind),parameter     ::Zero=0.D0,One=1.D0,Two=2.D0,Three=3.D0,Four=4.D0, &
                               Five=5.D0,Six=6.D0,Seven=7.D0,Eight=8.D0, &
							   Nine=9.D0,Ten=10.D0,Eleven=11.D0,Twelve=12.D0,pi=3.1416
end module NumKind

!*************************************************************************
module TypeDef
!*************************************************************************
   use Numkind
   implicit none
   type,public:: Typ_Kcol
      real(kind=rkind),pointer::row(:)
   end type Typ_kcol   
   
   type,public::Tptload
       integer(kind=ikind)::Jno
	   real(kind=rkind)::xload,yload
	end type Tptload
end module TypeDef

!*************************************************************************
module CoreCalculation
!*************************************************************************
   use TypeDef
   implicit none
   contains
   !============================================================================
   subroutine confirm_parameter(type3,Nnode,Ngaus,Nevab,Ndofn)
   !============================================================================      
      integer(kind=skind),intent(in)::type3
	  integer(kind=skind),intent(out)::Nnode,Ngaus,Nevab,Ndofn
      select case(type3)      
      case(1)  	     
	     nnode=3;Ngaus=1;Nevab=6;ndofn=2
	  case(2)  
	     nnode=4;Ngaus=2;Nevab=8;ndofn=2	     
	  case(3)  	     
	     nnode=4;Ngaus=2;Nevab=8;ndofn=2
      case(4) 	     
	     nnode=4;Ngaus=2;Nevab=12;ndofn=3	  
	  case default          
		 nnode=4;Ngaus=2;Nevab=12;ndofn=3
	  	    
      end select
    end subroutine confirm_parameter
   
    !============================================================================
   subroutine form_G_Pdof(G_Pdof) ! reform G_Pdof 节点定位向量     
   !============================================================================
      integer,intent(in out)::G_Pdof(:,:)
      integer(kind=ikind):: i,j,m
      m=0
      do j=1,size(G_Pdof,2)
         do i=1,ubound(G_Pdof,1)
            if(G_Pdof(i,j)/=0) then
               m=m+1
			   G_Pdof(i,j)=m
            end if
         end do
      end do
   end subroutine form_G_Pdof

   !============================================================================
   subroutine EJointNo_to_Edof(EJointNo,G_Pdof,Edof) !节点定位向量到单元定位向量
   !============================================================================
      integer(kind=ikind),intent(in)::EJointNo(:),G_Pdof(:,:)  
      integer(kind=ikind),intent(out):: Edof(:)
      integer(kind=ikind)::i,k,Nnode,Ndofn 
      Nnode=ubound(EJointNo,1) 
      Ndofn=ubound(G_Pdof,1)
      do i = 1 , Nnode
         k = i*Ndofn  
	     Edof(k-Ndofn+1:k) = G_Pdof( : , EJointNo(i) )
      end do
   end subroutine EJointNo_to_Edof 

   !============================================================================
   subroutine setmatband(Kcol,Ntotv,NCelem,NSelem,NFelem,G_Ecdof,G_ESdof,G_EFdof)   !确定变带宽
   !============================================================================
      type(typ_Kcol),intent(out)::Kcol(:)
      integer(kind=ikind),intent(in)::G_Ecdof(:,:),G_ESdof(:,:),G_EFdof(:,:)
	  integer(kind=ikind),intent(in)::Ntotv,NCelem,NSelem,NFelem
      integer(kind=ikind)::mindof,ielem,j
      integer(kind=ikind)::row1(Ntotv)  
      row1=Ntotv   !设置最大数

      do ielem=1,NCelem
         mindof=minval(G_Ecdof(:,ielem),mask=G_Ecdof(:,ielem)>0)
	     where (G_Ecdof(:,ielem)>0)
	        row1(G_Ecdof(:,ielem))=min(row1(G_Ecdof(:,ielem)),mindof)
	     end where
      end do

	  do ielem=1,NSelem
         mindof=minval(G_ESdof(:,ielem),mask=G_ESdof(:,ielem)>0)
	     where (G_ESdof(:,ielem)>0)
	        row1(G_ESdof(:,ielem))=min(row1(G_ESdof(:,ielem)),mindof)
	     end where
      end do

	  do ielem=1,NFelem
         mindof=minval(G_EFdof(:,ielem),mask=G_EFdof(:,ielem)>0)
	     where (G_EFdof(:,ielem)>0)
	        row1(G_EFdof(:,ielem))=min(row1(G_EFdof(:,ielem)),mindof)
	     end where
      end do



      do j=1,Ntotv
         allocate(kcol(j)%row(row1(j):j))
	     kcol(j)%row=0.0
      end do
	 
   end subroutine setmatband

   !============================================================================
   subroutine Delmatband(Kcol)   !清空
   !============================================================================
      type(typ_Kcol),intent(in out)::Kcol(:)
	  integer(kind=ikind)::j,nglbdof
      nglbdof=size(kcol,1)
	  do j=1,nglbdof
         deallocate(kcol(j)%row)
    end do
   end subroutine Delmatband
 
   !============================================================================
   subroutine shape_YL_translate_matrix(thera,T) !形成应力转换矩阵T，thera为角度制
   !============================================================================
      real(kind=rkind),intent(in)::thera
	  real(kind=rkind),intent(out)::T(3,3)

	  real(kind=rkind)::c,s
	  c=cos(thera*pi/180)
	  s=sin(thera*pi/180)
	  
	  T(1,1)=c*c;T(1,2)=s*s;T(1,3)=-2.0*s*c
	  T(2,1)=s*s;T(2,2)=c*c;T(2,3)=2.0*s*c
	  T(3,1)=s*c;T(3,2)=-1.0*s*c;T(3,3)=c*c-s*s
   end subroutine shape_YL_translate_matrix
   
   !============================================================================
   subroutine shape_T3_B(Eccoord,Bmatx,s) !三角形单元的B矩阵
   !============================================================================
      real(kind=rkind),intent(in)::Eccoord(:,:)
	  real(kind=rkind),intent(out)::bmatx(:,:),s 
	
	  integer(kind=ikind)::i,j
	  bmatx=0.0
	  bmatx(1,1)=Eccoord(2,2)-Eccoord(2,3)
	  bmatx(3,2)=bmatx(1,1)

	  bmatx(1,3)=Eccoord(2,3)-Eccoord(2,1)  
	  bmatx(3,4)=bmatx(1,3)

	  bmatx(1,5)=Eccoord(2,1)-Eccoord(2,2)
	  bmatx(3,6)=bmatx(1,5)
	  
	  bmatx(2,2)=Eccoord(1,3)-Eccoord(1,2)
	  bmatx(3,1)=bmatx(2,2)
	  
	  bmatx(2,4)=Eccoord(1,1)-Eccoord(1,3)
	  bmatx(3,3)=bmatx(2,4)
	  
	  bmatx(2,6)=Eccoord(1,2)-Eccoord(1,1)
	  bmatx(3,5)=bmatx(2,6)
      
	  s=abs((bmatx(1,3)*bmatx(2,6)-bmatx(1,5)*bmatx(2,4))*0.5)
	  
	  do i=1,3
	  
	     do j=1,6
		    bmatx(i,j)=0.5*bmatx(i,j)/s
	     end do
	  end do
   end subroutine shape_T3_B
   
   !============================================================================
   subroutine D_matrix_2d_plane(type2,Dmatx,e1,e2,pu)
   !============================================================================
      integer(kind=skind),intent(in)::type2
      real(kind=rkind),intent(in)::e1,e2,pu
      real(kind=rkind),intent(out)::Dmatx(:,:)
      real(kind=rkind)::tmpe1,tmpe2,tmppu  
      tmpe1=e1;tmpe2=e2;tmppu=pu
      Dmatx=0.0  
	  if (type2==2) then
	     tmpe1=tmpe1/(1-tmppu*tmppu)
		 tmpe2=tmpe2/(1-tmppu*tmppu)
		 tmppu=tmppu/(1-tmppu)
	  end if   
	  dmatx(1,1)=tmpe1/(1-tmppu*tmppu)
	  dmatx(2,2)=tmpe2/(1-tmppu*tmppu)
	  dmatx(1,2)=tmppu*sqrt(tmpe1*tmpe2)/(1-tmppu*tmppu)
	  dmatx(2,1)=dmatx(1,2)
      dmatx(3,3)=0.25*(tmpe1+tmpe2-2*tmppu*sqrt(tmpe1*tmpe2))/(1-tmppu*tmppu)	 
	end subroutine D_matrix_2d_plane   

   !============================================================================
   subroutine shape_Astif(kcol,Estif,Edof,Nevab)
   !============================================================================
      type(typ_Kcol),intent(out)::Kcol(:)
      real(kind=rkind),intent(in)::Estif(:,:)
      integer(kind=ikind),intent(in)::Edof(:)
	  integer(kind=skind),intent(in)::Nevab
      integer(kind=ikind)::j,jgdof
      do j=1,Nevab
	     jgdof=Edof(j)
		 if(jgdof==0) cycle
		 where (Edof<=jgdof .and. Edof>0)
		    kcol(jgdof)%row(Edof)=kcol(jgdof)%row(Edof)+Estif(:,j)
		 end where
	  end do
   end subroutine shape_Astif
   
   !============================================================================
   subroutine shape_ste_stif(ESCoord,ESStif)
   !============================================================================
      real(kind=rkind),intent(in)::ESCoord(:,:)
	  real(kind=rkind),intent(out)::ESStif(:,:)

	  real(kind=rkind)::L,x1,x2,y1,y2,c,s
	  x1=ESCoord(1,1)
	  y1=ESCoord(2,1)
	  x2=ESCoord(1,2)
	  y2=ESCoord(2,2)
	  L=sqrt((x1-x2)**2+(y1-y2)**2)
      c=(x2-x1)/l
	  s=(y2-y1)/l
	  ESStif(1,1)=c*c/l
	  ESStif(1,2)=c*s/l
	  ESStif(1,3)=-c*c/l
	  ESStif(1,4)=-c*s/l
	  ESStif(2,1)=ESStif(1,2)
      ESStif(3,1)=ESStif(1,3)
      ESStif(4,1)=ESStif(1,4)

	  ESStif(2,2)=s*s/l
	  ESStif(2,3)=-c*s/l
	  ESStif(2,4)=-s*s/l
      ESStif(3,2)=ESStif(2,3)
	  ESStif(4,2)=ESStif(2,4)

	  ESStif(3,3)=c*c/l
	  ESStif(3,4)=c*s/l
	  ESStif(4,3)=ESStif(3,4)

	  ESStif(4,4)=s*s/l
   end subroutine 
   
   !============================================================================
   subroutine shape_flet_stif(c,s,kh,kv,EFStif)
   !============================================================================
      real(kind=rkind),intent(in)::c,s,kh,kv
	  real(kind=rkind),intent(out)::EFStif(:,:)

	  EFStif(1,1)=kh*c*c+kv*s*s
	  EFStif(1,2)=kh*c*s-kv*c*s
	  EFStif(1,3)=-kh*c*c-kv*s*s	  
	  EFStif(1,4)=-kh*c*s+kv*c*s

	  EFStif(2,1)=EFStif(1,2)			
	  EFStif(3,1)=EFStif(1,3)
   	  EFStif(4,1)=EFStif(1,4)
	  
	  EFStif(2,2)=kh*s*s+kv*c*c
	  EFStif(2,3)=-kh*c*s+kv*c*s
	  EFStif(2,4)=-kh*s*s-kv*c*c
      
	  EFStif(3,2)=EFStif(2,3)
      EFStif(4,2)=EFStif(2,4)
	  
	  EFStif(3,3)=kh*c*c+kv*s*s
	  EFStif(3,4)=kh*c*s-kv*c*s
	  EFStif(4,3)=EFStif(3,4)
	  EFStif(4,4)=kh*s*s+kv*c*c
   end subroutine  shape_flet_stif
    
   !============================================================================
   subroutine kupf( Fcu,CYL,Fz)
   !============================================================================
      real(kind=rkind),intent(in):: Fcu,CYL(:,:)
	  real(kind=rkind),intent(out)::Fz(:,:)
	  real(kind=rkind)::fc,ft,s1,s2,alf
	  
	  fc=-0.67*fcu
	  ft=0.23*fcu**0.6667
	  s1=CYL(4,1)
	  s2=CYL(5,1)
	  alf=s1/s2
	  if (s1<0 .and. s2<0 .and. alf<=1 .and. alf>=0) then
	     fz(2,1)=(1+3.65*alf)*fc/(1+alf)**2
		 fz(1,1)=alf*fz(2,1)
	  elseif(s1>0 .and. s2<=0) then
	     if (alf<=0 .and. alf>=-0.17) then
		    fz(2,1)=(1+3.28*alf)*fc/(1+alf)**2
			fz(1,1)=alf*fz(2,1)
		 else
		    fz(1,1)=ft
			fz(2,1)=0.65*fc
		 end if
	  elseif(s1>0 .and.s2>0 .and. alf>=1) then
	     fz(1,1)=ft
		 fz(2,1)=ft
	  end if
   end subroutine

   !============================================================================
   subroutine saenz(Ki,Fz,Fc,sic,siu,E,R_CYL,Eo,YL)
   !============================================================================
      real(kind=rkind),intent(in out)::Fz(:,:),sic(:,:),siu(:,:),e(:,:),R_CYL(:,:),fc,Eo,YL(:)
	  integer(kind=skind),intent(in)::ki
	  real(kind=rkind)::T,ep,es,q
      integer(kind=skind)::i
!fc为负sic也为负
	  do i=2,ki,-1
	     t=Fz(i,1)/fc
		 ep=-0.0033
		 if (abs(fz(i,1))>abs(fc)) then
		    sic(i,1)=ep*(3*t-2)
		 else
		    sic(i,1)=ep*(-1.6*t**3+2.25*t**2+0.35*t)
		 end if  !得到对于极限应力的的等效单向极限应变
		 if (abs(siu(i,1))<abs(sic(i,1))) then  !如果累加应变小于极限应变
		    Es=abs(fz(i,1)/sic(i,1))
            q=siu(i,1)/sic(i,1)
			if (q<0) then
			  continue
			  !write(*,*) "Fine error"
			   siu(i,1)=2*yl(i+3)/(eo+e(i,1))
			   !if (siu(i,1)<0.0) siu(i,1)=0.0
			   q=siu(i,1)/sic(i,1)
			  
			   
			end if
			E(i,1)=Eo*(1.0-q**2)/((1.0+(Eo/es-2.0)*q+q**2)**2)  ! 得到E
	        if (e(i,1)<0.0) then 
			   e(i,1)=eo
			end if
			R_CYL(i,1)=Eo*siu(i,1)/(1.0+(Eo/es-2.0)*q+q**2)  !得到真正累加应力
		 end if
	  end do   
   end subroutine saenz
  
   !============================================================================
   subroutine Varbandsolv(Kcol,G_load,Ntotv)   !解方程
   !============================================================================
      integer(kind=ikind):: i,j,row1j,row_1,Ntotv,row1
      type(typ_Kcol),intent(in out)::Kcol(:)
      real(kind=rkind),intent(in out)::G_load(0:Ntotv)
      real(kind=rkind)::   s,diag(Ntotv)  

	  diag(:)=(/(kcol(j)%row(j),j=1,Ntotv)/) 
	  do j=2,Ntotv
	     row1j=lbound(kcol(j)%row,1)
	     do i=row1j,j-1
	        row_1=max(row1j,lbound(kcol(i)%row,1))
		    s=sum(diag(row_1:i-1)*kcol(i)%row(row_1:i-1)*kcol(j)%row(row_1:i-1))
            kcol(j)%row(i)=(kcol(j)%row(i)-s)/diag(i)
	     end do
	     s=sum(diag(row1j:j-1)*kcol(j)%row(row1j:j-1)**2)
	     diag(j)=diag(j)-s
      end do 
	  do j=2,Ntotv
	     row1=lbound(kcol(j)%row,1)
	     G_load(j)=G_load(j)-sum(kcol(j)%row(row1:j-1)*G_load(row1:j-1))
	  end do
	  do j=1,Ntotv
	     G_load(j)=G_load(j)/diag(j)
	  end do
	  do j=Ntotv,1,-1
	     row1=lbound(kcol(j)%row,1)
	     G_load(row1:j-1)=G_load(row1:j-1)-G_load(j)*kcol(j)%row(row1:j-1)
	  end do
   end subroutine Varbandsolv
   

end module CoreCalculation