⭐ 欢迎来到虫虫下载站！ | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站
📤 上传资源

📄 atf1.f90

📁 薄厚通用的板单元
💻 F90
📖 第 1 页 / 共 5 页
字号:
🔍
12345 
  DV=W(igaus)*dl*0.5
  CALL Shape_Function_T(nel,L,coren)
  DO inode=1,NNODE(IELES(nel,MNODE+1))
   SELECT CASE(nconp)
    CASE(1)  
	 rvece(3*inode-2)=rvece(3*inode-2)                        &
	  +(Rx(1)*L(1)+Rx(2)*L(2))*SHAPE_N_matrix(1,3*inode-2)*DV &
	  +(Ry(1)*L(1)+Ry(2)*L(2))*SHAPE_N_matrix(2,3*inode-2)*DV 
	 rvece(3*inode-1)=rvece(3*inode-1)            &
	  +(Rx(1)*L(1)+Rx(2)*L(2))*SHAPE_N_matrix(1,3*inode-1)*DV &
	  +(Ry(1)*L(1)+Ry(2)*L(2))*SHAPE_N_matrix(2,3*inode-1)*DV 
	 rvece(3*inode)=rvece(3*inode)                &
	  +(Rx(1)*L(1)+Rx(2)*L(2))*SHAPE_N_matrix(1,3*inode)*DV   &
	  +(Ry(1)*L(1)+Ry(2)*L(2))*SHAPE_N_matrix(2,3*inode)*DV 
    CASE(2)
	 rvece(3*inode-2)=rvece(3*inode-2)                        &
	  +(Rx(1)*L(2)+Rx(2)*L(3))*SHAPE_N_matrix(1,3*inode-2)*DV &
	  +(Ry(1)*L(2)+Ry(2)*L(3))*SHAPE_N_matrix(2,3*inode-2)*DV 
	 rvece(3*inode-1)=rvece(3*inode-1)                        &
	  +(Rx(1)*L(2)+Rx(2)*L(3))*SHAPE_N_matrix(1,3*inode-1)*DV &
	  +(Ry(1)*L(2)+Ry(2)*L(3))*SHAPE_N_matrix(2,3*inode-1)*DV 
	 rvece(3*inode)=rvece(3*inode)                            &
	  +(Rx(1)*L(2)+Rx(2)*L(3))*SHAPE_N_matrix(1,3*inode)*DV   &
	  +(Ry(1)*L(2)+Ry(2)*L(3))*SHAPE_N_matrix(2,3*inode)*DV 
    CASE(3)  
	 rvece(3*inode-2)=rvece(3*inode-2)                        &
	  +(Rx(1)*L(3)+Rx(2)*L(1))*SHAPE_N_matrix(1,3*inode-2)*DV &
	  +(Ry(1)*L(3)+Ry(2)*L(1))*SHAPE_N_matrix(2,3*inode-2)*DV 
	 rvece(3*inode-1)=rvece(3*inode-1)                        &
	  +(Rx(1)*L(3)+Rx(2)*L(1))*SHAPE_N_matrix(1,3*inode-1)*DV &
	  +(Ry(1)*L(3)+Ry(2)*L(1))*SHAPE_N_matrix(2,3*inode-1)*DV 
     rvece(3*inode)=rvece(3*inode)                            &
      +(Rx(1)*L(3)+Rx(2)*L(1))*SHAPE_N_matrix(1,3*inode)*DV   &
      +(Ry(1)*L(3)+Ry(2)*L(1))*SHAPE_N_matrix(2,3*inode)*DV 
    END SELECT
  END DO
  DEALLOCATE(SHAPE_N_matrix)
 END DO
   
 DO ievab=1,MEVAB
  eload(nel,ievab)=eload(nel,ievab)+rvece(ievab)
 END DO
  
END SUBROUTINE LOADD_2D_PLANE_METHOD_3

!*============================================================================*!
!*S16. LOADG   功能: 二维平面单元（无旋转自由度）重力荷载                     *!
!*============================================================================*!

SUBROUTINE LOADG(eload)
 USE M0_PROBLEM_TYPE
 USE M1_Gobal_information
 USE M2_Local_information
 USE M3_Stru_information
 USE M4_Support_information
 USE M5_Para_Gauss_Integral
 USE M6_Matrix_Element
 USE M8_Shape_Function
 REAL*8, DIMENSION(NELEM,MEVAB) :: eload
 REAL*8,DIMENSION(:),ALLOCATABLE:: rvece
 REAL*8, DIMENSION(:,:),ALLOCATABLE:: coren
 REAL*8, DIMENSION(3) :: L
 REAL*8:: dense,djacb,DV,gx,gy,thick,theta,gravy,XX,YY
 REAL*8:: GS1,GS2,GS3,GS4,A1,B3,C3
 INTEGER:: nel,nty,nic,kg,j,m
 INTEGER:: idime,ielem,igaus,ievab,inode,jgaus

 ALLOCATE(rvece(MEVAB),coren(NDIME,MNODE))
 WRITE(3,'(A)') ' -----------------------------------------------------------&
