atf1.f90

薄厚通用的板单元

 ALLOCATE(FIXED(NTOTV))
 FIXED=0.0
 
 DO ivfix=1, NVFIX
  nloca=(nofix(ivfix)-1)*NDOFN

  DO idofn=1,NDOFN
   ngash=nloca+idofn
   fixed(ngash)=PRESC(ivfix,idofn)
  END DO

 END DO
 
 DEALLOCATE(PRESC)
 WRITE(3,*)

!(6).读写单元材料参数

 WRITE(3,*) '*<5>.ELEMENT MATERIAL PROPERTIES'
 IF(TYPE1==1) THEN
  ALLOCATE(PROPS(NMATS, NPROP))
  PROPS=0.0

  DO imats=1,NMATS
   READ(1,*) nty,(PROPS(imats,jnmat),jnmat=1,NPROP)
   IF(TYPE2/=4) THEN
    WRITE(3,FMT='(1X,I3,E11.3,4F10.3,3F10.3)') nty, &
      (PROPS(imats,jnmat),jnmat=1,NPROP)
   ELSE
    WRITE(3,FMT='(1X,I4,2E11.3,7F10.3)')nty, &
	  (PROPS(imats,jnmat),jnmat=1,NPROP)
   END IF
  END DO
  IF(TYPE2==4) THEN
   READ(1,*) Num_Layer, K1,K2
   ALLOCATE(Layer_type(Num_Layer),Inform_Layer(Num_Layer,2))
   WRITE(3,*) '** THE INFORMATION FOR EACH LAYER OF COMPOSITE PLATE'
   DO ilayer=1,Num_Layer
    READ(1,*) Layer_type(ilayer),(Inform_Layer(ilayer, jnmat),jnmat=1,2) 
    WRITE(3,FMT='(1X,2I4, F10.3, E11.3)') ilayer,Layer_type(ilayer),&
			  (Inform_Layer(ilayer, jnmat),jnmat=1,2)
   END DO
   CALL COMPOSITE_PLATE_MATRIX
  END IF
 END IF
 WRITE(3,*)

END SUBROUTINE Input_Data

!*============================================================================*!
!*S3. GCVAR_A_1    功能: 给出平面膜元单元常数                                 *!
!*============================================================================*!

SUBROUTINE GCVAR_A_1(etp,n_node,n_gaus,n_evab,n_elsp,n_dofn)
 INTEGER::etp,n_node,n_gaus,n_evab,n_elsp,n_dofn

 SELECT CASE(etp)
  CASE(1)
   WRITE(5,*)' 4-node isoparametric element Q4 is used'
   WRITE(3,*)' 4-node isoparametric element Q4 is used, ELEMENT	MODEL NO.=',etp
   n_node=4; n_dofn=2; n_gaus=2; n_elsp=4; n_evab=n_dofn*n_node
  CASE(2)
   WRITE(5,*)' 8-node isoparametric element Q8 is used'
   WRITE(3,*)' 8-node isoparametric element Q8 is used, ELEMENT MODEL NO.=',etp
   n_node=8; n_dofn=2; n_gaus=3; n_elsp=4; n_evab=n_dofn*n_node
  CASE(3)
   WRITE(5,*)' 8-node Area coodinate element MQ8-I is used'
   WRITE(3,*)' 8-node Area coodinate element MQ8-I is used, ELEMENT MODEL NO.=',etp
   n_node=8; n_dofn=2; n_gaus=3; n_elsp=4; n_evab=n_dofn*n_node
  CASE(4)
   WRITE(5,*)' 8-node Area coodinate element MQ8-II is used'
   WRITE(3,*)' 8-node Area coodinate element MQ8-II is used, ELEMENT MODEL NO.=',etp
   n_node=8; n_dofn=2; n_gaus=3; n_elsp=4; n_evab=n_dofn*n_node
  CASE(61)
   WRITE(5,*)' 3-node triangular element T3 is used'
   WRITE(3,*)' 3-node triangular element T3 is used, ELEMENT MODEL NO.=',etp
   n_node=3; n_dofn=2; n_gaus=1; n_elsp=3; n_evab=n_dofn*n_node
  CASE(62)
   WRITE(5,*)' 6-node triangular element T6 is used'
   WRITE(3,*)' 6-node triangular element T6 is used, ELEMENT MODEL NO.=',etp
   n_node=6; n_dofn=2; n_gaus=3; n_elsp=3; n_evab=n_dofn*n_node
  CASE(81)
   WRITE(5,*)' Triangular element GT9 with Drilling DOF is used'
   WRITE(3,*)' Triangular element GT9 with Drilling DOF is used, ELEMENT MODEL NO.=',etp
   n_node=3; n_dofn=3; n_gaus=3; n_elsp=3; n_evab=n_dofn*n_node
  CASE(82)
   WRITE(5,*)' Triangular element GT9M8 with Drilling DOF is used'
   WRITE(3,*)' Triangular element GT9M8 with Drilling DOF is used, ELEMENT MODEL NO.=',etp
   n_node=3; n_dofn=3; n_gaus=7; n_elsp=3; n_evab=n_dofn*n_node
   !+++++++  Here you can add constants of new elements  +++++++++++
  CASE DEFAULT
   WRITE(*,*)' Sorry, we haven"t established the database of this element!'
   WRITE(5,*)' Sorry, we haven"t established the database of this element!'
   WRITE(3,*)' Sorry, we haven"t established the database of this element!'
   STOP
 END SELECT
 WRITE(5,*)

