anybody has experience on VUMAT of hyperelastic mat'l?

Discussion:

M. Qian

2004-02-06 07:33:34 UTC

hey, guyz

i am modeling large deformation with VUMAT, for 1 element tension test, VUMAT looks fine. but for 1 element shearing, solution has large variance (>10%) from built-in model. for multi-element tests, most cases can not be completed. Also when stepTime=0, I always got negative strain even it is tension test, then strain came back to positive in the next steps.
any reason for that?
kind wondering if i need add *orientation in my inp file due to large deformation? or why I am having such results?

any advices are appreciated,

john

---------------------------------
Do you Yahoo!?
Yahoo! Finance: Get your refund fast by filing online

[Non-text portions of this message have been removed]

------------------------ Yahoo! Groups Sponsor ---------------------~-->
Buy Ink Cartridges or Refill Kits for your HP, Epson, Canon or Lexmark
Printer at MyInks.com. Free s/h on orders $50 or more to the US & Canada.
http://www.c1tracking.com/l.asp?cid=5511
http://us.click.yahoo.com/mOAaAA/3exGAA/qnsNAA/PMYolB/TM
---------------------------------------------------------------------~->

M. Qian

2004-02-08 22:53:06 UTC

Permalink

hi,

I am attaching the VUMAT subroutine of Neo-hookean model.

Outputs from it are different from built-in neo-hookean model even for same C10 and D1, especially under large shear strains.

Could not figure it out by myself. hope someone can help me take a look.

Really thanks.

john

+++++++++++++++++

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
dimension devia(nblock,ndir+nshr),
1 BBar(nblock,4), defgradNewBar(nblock,5), intv(2)
C
character*80 cmname
parameter (zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0,
* four = 4.d0)
C
intv(1) = ndir
intv(2) = nshr
C
C if (ndir .ne. 3 .or. nshr .ne. 1) then
C call xplb_abqerr(1,'Subroutine VUMAT is implemented '//
C * 'only for plane strain and axisymmetric cases '//
C * '(ndir=3 and nshr=1)',0,zero,' ')
C call xplb_abqerr(-2,'Subroutine VUMAT has been called '//
C * 'with ndir=%I and nshr=%I',intv,zero,' ')
C call xplb_exit
C end if
C
C10 = props(1)
D1 = props(2)
C C10=1.11619E6 D1=4.48E-8
C
ak=two/D1
twomu=four*C10
amu=two*C10
alamda=(three*ak-twomu)/three
C
C if stepTime equals zero, assume pure elastic material and use initial elastic modulus
C
if(stepTime .EQ. zero) then
do k=1,nblock
trace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)
stressNew(k,1) = stressOld(k,1)
* + twomu * strainInc(k,1) + alamda * trace1
stressNew(k,2) = stressOld(k,2)
* + twomu * strainInc(k,2) + alamda * trace1
stressNew(k,3) = stressOld(k,3)
* + twomu * strainInc(k,3) + alamda * trace1
stressNew(k,4) = stressOld(k,4)
* + twomu * strainInc(k,4)
write(6,*) stepTime,k,stressNew(k,1),trace1
end do
else
C
C JACOBIAN AND DISTORTION TENSOR (F is non-symmetric)
C
do k=1,nblock
det=defgradNew(k, 1)*defgradNew(k, 2)*defgradNew(k, 3) -
1 defgradNew(k, 3)*defgradNew(k, 4)*defgradNew(k, 5)
scale=det**(-ONE/THREE)
defgradNewBar(k, 1)=defgradNew(k, 1)*scale
defgradNewBar(k, 2)=defgradNew(k, 2)*scale
defgradNewBar(k, 3)=defgradNew(k, 3)*scale
defgradNewBar(k, 4)=defgradNew(k, 4)*scale
defgradNewBar(k, 5)=defgradNew(k, 5)*scale
C
C CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)
C
BBar(k,1)=defgradNewBar(k, 1)**2.+defgradNewBar(k, 4)**2.
BBar(k,2)=defgradNewBar(k, 2)**2.+defgradNewBar(k, 5)**2.
BBar(k,3)=defgradNewBar(k, 3)**2.
BBar(k,4)=defgradNewBar(k, 1)*defgradNewBar(k, 5) +
1 defgradNewBar(k, 2)*defgradNewBar(k, 4)
C
C CALCULATE STRESS tensor
C
TRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)
EG=two*C10/det
PR=TWO/D1*(det-one)
stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PR
stressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PR
stressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PR
stressNew(k,4)=EG *BBar(k,4)
C
write(6,*) stepTime,k,stressNew(k,1),stressNew(k,2),
1 stressNew(k,3),stressNew(k,4)
end do
end if
return
end