-----------'

 READ(1,*) theta, gravy
 WRITE(3,*) 'theta=', theta, '       gravy=',gravy
 WRITE(3,'(A)') ' -----------------------------------------------------------&
-----------'
 DO ielem=1,NELEM
  nel=IELES(ielem,MNODE+1)
  nty=IELES(ielem,MNODE+2)
  thick=PROPS(nty,3)
  dense=PROPS(nty,4)
  IF(dense==0.0) CYCLE
  gx= dense*gravy*SIN(theta)
  gy=-dense*gravy*COS(theta)
  rvece=0.0
  DO inode=1,MNODE
   nic=IELES(ielem,inode)
   DO idime=1,NDIME
    coren(idime,inode)=COORD(nic,idime)
   END DO
  END DO

  SELECT CASE(ELESP(nel))
   CASE(4)
    CALL GAUSS_quadrilateral(NGAUS(nel))
    DO igaus=1,NGAUS(nel)
     DO jgaus=1,NGAUS(nel)
	  XX=POSGP(igaus)
	  YY=POSGP(jgaus)
      CALL Shape_Function_Q(ielem,XX,YY,coren)
	  kg=1
	  SELECT CASE(ELMOD(nel))
       CASE(1:2)
	    CALL JACOBI_2or3_DIMENSION(ielem,djacb,kg,coren)
	   CASE(3:4)
	    CALL JACOBI_FOR_AREA_COODINATE(ielem,djacb,kg,coren)
	  END SELECT
	  DV=djacb*WEIGP(igaus)*WEIGP(jgaus)
	  IF(thick/=0.0) DV=DV*thick
	  DO inode=1,NNODE(nel)
	   rvece(2*inode-1)=rvece(2*inode-1)+gx*SHAPE_N(inode)*DV
	   rvece(2*inode)  =rvece(2*inode)  +gy*SHAPE_N(inode)*DV
	  END DO
      DEALLOCATE(SHAPE_N,DERIV_N,CORGP,CARTD) 
     END DO
    END DO
   CASE(3)
    CALL GAUSS_triangular(NGAUS(nel))
    GS1=POSGP(1);GS2=POSGP(2);GS3=POSGP(3);GS4=POSGP(4)
    A1=WEIGP(1) ; B3=WEIGP(2) ; C3=WEIGP(3)
    kg=1

    IF(A1/=0.0) THEN
     L=1.0/3.0
     CALL Shape_Function_T(ielem,L,coren)
     CALL JACOBI_2or3_DIMENSION(ielem,djacb,kg,coren)
     DV=djacb*A1
     IF(thick/=0.0) DV=DV*thick
     DO inode=1,NNODE(nel)
      rvece(2*inode-1)=rvece(2*inode-1)+gx*SHAPE_N(inode)*DV
      rvece(2*inode)  =rvece(2*inode)  +gy*SHAPE_N(inode)*DV
     END DO
     DEALLOCATE(SHAPE_N,DERIV_N,CORGP,CARTD) 
	END IF

	IF(B3/=0.0) THEN
	 DO igaus=1, 3
      SELECT CASE(igaus)
       CASE(1)
	    L(1)=GS1; L(2)=GS1; L(3)=GS2
       CASE(2)
	    L(1)=GS1; L(2)=GS2; L(3)=GS1
       CASE(3)
	    L(1)=GS2; L(2)=GS1; L(3)=GS1
      END SELECT
      CALL Shape_Function_T(ielem,L,coren)
      CALL JACOBI_2or3_DIMENSION(ielem,djacb,kg,coren)
      DV=djacb*B3
      IF(thick/=0.0) DV=DV*thick
      DO inode=1,NNODE(nel)
       rvece(2*inode-1)=rvece(2*inode-1)+gx*SHAPE_N(inode)*DV
       rvece(2*inode)  =rvece(2*inode)  +gy*SHAPE_N(inode)*DV
      END DO
      DEALLOCATE(SHAPE_N,DERIV_N,CORGP,CARTD)
	 END DO 
	END IF

	IF(C3/=0.0) THEN
	 DO igaus=1, 3
      SELECT CASE(igaus)
       CASE(1)
        L(1)=GS3; L(2)=GS3; L(3)=GS4
       CASE(2)
	    L(1)=GS3; L(2)=GS4; L(3)=GS3
	   CASE(3)
	    L(1)=GS4; L(2)=GS3; L(3)=GS3
      END SELECT
      CALL Shape_Function_T(ielem,L,coren)
      CALL JACOBI_2or3_DIMENSION(ielem,djacb,kg,coren)
      DV=djacb*C3
      IF(thick/=0.0) DV=DV*thick
      DO inode=1,NNODE(nel)
       rvece(2*inode-1)=rvece(2*inode-1)+gx*SHAPE_N(inode)*DV
       rvece(2*inode)  =rvece(2*inode)  +gy*SHAPE_N(inode)*DV
      END DO
      DEALLOCATE(SHAPE_N,DERIV_N,CORGP,CARTD)
	 END DO 
	END IF
    
  END SELECT

  DO ievab=1,MEVAB
   eload(ielem,ievab)=eload(ielem,ievab)+rvece(ievab)
  END DO
  DEALLOCATE(POSGP,WEIGP)
 END DO

 DEALLOCATE(rvece,coren)

