SUBROUTINE UHYPER(BI1,BI2,AJ,U,UI1,UI2,UI3,TEMP,NOEL,CMNAME,
     $                  INCMPFLAG,NUMSTATEV,STATEV,NUMFIELDV,
     $                  FIELDV,FIELDVINC,NUMPROPS,PROPS)
C===============================================================
C     USER DEFINED HYPERELASTIC MATERIAL SUBROUTINE 
C       FOR MECHANOCHEMICALLY RESPONSIVE ELASTOMER
C       TO BE USED IN ABAQUS STANDARD
C     FOMULATED AND WRITTEN BY QIMING WANG, JUNE 17, 2014
C---------------------------------------------------------------
C     MATERIAL PROPERTIES TO BE PASSED TO THE SUBROUTINE:
C       PROPS(1) - TOTAL CHAIN NUMBER NCHAIN
C       PROPS(2) - CHAIN DISTRIBUTION PARAMETER DETA
C       PROPS(3) - LOG (PROPS3) IS MEAN OF LOG (CHAIN LENGTH)
C       PROPS(4) - PARAMETER TO REIFORCE THE INCOMPRESSIBILITY D1
C     PARAMETERS USED IN THE SUBROUTINE
C       LAM - CHAIN STRETCH
C       BETA - BETA=L^(-1)(LAM/SQRT(NI)),ESTIMATED BY X(3-X^2)/(1-X^2)
C       NI - CHAIN NUMBER PROPERBILITY 
C       KI - CHAIN LENGTH
C===============================================================
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*80 CMNAME
      DIMENSION U(2),UI1(3),UI2(6),UI3(6),STATEV(*),FIELDV(*),
     $          FIELDVINC(*),PROPS(*)
      REAL(8) LAM,BETA,BETAI,NI,KI
      REAL(8) NCHAIN,DETA,NBAR,D1
C    
      NCHAIN = PROPS(1)
      DETA = PROPS(2)
      NBAR = PROPS(3)
      D1 = PROPS(4)
C
      U(1)=0.0
      U(2)=0.0
      UI1(1)=0.0
      UI1(2)=0.0
      UI1(3)=0.0      
      UI2(1)=0.0
      UI2(2)=0.0
      UI2(3)=0.0
      UI2(4)=0.0
      UI2(5)=0.0
      UI2(6)=0.0
      UI3(1)=0.0
      UI3(2)=0.0
      UI3(3)=0.0
      UI3(4)=0.0
      UI3(5)=0.0
      UI3(6)=0.0
