Description  ：

    This is a self-contained computer code—Fp_DBIT.F90 (Fracture Parameters evaluation based on Dual BEM and I-Integral method) developed for evaluating stress intensity factors and T-stress for mixed-mode cracks of elastic solids.   

    This code is written in Fortran 90 (keeping to Fortran77 style make it easy to understand) and consists of the main program and six subroutines. This code can evaluate the fracture parameters given that the elastic variables, displacement gradients and stresses, has been solved by the Dual BEM. Fig.1 shows the computation model of combined Dual BEM and I-integral method.

Input File: (FPPOST.dat)

The Main program Fp_DBIT inputs the coordinates of nodes which constructs the path of I-integral, as shown in Fig.2. The detailed description of the input data are listed below:

Card Set 1: Information for the No. of crack-tip and the No. of the I-integral path;

Card Set 2:

G: shear modulus;

PR: Poisson’s ratio;

NSIG: 3 for plane stress, or 4 for plane strain;

CDTIP(i), i=x or y: The i-coordinate of the crack-tip in global coordinate system;

COSL(i,1:2), i =x or y : The i-th coordinate axis of the local Cartesian coordinate system.

Card Set 3:

NJP: The number of nodes those construct the I-integral path as shown in Fig.2;

NJE: The number of boundary elements those construct the I-integral path as shown in Fig.2;

NODJ: The number of nodes within one I-integral elements;

NGSSJ: The number of Gauss points used for do Gauss quadrature of the I-integral elements.

Card Set 4:

ID: The No. of points used in the Fp_DBIT;

IP: The original ID of the point which is input from Dual_BEM program;

CD(i, ID), i=x or y: The i-coordinate for the location of point ID in global coordinates;

UGD(i, j, ID), i, j=x, y: The displacement gradient of the point with the No. of ID.

SIG(i, j, ID), i, j=x, y: The stress for the point with the No. of ID.

Card Set 5:

LNDJ(1, IE): Global node number for the 1st node of element IE;

LNDJ(2, IE): Global node number for the 2nd node of element IE;

LNDJ(2, IE): Global node number for the 3rd node of element IE;

Output File: (RESULTS.out)

FP(1): Stress intensity factor KI for pure I mode;

FP(2): Stress intensity factor KII for pure II mode;

FP(3): T-stress.

Source File：

PROGRAM Fp_DBIT
      IMPLICIT REAL*8 (A-H,O-Z)
      ALLOCATABLE:: IDP(:),CD(:,:),UGD(:,:,:),SIG(:,:,:),LND(:,:)
      DIMENSION CDTIP(2),COSL(2,2),CKL(2,3),UGDL(2,2,3),SIGL(2,2,3),     &
     &          POSGJ(4),WEITJ(4),VI(3),VI0(3),FP(3)
      OPEN(5,FILE='FPPOST.DAT',STATUS='OLD')
      OPEN(7,FILE='FPRESULTS.OUT',STATUS='UNKNOWN')
      READ(5,*); READ(5,*); READ(5,*)
      READ(5,*) G,PR,NS,(CDTIP(I),I=1,2),((COSL(I,J),J=1,2),I=1,2)
      READ(5,*); READ(5,*)
      READ(5,*) NJP,NJE,NODJ,NGSSJ; READ(5,*); READ(5,*)
      ALLOCATE(IDP(NJP),CD(2,NJP),UGD(2,2,NJP),SIG(2,2,NJP),LND(NODJ,NJE))
      DO IP=1,NJP
        READ(5,*)ID,IDP(IP),(CD(I,IP),I=1,2),((UGD(I,J,IP),J=1,2),I=1,2),&
     &         (( SIG(I,J,IP),J=1,2),I=1,2)
      ENDDO
      READ(5,*); DO IE=1,NJE; READ(5,*) (LND(ID,IE),ID=1,NODJ); ENDDO
      IF(NS.EQ.4)THEN; ALFA=2.*G/(1.-PR); ELSEIF(NS.EQ.3)THEN
      ALFA=2.*G*(1+PR); ENDIF
      CALL GAUSSVI(NGSSJ,POSGJ,WEITJ)   ! FOR INTEGRATION FOR r
      VI(1:3)=0.
      DO 100 IE=1,NJE
      CALL LOCAL_VAB(NODJ,CDTIP,COSL,CD,UGD,SIG,LND(:,IE),CKL,UGDL,SIGL)
      DO 100 IG = 1 , NGSSJ
      CALL EVAL_IKE(NODJ,CKL,UGDL,SIGL,POSGJ(IG),G,PR,NS,FJCB,VI0)