"M. Qian" <***@yahoo.com> wrote:
hey, guyz

i am modeling large deformation with VUMAT, for 1 element tension test, VUMAT looks fine. but for 1 element shearing, solution has large variance (>10%) from built-in model. for multi-element tests, most cases can not be completed. Also when stepTime=0, I always got negative strain even it is tension test, then strain came back to positive in the next steps.
any reason for that?
kind wondering if i need add *orientation in my inp file due to large deformation? or why I am having such results?

any advices are appreciated,

john

---------------------------------
Do you Yahoo!?
Yahoo! Finance: Get your refund fast by filing online

[Non-text portions of this message have been removed]

benjamin_hagege

2004-02-09 07:16:25 UTC

Permalink

If you use solid elements, the formulation is a Green-Naghdi one (so
it is hypoelastic). This is why you have not the same values when
compared to your hyperelastic model !

BenZ.

Post by M. Qian
hi,
I am attaching the VUMAT subroutine of Neo-hookean model.
Outputs from it are different from built-in neo-hookean model even

for same C10 and D1, especially under large shear strains.

Post by M. Qian
Could not figure it out by myself. hope someone can help me take a look.
Really thanks.
john
+++++++++++++++++
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),

Post by M. Qian
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
dimension devia(nblock,ndir+nshr),
1 BBar(nblock,4), defgradNewBar(nblock,5), intv(2)
C
character*80 cmname
parameter (zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0,
* four = 4.d0)
C
intv(1) = ndir
intv(2) = nshr
C
C if (ndir .ne. 3 .or. nshr .ne. 1) then
C call xplb_abqerr(1,'Subroutine VUMAT is implemented '//
C * 'only for plane strain and axisymmetric cases '//
C * '(ndir=3 and nshr=1)',0,zero,' ')
C call xplb_abqerr(-2,'Subroutine VUMAT has been called '//
C * 'with ndir=%I and nshr=%I',intv,zero,' ')
C call xplb_exit
C end if
C
C10 = props(1)
D1 = props(2)
C C10=1.11619E6 D1=4.48E-8
C
ak=two/D1
twomu=four*C10
amu=two*C10
alamda=(three*ak-twomu)/three
C
C if stepTime equals zero, assume pure elastic material and use initial elastic modulus
C
if(stepTime .EQ. zero) then
do k=1,nblock
trace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)
stressNew(k,1) = stressOld(k,1)
* + twomu * strainInc(k,1) + alamda * trace1
stressNew(k,2) = stressOld(k,2)
* + twomu * strainInc(k,2) + alamda * trace1
stressNew(k,3) = stressOld(k,3)
* + twomu * strainInc(k,3) + alamda * trace1
stressNew(k,4) = stressOld(k,4)
* + twomu * strainInc(k,4)
write(6,*) stepTime,k,stressNew(k,1),trace1
end do
else
C
C JACOBIAN AND DISTORTION TENSOR (F is non-symmetric)
C
do k=1,nblock
det=defgradNew(k, 1)*defgradNew(k, 2)*defgradNew(k, 3) -
1 defgradNew(k, 3)*defgradNew(k, 4)*defgradNew(k, 5)
scale=det**(-ONE/THREE)
defgradNewBar(k, 1)=defgradNew(k, 1)*scale
defgradNewBar(k, 2)=defgradNew(k, 2)*scale
defgradNewBar(k, 3)=defgradNew(k, 3)*scale
defgradNewBar(k, 4)=defgradNew(k, 4)*scale
defgradNewBar(k, 5)=defgradNew(k, 5)*scale
C
C CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)
C
BBar(k,1)=defgradNewBar(k, 1)**2.+defgradNewBar(k, 4)**2.
BBar(k,2)=defgradNewBar(k, 2)**2.+defgradNewBar(k, 5)**2.
BBar(k,3)=defgradNewBar(k, 3)**2.
BBar(k,4)=defgradNewBar(k, 1)*defgradNewBar(k, 5) +
1 defgradNewBar(k, 2)*defgradNewBar(k, 4)
C
C CALCULATE STRESS tensor
C
TRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)
EG=two*C10/det
PR=TWO/D1*(det-one)
stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PR
stressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PR
stressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PR
stressNew(k,4)=EG *BBar(k,4)
C
write(6,*) stepTime,k,stressNew(k,1),stressNew(k,2),
1 stressNew(k,3),stressNew(k,4)
end do
end if
return
end
hey, guyz
i am modeling large deformation with VUMAT, for 1 element tension

