C C ABAQUS user subroutine for a 2D cohesive zone model, based on the paper by C Y.F. Gao and A.F. Bower, "A simple technique for avoiding convergence C problems in finite element simulations of crack nucleation and growth on C cohesive interfaces," Modelling Simul. Mater. Sci. Eng. 12, 453-463, 2004. C C The Fortran code and the user's manual are last modified on December 9, 2004. C C Y.F. Gao & A.F. Bower, Divsion of Engineering, Brown University C C=========================== SUBROUTINE UEL =================== C SUBROUTINE UEL(RHS,STIF,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, & PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME, & KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF, & LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD) C C INCLUDE 'ABA_PARAM.INC' IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION RHS(MLVARX,*),STIF(NDOFEL,NDOFEL),PROPS(*) DIMENSION SVARS(*),ENERGY(8),COORDS(MCRD,NNODE),U(NDOFEL) DIMENSION DU(MLVARX,*),V(NDOFEL),A(NDOFEL),TIME(2),PARAMS(*) DIMENSION JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*) DIMENSION PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*) C C Only relevant variables are commented. Refer to ABAQUS manual for others. C C Variables with intent(out) C RHS Residual Vector C STIF Stiffness matrix C C Variables with intent(in) C PROPS Element property array for the Xu-Needleman model C PROPS(1) SIGMA_max C PROPS(2) Delta_n C PROPS(3) Delta_t C PROPS(4) Q C PROPS(5) R C COORDS Nodal coordinate array: COORD(J,N) is jth coord of nth node C U Total accumulated DOF array. Contains accumulated displacements, C ordered as (u_i^1, u_i^2) C DU Incremental displacements, ordered as C DU(2*(N-1)+I,1) = u_i^n C C NNODE No. nodes on element. C NDOFEL No. degrees of freedom for element C NPROPS No. real valued element properties C MLVARX Dimensioning variable. C NRHS No. RHS vectors. C MCRD Largest of max value of COORDINATES parameter or active DOF <3. C C NNODE=4 Four-point element (4,3 should collapse onto 1,2 respectively) C C 4----3 C | | C 1----2 C C NDOFEL=8 Two-dimensional, d.o.f=4*2 C NPROPS=5 Five parameters passed in C MLVARX=1 C NRHS=1 C MCRD=2 (X,Y) C DIMENSION W(4),PQ(4) DIMENSION F(10),DF(10) DIMENSION GAP(2),TRACT(2),DTDG(2,2) DIMENSION RNM(2),TANGENT(2),UREL(2) DIMENSION VGAP(2),VREL(2) NINTP=2 C C The above arrays are used inside the subroutine C NINTP No. integration points C PQ(NINTP) Local coordinates of integration points C W(NINTP) Weights of integration points C F(NNODE) Shape function C DF(NNODE) Derivative of shape function w.r.t local coordinate C DO I = 1,NDOFEL RHS(I,1)=0.D0 DO J = 1,NDOFEL STIF(J,I) = 0.D0 END DO END DO C C Set up integration points and weights CALL KINTPT(PQ,W,NNODE,NINTP) C DO I=1,NINTP C C Shape functions and derivatives CALL KSHAPE(NNODE,PQ(I),F) CALL DSHAPE(NNODE,PQ(I),DF) C C Compute Normal and tangent vectors to boundary at int pt. C Plane of interface defined as average of upper and lower C surfaces (hopefully coincident, but you never know...) TANGENT(1)=.5D0*(COORDS(1,2)-COORDS(1,1) & +COORDS(1,3)-COORDS(1,4)) TANGENT(2)=.5D0*(COORDS(2,2)-COORDS(2,1) & +COORDS(2,3)-COORDS(2,4)) CALL KUNITV(TANGENT,DET) DET=.5D0*DET RNM(1)=-TANGENT(2) RNM(2)=TANGENT(1) C C Relative displacement, in global coords UREL(1) = 0.D0 UREL(2) = 0.D0 VREL(1)=0.D0 VREL(2)=0.D0 DO N = 1,NNODE/2 N2 = N+NNODE/2 DO J = 1,2 UREL(J) = UREL(J) + F(N2)*(U(2*(N2-1)+J)) & - F(N)*(U(2*(N-1)+J)) VREL(J) = VREL(J) + F(N2)*(V(2*(N2-1)+J)) & - F(N)*(V(2*(N-1)+J)) END DO END DO C C GAP(1) is normal separation, GAP(2) is separation in tan dirn GAP(1) = RNM(1)*UREL(1)+RNM(2)*UREL(2) GAP(2) = TANGENT(1)*UREL(1)+TANGENT(2)*UREL(2) VGAP(1)=RNM(1)*VREL(1)+RNM(2)*VREL(2) VGAP(2)=TANGENT(1)*VREL(1)+TANGENT(2)*VREL(2) C CALL SEPLAW(PROPS,GAP,VGAP,TRACT,DTDG,DTIME) C SIG = 1.D0 DO N = 1,NNODE IF (N.GT.NNODE/2) SIG = -1.D0 DO K = 1,2 RHS(2*(N-1)+K,1) = RHS(2*(N-1)+K,1) & + SIG*F(N)*( TRACT(1)*RNM(K) & + TRACT(2)*TANGENT(K))*W(I)*DET END DO END DO SIGN = 1.D0 DO N = 1,NNODE IF (N.GT.NNODE/2) SIGN = -1.D0 SIGM = 1.D0 DO M=1,NNODE IF (M.GT.NNODE/2) SIGM = -1.D0 DO KN = 1,2 DO KM = 1,2 ICOL = 2*(M-1) + KM IROW = 2*(N-1) + KN C STIF(ICOL,IROW) = STIF(ICOL,IROW) & +( (DTDG(1,1)*RNM(KN)+DTDG(2,1)*TANGENT(KN))*RNM(KM) & +(DTDG(1,2)*RNM(KN)+DTDG(2,2)*TANGENT(KN))*TANGENT(KM) ) & *SIGM*SIGN*F(M)*F(N)*W(I)*DET END DO END DO C END DO END DO C END DO C RETURN END C C===================== SUBROUTINE SEPLAW ===================== C SUBROUTINE SEPLAW(PROPS,GAP,VGAP,TRACT,DTDGAP,DTIME) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DIMENSION PROPS(5) DIMENSION GAP(2),VGAP(2),TRACT(2),DTDGAP(2,2) C C Subroutine to specify traction-separation law for C a debonding interface. C C Currently coded for Xu-Needleman constitutive law, i.e. C X.P. Xu and A. Needleman, J. Mech. Phys. Solids, 42, 1397-1434. C C PROPS(1) SIGMA_max C PROPS(2) Delta_n C PROPS(3) Delta_t C PROPS(4) Q C PROPS(5) R C C GAP(1) Normal separation C GAP(2) Tangential separation C TRACT(1) Normal traction C TRACT(2) Tangential traction C C C -- Work of separation SEPWRK = DEXP(1.D0)*PROPS(1)*PROPS(2) DN = PROPS(2) DT = PROPS(3) Q = PROPS(4) R = PROPS(5) C C C1 = ( 1.D0-DEXP(-GAP(2)*GAP(2)/(DT*DT)) ) C1 = C1 * (1.D0-Q)/(R-1.D0) * (R - GAP(1)/DN) C2 = (GAP(1)/DN)*DEXP(-GAP(2)*GAP(2)/(DT*DT)) TRACT(1) = (SEPWRK/DN)*DEXP(-GAP(1)/DN)*(C2+C1) C C1 = Q + (R-Q)/(R-1.D0)*(GAP(1)/DN) C1 = C1 * DEXP(-GAP(1)/DN) * DEXP(-GAP(2)*GAP(2)/(DT*DT)) C1 = C1 * 2.D0*(DN/DT)*(SEPWRK/DN) TRACT(2) = C1*GAP(2)/DT C C1 = (1.D0-Q)/(R-1)*(1.D0-DEXP(-GAP(2)*GAP(2)/(DT*DT))) C1 = C1 * (R+1.D0-GAP(1)/DN) C1 = (1.D0-GAP(1)/DN)*DEXP(-GAP(2)*GAP(2)/(DT*DT)) - C1 DTDGAP(1,1) = (SEPWRK/(DN*DN))*DEXP(-GAP(1)/DN) * C1 C C1 = Q + (GAP(1)/DN) * (R-Q)/(R-1.D0) C1 = C1 * DEXP(-GAP(1)/DN)*DEXP(-GAP(2)*GAP(2)/(DT*DT)) C1 = 2.D0*(SEPWRK/(DT*DT)) * C1 DTDGAP(2,2) = C1*(1.D0-2.D0*(GAP(2)*GAP(2)/(DT*DT))) C C1 = -GAP(1)/DN + (1.D0-Q)/(R-1.D0)*(R-GAP(1)/DN) C1 = C1 * DEXP(-GAP(1)/DN)*DEXP(-GAP(2)*GAP(2)/(DT*DT)) C1 = 2.D0*(SEPWRK/(DT*DN))*C1 DTDGAP(1,2) = (GAP(2)/DT)*C1 DTDGAP(2,1) = DTDGAP(1,2) C C ZETA is the fictitious viscosity used to regularize the instability C problem. In our papere, \zeta_n is ZETA*PROPS(1) here. Simply C setting ZETA=0.D0 (or delete the following three lines C gives the standard Xu-Needleman model. ZETA=0.0001D0 TRACT(1)=TRACT(1)+ZETA*PROPS(1)*VGAP(1)/DN DTDGAP(1,1)=DTDGAP(1,1)+ZETA*PROPS(1)/DN/DTIME C RETURN END C C C========================== SUBROUTINE INTPTS ================================= C SUBROUTINE KINTPT(PQ,W,KNODE,NINTP) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C DIMENSION PQ(*),W(*) C C Subroutine to initialize integration points PQ and weights W C NINTP = No. integration points (DEFAULT=2) C KNODE = No. of nodes on element (DEFAULT=4) C C Integration points for 2D interface elements C C 4----3 C | | C | | C 1----2 c IF (NINTP.EQ.2) THEN CN=0.5773502691896260D0 PQ(1)=-CN PQ(2)=CN W(1)=1.D0 W(2)=1.D0 ENDIF C RETURN END C C=========================== SUBROUTINE SHAPE =============================== C SUBROUTINE KSHAPE(KNODE,PQ,F) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C DIMENSION F(*) C C IF(KNODE.EQ.4) THEN F(1)=.5D0*(1-PQ) F(2)=.5D0*(1+PQ) F(3)=F(2) F(4)=F(1) ENDIF RETURN END C C============================== SUBROUTINE DSHAPE ==================== C SUBROUTINE DSHAPE(KNODE,PQ,DF) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C DIMENSION DF(*) C C IF(KNODE.EQ.4) THEN DF(1)=-.5D0 DF(2)=.5D0 DF(3)=DF(2) DF(4)=DF(1) ENDIF RETURN END C C========================== SUBROUTINE KUNITV =================== C SUBROUTINE KUNITV(A,AMAG) DOUBLE PRECISION A(2), AMAG C C Normalize vector A and return its magnitude as AMAG C AMAG = DSQRT(A(1)*A(1)+A(2)*A(2)) A(1) = A(1)/AMAG A(2) = A(2)/AMAG C RETURN END