100   VI(1:3)=VI(1:3)+VI0(1:3)*FJCB*WEITJ(IG)
      FP(1:2)=0.5*ALFA*VI(1:2); FP(3)=ALFA*VI(3)
      WRITE(*,'(3(A8,F12.6))') 'KI:',FP(1),'KII:',FP(2),'T:',FP(3)
      WRITE(7,'(3(A8,F12.6))') 'KI:',FP(1),'KII:',FP(2),'T:',FP(3)
      END      

      SUBROUTINE LOCAL_VAB(NODJ,CDTIP,COSL,CD,UGD,SIG,LND,CKL,UGDL,SIGL)
      IMPLICIT REAL*8 (A-H,O-Z)      
      DIMENSION CDTIP(2),COSL(2,2),LND(NODJ),CD(2,*),UGD(2,2,*),SIG(2,2,*)
      DIMENSION CKL(2,3),UGDL(2,2,3),SIGL(2,2,3)
      DO 100 ID = 1, NODJ
        IP=LND(ID); CKL(:,ID)=0.; UGDL(:,:,ID)=0.; SIGL(:,:,ID)=0.
        DO K=1,2
        DO L=1,2
          CKL(K,ID)=CKL(K,ID)+COSL(K,L)*(CD(L,IP)-CDTIP(L))
          DO I=1,2
          DO J=1,2
            UGDL(K,L,ID)=UGDL(K,L,ID)+COSL(K,I)*UGD(I,J,IP)*COSL(L,J)
            SIGL(K,L,ID)=SIGL(K,L,ID)+COSL(K,I)*SIG(I,J,IP)*COSL(L,J)
          ENDDO
          ENDDO
        ENDDO
        ENDDO
100   CONTINUE
      ENDSUBROUTINE
!      
      SUBROUTINE GAUSSVI(NGAUS,GP,GW)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION GP(NGAUS),GW(NGAUS)
      SELECT CASE(NGAUS)
        CASE(2)
        GP(1)=-0.57735026918962D0; GW(1)=1.D0
        GP(2)= 0.57735026918962D0; GW(2)=1.D0
        CASE(3)
        GP(1)=-0.77459666924148D0; GP(2)=0.D0
        GP(3)= 0.77459666924148D0
        GW(1)= 0.55555555555556D0; GW(2)=0.88888888888888D0
        GW(3)= 0.55555555555556D0
        CASE(4)
        GP(1)=-0.86113631159405; GP(2)=-0.33998104358486
        GP(3)= 0.33998104358486; GP(4)= 0.86113631159405
        GW(1)= 0.34785484513745; GW(2)= 0.65214515486255
        GW(3)= 0.65214515486255; GW(4)= 0.34785484513745
      ENDSELECT
      ENDSUBROUTINE
      
      
      SUBROUTINE SHAPI(NODE,SP,X,CK,CP,RI,R2)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION SP(NODE),CK(2,*),CP(2),RI(2),X(*)
        ! 2 noded-element (line element)
      SP(1)=0.5D0*(1.D0-X(1)); SP(2)=0.5D0*(1.D0+X(1))
      IF(NODE.EQ.2) GOTO 50
        ! 3 noded-element (line element)
      SP(1)=-X(1)*SP(1); SP(2)=X(1)*SP(2); SP(3)=1.D0-X(1)*X(1)
      GOTO 50
        ! Calculate r and its vector components
  50  R2=0.D0; DO 70 I=1,2; RI(I)=-CP(I)
      DO ID=1,NODE; RI(I)=RI(I)+SP(ID)*CK(I,ID); ENDDO
  70  R2=R2+RI(I)*RI(I)
      ENDSUBROUTINE
      
      SUBROUTINE DSHAPI(NODE,X,CK,COSB,FJCB,DN,GD,IFCOS)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION X(*),CK(2,*),DN(2,NODE),GD(2,2),COSB(2),GR(2)
      DN(1,1)=-0.5D0; DN(1,2)=0.5D0       ! For 2D, 2 noded-element
      IF(NODE.EQ.2) GOTO 30
      DN(1,1)=-0.5D0*(1.D0-2.D0*X(1))     ! For 2D, 3 noded-element
      DN(1,2)=0.5D0*(1.D0+2.D0*X(1)); DN(1,3)=-2.D0*X(1)
  30  DO 50 I=1,2; DO 50 J=1,1; GD(I,J)=0.D0; DO 50 ID=1,NODE
  50  GD(I,J)=GD(I,J)+DN(J,ID)*CK(I,ID)
      GR(1)=GD(2,1); GR(2)=-GD(1,1)                  ! For 2D normals
      FJCB=DSQRT(DOT_PRODUCT(GD(1:2,1),GD(1:2,1))) ! LINE JACOBIAN
      IF(IFCOS.EQ.0) RETURN
      COSB(1:2)=GR(1:2)/FJCB                   ! Normal Cosine
      ENDSUBROUTINE