test, VUMAT looks fine. but for 1 element shearing, solution has
large variance (>10%) from built-in model. for multi-element tests,
most cases can not be completed. Also when stepTime=0, I always got
negative strain even it is tension test, then strain came back to
positive in the next steps.

Post by M. Qian
any reason for that?
kind wondering if i need add *orientation in my inp file due to

large deformation? or why I am having such results?

Post by M. Qian
any advices are appreciated,
john
---------------------------------
Do you Yahoo!?
Yahoo! Finance: Get your refund fast by filing online
[Non-text portions of this message have been removed]
http://groups.yahoo.com/group/abaqus
Yahoo! Groups SponsorADVERTISEMENT
---------------------------------
Yahoo! Groups Links
http://groups.yahoo.com/group/ABAQUS/
Your use of Yahoo! Groups is subject to the Yahoo! Terms of

Service.

Post by M. Qian
---------------------------------
Do you Yahoo!?
Yahoo! Finance: Get your refund fast by filing online
[Non-text portions of this message have been removed]

------------------------ Yahoo! Groups Sponsor ---------------------~-->
Buy Ink Cartridges or Refill Kits for your HP, Epson, Canon or Lexmark
Printer at MyInks.com. Free s/h on orders $50 or more to the US & Canada.
http://www.c1tracking.com/l.asp?cid=5511
http://us.click.yahoo.com/mOAaAA/3exGAA/qnsNAA/PMYolB/TM
---------------------------------------------------------------------~->

M. Qian

2004-02-12 04:51:07 UTC

Permalink

Hi, Benz and other workmates,

I am looking at ABAQUS manuals carefully these days. Only found 2 locations with Green-Naghdi, no one explained it clearly.
In Explicit, the example after VUMAT does not apply any rotation in it, I assume the example is good for large deformations of solid elements,including rotations. So where does the rotation go? sth. need do with relSpinInc(*)=w-omega? How?

Thanks for any lights,

John

Btw, Benz, why G-N strain rate is related to hypoelasticity ? where I can find these concepts? Thanks.

benjamin_hagege <***@yahoo.fr> wrote:
If you use solid elements, the formulation is a Green-Naghdi one (so
it is hypoelastic). This is why you have not the same values when
compared to your hyperelastic model !

BenZ.

Post by M. Qian
hi,
I am attaching the VUMAT subroutine of Neo-hookean model.
Outputs from it are different from built-in neo-hookean model even

for same C10 and D1, especially under large shear strains.

(nprops),

Post by M. Qian
any reason for that?
kind wondering if i need add *orientation in my inp file due to

large deformation? or why I am having such results?

Service.

Post by M. Qian
---------------------------------
Do you Yahoo!?
Yahoo! Finance: Get your refund fast by filing online
[Non-text portions of this message have been removed]

benjamin_hagege

2004-02-12 13:33:52 UTC

Permalink

Heyhey, here you are in complex issues !!!

1/the rotation doesn't appear in the vumat because the tensorial
equations are projected in the rotating frame ! So they takes the
same form as if you were in small perturbations. You can not choose
the rotation, it is the Green-Naghdi's one, I.E. R=F*Inv(U).

2/You do not need to use the realspin, just apply the last formula if
you wanna get the rotation of the Green-Naghdi's frame.

3/You can find these concepts in my PhD thesis... that is still not
written sorry ! :) Try to get the papers of Dafalias or Nagtegaal
that compares the Green-Naghdi formulations with the Jaumann one, and
get also the fundamental paper of Hughes and Winget in IJNME (???
when large rotations arising).

Post by M. Qian
Hi, Benz and other workmates,
I am looking at ABAQUS manuals carefully these days. Only found 2

locations with Green-Naghdi, no one explained it clearly.

Post by M. Qian
In Explicit, the example after VUMAT does not apply any rotation in

