Hi Guys,

First time poster here. Im writing a UMAT to implement the fung elastic model into abaqus. The equation is U=c/2(e^^Q -1) and Q is int terms of the material coefficients of the strain components Eij.

Im using the UANISOHYPER Strain in abaqus, and what i need to define is U - the strain energy density function, UA(1)- Derivatives of strain energy potential with respect to the components of the modified Green strain tensor, UA(2)-2nd Derivatives of strain energy potential with respect to the components of the modified Green strain tensor.

Im modelling an incompressible soft tissue so the third derivative and the derivatives with respect to J ( volume ratio) are ignored.

I'm running a very simple job, 1 element, fixed at one end and being pulled at the opposite end. The umat compiles properly but i get an unexpected error 193 running the job. I have attached the UMAT if anyone has any experience writing something similar. Any help guys would be greatly appreciated.

    SUBROUTINE UANISOHYPER_STRAIN (EBAR, AJ, UA, DU1, DU2, DU3,
     1 TEMP, NOEL, CMNAME, INCMPFLAG, IHYBFLAG, NDI, NSHR, NTENS,
     2 NUMSTATEV, STATEV, NUMFIELDV, FIELDV, FIELDVINC,
     3 NUMPROPS, PROPS)
C
     INCLUDE 'ABA_PARAM.INC'  
C
      DIMENSION EBAR(NTENS), UA(2), DU1(NTENS+1),
     2     DU2((NTENS+1)*(NTENS+2)/2),
     3     DU3((NTENS+1)*(NTENS+2)/2),
     4     STATEV(NUMSTATEV), FIELDV(NUMFIELDV),
     5     FIELDVINC(NUMFIELDV), PROPS(NUMPROPS)
C     
      DIMENSION D(3,3)
      DIMENSION ET(3,3)
C where ET = elasticity tensor and D = term1*ET      
      CHARACTER*80 FUNG_HYPERELASTIC  
C User Code to define
C
C UA(1)=Strain energy density function
C UA(2)=The deviatoric part of the strain energy desnity function
C DU1(NTENS+1)=Derivatives of the strain enery potential wrt to modified green strain
C DU2=Second Derivatives of the strain enery potential wrt to modified green strain
C     Parameters
C     PROPS(X) = Used to determine constant value and location in Abaqus
C     Material Property Interface
        zero = 0.d0
        half = 0.5d0
        one = 1.d0
        two = 2.d0
        twothirds = 2.d0/3.d0
        fourthirds = 4.d0/3.d0
        oneninth = 1.d0/9.d0
        twoninths = 2.d0/9.d0
        fourninths = 4.d0/9.d0
        print*, 'in'                                        
C    Constants
       C=PROPS(1)
       A_1=PROPS(2)
       A_2=PROPS(3)
       A_3=PROPS(4)
       A_4=PROPS(5)
       A_5=PROPS(6)
       A_6=PROPS(7)     
C       converting Green strain to modified green strain        
C       xpow= j**2/3= determination of deformation gradient (volume ratio)
       xpow=exp(log(aj)*twothirds)
       E11 = xpow * ebar(1) + half * ( xpow - one )
       E22 = xpow * ebar(2) + half * ( xpow - one )       
       E33 = xpow * ebar(3) + half * ( xpow - one )
       E12 = xpow * ebar(4)
       E13 = xpow * ebar(5)
       E23 = xpow * ebar(6)
C      term 1 is c/2*expQ and term2 is expQ  
       term1=half*C*(exp((A_1*(E11)**2)+(A_2*(E22)**2)
     *          +(two*A_3*E11*E22)+(A_4*(E12)**2)
     *          +(two*A_5*E12*E11)+(two*A_6*E12*E22)))
       term2=(exp((A_1*(E11)**2)+(A_2*(E22)**2)
     *          +(two*A_3*E11*E22)+(A_4*(E12)**2)
     *          +(two*A_5*E12*E11)+(two*A_6*E12*E22)))     
C      let z1 - z3 be the derivative of the u substition of the first derivative.
C      i.e z1 wrt to E11 z2 wrt E22 etc    
        z1=two*A_1*E11+two*A_3*E22+two*A_5*E12
        z2=two*A_2*E22+two*A_3*E11+two*A_6*E12
        z3=two*A_4*E12+two*A_5*E11+two*A_6*E22
C        
        ET(1,1)=2*A_1+(z1**2)
        ET(1,2)=2*A_3+(z1*z2)
        ET(2,2)=2*A_2+(z2**2)
        ET(1,3)=2*A_5+(z1*z2)
        ET(2,3)=2*A_6+(z2*z3)
        ET(3,3)=2*A_4+(z3**2)
C Symmetry
        ET(2,1)=ET(1,2)
        ET(3,1)=ET(1,3)
        ET(3,2)=ET(2,3)
C        
        D1111=term1*ET(1,1)
        D1122=term1*ET(1,2)
        D2222=term1*ET(2,2)
        D1133=term1*ET(1,3)
        D2233=term1*ET(2,3)
        D3333=term1*ET(3,3)
        D2211=D1122
        D3311=D1133
        D3322=D2233
  