C      
      DO KI=3, 100, 1
            LAM=SQRT(BI1*(AJ**(2.0/3.0))/3.0)
            BETA=(LAM/SQRT(KI))*(3.0-(LAM/SQRT(KI))**2)
     $           /(1.0-(LAM/SQRT(KI))**2)
     				BETAI=(1.0/SQRT(KI))*(3.0-(1.0/SQRT(KI))**2)
     $           /(1.0-(1.0/SQRT(KI))**2)
            NI=1.0/(KI*SQRT(2.0*3.14)*DETA)
     $       *EXP(-1.0*(LOG(KI)-LOG(NBAR))**2/(2.0*DETA**2))
      U(1)=U(1)+NCHAIN*NI*KI*(BETA/TANH(BETA)+LOG(BETA/SINH(BETA)))
     $		-NCHAIN*NI*KI*(BETAI/TANH(BETAI)+LOG(BETAI/SINH(BETAI)))
            UI1(1)=UI1(1)+NCHAIN*NI*SQRT(KI)*BETA/SQRT(12.0*BI1)
     $              *(AJ**(1.0/3.0))
			UI1(3)=UI1(3)+NCHAIN*NI*SQRT(KI)*BETA/SQRT(27.0*BI1)
     $				*(AJ**(-2.0/3.0))
            UI2(1)=UI2(1)+NCHAIN*NI*SQRT(KI)/SQRT(12.0)
     $       *( 1.0/BI1/SQRT(12.0*KI)*AJ**(2.0/3.0)
     $       /(1.0/BETA**2-1.0/(SINH(BETA)**2))
     $       -BETA/2.0/SQRT(BI1**3)*AJ**(1.0/3.0) )
			UI2(3)=UI2(3)+NCHAIN*NI*SQRT(KI)/SQRT(27.0)
     $       *(1.0/SQRT(27.0*KI)*AJ**(-4.0/3.0)
     $       /(1.0/BETA**2-1.0/(SINH(BETA)**2))
     $       -BETA*2.0/3.0*BI1**(-1.0/2.0)*AJ**(-5.0/3.0))
			UI2(5)=UI2(5)+NCHAIN*NI*SQRT(KI)/SQRT(12.0)
     $       *(1.0/SQRT(27.0*KI)*AJ**(-1.0/3.0)
     $       /(1.0/BETA**2-1.0/(SINH(BETA)**2))
     $       -BETA/3.0/SQRT(BI1)*AJ**(-2.0/3.0))
            UI3(1)=UI3(1)+NCHAIN*NI*SQRT(KI)/SQRT(12.0)
     $       *( 1.0/BI1/SQRT(12.0*KI)* 
     $          ( 2.0/3.0*AJ**(-1.0/3.0)
     $              /(1.0/BETA**2-1.0/(SINH(BETA)**2))
     $            - SQRT(BI1/27.0/KI)*
     $              (2.0*COSH(BETA)/SINH(BETA)**3-2.0*BETA**(-3))
     $              /(1.0/BETA**2-1.0/(SINH(BETA)**2))**3  ) 
     $          -BI1**(-3.0/2.0)/2*
     $            ( SQRT(BI1/27.0/KI)*AJ**(-1.0/3.0)
     $                /(1.0/BETA**2-1.0/(SINH(BETA)**2))
     $              + BETA*AJ**(-2.0/3.0)/3.0 )  )  
			UI3(4)=UI3(4)+NCHAIN*NI*SQRT(KI)/SQRT(12.0)
     $       *( 1.0/SQRT(27.0*KI)* 
     $          ( -1.0/3.0*AJ**(-4.0/3.0)
     $              /(1.0/BETA**2-1.0/(SINH(BETA)**2))
     $            - SQRT(BI1/27.0/KI)*AJ**(-1)*
     $              (2.0*COSH(BETA)/SINH(BETA)**3-2.0*BETA**(-3))
     $              /(1.0/BETA**2-1.0/(SINH(BETA)**2))**3  ) 
     $          -BI1**(-1.0/2.0)/3.0*
     $            ( SQRT(BI1/27.0/KI)*AJ**(-4.0/3.0)
     $                /(1.0/BETA**2-1.0/(SINH(BETA)**2))
     $              -2.0/3.0*BETA*AJ**(-5.0/3.0) )  )
            UI3(6)=UI3(6)+NCHAIN*NI*SQRT(KI)/SQRT(27.0)
     $       *( 1.0/SQRT(27.0*KI)* 
     $          ( -4.0/3.0*AJ**(-7.0/3.0)
     $              /(1.0/BETA**2-1.0/(SINH(BETA)**2))
     $            - SQRT(BI1/27.0/KI)*AJ**(-2)*
     $              (2.0*COSH(BETA)/SINH(BETA)**3-2.0*BETA**(-3))
     $              /(1.0/BETA**2-1.0/(SINH(BETA)**2))**3  ) 
     $          -2.0/3.0*BI1**(-1.0/2.0)*
     $            ( SQRT(BI1/27.0/KI)*AJ**(-7.0/3.0)
     $                /(1.0/BETA**2-1.0/(SINH(BETA)**2))
     $              -5.0/3.0*BETA*AJ**(-8.0/3.0) )  )
      END DO
C
	  U(1)=U(1)+(AJ-1.0)**2/D1
	  UI1(3)=UI1(3)+(AJ-1.0)*2.0/D1
	  UI2(3)=UI2(3)+2.0/D1
      RETURN
      END