it, I assume the example is good for large deformations of solid
elements,including rotations. So where does the rotation go? sth.
need do with relSpinInc(*)=w-omega? How?

Post by M. Qian
Thanks for any lights,
John
Btw, Benz, why G-N strain rate is related to hypoelasticity ? where

I can find these concepts? Thanks.

Post by M. Qian
If you use solid elements, the formulation is a Green-Naghdi one (so
it is hypoelastic). This is why you have not the same values when
compared to your hyperelastic model !
BenZ.

Post by M. Qian
hi,
I am attaching the VUMAT subroutine of Neo-hookean model.
Outputs from it are different from built-in neo-hookean model even

for same C10 and D1, especially under large shear strains.

Post by M. Qian
Could not figure it out by myself. hope someone can help me take

Post by M. Qian
look.

Post by M. Qian
Really thanks.
john
+++++++++++++++++
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),

called '//

Post by M. Qian

Post by M. Qian
C * 'with ndir=%I and nshr=%I',intv,zero,' ')
C call xplb_exit
C end if
C
C10 = props(1)
D1 = props(2)
C C10=1.11619E6 D1=4.48E-8
C
ak=two/D1
twomu=four*C10
amu=two*C10
alamda=(three*ak-twomu)/three
C
C if stepTime equals zero, assume pure elastic material and use

initial elastic modulus

Post by M. Qian
C
if(stepTime .EQ. zero) then
do k=1,nblock
trace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)
stressNew(k,1) = stressOld(k,1)
* + twomu * strainInc(k,1) + alamda * trace1
stressNew(k,2) = stressOld(k,2)
* + twomu * strainInc(k,2) + alamda * trace1
stressNew(k,3) = stressOld(k,3)
* + twomu * strainInc(k,3) + alamda * trace1
stressNew(k,4) = stressOld(k,4)
* + twomu * strainInc(k,4)
write(6,*) stepTime,k,stressNew(k,1),trace1
end do
else
C
C JACOBIAN AND DISTORTION TENSOR (F is non-symmetric)
C
do k=1,nblock
det=defgradNew(k, 1)*defgradNew(k, 2)*defgradNew(k, 3) -
1 defgradNew(k, 3)*defgradNew(k, 4)*defgradNew(k, 5)
scale=det**(-ONE/THREE)
defgradNewBar(k, 1)=defgradNew(k, 1)*scale
defgradNewBar(k, 2)=defgradNew(k, 2)*scale
defgradNewBar(k, 3)=defgradNew(k, 3)*scale
defgradNewBar(k, 4)=defgradNew(k, 4)*scale
defgradNewBar(k, 5)=defgradNew(k, 5)*scale
C
C CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)
C
BBar(k,1)=defgradNewBar(k, 1)**2.+defgradNewBar(k, 4)**2.
BBar(k,2)=defgradNewBar(k, 2)**2.+defgradNewBar(k, 5)**2.
BBar(k,3)=defgradNewBar(k, 3)**2.
BBar(k,4)=defgradNewBar(k, 1)*defgradNewBar(k, 5) +
1 defgradNewBar(k, 2)*defgradNewBar(k, 4)
C
C CALCULATE STRESS tensor
C
TRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)
EG=two*C10/det
PR=TWO/D1*(det-one)
stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PR
stressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PR
stressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PR
stressNew(k,4)=EG *BBar(k,4)
C
write(6,*) stepTime,k,stressNew(k,1),stressNew(k,2),
1 stressNew(k,3),stressNew(k,4)
end do
end if
return
end
hey, guyz
i am modeling large deformation with VUMAT, for 1 element tension

Post by M. Qian
any reason for that?
kind wondering if i need add *orientation in my inp file due to

large deformation? or why I am having such results?

Service.

Post by M. Qian
---------------------------------
Do you Yahoo!?
Yahoo! Finance: Get your refund fast by filing online
[Non-text portions of this message have been removed]

http://groups.yahoo.com/group/abaqus
Yahoo! Groups SponsorADVERTISEMENT
---------------------------------
Yahoo! Groups Links
http://groups.yahoo.com/group/ABAQUS/
Your use of Yahoo! Groups is subject to the Yahoo! Terms of

Service.

Post by M. Qian
---------------------------------
Do you Yahoo!?
Yahoo! Finance: Get your refund fast by filing online
[Non-text portions of this message have been removed]