END SUBROUTINE GCVAR_A_1

!*============================================================================*!
!*S4. GCVAR_A_2    功能: 给出平板弯曲单元常数                                 *!
!*============================================================================*!

SUBROUTINE GCVAR_A_2(etp,n_node,n_gaus,n_evab,n_elsp,n_dofn)
 INTEGER::etp,n_node,n_gaus,n_evab,n_elsp,n_dofn

 SELECT CASE(etp)
  CASE(1)
   WRITE(5,*)' 3-node Triangular Mindlin-plate element TCGC-T9 is used'
   WRITE(3,*)' 3-node Triangular Mindlin-plate element TCGC-T9 is used,'
   WRITE(3,*)'  ELEMENT MODEL NO.=',etp
   n_node=3; n_dofn=3; n_gaus=4; n_elsp=3; n_evab=n_dofn*n_node
  CASE(2)
   WRITE(5,*)' 3-node Triangular Mindlin-plate element TMT is used'
   WRITE(3,*)' 3-node Triangular Mindlin-plate element TMT is used,'
   WRITE(3,*)'  ELEMENT MODEL NO.=',etp
   n_node=3; n_dofn=3; n_gaus=4; n_elsp=3; n_evab=n_dofn*n_node
  CASE(3)
   WRITE(5,*)' 3-node Triangular Mindlin-plate element TLSL-T9 is used'
   WRITE(3,*)' 3-node Triangular Mindlin-plate element TLSL-T9 is used,'
   WRITE(3,*)'  ELEMENT MODEL NO.=',etp
   n_node=3; n_dofn=3; n_gaus=4; n_elsp=3; n_evab=n_dofn*n_node
  CASE(4)
   WRITE(5,*)' 3-node Triangular Mindlin-plate element ATF-T9 is used'
   WRITE(3,*)' 3-node Triangular Mindlin-plate element ATF-T9 is used,'
   WRITE(3,*)'  ELEMENT MODEL NO.=',etp
   n_node=3; n_dofn=3; n_gaus=4; n_elsp=3; n_evab=n_dofn*n_node
   !!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!Input
  CASE(5)
   WRITE(5,*)' 3-node Triangular Mindlin-plate element ATF1-T9 is used'
   WRITE(3,*)' 3-node Triangular Mindlin-plate element ATF1-T9 is used,'
   WRITE(3,*)'  ELEMENT MODEL NO.=',etp
	n_node=3; n_dofn=3; n_gaus=4; n_elsp=3; n_evab=n_dofn*n_node
  CASE(21)
   WRITE(5,*)' 4-node Quadrilateral Mindlin-plate element TACQ is used'
   WRITE(3,*)' 4-node Quadrilateral Mindlin-plate element TACQ is used'
   WRITE(3,*)'  ELEMENT MODEL NO.=',etp
   n_node=4; n_dofn=3; n_gaus=3; n_elsp=4; n_evab=n_dofn*n_node
  CASE(22)
   WRITE(5,*)' 4-node Quadrilateral Mindlin-plate element TMQ is used'
   WRITE(3,*)' 4-node Quadrilateral Mindlin-plate element TMQ is used'
   WRITE(3,*)'  ELEMENT MODEL NO.=',etp
   n_node=4; n_dofn=3; n_gaus=3; n_elsp=4; n_evab=n_dofn*n_node
 ! CASE(11)                                                       !11板弯曲单元为自己所加