C                            
C        write(6,*)aj,twothirds,xpow
C                     
C        write(6,*)NOEL,E11  
C     Strain Energy Density Function
       UA(1)=U
       U=half*C*(exp((A_1*(E11)**2)+(A_2*(E22)**2)
     *          +(two*A_3*E11*E22)+(A_4*(E12)**2)
     *          +(two*A_5*E12*E11)+(two*A_6*E12*E22))-1)
C
C      UA(2)=not needed             
C
C      let z1 - z3 be the derivative of the u substition of the first derivative.
C      i.e z1 wrt to E11 z2 wrt E22 etc
       z1=two*A_1*E11+two*A_3*E22+two*A_5*E12
       z2=two*A_2*E22+two*A_3*E11+two*A_6*E12
       z3=two*A_4*E12+two*A_5*E11+two*A_6*E22
C        
       dUde11=half*C*(exp((A_1*(E11)**2)+(A_2*(E22)**2)
     *          +(two*A_3*E11*E22)+(A_4*(E12)**2)
     *          +(two*A_5*E12*E11)+(two*A_6*E12*E22))*(z1))
       dUde22=half*C*(exp((A_1*(E11)**2)+(A_2*(E22)**2)
     *          +(two*A_3*E11*E22)+(A_4*(E12)**2)
     *          +(two*A_5*E12*E11)+(two*A_6*E12*E22))*(z2))
       dUde33=zero
       dUde12=half*C*(exp((A_1*(E11)**2)+(A_2*(E22)**2)
     *          +(two*A_3*E11*E22)+(A_4*(E12)**2)
     *          +(two*A_5*E12*E11)+(two*A_6*E12*E22))*(z3))
       dUde13=zero
       dUde23=zero
       dUdeJ=zero         
C       
       dUde11=du1(1)
       dUde22=du1(2)
       dUde33=du1(3)
       dUde12=du1(4)
       dUde13=du1(5)
       dUde23=du1(6)
       dUdeJ =du1(7)                              
C      Second derivatives with respect to modified green strain tensors               
C      du2(X)   
       d2Ude11de11=dUde11*z1
       d2Ude11de22=dUde11*z2
       d2Ude22de22=dUde22*z2
       d2Ude11de12=dUde11*z3
       d2Ude22de12=dUde22*z3
       d2Ude12de12=dUde12*z3  
C
       du2(indx(1,1)) = d2Ude11de11
       du2(indx(1,2)) = d2Ude11de22
       du2(indx(2,2)) = d2Ude22de22
       du2(indx(1,3)) = zero
       du2(indx(2,3)) = zero
       du2(indx(3,3)) = zero
       du2(indx(1,4)) = d2Ude11de12
       du2(indx(2,4)) = d2Ude22de12
       du2(indx(3,4)) = zero
       du2(indx(4,4)) = d2Ude12de12
       du2(indx(1,5)) = zero
       du2(indx(2,5)) = zero
       du2(indx(3,5)) = zero
       du2(indx(4,5)) = zero
       du2(indx(5,5)) = zero
      return
      end
* Maps index from Square to Triangular storage
* of symmetric matrix
*
      integer function indx(i,j)
*
      include 'aba_param.inc'
C      write(6,*)'in'
*
      ii = min(i,j)
      jj = max(i,j)
*
      indx = ii + jj*(jj-1)/2
C            write(6,*)'out'
*
       
      return


Thanks Guys

Any ideas Guys?

Hi slimjimmy ,

I suppose that your problems have been resolved,
however, I am doing the same subroutine in Abaqus for the first time (I am a Ansys user)

-) I am not sure if this is correct:

dUde11=du1(1)
dUde22=du1(2)
dUde33=du1(3)
dUde12=du1(4)
dUde13=du1(5)
dUde23=du1(6)
dUdeJ =du1(7)

why not the inverse?...
du1(1)=dUde11
...
...

-)  additionally, did you try to define UA2 and don't use it?

These my two hypotheses arise from the example "UANISOHYPER_STRAIN" reported in the Abaqus manual

best
 

Hi Biobio,

Cheers for the reply, I'm still stuck on this! I noticed your first correction straight away and changed it as you said, but still to no avail.

I didnt use UA2. My reason being the documentation says this is the deviatoric part of the strain energy function,I however am just differentiating the whole strain energy function. Perhaps this could be a source error?

My umat is running and I'm using print statements "in" at the start of the code and "out" at the end of the code  and I can see that the umat is being called 10  times after which an "unexpected error 193" appears. I am also printing my arrays and none of them are filling.

How are you getting on implementing it into Ansys? There is an existing fung model in abaqus if you dont already know, may be of use to you. I need to changed it slightly hence why im doing a umat.

Cheers for the tips, any other help would be great


 

I have just noticed that the ebar variable (modified green strain -  variable to be passed in by abaqus) is not updating from one iteration to the next. If i assign an arbitrary value to ebar the umat runs but still errors except this time without the "unexpected error 193"

Is there any particular reason why ebar is not updating from its original value of 0?  

Try debugging it by writing a Fortran 'wrapper script' to call the subroutine and pass in fictitious strain values. Or else attach a proper debugger to the running analysis as shown here (http://www.ichec.ie/support/tutorials/abaqusdebug.pdf). Should help ya figure out what's going wrong.

Ebar probably isn't updating because you're probably not returning realistic strain energy density values to Abaqus, they might be very small etc.