M. Qian

2004-02-15 21:31:56 UTC

Permalink

Hi Benz,

in VUMAT I think total formulation is used instead rate formulation.
So we need not consider either Jaumann or Green-Naghdi rate, should we?
Similar to the neo-hookean example in UMAT, I used the same formula which is
sigma=2*C10/J*(B_bar-1/3*trace(B_bar)-I)+2/D1(J-1)*I in VUMAT, but my results using this VUMAT is different from ABAQUS built-in neo-hookean model for same C10 and D1 values.
We should not consider "stress rate" or "rotation" effects, because using total formulation and corotational coordinated used in VUMAT.
Could the formula itself be wrong? or VUMAT has different format?

Pls give lights!

John

benjamin_hagege <***@yahoo.fr> wrote:
Heyhey, here you are in complex issues !!!

1/the rotation doesn't appear in the vumat because the tensorial
equations are projected in the rotating frame ! So they takes the
same form as if you were in small perturbations. You can not choose
the rotation, it is the Green-Naghdi's one, I.E. R=F*Inv(U).

2/You do not need to use the realspin, just apply the last formula if
you wanna get the rotation of the Green-Naghdi's frame.

3/You can find these concepts in my PhD thesis... that is still not
written sorry ! :) Try to get the papers of Dafalias or Nagtegaal
that compares the Green-Naghdi formulations with the Jaumann one, and
get also the fundamental paper of Hughes and Winget in IJNME (???
when large rotations arising).

Post by M. Qian
Hi, Benz and other workmates,
I am looking at ABAQUS manuals carefully these days. Only found 2

locations with Green-Naghdi, no one explained it clearly.

Post by M. Qian
In Explicit, the example after VUMAT does not apply any rotation in

it, I assume the example is good for large deformations of solid
elements,including rotations. So where does the rotation go? sth.
need do with relSpinInc(*)=w-omega? How?

Post by M. Qian
Thanks for any lights,
John
Btw, Benz, why G-N strain rate is related to hypoelasticity ? where

I can find these concepts? Thanks.

Post by M. Qian
If you use solid elements, the formulation is a Green-Naghdi one (so
it is hypoelastic). This is why you have not the same values when
compared to your hyperelastic model !
BenZ.

Post by M. Qian
hi,
I am attaching the VUMAT subroutine of Neo-hookean model.
Outputs from it are different from built-in neo-hookean model even

for same C10 and D1, especially under large shear strains.

Post by M. Qian
Could not figure it out by myself. hope someone can help me take

Post by M. Qian
look.

(nprops),

called '//

Post by M. Qian

initial elastic modulus

Post by M. Qian
any reason for that?
kind wondering if i need add *orientation in my inp file due to

large deformation? or why I am having such results?

Service.

Post by M. Qian
---------------------------------
Do you Yahoo!?
Yahoo! Finance: Get your refund fast by filing online
[Non-text portions of this message have been removed]

Service.

Post by M. Qian
---------------------------------
Do you Yahoo!?
Yahoo! Finance: Get your refund fast by filing online
[Non-text portions of this message have been removed]

hariton ilia

2004-02-10 06:44:17 UTC

Permalink

Hi,

Had no time to look at your subroutine but for isotropic hyperelasticity (invariants formulation) you can use UHYPER it is much easier. All you have to do is to insert the strain potential and its derivatives.

Ilia
----- Original Message -----
From: M. Qian
To: ***@yahoogroups.com
Sent: Monday, February 09, 2004 00:53
Subject: [ABAQUS] Who can help me on VUMAT of hyperelasticity ?

hi,

I am attaching the VUMAT subroutine of Neo-hookean model.

Outputs from it are different from built-in neo-hookean model even for same C10 and D1, especially under large shear strains.

Could not figure it out by myself. hope someone can help me take a look.

Really thanks.

john

