kkbalaji
2007-12-12, 08:52
Hello Dr. Jorgen,
I deal with thin sheet-like structures that are non-linearly
orthotropic in 1 and 2 directions. I have been unsuccessful at
creating a material user subroutine VUMAT to simulate this behavior.
I created an input file using a 2D planar solid that is 3 X 6 inches.
The problem involves stretching the sheet in 1-direction by fixing
the sheet at one edge and pulling the opposite edge to a fixed
deformation (~30% Strain). The other 2 edges are free to move inward
as the sheet is pulled.
I have attached the material subroutine I created. Once a yield
strain (say ~0.03) is reached, E1 & E2 become dependent on total
strain in 1-direction and decrease in a logarithmic fashion. Material
tests show that the material fails at ~50% strain and in the
subroutine I use a break strain criteria and set all the values as
constant once that strain is reached. I have already checked that the
orthotropic material laws are not violated (like sqrt(E1/E2)>v12).
So long as the strain stays below the yield strain, the analysis
computes the first step. Once it exceeds the yield strain, it gives
erroneous deformation and the analysis exits with severe element
distortion error. Although the overall strain seems less (like 5%) some of the regions exceed strain levels of 50%). I appreaciate your help and any suggestions are
welcome.
Thanks,
Balaji.
Input properties: E1=10000, E2=3100, E3=200, V12=0.4, V13=0.1,
V23=0.1, yStr=0.03, bStr=0.5
C User subroutine VUMAT
subroutine vumat (
C Read only -
* nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
* stepTime, totalTime, dt, cmname, coordMp, charLength,
* props, density, strainInc, relSpinInc,
* tempOld, stretchOld, defgradOld, fieldOld,
* stressOld, stateOld, enerInternOld, enerInelasOld,
* tempNew, stretchNew, defgradNew, fieldNew,
C Write only -
* stressNew, stateNew, enerInternNew, enerInelasNew )
C
include 'vaba_param.inc'
C
dimension coordMp(nblock,*), charLength(nblock), props(nprops),
1 density(nblock), strainInc(nblock,ndir+nshr),
2 relSpinInc(nblock,nshr), tempOld(nblock),
3 stretchOld(nblock,ndir+nshr),
4 defgradOld(nblock,ndir+nshr+nshr),
5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(nblock),
8 stretchNew(nblock,ndir+nshr),
9 defgradNew(nblock,ndir+nshr+nshr),
1 fieldNew(nblock,nfieldv),
2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
3 enerInternNew(nblock), enerInelasNew(nblock)
C
character*80 cmname
dimension intv(2)
parameter ( zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0)
C MD, CD Modulus use equation Modulus = C . Ln(X) + B
parameter (MdModB=-793.85d0, MdModC=-1721.3d0, CdModB=-
73.746d0,
* CdModC=-629.97d0)
!C
!C Check that ndir=3 and nshr=1. If not, exit.
!C
intv(1) = ndir
intv(2) = nshr
if (ndir .ne. 3 .or. nshr .ne. 1) then
call xplb_abqerr(1,'Subroutine VUMAT is implemented '//
* 'only for plane strain and axisymmetric cases '//
* '(ndir=3 and nshr=1)',0,zero,' ')
call xplb_abqerr(-2,'Subroutine VUMAT has been called '//
* 'with ndir=%I and nshr=%I',intv,zero,' ')
call xplb_exit
end if
! Input Parameters
E1 = props(1)
E2 = props(2)
E3 = props(3)
Pr12 = props(4)
Pr13 = props(5)
Pr23 = props(6)
yStr = props(7)
bStr = props(8)
Pr21 = Pr12*E2/E1
Pr31 = Pr13*E3/E1
Pr32 = Pr23*E3/E2
G12 = sqrt(E1*E2)/(2*(1+sqrt(Pr12*Pr21)))
G13 = sqrt(E1*E3)/(2*(1+sqrt(Pr13*Pr31)))
G23 = sqrt(E2*E3)/(2*(1+sqrt(Pr23*Pr32)))
C If stepTime equals to zero, assume the material pure elastic
C and use initial elastic modulus
do k = 1, nblock
curStr = stateOld(k,1)
if ( abs(curStr) .gt. yStr ) then
if ( abs(curStr) .gt. bStr ) then
E1 = MdModC * log(abs(bStr)) + MdModB
E2 = CdModC * log(abs(bStr)) + CdModB
E3 = props(3)
Pr12 = props(4)
Pr13 = props(5)
Pr23 = props(6)
Pr21 = Pr12*E2/E1
Pr31 = Pr13*E3/E1
Pr32 = Pr23*E3/E2
G12 = sqrt(E1*E2)/(2*(1+sqrt(Pr12*Pr21)))
G13 = sqrt(E1*E3)/(2*(1+sqrt(Pr13*Pr31)))
G23 = sqrt(E2*E3)/(2*(1+sqrt(Pr23*Pr32)))
else
E1 = MdModC * log(abs(curStr)) + MdModB
E2 = CdModC * log(abs(curStr)) + CdModB
E3 = props(3)
Pr12 = props(4)
Pr13 = props(5)
Pr23 = props(6)
Pr21 = Pr12*E2/E1
Pr31 = Pr13*E3/E1
Pr32 = Pr23*E3/E2
G12 = sqrt(E1*E2)/(2*(1+sqrt(Pr12*Pr21)))
G13 = sqrt(E1*E3)/(2*(1+sqrt(Pr13*Pr31)))
G23 = sqrt(E2*E3)/(2*(1+sqrt(Pr23*Pr32)))
end if
end if
C Calculate the Stiffness Matrix
Gm = 1/(1-Pr12*Pr21-Pr23*Pr32-Pr31*Pr13-2*Pr21
1 *Pr32*Pr13)
D1111 = E1*(1-Pr23*Pr32)*Gm
D2222 = E2*(1-Pr13*Pr31)*Gm
D3333 = E3*(1-Pr12*Pr21)*Gm
D1122 = E1*(Pr21+Pr31*Pr23)*Gm
D1133 = E1*(Pr31+Pr21*Pr32)*Gm
D2233 = E2*(Pr32+Pr12*Pr31)*Gm
D1212 = G12
D1313 = G13
D2323 = G23
C Copy the values of strains into temp variables
ES1 = strainInc(k,1)
ES2 = strainInc(k,2)
ES3 = strainInc(k,3)
ES4 = strainInc(k,4)
C Update the stress
stressNew(k,1) = stressOld(k,1)
+D1111*ES1+D1122*ES2+D1133*ES3
stressNew(k,2) = stressOld(k,2)
+D1122*ES1+D2222*ES2+D2233*ES3
stressNew(k,3) = stressOld(k,3)
+D1133*ES1+D2233*ES2+D3333*ES3
stressNew(k,4) = stressOld(k,4)+D1212*ES4*2
C Calculate the total strain in 11 direction based on new stress
values
Str11 = stressOld(k,1)
+D1111*ES1+D1122*ES2+D1133*ES3
Str22 = stressOld(k,2)
+D1122*ES1+D2222*ES2+D2233*ES3
Str33 = stressOld(k,3)
+D1133*ES1+D2233*ES2+D3333*ES3
curStrNew = Str11/E1-Str22/E2*Pr21-
Str33/E3*Pr31
C Update the state variable with current value of strain in 11
direction
stateNew(k,1) = curStrNew
C Update the specific internal energy -
stressPower = half * (
* ( stressOld(k,1)+stressNew(k,1) ) * strainInc(k,1) +
* ( stressOld(k,2)+stressNew(k,2) ) * strainInc(k,2) +
* ( stressOld(k,3)+stressNew(k,3) ) * strainInc(k,3) )
+
* ( stressOld(k,4)+stressNew(k,4) ) * strainInc
(k,4)
enerInternNew(k) = enerInternOld(k)
* + stressPower / density(k)
end do
return
end
I deal with thin sheet-like structures that are non-linearly
orthotropic in 1 and 2 directions. I have been unsuccessful at
creating a material user subroutine VUMAT to simulate this behavior.
I created an input file using a 2D planar solid that is 3 X 6 inches.
The problem involves stretching the sheet in 1-direction by fixing
the sheet at one edge and pulling the opposite edge to a fixed
deformation (~30% Strain). The other 2 edges are free to move inward
as the sheet is pulled.
I have attached the material subroutine I created. Once a yield
strain (say ~0.03) is reached, E1 & E2 become dependent on total
strain in 1-direction and decrease in a logarithmic fashion. Material
tests show that the material fails at ~50% strain and in the
subroutine I use a break strain criteria and set all the values as
constant once that strain is reached. I have already checked that the
orthotropic material laws are not violated (like sqrt(E1/E2)>v12).
So long as the strain stays below the yield strain, the analysis
computes the first step. Once it exceeds the yield strain, it gives
erroneous deformation and the analysis exits with severe element
distortion error. Although the overall strain seems less (like 5%) some of the regions exceed strain levels of 50%). I appreaciate your help and any suggestions are
welcome.
Thanks,
Balaji.
Input properties: E1=10000, E2=3100, E3=200, V12=0.4, V13=0.1,
V23=0.1, yStr=0.03, bStr=0.5
C User subroutine VUMAT
subroutine vumat (
C Read only -
* nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
* stepTime, totalTime, dt, cmname, coordMp, charLength,
* props, density, strainInc, relSpinInc,
* tempOld, stretchOld, defgradOld, fieldOld,
* stressOld, stateOld, enerInternOld, enerInelasOld,
* tempNew, stretchNew, defgradNew, fieldNew,
C Write only -
* stressNew, stateNew, enerInternNew, enerInelasNew )
C
include 'vaba_param.inc'
C
dimension coordMp(nblock,*), charLength(nblock), props(nprops),
1 density(nblock), strainInc(nblock,ndir+nshr),
2 relSpinInc(nblock,nshr), tempOld(nblock),
3 stretchOld(nblock,ndir+nshr),
4 defgradOld(nblock,ndir+nshr+nshr),
5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(nblock),
8 stretchNew(nblock,ndir+nshr),
9 defgradNew(nblock,ndir+nshr+nshr),
1 fieldNew(nblock,nfieldv),
2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
3 enerInternNew(nblock), enerInelasNew(nblock)
C
character*80 cmname
dimension intv(2)
parameter ( zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0)
C MD, CD Modulus use equation Modulus = C . Ln(X) + B
parameter (MdModB=-793.85d0, MdModC=-1721.3d0, CdModB=-
73.746d0,
* CdModC=-629.97d0)
!C
!C Check that ndir=3 and nshr=1. If not, exit.
!C
intv(1) = ndir
intv(2) = nshr
if (ndir .ne. 3 .or. nshr .ne. 1) then
call xplb_abqerr(1,'Subroutine VUMAT is implemented '//
* 'only for plane strain and axisymmetric cases '//
* '(ndir=3 and nshr=1)',0,zero,' ')
call xplb_abqerr(-2,'Subroutine VUMAT has been called '//
* 'with ndir=%I and nshr=%I',intv,zero,' ')
call xplb_exit
end if
! Input Parameters
E1 = props(1)
E2 = props(2)
E3 = props(3)
Pr12 = props(4)
Pr13 = props(5)
Pr23 = props(6)
yStr = props(7)
bStr = props(8)
Pr21 = Pr12*E2/E1
Pr31 = Pr13*E3/E1
Pr32 = Pr23*E3/E2
G12 = sqrt(E1*E2)/(2*(1+sqrt(Pr12*Pr21)))
G13 = sqrt(E1*E3)/(2*(1+sqrt(Pr13*Pr31)))
G23 = sqrt(E2*E3)/(2*(1+sqrt(Pr23*Pr32)))
C If stepTime equals to zero, assume the material pure elastic
C and use initial elastic modulus
do k = 1, nblock
curStr = stateOld(k,1)
if ( abs(curStr) .gt. yStr ) then
if ( abs(curStr) .gt. bStr ) then
E1 = MdModC * log(abs(bStr)) + MdModB
E2 = CdModC * log(abs(bStr)) + CdModB
E3 = props(3)
Pr12 = props(4)
Pr13 = props(5)
Pr23 = props(6)
Pr21 = Pr12*E2/E1
Pr31 = Pr13*E3/E1
Pr32 = Pr23*E3/E2
G12 = sqrt(E1*E2)/(2*(1+sqrt(Pr12*Pr21)))
G13 = sqrt(E1*E3)/(2*(1+sqrt(Pr13*Pr31)))
G23 = sqrt(E2*E3)/(2*(1+sqrt(Pr23*Pr32)))
else
E1 = MdModC * log(abs(curStr)) + MdModB
E2 = CdModC * log(abs(curStr)) + CdModB
E3 = props(3)
Pr12 = props(4)
Pr13 = props(5)
Pr23 = props(6)
Pr21 = Pr12*E2/E1
Pr31 = Pr13*E3/E1
Pr32 = Pr23*E3/E2
G12 = sqrt(E1*E2)/(2*(1+sqrt(Pr12*Pr21)))
G13 = sqrt(E1*E3)/(2*(1+sqrt(Pr13*Pr31)))
G23 = sqrt(E2*E3)/(2*(1+sqrt(Pr23*Pr32)))
end if
end if
C Calculate the Stiffness Matrix
Gm = 1/(1-Pr12*Pr21-Pr23*Pr32-Pr31*Pr13-2*Pr21
1 *Pr32*Pr13)
D1111 = E1*(1-Pr23*Pr32)*Gm
D2222 = E2*(1-Pr13*Pr31)*Gm
D3333 = E3*(1-Pr12*Pr21)*Gm
D1122 = E1*(Pr21+Pr31*Pr23)*Gm
D1133 = E1*(Pr31+Pr21*Pr32)*Gm
D2233 = E2*(Pr32+Pr12*Pr31)*Gm
D1212 = G12
D1313 = G13
D2323 = G23
C Copy the values of strains into temp variables
ES1 = strainInc(k,1)
ES2 = strainInc(k,2)
ES3 = strainInc(k,3)
ES4 = strainInc(k,4)
C Update the stress
stressNew(k,1) = stressOld(k,1)
+D1111*ES1+D1122*ES2+D1133*ES3
stressNew(k,2) = stressOld(k,2)
+D1122*ES1+D2222*ES2+D2233*ES3
stressNew(k,3) = stressOld(k,3)
+D1133*ES1+D2233*ES2+D3333*ES3
stressNew(k,4) = stressOld(k,4)+D1212*ES4*2
C Calculate the total strain in 11 direction based on new stress
values
Str11 = stressOld(k,1)
+D1111*ES1+D1122*ES2+D1133*ES3
Str22 = stressOld(k,2)
+D1122*ES1+D2222*ES2+D2233*ES3
Str33 = stressOld(k,3)
+D1133*ES1+D2233*ES2+D3333*ES3
curStrNew = Str11/E1-Str22/E2*Pr21-
Str33/E3*Pr31
C Update the state variable with current value of strain in 11
direction
stateNew(k,1) = curStrNew
C Update the specific internal energy -
stressPower = half * (
* ( stressOld(k,1)+stressNew(k,1) ) * strainInc(k,1) +
* ( stressOld(k,2)+stressNew(k,2) ) * strainInc(k,2) +
* ( stressOld(k,3)+stressNew(k,3) ) * strainInc(k,3) )
+
* ( stressOld(k,4)+stressNew(k,4) ) * strainInc
(k,4)
enerInternNew(k) = enerInternOld(k)
* + stressPower / density(k)
end do
return
end