END SUBROUTINE LOADG

!*============================================================================*!
!*S17. LOADD_2D_Plate   功能: 二维板弯曲单元均布荷载                          *!
!*============================================================================*!

SUBROUTINE LOADD_2D_Plate(eload)
 USE M0_PROBLEM_TYPE
 USE M1_Gobal_information
 USE M2_Local_information
 USE M3_Stru_information
 USE M5_Para_Gauss_Integral
 USE M6_Matrix_Element
 USE M8_Shape_Function
 REAL*8, DIMENSION(NELEM,MEVAB):: eload
 INTEGER,DIMENSION(:),   ALLOCATABLE:: noprs 
 INTEGER::ipart,ielem 
 INTEGER::ntp,nelef,nelel,nqpart
 REAL*8 ::q

 READ(1,*) nqpart
 WRITE(3,*) ' *There are ',nqpart,'   element groups under uniform loadings.'
 IF(nqpart>0) THEN
 WRITE(3,'(A)') ' -----------------------------------------------------------&
-----------'
  DO ipart=1,nqpart 
   READ(1,*) nelef,nelel,q
   WRITE(3,*) '  *The No.', ipart, '   uniform loading from element',nelef,'   to',nelel
   WRITE(3,*) '    and the loading q=',q
   DO ielem=nelef, nelel
    ntp=IELES(ielem,MNODE+1)
     SELECT CASE(ELESP(ntp))
      CASE(3)
       CALL LOADD_2D_PLATE_T(ielem,q,eload)
	  CASE(4)
       CALL LOADD_2D_PLATE_Q(ielem,q,eload)
	 END SELECT
   END DO
  END DO
 END IF
 WRITE(3,'(A)') ' -----------------------------------------------------------&
-----------'

END SUBROUTINE LOADD_2D_Plate

!*============================================================================*!
!*S18. LOADD_2D_PLATE_T   功能: 二维三角形板单元均布荷载                      *!
!*============================================================================*!