!   WRITE(5,*)' 4-node Quadrilateral Mindlin-plate element ATF_PQ4a is used'
 !  WRITE(3,*)' 4-node Quadrilateral Mindlin-plate element ATF_PQ4a is used'
 !  WRITE(3,*)'  ELEMENT MODEL NO.=',etp
 !  n_node=4; n_dofn=3; n_gaus=3; n_elsp=4; n_evab=n_dofn*n_node
!!+++++++  Here you can add constants of new elements  +++++++++++
  CASE DEFAULT
   WRITE(*,*)' Sorry, we haven"t established the database of this element!'
   WRITE(5,*)' Sorry, we haven"t established the database of this element!'
   WRITE(3,*)' Sorry, we haven"t established the database of this element!'
   STOP
 END SELECT
 WRITE(5,*)

END SUBROUTINE GCVAR_A_2

!*============================================================================*!
!*S5. GCVAR_A_3    功能: 给出层合板弯曲单元常数                               *!
!*============================================================================*!

SUBROUTINE GCVAR_A_3(etp,n_node,n_gaus,n_evab,n_elsp,n_dofn)
 INTEGER::etp,n_node,n_gaus,n_evab,n_elsp,n_dofn

 SELECT CASE(etp)
!  CASE(1)
!   WRITE(5,*)' 3-node Triangular Mindlin-plate element TCGC-T9 is used'
!   WRITE(3,*)' 3-node Triangular Mindlin-plate element TCGC-T9 is used,'
!   WRITE(3,*)'  ELEMENT MODEL NO.=',etp
!   n_node=3; n_dofn=3; n_gaus=4; n_elsp=3; n_evab=n_dofn*n_node
!  CASE(2)
!   WRITE(5,*)' 3-node Triangular Mindlin-plate element TMT is used'
!   WRITE(3,*)' 3-node Triangular Mindlin-plate element TMT is used,'
!   WRITE(3,*)'  ELEMENT MODEL NO.=',etp
!   n_node=3; n_dofn=3; n_gaus=4; n_elsp=3; n_evab=n_dofn*n_node
  CASE(1)
   WRITE(5,*)' 4-node Quadrilateral Mindlin Composite plate element CTACQH is used'
   WRITE(3,*)' 4-node Quadrilateral Mindlin Composite plate element CTACQH is used'
   WRITE(3,*)'  ELEMENT MODEL NO.=',etp
   n_node=4; n_dofn=5; n_gaus=2; n_elsp=4; n_evab=n_dofn*n_node
  CASE(2)
   WRITE(5,*)' 4-node Quadrilateral Mindlin Composite plate element CTMQH is used'
   WRITE(3,*)' 4-node Quadrilateral Mindlin Composite plate element CTMQH is used'
   WRITE(3,*)'  ELEMENT MODEL NO.=',etp
   n_node=4; n_dofn=5; n_gaus=2; n_elsp=4; n_evab=n_dofn*n_node
  CASE(101)
   WRITE(5,*)' 3-node Triangular Mindlin Composite plate element CTLSL-T9 is used'
   WRITE(3,*)' 3-node Triangular Mindlin Composite plate element CTLSL-T9 is used,'
   WRITE(3,*)'  ELEMENT MODEL NO.=',etp
   n_node=3; n_dofn=5; n_gaus=4; n_elsp=3; n_evab=n_dofn*n_node
!!+++++++  Here you can add constants of new elements  +++++++++++
  CASE DEFAULT
   WRITE(*,*)' Sorry, we haven"t established the database of this element!'
   WRITE(5,*)' Sorry, we haven"t established the database of this element!'
   WRITE(3,*)' Sorry, we haven"t established the database of this element!'
   STOP
 END SELECT
 WRITE(5,*)

END SUBROUTINE GCVAR_A_3

!*============================================================================*!
!*S6. CORMN      功能: 计算边中结点坐标                                       *!
!*============================================================================*!