+++++++++++++++++

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
dimension devia(nblock,ndir+nshr),
1 BBar(nblock,4), defgradNewBar(nblock,5), intv(2)
C
character*80 cmname
parameter (zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0,
* four = 4.d0)
C
intv(1) = ndir
intv(2) = nshr
C
C if (ndir .ne. 3 .or. nshr .ne. 1) then
C call xplb_abqerr(1,'Subroutine VUMAT is implemented '//
C * 'only for plane strain and axisymmetric cases '//
C * '(ndir=3 and nshr=1)',0,zero,' ')
C call xplb_abqerr(-2,'Subroutine VUMAT has been called '//
C * 'with ndir=%I and nshr=%I',intv,zero,' ')
C call xplb_exit
C end if
C
C10 = props(1)
D1 = props(2)
C C10=1.11619E6 D1=4.48E-8
C
ak=two/D1
twomu=four*C10
amu=two*C10
alamda=(three*ak-twomu)/three
C
C if stepTime equals zero, assume pure elastic material and use initial elastic modulus
C
if(stepTime .EQ. zero) then
do k=1,nblock
trace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)
stressNew(k,1) = stressOld(k,1)
* + twomu * strainInc(k,1) + alamda * trace1
stressNew(k,2) = stressOld(k,2)
* + twomu * strainInc(k,2) + alamda * trace1
stressNew(k,3) = stressOld(k,3)
* + twomu * strainInc(k,3) + alamda * trace1
stressNew(k,4) = stressOld(k,4)
* + twomu * strainInc(k,4)
write(6,*) stepTime,k,stressNew(k,1),trace1
end do
else
C
C JACOBIAN AND DISTORTION TENSOR (F is non-symmetric)
C
do k=1,nblock
det=defgradNew(k, 1)*defgradNew(k, 2)*defgradNew(k, 3) -
1 defgradNew(k, 3)*defgradNew(k, 4)*defgradNew(k, 5)
scale=det**(-ONE/THREE)
defgradNewBar(k, 1)=defgradNew(k, 1)*scale
defgradNewBar(k, 2)=defgradNew(k, 2)*scale
defgradNewBar(k, 3)=defgradNew(k, 3)*scale
defgradNewBar(k, 4)=defgradNew(k, 4)*scale
defgradNewBar(k, 5)=defgradNew(k, 5)*scale
C
C CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)
C
BBar(k,1)=defgradNewBar(k, 1)**2.+defgradNewBar(k, 4)**2.
BBar(k,2)=defgradNewBar(k, 2)**2.+defgradNewBar(k, 5)**2.
BBar(k,3)=defgradNewBar(k, 3)**2.
BBar(k,4)=defgradNewBar(k, 1)*defgradNewBar(k, 5) +
1 defgradNewBar(k, 2)*defgradNewBar(k, 4)
C
C CALCULATE STRESS tensor
C
TRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)
EG=two*C10/det
PR=TWO/D1*(det-one)
stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PR
stressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PR
stressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PR
stressNew(k,4)=EG *BBar(k,4)
C
write(6,*) stepTime,k,stressNew(k,1),stressNew(k,2),
1 stressNew(k,3),stressNew(k,4)
end do
end if
return
end

"M. Qian" <***@yahoo.com> wrote:
hey, guyz

i am modeling large deformation with VUMAT, for 1 element tension test, VUMAT looks fine. but for 1 element shearing, solution has large variance (>10%) from built-in model. for multi-element tests, most cases can not be completed. Also when stepTime=0, I always got negative strain even it is tension test, then strain came back to positive in the next steps.
any reason for that?
kind wondering if i need add *orientation in my inp file due to large deformation? or why I am having such results?

any advices are appreciated,

john

---------------------------------
Do you Yahoo!?
Yahoo! Finance: Get your refund fast by filing online

[Non-text portions of this message have been removed]

Community email addresses:
Post message: ***@yahoogroups.com
Subscribe: ABAQUS-***@yahoogroups.com
Unsubscribe: ABAQUS-***@yahoogroups.com
List owner: ABAQUS-***@yahoogroups.com

Shortcut URL to this page:
http://groups.yahoo.com/group/abaqus

Yahoo! Groups SponsorADVERTISEMENT

---------------------------------
Yahoo! Groups Links

To visit your group on the web, go to:
http://groups.yahoo.com/group/ABAQUS/

To unsubscribe from this group, send an email to:
ABAQUS-***@yahoogroups.com

Your use of Yahoo! Groups is subject to the Yahoo! Terms of Service.

---------------------------------
Do you Yahoo!?
Yahoo! Finance: Get your refund fast by filing online

[Non-text portions of this message have been removed]

Community email addresses:
Post message: ***@yahoogroups.com
Subscribe: ABAQUS-***@yahoogroups.com
Unsubscribe: ABAQUS-***@yahoogroups.com
List owner: ABAQUS-***@yahoogroups.com

