BBarboza
Geotechnical
- Oct 8, 2020
- 1
Hi folks,
I'm trying to delete elements from my mesh when a critical strain is reached. To be able to do that I wrote the following subroutine:
Do you see any mistake in that? I can't have it compiled and also can't find the error.
PS: I must use vusdfld cause I'm in the explicit solver.
PS: I'm dealing with a elastoplastic material (Mohr-Coulomb Plasticity)
Thanks very much!
I'm trying to delete elements from my mesh when a critical strain is reached. To be able to do that I wrote the following subroutine:
Code:
subroutine vusdfld(nblock,nstatev,nfieldv,nprops,ndir,
& nshr,jElem,kIntPt,kLayer,kSecPt,
& stepTime,totalTime,dt,cmname,
& coordMp,direct,T,charLength,props,
& stateOld,stateNew,field )
C
C include 'ABA_PARAM.INC'
C
include 'vaba_param.inc'
C
dimension jElem(nblock), coordMp(nblock,*),
1 direct(nblock,3,3), T(nblock,3,3),
2 charLength(nblock), props(nprops),
3 stateOld(nblock,nstatev),
4 stateNew(nblock,nstatev),
5 field(nblock,nfieldv)
character*80 cmname
parameter( nrData=6 )
character*3 cData(maxblk*nrData)
dimension rData(maxblk*nrData), jData(maxblk*nrData)
C
real*8,dimension(6)::sigmaNew
real*8,dimension(1)::Critical_Strain
real*8,dimension(1)::MaxStress
real*8,dimension(1)::i1
real*8,dimension(1)::i2
real*8,dimension(1)::i3
real*8,dimension(1)::term1
real*8,dimension(1)::term2
real*8,dimension(1)::term3
real*8,dimension(1)::termaux
real*8,dimension(1)::phiaux
real*8,dimension(1)::phi
real*8,dimension(1)::maxprincipalE
C
jStatus = 1
call vgetvrm('LE',rdata,jdata,cdata,jStatus)
if( jStatus .ne. 0 ) then
call xplb_abqerr(-2,'Utility routine VGETVRM '//
. 'failed to get variable.',0,zero,' ')
call xplb_exit
end if
C
Critical_Strain(1)=0.3000
C
do k = 1, nblock
StrainE11 = rData(k)
StrainE22 = rData(nblock+k)
StrainE33 = rData(2*nblock+k)
StrainE12 = rData(3*nblock+k)
StrainE23 = rData(4*nblock+k)
StrainE31 = rData(5*nblock+k)
i1(1) = StrainE11 + StrainE22 + StrainE33
i2(1) = StrainE11*StrainE22 + StrainE22*StrainE33 + StrainE33*StrainE11 - StrainE12*StrainE12 - StrainE23*StrainE23 - StrainE31*StrainE31
i3(1) = StrainE11*StrainE22*StrainE33 - StrainE11*StrainE23*StrainE23 - StrainE22*StrainE31*StrainE31 - StrainE33*StrainE12*StrainE12 + 2*StrainE12*StrainE23*StrainE31
phiaux(1) = (2*i1(1)*i1(1)*i1(1) - 9*i1(1)*i2(1) + 27*i3(1)*i3(1)*i3(1))/(2*(i1(1)*i1(1)-3*i2(1))**(3./2.))
phi(1) = 1/3 * 1/cos(phiaux(1))
termaux(1) = i1(1)*i1(1) - 3*i2(1)*i2(1)
term1(1) = i1(1)/3 + 2/3*SQRT(termaux(1))*cos(phi(1))
term2(1) = i1(1)/3 + 2/3*SQRT(termaux(1))*cos(phi(1) - 2*3.141592/3.)
term3(1) = i1(1)/3 + 2/3*SQRT(termaux(1))*cos(phi(1) - 4*3.141592/3.)
maxprincipalE(1) = term1(1) + term2(1) + term3(1)
if (maxprincipalE(1) .ge. Critical_Strain(1)) then
stateNew(k,1) = 0
end if
end do
RETURN
END
Do you see any mistake in that? I can't have it compiled and also can't find the error.
PS: I must use vusdfld cause I'm in the explicit solver.
PS: I'm dealing with a elastoplastic material (Mohr-Coulomb Plasticity)
Thanks very much!