!
      
      SUBROUTINE EVAL_IKE(NODJ,CKL,UGDL,SIGL,POSGJ,G,PR,NS,FJCB,VI0)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION CKL(2,3),UGDL(2,2,3),SIGL(2,2,3),POSGJ(4),SHAP(NODJ),      &
     &  CP0(2),RI(2),COSB(2),DN(2,NODJ),GD(2,2),ACTDU(2,2),ACTSIG(2,2),    &
     &  auxdu(2,2),auxsig(2,2)
      DIMENSION TRAC1(2),TRAC2(2),VI0(3)
      CP0=0.; PI=3.141592654d0
      CALL SHAPI(NODJ,SHAP,POSGJ,CKL,CP0,RI,R2)
      CALL DSHAPI(NODJ,POSGJ,CKL,COSB,FJCB,DN,GD,1)
      IF(RI(1).LT.0.)THEN; IF(RI(2).LT.0.)THEN
      THETA=-(PI-DATAN(RI(2)/RI(1)))
      ELSEIF(RI(2).GT.0.)THEN; THETA=PI+DATAN(RI(2)/RI(1)); ENDIF
      ELSE; THETA=DATAN(RI(2)/RI(1))
      ENDIF; R=DSQRT(R2)
      IF(NS.EQ.3)THEN; AK=(3.D0-PR)/(1.D0+PR); ELSEIF(NS.EQ.4)THEN
      AK=3.D0-4.D0*PR; ENDIF
      ACTDU(:,:)=0.; ACTSIG(:,:)=0.
      DO ID=1,NODJ
        ACTSIG(:,:) = ACTSIG(:,:) + SHAP(ID)*SIGL(:,:,ID) ! Stress
        ACTDU(:,:) = ACTDU(:,:) + SHAP(ID)*UGDL(:,:,ID) ! Displacement gradient
      ENDDO
      DO 50 IMOD=1,3
      CALL AUXI_SOLUTION(G,PI,PR,AK,R,THETA,auxdu,auxsig,IMOD)        ! Auxiliary solution for mixed-mode SIF
      SEGY = 0.
      DO I = 1,2
        DO J=1,2; SEGY = SEGY + ACTSIG(I,J) * auxdu(I,J); ENDDO
        TRAC1(I)=DOT_PRODUCT(ACTSIG(I,1:2),COSB(1:2))
        TRAC2(I)=DOT_PRODUCT(auxsig(I,1:2),COSB(1:2))
      ENDDO
      TMP1=DOT_PRODUCT(TRAC1(1:2),auxdu(1:2,1))
      TMP2=DOT_PRODUCT(TRAC2(1:2),ACTDU(1:2,1))
50    VI0(IMOD)=SEGY*COSB(1)-(TMP1+TMP2)     
      ENDSUBROUTINE

Input File for Example:

Computed Results：

     KI:    2.615375    KII:    2.403349      T:   -0.352381

File package：Fp_DBIT.zip