SUBROUTINE CORMN
 USE M1_Gobal_information
 USE M2_Local_information
 USE M3_Stru_information
 INTEGER:: ielem,ielet,inode,jnode,k,midpt,nodfn,nodmd,nodst
 REAL*8:: total
 DO ielem=1,NELEM
   ielet=IELES(ielem, MNODE+1)
   IF(NNODE(ielet)==2*ELESP(ielet)) THEN
   DO inode=1,ELESP(ielet)
    nodst=IELES(ielem,inode)
	jnode=inode+1
	IF(jnode>ELESP(ielet)) jnode=1 
	nodfn=IELES(ielem,jnode)
	midpt=inode+ELESP(ielet)
	nodmd=IELES(ielem,midpt)
	total=ABS(COORD(nodmd,1))+ABS(COORD(nodmd,2))
	IF(total<=0.0) THEN
	  k=1
	  COORD(nodmd,k)=(COORD(nodst,k)+COORD(nodfn,k))/2.0
      K=2
	  COORD(nodmd,k)=(COORD(nodst,k)+COORD(nodfn,k))/2.0
	ENDIF
   ENDDO
   END IF
ENDDO

END SUBROUTINE CORMN

!*============================================================================*!
!*S7. COMPOSITE_PLATE_MATRIX      功能: 计算层合板偏轴弹性矩阵                *!
!*============================================================================*!

SUBROUTINE COMPOSITE_PLATE_MATRIX
 USE M3_Stru_information
 USE M21_INFORMATION_COMPOSITE_Layer
 REAL*8 :: Pi=3.141592653589793,lk,mk
 REAL*8 :: E1,E2,PU12,PU21,G12,G23,G13
 REAL*8 :: Q11,Q12,Q22,Q66,Q55,Q44
 INTEGER :: ilayer

 ALLOCATE(Q_MATRIX_PZ(Num_Layer,3,3),Q_MATRIX_ZZ(Num_Layer,3,3))
 ALLOCATE(C_MATRIX_PZ(Num_Layer,2,2))
 Q_MATRIX_PZ=0.0;Q_MATRIX_ZZ=0.0

 DO ilayer=1,Num_Layer
  E1=PROPS(Layer_type(ilayer),1); E2=PROPS(Layer_type(ilayer),2)
  PU12=PROPS(Layer_type(ilayer),3);PU21=PROPS(Layer_type(ilayer),4)
  G12=PROPS(Layer_type(ilayer),5)
  G23=PROPS(Layer_type(ilayer),6); G13=PROPS(Layer_type(ilayer),7)

  IF(PU21==0.0)   PU21=E2*PU12/E1
  Q_MATRIX_ZZ(ilayer,1,1)=E1/(1-PU12*PU21)
  Q_MATRIX_ZZ(ilayer,1,2)=PU21*E1/(1-PU12*PU21)
  Q_MATRIX_ZZ(ilayer,2,2)=E2/(1-PU12*PU21)
  Q_MATRIX_ZZ(ilayer,2,1)=Q_MATRIX_ZZ(ilayer,1,2)
  Q_MATRIX_ZZ(ilayer,3,3)=G12

  Q11=Q_MATRIX_ZZ(ilayer,1,1); Q12=Q_MATRIX_ZZ(ilayer,1,2)
  Q22=Q_MATRIX_ZZ(ilayer,2,2); Q66=Q_MATRIX_ZZ(ilayer,3,3)
  Q55=G13; Q44=G23

  lk=DCOS(Inform_Layer(ilayer,1)*PI/180.0)
  mk=DSIN(Inform_Layer(ilayer,1)*PI/180.0)

  Q_MATRIX_PZ(ilayer,1,1)=Q11*lk**4+2.*(Q12+2*Q66)*lk**2*mk**2+Q22*mk**4
  Q_MATRIX_PZ(ilayer,1,2)=(Q11+Q22-4*Q66)*lk**2*mk**2+Q12*(lk**4+mk**4)
  Q_MATRIX_PZ(ilayer,2,1)=Q_MATRIX_PZ(ilayer,1,2)
  Q_MATRIX_PZ(ilayer,2,2)=Q11*mk**4+2.*(Q12+2*Q66)*lk**2*mk**2+Q22*lk**4
  Q_MATRIX_PZ(ilayer,1,3)=(Q11-Q12-2*Q66)*lk**3*mk+(Q12-Q22+2*Q66)*lk*mk**3
  Q_MATRIX_PZ(ilayer,3,1)=Q_MATRIX_PZ(ilayer,1,3)
  Q_MATRIX_PZ(ilayer,2,3)=(Q11-Q12-2*Q66)*lk*mk**3+(Q12-Q22+2*Q66)*lk**3*mk
  Q_MATRIX_PZ(ilayer,3,2)=Q_MATRIX_PZ(ilayer,2,3)
  Q_MATRIX_PZ(ilayer,3,3)=(Q11+Q22-2*Q12-2*Q66)*lk**2*mk**2+Q66*(lk**4+mk**4)
  C_MATRIX_PZ(ilayer,2,2)=Q44*lk**2+Q55*mk**2
  C_MATRIX_PZ(ilayer,1,2)=(Q55-Q44)*lk*mk
  C_MATRIX_PZ(ilayer,2,1)=C_MATRIX_PZ(ilayer,1,2)
  C_MATRIX_PZ(ilayer,1,1)=Q44*mk**2+Q55*lk**2
 END DO  