Shortcut URL to this page:
http://groups.yahoo.com/group/abaqus

Yahoo! Groups Sponsor
ADVERTISEMENT

------------------------------------------------------------------------------
Yahoo! Groups Links

a.. To visit your group on the web, go to:
http://groups.yahoo.com/group/ABAQUS/

b.. To unsubscribe from this group, send an email to:
ABAQUS-***@yahoogroups.com

c.. Your use of Yahoo! Groups is subject to the Yahoo! Terms of Service.

[Non-text portions of this message have been removed]

------------------------ Yahoo! Groups Sponsor ---------------------~-->
Buy Ink Cartridges or Refill Kits for your HP, Epson, Canon or Lexmark
Printer at MyInks.com. Free s/h on orders $50 or more to the US & Canada.
http://www.c1tracking.com/l.asp?cid=5511
http://us.click.yahoo.com/mOAaAA/3exGAA/qnsNAA/PMYolB/TM
---------------------------------------------------------------------~->

M. Qian

2004-02-10 16:26:42 UTC

Permalink

Thanks for reply, UHYPER is only available in standard.
I want to use VUMAT which is explicit version.

hariton ilia <***@bgumail.bgu.ac.il> wrote:
Hi,

Had no time to look at your subroutine but for isotropic hyperelasticity (invariants formulation) you can use UHYPER it is much easier. All you have to do is to insert the strain potential and its derivatives.

Ilia
----- Original Message -----
From: M. Qian
To: ***@yahoogroups.com
Sent: Monday, February 09, 2004 00:53
Subject: [ABAQUS] Who can help me on VUMAT of hyperelasticity ?

hi,

I am attaching the VUMAT subroutine of Neo-hookean model.

Outputs from it are different from built-in neo-hookean model even for same C10 and D1, especially under large shear strains.

Could not figure it out by myself. hope someone can help me take a look.

Really thanks.

john

+++++++++++++++++

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
dimension devia(nblock,ndir+nshr),
1 BBar(nblock,4), defgradNewBar(nblock,5), intv(2)
C
character*80 cmname
parameter (zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0,
* four = 4.d0)
C
intv(1) = ndir
intv(2) = nshr
C
C if (ndir .ne. 3 .or. nshr .ne. 1) then
C call xplb_abqerr(1,'Subroutine VUMAT is implemented '//
C * 'only for plane strain and axisymmetric cases '//
C * '(ndir=3 and nshr=1)',0,zero,' ')
C call xplb_abqerr(-2,'Subroutine VUMAT has been called '//
C * 'with ndir=%I and nshr=%I',intv,zero,' ')
C call xplb_exit
C end if
C
C10 = props(1)
D1 = props(2)
C C10=1.11619E6 D1=4.48E-8
C
ak=two/D1
twomu=four*C10
amu=two*C10
alamda=(three*ak-twomu)/three
C
C if stepTime equals zero, assume pure elastic material and use initial elastic modulus
C
if(stepTime .EQ. zero) then
do k=1,nblock
trace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)
stressNew(k,1) = stressOld(k,1)
* + twomu * strainInc(k,1) + alamda * trace1
stressNew(k,2) = stressOld(k,2)
* + twomu * strainInc(k,2) + alamda * trace1
stressNew(k,3) = stressOld(k,3)
* + twomu * strainInc(k,3) + alamda * trace1
stressNew(k,4) = stressOld(k,4)
* + twomu * strainInc(k,4)
write(6,*) stepTime,k,stressNew(k,1),trace1
end do
else
C
C JACOBIAN AND DISTORTION TENSOR (F is non-symmetric)
C
do k=1,nblock
det=defgradNew(k, 1)*defgradNew(k, 2)*defgradNew(k, 3) -
1 defgradNew(k, 3)*defgradNew(k, 4)*defgradNew(k, 5)
scale=det**(-ONE/THREE)
defgradNewBar(k, 1)=defgradNew(k, 1)*scale
defgradNewBar(k, 2)=defgradNew(k, 2)*scale
defgradNewBar(k, 3)=defgradNew(k, 3)*scale
defgradNewBar(k, 4)=defgradNew(k, 4)*scale
defgradNewBar(k, 5)=defgradNew(k, 5)*scale
C
C CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)
C
BBar(k,1)=defgradNewBar(k, 1)**2.+defgradNewBar(k, 4)**2.
BBar(k,2)=defgradNewBar(k, 2)**2.+defgradNewBar(k, 5)**2.
BBar(k,3)=defgradNewBar(k, 3)**2.
BBar(k,4)=defgradNewBar(k, 1)*defgradNewBar(k, 5) +
1 defgradNewBar(k, 2)*defgradNewBar(k, 4)
C
C CALCULATE STRESS tensor
C
TRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)
EG=two*C10/det
PR=TWO/D1*(det-one)
stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PR
stressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PR
stressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PR
stressNew(k,4)=EG *BBar(k,4)
C
write(6,*) stepTime,k,stressNew(k,1),stressNew(k,2),
1 stressNew(k,3),stressNew(k,4)
end do
end if
return
end

