Hello dear all,
for my master's thesis I'd like to write a simple VUMAT for a material with nearly zero compression stiffness.
So for a better understanding of VUMAT subroutines I did the ABQ tutorial for a VUMAT with isotropic, isothermal elastic material behaviour and compared the results with the built-in material definition.
There is a discrepancy of almost 35% of the S11 stress. ABQ Manual says that both times the Green-Naghdi stress-rate formulation is used.
Why are the results not the same? After hours of trying I still could not find the solution!
The files are attached below.
I am thankful for any advice!
Greetings,
Johannes
VUMAT.FOR*****************************************************************************
**************************************************************************************
!F90 freecode
!VUMAT with linear elastic material behaviour for steel
subroutine vumat( &
! Read only (unmodifiable)variables -
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, &
! Write only (modifiable) variables -
stressNew, stateNew, enerInternNew, enerInelasNew )
include 'vaba_param.inc'
dimension props(nprops), density(nblock), coordMp(nblock,*), &
charLength(nblock), strainInc(nblock,ndir+nshr), &
relSpinInc(nblock,nshr), tempOld(nblock), &
stretchOld(nblock,ndir+nshr), &
defgradOld(nblock,ndir+nshr+nshr), &
fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr), &
stateOld(nblock,nstatev), enerInternOld(nblock), &
enerInelasOld(nblock), tempNew(nblock), &
stretchNew(nblock,ndir+nshr), &
defgradNew(nblock,ndir+nshr+nshr), &
fieldNew(nblock,nfieldv), &
stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev), &
enerInternNew(nblock), enerInelasNew(nblock)
character*80 cmname
!************************************EXPLANATIONS*************************************
! nblock: Number of material points to be processed in this call to VUMAT
! ndir: Number of direct components in a symmetric tensor.
! nshr: Number of indirect components in a symmetric tensor.
! strainInc (nblock, ndir+nshr): Strain increment tensor at each material point.
! stressOld (nblock, ndir+nshr): Stress tensor at each material point at the
! beginning of the increment.
!
!*************************************USER CODING*************************************
double precision, parameter :: one=1.d0
double precision, parameter :: two=2.d0
double precision, parameter :: three=3.d0
integer :: k
double precision :: e, xnu
! Definition of elastic constants
e=210000
xnu=0.3
twomu=e/(one+xnu)
sixmu=three*twomu
alambda=twomu*(e-twomu)/(sixmu-two*e)
! if(nblock == 1) then
! print *, ndir
! end if
do k=1,nblock
print *,stretchOld(k,:)
! print *, 's11',stressOld(k,1)
! print *, 's22',stressOld(k,2)
! print *, 's33',stressOld(k,3)
! print *, 's12',stressOld(k,4)
trace = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)
stressNew(k,1) = stressOld(k,1) + &
twomu * strainInc(k,1) + alambda * trace
stressNew(k,2) = stressOld(k,2) + &
twomu * strainInc(k,2) + alambda * trace
stressNew(k,3) = stressOld(k,3) + &
twomu * strainInc(k,3) + alambda * trace
stressNew(k,4) = stressOld(k,4) + &
twomu * strainInc(k,4)
if (nshr .gt. 1) then
stressNew(k,5) = stressOld(k,5) + &
twomu * strainInc(k,5)
stressNew(k,6) = stressOld(k,6) + &
twomu * strainInc(k,6)
end if
print *, totalTime
end do
return
end
BAR-ELASTIC.INP*******************************************************************
**********************************************************************************
*Heading
** Job name: bar-elastic Model name: bar-elastic
** Generated by: Abaqus/CAE 6.12-3
*Preprint, echo=NO, model=NO, history=NO, contact=NO
**
** PARTS
**
*Part, name=Part-1
*End Part
**
**
** ASSEMBLY
**
*Assembly, name=Assembly
**
*Instance, name=Part-1-1, part=Part-1
*Node
1, 0., 0., 0.
2, 10., 0., 0.
3, 20., 0., 0.
4, 30., 0., 0.
5, 40., 0., 0.
6, 50., 0., 0.
7, 60., 0., 0.
8, 70., 0., 0.
9, 80., 0., 0.
10, 90., 0., 0.
11, 100., 0., 0.
*Element, type=T3D2
1, 1, 2
2, 2, 3
3, 3, 4
4, 4, 5
5, 5, 6
6, 6, 7
7, 7, 8
8, 8, 9
9, 9, 10
10, 10, 11
*Nset, nset=Encastre-Set
1,
*Nset, nset=Set-2, generate
1, 11, 1
*Elset, elset=Set-2, generate
1, 10, 1
*Nset, nset=Set-3
11,
*Nset, nset=RP-Set
11,
*Nset, nset=Set-5, generate
1, 11, 1
*Elset, elset=Set-5, generate
1, 10, 1
*Nset, nset=Set-6, generate
1, 11, 1
*Elset, elset=Set-6, generate
1, 10, 1
*Nset, nset=Set-7, generate
1, 11, 1
*Elset, elset=Set-7, generate
1, 10, 1
** Section: Stitch-Section
*Solid Section, elset=Set-2, material=elastic
5.,
*End Instance
**
*Node
1, 100., 0., 0.
*Nset, nset=Set-1, instance=Part-1-1
1,
*Nset, nset=Set-2
1,
*Nset, nset=m_Set-4
1,
*Nset, nset=Set-5
1,
*Nset, nset=Set-6
1,
*Surface, type=NODE, name=Part-1-1_Set-3_CNS_, internal
Part-1-1.Set-3, 1.
** Constraint: Constraint-1
*Coupling, constraint name=Constraint-1, ref node=m_Set-4, surface=Part-1-1_Set-3_CNS_
*Kinematic
*End Assembly
*Amplitude, name=Amp-1
0., 0., 2., 1.
**
** MATERIALS
**
*Material, name=elastic
*Density
7.86e-06,
*Elastic
210000., 0.3
**
** BOUNDARY CONDITIONS
**
** Name: Encastre Type: Symmetry/Antisymmetry/Encastre
*Boundary
Part-1-1.Encastre-Set, ENCASTRE
** ----------------------------------------------------------------
**
** STEP: Tensile
**
*Step, name=Tensile
*Dynamic, Explicit
, 2.
*Bulk Viscosity
0.06, 1.2
**
** BOUNDARY CONDITIONS
**
** Name: tension-u1 Type: Displacement/Rotation
*Boundary, amplitude=Amp-1
Part-1-1.RP-Set, 1, 1, 1.
**
** OUTPUT REQUESTS
**
*Restart, write, number interval=1, time marks=NO
**
** FIELD OUTPUT: F-Output-1
**
*Output, field
*Node Output
A, RF, U, V
*Element Output, directions=YES
EVF, LE, NE, PE, PEEQ, PEEQVAVG, PEVAVG, S, SVAVG
*Contact Output
CSTRESS,
**
** HISTORY OUTPUT: RP-Output
**
*Output, history
*Node Output, nset=Part-1-1.RP-Set
RF1, U1
**
** HISTORY OUTPUT: H-Output-1
**
*Output, history, variable=PRESELECT
*End Step
BAR-VUMAT.INP************************************************************
*************************************************************************
*Heading
** Job name: bar-vumat Model name: bar-vumat
** Generated by: Abaqus/CAE 6.12-3
*Preprint, echo=NO, model=NO, history=NO, contact=NO
**
** PARTS
**
*Part, name=Part-1
*End Part
**
**
** ASSEMBLY
**
*Assembly, name=Assembly
**
*Instance, name=Part-1-1, part=Part-1
*Node
1, 0., 0., 0.
2, 100., 0., 0.
*Element, type=T3D2
1, 1, 2
*Nset, nset=Encastre-Set
1,
*Nset, nset=Set-2
1, 2
*Elset, elset=Set-2
1,
*Nset, nset=Set-3
2,
*Nset, nset=RP-Set
2,
*Nset, nset=Set-5
1, 2
*Elset, elset=Set-5
1,
*Nset, nset=Set-6
1, 2
*Elset, elset=Set-6
1,
*Nset, nset=Set-7
1, 2
*Elset, elset=Set-7
1,
** Section: Stitch-Section
*Solid Section, elset=Set-2, material=vumat-el
5.,
*End Instance
**
*Node
1, 100., 0., 0.
*Nset, nset=Set-1, instance=Part-1-1
1,
*Nset, nset=Set-2
1,
*Nset, nset=m_Set-4
1,
*Nset, nset=Set-5
1,
*Nset, nset=Set-6
1,
*Surface, type=NODE, name=Part-1-1_Set-3_CNS_, internal
Part-1-1.Set-3, 1.
** Constraint: Constraint-1
*Coupling, constraint name=Constraint-1, ref node=m_Set-4, surface=Part-1-1_Set-3_CNS_
*Kinematic
*End Assembly
*Amplitude, name=Amp-1
0., 0., 2., 1.
**
** MATERIALS
**
*Material, name=vumat-el
*Density
7.86e-06,
*User Material, constants=2
210000., 0.3
**
** BOUNDARY CONDITIONS
**
** Name: Encastre Type: Symmetry/Antisymmetry/Encastre
*Boundary
Part-1-1.Encastre-Set, ENCASTRE
** ----------------------------------------------------------------
**
** STEP: Tensile
**
*Step, name=Tensile
*Dynamic, Explicit
, 2.
*Bulk Viscosity
0.06, 1.2
**
** BOUNDARY CONDITIONS
**
** Name: tension-u1 Type: Displacement/Rotation
*Boundary, amplitude=Amp-1
Part-1-1.RP-Set, 1, 1, 1.
**
** OUTPUT REQUESTS
**
*Restart, write, number interval=1, time marks=NO
**
** FIELD OUTPUT: F-Output-1
**
*Output, field
*Node Output
A, RF, U, V
*Element Output, directions=YES
EVF, LE, NE, PE, PEEQ, PEEQVAVG, PEVAVG, S, SVAVG
*Contact Output
CSTRESS,
**
** HISTORY OUTPUT: RP-Output
**
*Output, history
*Node Output, nset=Part-1-1.RP-Set
RF1, U1
**
** HISTORY OUTPUT: H-Output-1
**
*Output, history, variable=PRESELECT
*End Step
Mismatch between built-in linear-elastic material behaviour and tutorial VUMAT in ABQ/Explicit using truss elements
Forums