END SUBROUTINE COMPOSITE_PLATE_MATRIX

!*============================================================================*!
!*S8. GAUSS_quadrilateral    功能: 给出四边形单元高斯积分常数                 *!
!*============================================================================*!
 
SUBROUTINE GAUSS_quadrilateral(n_gaus)
 USE M5_Para_Gauss_Integral

 INTEGER::n_gaus

 ALLOCATE(POSGP(n_gaus),WEIGP(n_gaus))
 POSGP=0.0
 WEIGP=0.0
 SELECT CASE(n_gaus)
  CASE(1)
   POSGP(1)=0.0; WEIGP(1)=2.0
  CASE(2)
   POSGP(1)=-0.577350269189626
   POSGP(2)=-POSGP(1)
   WEIGP(1)=1.0
   WEIGP(2)=1.0
  CASE(3)
   POSGP(1)=-0.774596669241483
   POSGP(2)=0.0
   POSGP(3)=-POSGP(1)
   WEIGP(1)=0.555555555555556
   WEIGP(2)=0.888888888888889
   WEIGP(3)=WEIGP(1)     
  CASE(4)
   POSGP(1)=-0.861136311594053  
   POSGP(2)=-0.339981043584856
   POSGP(3)=-POSGP(2)
   POSGP(4)=-POSGP(1)
   WEIGP(1)=0.347854845137454
   WEIGP(2)=0.652145154862546
   WEIGP(3)=WEIGP(2)     
   WEIGP(4)=WEIGP(1)
  CASE(5)
   POSGP(1)=-0.906179845938664  
   POSGP(2)=-0.538469310105683
   POSGP(3)=0.0
   POSGP(4)=-POSGP(2)
   POSGP(5)=-POSGP(1)
   WEIGP(1)=0.236926885056189
   WEIGP(2)=0.478628670499366
   WEIGP(3)=0.568888888888889
   WEIGP(4)=WEIGP(2)     
   WEIGP(5)=WEIGP(1)
  CASE(6)
   POSGP(1)=-0.932469514203152  
   POSGP(2)=-0.661209386466265
   POSGP(3)=-0.238619186083197
   POSGP(4)=-POSGP(3)
   POSGP(5)=-POSGP(2)
   POSGP(6)=-POSGP(1)
   WEIGP(1)=0.171324492379170
   WEIGP(2)=0.360761573048139
   WEIGP(3)=0.467913934572691
   WEIGP(4)=WEIGP(3)     
   WEIGP(5)=WEIGP(2)
   WEIGP(6)=WEIGP(1)
  CASE DEFAULT
   WRITE(*,*)' Error: you have type a wrong gauss point number!'
   WRITE(3,*)' Error: you have type a wrong gauss point number!'
   STOP
 END SELECT
   
END SUBROUTINE GAUSS_quadrilateral

!*============================================================================*!
!*S9. GAUSS_triangular    功能: 给出三角形单元高斯积分常数                    *!
!*============================================================================*!

SUBROUTINE GAUSS_triangular(n_gaus)
 USE M5_Para_Gauss_Integral
 INTEGER:: n_gaus

 ALLOCATE(POSGP(4),WEIGP(3))
 POSGP=0.0
 WEIGP=0.0
 SELECT CASE(n_gaus)
  CASE(1)
   WEIGP(1)=1.0
  CASE(3)
   POSGP(1)=0.5
   WEIGP(2)=1.0/3.0
  CASE(4)
   POSGP(1)=0.2
   POSGP(2)=0.6
   WEIGP(1)=-27.0/48.0
   WEIGP(2)=25.0/48.0
  CASE DEFAULT
   POSGP(1)=0.1012865037