"M. Qian" <***@yahoo.com> wrote:
hey, guyz

i am modeling large deformation with VUMAT, for 1 element tension test, VUMAT looks fine. but for 1 element shearing, solution has large variance (>10%) from built-in model. for multi-element tests, most cases can not be completed. Also when stepTime=0, I always got negative strain even it is tension test, then strain came back to positive in the next steps.
any reason for that?
kind wondering if i need add *orientation in my inp file due to large deformation? or why I am having such results?

any advices are appreciated,

john

---------------------------------
Do you Yahoo!?
Yahoo! Finance: Get your refund fast by filing online

[Non-text portions of this message have been removed]

Community email addresses:
Post message: ***@yahoogroups.com
Subscribe: ABAQUS-***@yahoogroups.com
Unsubscribe: ABAQUS-***@yahoogroups.com
List owner: ABAQUS-***@yahoogroups.com

Shortcut URL to this page:
http://groups.yahoo.com/group/abaqus

Yahoo! Groups SponsorADVERTISEMENT

---------------------------------
Yahoo! Groups Links

To visit your group on the web, go to:
http://groups.yahoo.com/group/ABAQUS/

To unsubscribe from this group, send an email to:
ABAQUS-***@yahoogroups.com

Your use of Yahoo! Groups is subject to the Yahoo! Terms of Service.

---------------------------------
Do you Yahoo!?
Yahoo! Finance: Get your refund fast by filing online

[Non-text portions of this message have been removed]

Community email addresses:
Post message: ***@yahoogroups.com
Subscribe: ABAQUS-***@yahoogroups.com
Unsubscribe: ABAQUS-***@yahoogroups.com
List owner: ABAQUS-***@yahoogroups.com

Shortcut URL to this page:
http://groups.yahoo.com/group/abaqus

Yahoo! Groups Sponsor
ADVERTISEMENT

------------------------------------------------------------------------------
Yahoo! Groups Links

a.. To visit your group on the web, go to:
http://groups.yahoo.com/group/ABAQUS/

b.. To unsubscribe from this group, send an email to:
ABAQUS-***@yahoogroups.com

c.. Your use of Yahoo! Groups is subject to the Yahoo! Terms of Service.

[Non-text portions of this message have been removed]

l_xianwu

2007-02-24 22:41:09 UTC

Permalink

Hi, Qian,

I have now met exactly the same problem as you did in using VUMAT for
Ogden's model. One of my friends also met the same problem before. He
suggested to use U as F. This suggestion is also made by the latest
Abaqus (v6.51) VUMAT document. I did. It kind of works, but does not
solve my problem totally.

Here is what I did. At the beginning of the VUMAT, I used U as F. At
the end of VUMAT, I wrote the stressnew to unit 6, which I can then
use to plot the stress data using Tecplot. This stressnew data is
similar to my FEAP results for exactly the same model. On the other
hand, the stress data (say, S12) that Abaqus calculates (which I plot
using Abaqus CAE) are quite different from the stressnew that I myself
output at the end of VUMAT.

Have you solved your problem and how? Thanks a lot!

Xianwu

Post by M. Qian
hey, guyz
i am modeling large deformation with VUMAT, for 1 element tension

test, VUMAT looks fine. but for 1 element shearing, solution has large
variance (>10%) from built-in model. for multi-element tests, most
cases can not be completed. Also when stepTime=0, I always got
negative strain even it is tension test, then strain came back to
positive in the next steps.

Post by M. Qian
any reason for that?
kind wondering if i need add *orientation in my inp file due to

large deformation? or why I am having such results?