UANISOHYPER STRAIN help

[slimjimmy]

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

 

[slimjimmy]

Any ideas Guys?

 

[biobio2011]

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

 

[slimjimmy]

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

 

[slimjimmy]

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?

 

[MechIrl]

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.