simple failure criterion in vuanisohyper

[Aisling]

Hi all,

I would like to include a simple element deletion into the vuanisohyper subroutine I got from the abaqus verification manual.

My plan is to set the element deletion flag to zero when some specific criterion is met, for example some max tensile stress. My material is defined as follows in the .inp, where property 6 is our critical condition (say 10MPa).

*Material, name=VUANISO_HGO
*density
<density>
*Depvar, delete=6
1,
*anisotropic hyperelastic,user,formulation=invariant, local directions=2,type=compressible,properties=6
<C10>, <D>, <K1>, <K2>, <Kappa>, <FT>

The subroutine looks like this...

c
c HGO model
c
subroutine vuanisohyper_invhgo (ainv, ua, zeta, nfibers, ninv,
$ ui1, ui2, nb, numprops, props, stateOld, stateNew)
C
include 'vaba_param.inc'
C
dimension ua(nb), ainv(nb,ninv), ui1(nb,ninv),
$ ui2(nb,ninv*(ninv+1)/2), props(numprops), stateOld(nb, nstatev), stateNew(nb, nstatev)
C
c ainv: invariants
c ua : udev
c ui1 : dUdI
c ui2 : d2U/dIdJ
C
parameter ( half = 0.5d0,
* zero = 0.d0,
* one = 1.d0,
* two = 2.d0,
* three= 3.d0,
* four = 4.d0,
* five = 5.d0,
* six = 6.d0,
c
* index_I1 = 1,
* index_J = 3,
* asmall = 2.d-16 )
C
C HGO model
C
C10 = props(1)
rk1 = props(3)
rk2 = props(4)
rkp = props(5)
FT = props(6)
c
do kb = 1,nb
ua(kb) = zero
om3kp = one - three * rkp
do k1 = 1, nfibers
index_i4 = indxInv4(k1,k1)
E_alpha1 = rkp * (ainv(kb,index_i1) - three)
* + om3kp * (ainv(kb,index_i4) - one )
E_alpha = max(E_alpha1, zero)
ht4a = half + sign(half,E_alpha1 + asmall)
aux = exp(rk2*E_alpha*E_alpha)
c energy
ua(kb) = ua(kb) + aux - one
c ui1
ui1(kb,index_i1) = ui1(kb,index_i1) + aux * E_alpha
ui1(kb,index_i4) = rk1 * om3kp * aux * E_alpha
c ui2
aux2 = ht4a + two * rk2 * E_alpha * E_alpha
ui2(kb,indx(index_I1,index_I1))=
* ui2(kb,indx(index_I1,index_I1))
* + aux * aux2
ui2(kb,indx(index_I1,index_i4)) = rk1*rkp*om3kp * aux * aux2
ui2(kb,indx(index_i4,index_i4)) = rk1*om3kp*om3kp*aux * aux2
end do
c
c deviatoric energy
c
ua(kb) = ua(kb) * rk1 / (two * rk2)
ua(kb) = ua(kb) + C10 * (ainv(kb,index_i1) - three)
c
c compute derivatives
c
ui1(kb,index_i1) = rk1 * rkp * ui1(kb,index_i1) + C10
ui2(kb,indx(index_I1,index_I1))= ui2(kb,indx(index_I1,index_I1))
* * rk1 * rkp * rkp
end do
c
c compressible case
if(props(2).gt.zero) then
do kb = 1,nb
Dinv = one / props(2)
det = ainv(kb,index_J)
ui1(kb,index_J) = Dinv * (det - one/det)
ui2(kb,indx(index_J,index_J))= Dinv * (one + one / det / det)
end do
end if

return
end
C-------------------------------------------------------------
C Function to map index from Square to Triangular storage
C of symmetric matrix
C
integer function indx( i, j )
include 'vaba_param.inc'
ii = min(i,j)
jj = max(i,j)
indx = ii + jj*(jj-1)/2
return
end
C-------------------------------------------------------------
C
C Function to generate enumeration of scalar
C Pseudo-Invariants of type 4

integer function indxInv4( i, j )
include 'vaba_param.inc'
ii = min(i,j)
jj = max(i,j)
indxInv4 = 4 + jj*(jj-1) + 2*(ii-1)
return
end
C-------------------------------------------------------------
C
C Function to generate enumeration of scalar
C Pseudo-Invariants of type 5
C
integer function indxInv5( i, j )
include 'vaba_param.inc'
ii = min(i,j)
jj = max(i,j)
indxInv5 = 5 + jj*(jj-1) + 2*(ii-1)
return
end
C------------------------------------------------------------

The various stresses and strains are not calculated in the subroutine, only the strain energy potential derivatives. In order to implement my condition I believe I need to first calculate /sigma_11. Can somebody help me in how I might define \sigma_11 within the subroutine? My question has more to do continuum mechanics I guess but I am thoroughly confused!

Aisling

[Aisling]

more specifically, how can I access the deformation gradient within this subroutine?

[Aisling]

I have now realised that I probably need to call vusdfld subroutine and then call vgtvrm in order to get /sigma_11 directly.

does anybody have any advice on how I should call 2 subroutines like this?

[lumpwood]

I think that calling vgetvrm is a good idea for what you are looking for. Just put the vusdfld subroutine in the same file as the vuanisohyper routine and remember to include the *field keyword in your input file, otherwise the subroutine won't be called. If you are looking for something like max principal stress then you can call the vsprinc routine from within the vusdfld subroutine. Store the value of stress in a state variable for use in the vuanisohyper routine.

[Aisling]

Thanks for the advice lumpwood. I'm using 6.9 at the moment which doesn't have the vsprinc routine. This is my approach at the moment...

do k = 1, nblock
call vgetvrm( 'S', rData, jData, cData, jStatus )
s11=S(1)
field (k,1) = S11
stateNew(k,1) = field(k,1)
end do

it compiles and runs fine. in the .inp then I have

*Element Output, directions=YES
S, SDV, FV

thing is the SDV and S(s11) don't match at all!! the SDV stays the same for the whole simulation. I'm sure it has to do with the way I'm saving the data. Perhaps the array dimensions are wrong? Or I am overwriting something?

[Aisling]

by putting in some print statement i found the problem is in calling vgetvrm. anybody have an explanation for this? The code is from an example! I don't get any detailed errors....

parameter( nrData=6 )
character*3 cData(maxblk*nrData)
dimension S(ndir+nshr), rData(maxblk*nrData), jData(maxblk*nrData)

jStatus = 1
call vgetvrm( 'S', rdata, jData, cData, jStatus)