SUBROUTINE LOADD_2D_PLATE_T(nel,q,eload)
 USE M0_PROBLEM_TYPE
 USE M1_Gobal_information
 USE M2_Local_information
 USE M3_Stru_information
 USE M5_Para_Gauss_Integral
 USE M6_Matrix_Element
 USE M8_Shape_Function
 REAL*8, DIMENSION(MEVAB)      :: rvece
 REAL*8, DIMENSION(NDIME,MNODE):: coren
 REAL*8, DIMENSION(NELEM,MEVAB):: eload
 REAL*8,DIMENSION(3):: L 
 REAL*8 :: q,A1,B3,C3,GS1,GS2,GS3,GS4,DV,djacb 
 INTEGER:: kg,nel,etp,nic
 INTEGER:: idime,inode,igaus,ievab

 etp=IELES(nel,MNODE+1)
 
 DO inode=1,MNODE
  nic=IELES(nel,inode)
  DO idime=1,NDIME
   IF(nic/=0) coren(idime,inode)=COORD(nic,idime)
  END DO
 END DO
 
 CALL GAUSS_triangular(NGAUS(etp))
 GS1=POSGP(1);GS2=POSGP(2);GS3=POSGP(3);GS4=POSGP(4)
 A1=WEIGP(1) ; B3=WEIGP(2) ; C3=WEIGP(3)
 rvece=0.0; kg=0  

 IF(A1/=0.0) THEN
  L=1.0/3.0
  kg=kg+1
  CALL Shape_Function_T(nel,L,coren)
  CALL JACOBI_2or3_DIMENSION(nel,djacb,kg,coren)
  DV=0.5*djacb*A1
  rvece=rvece+DV*q*SHAPE_N
  DEALLOCATE(SHAPE_N,DERIV_N,CORGP,CARTD) 
 END IF

 IF(B3/=0.0) THEN       
  DO igaus=1,3
   kg=kg+1
   SELECT CASE(igaus)
    CASE(1)
	 L(1)=GS1; L(2)=GS1; L(3)=GS2
    CASE(2)
	 L(1)=GS1; L(2)=GS2; L(3)=GS1
	CASE(3)
	 L(1)=GS2; L(2)=GS1; L(3)=GS1
   END SELECT
   CALL Shape_Function_T(nel,L,coren)
   CALL JACOBI_2or3_DIMENSION(nel,djacb,kg,coren)
   DV=0.5*djacb*B3
   rvece=rvece+DV*q*SHAPE_N
   DEALLOCATE(SHAPE_N,DERIV_N,CORGP,CARTD) 
  END DO
 END IF
 
 IF(C3/=0.0) THEN
  DO igaus=1,3
   kg=kg+1
   SELECT CASE(igaus)
    CASE(1)
     L(1)=GS3; L(2)=GS3; L(3)=GS4
	CASE(2)
     L(1)=GS3; L(2)=GS4; L(3)=GS3
	CASE(3)
     L(1)=GS4; L(2)=GS3; L(3)=GS3
   END SELECT
   CALL Shape_Function_T(nel,L,coren)
   CALL JACOBI_2or3_DIMENSION(nel,djacb,kg,coren)
   DV=0.5*djacb*C3
   rvece=rvece+DV*q*SHAPE_N
   DEALLOCATE(SHAPE_N,DERIV_N,CORGP,CARTD) 
  END DO
 END IF

 DO ievab=1,NEVAB(etp)
  eload(nel,ievab)=eload(nel,ievab)+rvece(ievab)
 END DO
	
 DEALLOCATE(POSGP,WEIGP)

END SUBROUTINE LOADD_2D_PLATE_T

!*============================================================================*!
!*S19. LOADD_2D_PLATE_Q   功能: 二维四边形板单元均布荷载                      *!
!*============================================================================*!

SUBROUTINE LOADD_2D_PLATE_Q(nel,q,eload)
 USE M0_PROBLEM_TYPE
 USE M1_Gobal_information
 USE M2_Local_information
 USE M3_Stru_information
 USE M5_Para_Gauss_Integral
 USE M6_Matrix_Element
 USE M8_Shape_Function
 REAL*8, DIMENSION(:),ALLOCATABLE:: rvece  !(MEVAB) 
 REAL*8, DIMENSION(NDIME,MNODE):: coren
 REAL*8, DIMENSION(NELEM,MEVAB):: eload
 REAL*8 :: q,S,T,djacb
 INTEGER:: kg,nel,nic,etp
 INTEGER:: inode,idime,igaus,jgaus,ievab,jdofn

 etp=IELES(nel,MNODE+1) 
 ALLOCATE(rvece(NEVAB(etp)))

 DO inode=1,MNODE
  nic=IELES(nel,inode)
  DO idime=1,NDIME
   IF(nic/=0) coren(idime,inode)=COORD(nic,idime)
  END DO
 END DO

 CALL GAUSS_quadrilateral(NGAUS(etp))
 rvece=0.0; kg=0  

 DO igaus=1,NGAUS(etp)
  DO jgaus=1,NGAUS(etp)
   kg=kg+1
   S=POSGP(igaus);T=POSGP(jgaus)
   CALL Shape_Function_Q(nel,S,T,coren)
   CALL JACOBI_2or3_DIMENSION(nel,djacb,kg,coren)
     rvece=rvece+djacb*WEIGP(igaus)*WEIGP(jgaus)*q*SHAPE_N
   DEALLOCATE(SHAPE_N,DERIV_N,CORGP,CARTD)
  END DO
 END DO

 DO ievab=1,NEVAB(etp)
  SELECT CASE(TYPE2)
   CASE(3)
    eload(nel,ievab)=eload(nel,ievab)+rvece(ievab)
   CASE(4)
    eload(nel,ievab)=eload(nel,ievab)+rvece(ievab)
  END SELECT
 END DO

 DEALLOCATE(POSGP,WEIGP,rvece)

END SUBROUTINE LOADD_2D_PLATE_Q

!*============================================================================*!
!*S20. LOADD_2D_Composite_Plate   功能: 层合板单元面内边界分布与横向均布荷载  *!
!*============================================================================*!
12345

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -