Sources Esope | Cast3M

incrb4
C INCRB4    SOURCE    CB215821  16/04/21    21:17:08     8920
C INCRB4    SOURCE    AC2       96/01/19    21:29:40     1995
       SUBROUTINE INCRB4(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 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(9) * VAR(     I)
       XXM2(I) = DETIER * XMAT(12) * VAR(2*NS+I)
       XXB1(I) = DETIER * XMAT(9) * VAR(  NS+I)
       XXB2(I) = DETIER * XMAT(12) * VAR(3*NS+I)
       EPMM(I) = VAR(4*NS+I)
       EPBB(I) = VAR(5*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
       RISOT = XMAT(17) -
     . ((XMAT(17) - XMAT(7)) * EXP(-XMAT(18) * VAR(3*NSTRS+1)))
       SK = (AJ2 - RISOT)/ XMAT(8)
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(ABS(T1).LT.40.0D0) PPT2= EXP(T1)
       PPT = PPT1 * PPT2
       VARPT(3*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
         VARPT(2*NSTRS+I)=EPS0(I)
 30    CONTINUE
C----------------------------------------------------------------|
C ******* CALCULATION OF PI(P)                                   |
C----------------------------------------------------------------I
       P  = VAR(3*NSTRS+1)
       T1 = XMAT(14)*P
       T2 = EXP(-T1)
       PIP= XMAT(11)+(1.0D0-XMAT(11)) * 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(15).NE.0.D0)   X1IISM = X1II/XMAT(15)
C
       T1=PROCON(XXM2,XXM2,NS)
       T1=3.D0*T1/2.D0
       T2=PROCON(XXB2,XXB2,NS)
       T2=3.D0*T2/2.D0
C
       COX1 = 0.D0
       IF ((X1II.NE.0.D0).AND.(X1IISM.GT.0.D0))
     . COX1 = XPUISS(X1IISM,XMAT(16),AMAX)  / X1II
C ----------------------------------------------------------------|
C*******  CALCULATION OF DX1 AND DX2                              |
C-----------------------------------------------------------------I
       DO 74 I=1,NSTRS
         T1  =  EPS0(I)
         T2  =  XMAT(10) * PIP *  VAR(I) * PPT
         T3  =  XMAT(9) * COX1 *  VAR(I)
         VARPT(I)= T1 - T2 -T3
         T2  =  XMAT(13) * PIP *  VAR(NSTRS+I) * PPT
         VARPT(NSTRS+I)= T1 - T2
  74   CONTINUE
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