C INCRA4 SOURCE CHAT 05/01/13 00:34:59 5004 C INCRA4 SOURCE AC2 96/01/19 21:29:40 1995 SUBROUTINE INCRA4(SIG,VAR,EPSVPT,VARPT,XMAT,ALFA,NSTRS0, & NVARI,NCOMAT) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION SIG(*),VAR(*),EPSVPT(*),VARPT(*),XMAT(*) DIMENSION SM(4),SB(4),XXM1(4),XXB1(4),XXM2(4),XXB2(4) DIMENSION GAMM(4),GABB(4),EPMM(4),EPBB(4),XXM(4),XXB(4) DIMENSION ANSM(4),ANSB(4),ANM(4),ANB(4) DIMENSION EPS0(8) PARAMETER (AMAX = 1.0E20) DETIER=2.0D0/3.0D0 ROOT=SQRT(1.5D0) C-------------------------------------------------------------------| C STEP-1 SEPARATING MEMBRANE & BENDING PARTS | C-------------------------------------------------------------------I NSTRS = NSTRS0+2 NS = NSTRS/2 SM(1)=SIG(1) SM(2)=SIG(2) SM(3)=0.0D0 IF (NS.EQ.4) SM(4)=SIG(3) SB(1)=SIG(NS) SB(2)=SIG(NS+1) SB(3)=0.0D0 IF (NS.EQ.4) SB(4)=SIG(NS+2) DO 20 I=1,NS XXM1(I) = DETIER * XMAT(11) * VAR( I) XXM2(I) = DETIER * XMAT(14) * VAR(2*NS+I) XXB1(I) = DETIER * XMAT(11) * VAR( NS+I) XXB2(I) = DETIER * XMAT(14) * VAR(3*NS+I) GAMM(I) = VAR(4*NS+I) GABB(I) = VAR(5*NS+I) EPMM(I) = VAR(6*NS+I) EPBB(I) = VAR(7*NS+I) 20 CONTINUE C-------------------------------------------------------------------| C STEP-2: EVALUATION OF J2 SIGMA - X | C-------------------------------------------------------------------I XXM(1)= (2.D0*SM(1) -SM(2))/3.D0 -XXM1(1) -XXM2(1) XXM(2)= (2.D0*SM(2) -SM(1))/3.D0 -XXM1(2) -XXM2(2) XXM(3)= -1.D0*(SM(1)+SM(2))/3.D0 -XXM1(3) -XXM2(3) XXM(4)=0.0D0 IF(NS.EQ.4) XXM(4)= SM(4) -XXM1(4) -XXM2(4) XXB(1)= (2.D0*SB(1) -SB(2))/3.D0 -XXB1(1) -XXB2(1) XXB(2)= (2.D0*SB(2) -SB(1))/3.D0 -XXB1(2) -XXB2(2) XXB(3)= -1.D0*(SB(1)+SB(2))/3.D0 -XXB1(3) -XXB2(3) XXB(4)=0.0D0 IF(NS.EQ.4) XXB(4)= SB(4) -XXB1(4) -XXB2(4) T1 = PROCON(XXM,XXM,NS) AJM2 = 3.0D0*T1/2.D0 T1 = PROCON(XXB,XXB,NS) AJB2 = 3.0D0*T1/2.D0 AJ2 = SQRT( AJM2 + ALFA*ALFA*AJB2 ) C--------------------------------------------------------------------| C STEP-3: CALCUL OF SIGV/K | C -------------------------------------------------------------------I RR = VAR (4*NSTRS+2) RS = XMAT(10)*RR SK = (AJ2 - RS - XMAT(7))/(XMAT(8) + XMAT(9)*RR ) C--------------------------------------------------------------------| C STEP-4: CALCULATION OF EFFECTIVE INELASTIC STRAIN INCREMENT (P) | C -------------------------------------------------------------------I IF (SK.GT.0.0D0) THEN PPT1= XPUISS(SK,XMAT(6) ,AMAX) PPT2= XPUISS(SK,XMAT(6)+1,AMAX) ELSE PPT1= 0.0D0 PPT2= 0.0D0 ENDIF T1 = XMAT(5)*PPT2 PPT2=1.0D20 IF(DABS(T1).LT.40.0D0) PPT2= DEXP(T1) PPT = PPT1 * PPT2 VARPT(4*NSTRS+1) = PPT C---------------------------------------------------------------------| C******* CALCULATION OF INELASTIC STRAIN INCREMENTS (EPSV0) | C---------------------------------------------------------------------I IF(AJ2.EQ.0.0D0) AJ2=1.0D-20 T1 = 1.5D0*PPT/AJ2 EPS0(1) = T1*XXM(1) EPS0(2) = T1*XXM(2) EPS0(3) = T1*XXM(3) IF(NS.EQ.4) EPS0(4)= T1*XXM(4) EPS0(NS+1) = T1*XXB(1) *ALFA *ALFA EPS0(NS+2) = T1*XXB(2) *ALFA *ALFA EPS0(NS+3) = T1*XXB(3) *ALFA *ALFA IF(NS.EQ.4) EPS0(NS+4)= T1*XXB(4) *ALFA*ALFA DO 30 I=1,NSTRS 30 VARPT(3*NSTRS+I)=EPS0(I) C C----------------------------------------------------------------| C ******* CALCULATION OF PI(P) | C----------------------------------------------------------------I P = VAR(4*NSTRS+1) T1 = XMAT(16)*P T2 = DEXP(-T1) PIP= XMAT(13)+(1.0D0-XMAT(13)) * T2 C----------------------------------------------------------------| C******* CALCULATION OF XII | C----------------------------------------------------------------I T1=PROCON(XXM1,XXM1,NS) T1=3.D0*T1/2.D0 T2=PROCON(XXB1,XXB1,NS) T2=3.D0*T2/2.D0 C X1II = SQRT(T1+ALFA*ALFA*T2) X1IISM = 0.0D0 IF(XMAT(17).NE.0.D0) X1IISM = X1II/XMAT(17) C T1=PROCON(XXM2,XXM2,NS) T1=3.D0*T1/2.D0 T2=PROCON(XXB2,XXB2,NS) T2=3.D0*T2/2.D0 C X2II = SQRT(T1+ALFA*ALFA*T2) X2IISM = 0.0D0 IF(XMAT(18).NE.0.D0) X2IISM = X2II/XMAT(18) C COX1 = 0.D0 COX2 = 0.D0 IF ((X1II.NE.0.D0).AND.(X1IISM.GT.0.D0)) . COX1 = XPUISS(X1IISM,XMAT(19),AMAX) / X1II IF ((X2II.NE.0.D0).AND.(X2IISM.GT.0.D0)) . COX2 = XPUISS(X2IISM,XMAT(20),AMAX) / X2II C ----------------------------------------------------------------| C******* CALCULATION OF DX1 AND DX2 | C-----------------------------------------------------------------I DO 74 I=1,NSTRS T1 = EPS0(I) T2 = XMAT(12) * PIP * VAR(I) * PPT T3 = XMAT(11) * COX1 * VAR(I) VARPT(I)= T1 - T2 -T3 T2 = XMAT(15) * PIP * VAR(NSTRS+I) * PPT T3 = XMAT(14) * COX2 * VAR(NSTRS+I) VARPT(NSTRS+I)= T1 - T2 - T3 74 CONTINUE C-----------------------------------------------------------------| C******* CALCULATION OF DR | C-----------------------------------------------------------------I T1 = 1.0D0-(1.0D0-VAR(4*NSTRS+3)/XMAT(25))**2 QR = VAR(4*NSTRS+3)-XMAT(26)*T1 CO1= QR-VAR(4*NSTRS+2) CO0= ABS(CO1) CO = XPUISS(CO0,XMAT(22)-1.0D0,AMAX) T1= XMAT(16)*(VAR(4*NSTRS+3)-VAR(4*NSTRS+2))*PPT T2= XMAT(21)*CO*CO1 DR= T1 + T2 VARPT(4*NSTRS+2)=DR C----------------------------------------------------------------| C******* CALCULATION OF PROD (N X N*) | C----------------------------------------------------------------I TM1 = EPMM(1)-GAMM(1) TM2 = EPMM(2)-GAMM(2) TM3 = EPMM(3)-GAMM(3) TM4 = 0.0D0 IF(NS.EQ.4) TM4 = EPMM(4)-GAMM(4) AIM2 = TM1*TM1+TM2*TM2+TM3*TM3+2.D0*TM4*TM4 AIM2 = AIM2*3.D0/2.D0 TB1 = EPBB(1)-GABB(1) TB2 = EPBB(2)-GABB(2) TB3 = EPBB(3)-GABB(3) TB4 = 0.0D0 IF(NS.EQ.4) TB4 = EPBB(4)-GABB(4) AIB2 = TB1*TB1+TB2*TB2+TB3*TB3+2.D0*TB4*TB4 AIB2 = AIB2*3.D0/2.D0 AI2 = SQRT(AIM2+(AIB2/ALFA/ALFA)) IF (AI2.LT.1.D-10) THEN IF (AI2.EQ.0.D0) AI2=1.D-10 AI20=AI2*1.D20 TM1=TM1*1.D20 TM2=TM2*1.D20 TM3=TM3*1.D20 TM4=TM4*1.D20 TB1=TB1*1.D20 TB2=TB2*1.D20 TB3=TB3*1.D20 TB4=TB4*1.D20 ELSE AI20=AI2 ENDIF ANSM(1) = ROOT*TM1/AI20 ANSM(2) = ROOT*TM2/AI20 ANSM(3) = ROOT*TM3/AI20 ANSM(4) = 0.0D0 IF(NS.EQ.4) ANSM(4) = ROOT*TM4/AI20 ANSB(1) = ROOT*TB1/AI20 ANSB(2) = ROOT*TB2/AI20 ANSB(3) = ROOT*TB3/AI20 ANSB(4) = 0.0D0 IF(NS.EQ.4) ANSB(4) = ROOT*TB4/AI20 ANM (1) = ROOT*XXM(1)/AJ2 ANM (2) = ROOT*XXM(2)/AJ2 ANM (3) = ROOT*XXM(3)/AJ2 ANM (4) = 0.0D0 IF(NS.EQ.4) ANM (4) = ROOT*XXM(4)/AJ2 ANB (1) = ROOT*XXB(1)/AJ2 ANB (2) = ROOT*XXB(2)/AJ2 ANB (3) = ROOT*XXB(3)/AJ2 ANB (4) = 0.0D0 IF(NS.EQ.4) ANB (4) = ROOT*XXB(4)/AJ2 76 CONTINUE FF = DETIER*AI2-VAR(4*NSTRS+4) HF=0.0D0 IF(FF.GT.0.0D0) HF=1.0D0 ANAN1=PROCON(ANM,ANSM,NS) * ANAN1 = ANM(1)*ANSM(1)+ANM(2)*ANSM(2) * ANAN1 = ANAN1+2.D0*ANM(3)*ANSM(3) ANANSM=0.0D0 IF(ANAN1.GT.0.0D0) ANANSM=ANAN1 ANAN1=PROCON(ANB,ANSB,NS) * ANAN1 = ANB(1)*ANSB(1)+ANB(2)*ANSB(2) * ANAN1 = ANAN1+2.D0*ANB(3)*ANSB(3) ANANSB=0.0D0 IF(ANAN1.GT.0.0D0) ANANSB=ANAN1 C----------------------------------------------------------------| C******* CALCULATION OF DQ | C----------------------------------------------------------------I ANANS=ANANSM+ANANSB/ALFA/ALFA DQ = XMAT(23)*HF*ANANS*PPT VARPT(4*NSTRS+4)=DQ C----------------------------------------------------------------| C******* CALCULATION OF DGETA | C----------------------------------------------------------------I T1 = ROOT*(1.0D0-XMAT(23))*HF*ANANS*PPT DO 78 I=1,NS VARPT(2*NSTRS+I) = T1*ANSM(I) VARPT(2*NSTRS+NS+I)= T1*ANSB(I) 78 CONTINUE C----------------------------------------------------------------| C******* CALCULATION OF DQQ | C----------------------------------------------------------------I VARPT(4*NSTRS+3)=2.0D0*XMAT(24)*(XMAT(25)-VAR(4*NSTRS+3))*DQ C EPSVPT(1)=EPS0(1) EPSVPT(2)=EPS0(2) IF (NS.EQ.4) EPSVPT(3)=EPS0(4) EPSVPT(NS)=EPS0(NS+1) EPSVPT(NS+1)=EPS0(NS+2) IF (NS.EQ.4) EPSVPT(6)=EPS0(8) C RETURN END