Announcement

Collapse
No announcement yet.

VUMAT for maxwell zener model-validation

Collapse
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • VUMAT for maxwell zener model-validation

    Hello all ,

    Here is the VUMAT I wrote for a standard linear solid (maxwell-zener model obtained from Mr.Frank Richter's thesis) .When I tried to test it on a single element ,though the code runs fine giving accurate values of stresses etc. , the stress and strain curves are not smooth.Attached is my code ,the input file I used and the plots for stresses and strains .

    Please let me know if it is a valid response from a viscoelastic material.


    Best Regards,
    Sireesha

    VUMAT


    subroutine vumat(nblock,ndir,nshr,nstatev,nfieldv,nprops,lann eal,
    1 stepTime,totalTime,dt,cmname,coordMp,charLength,pr ops,
    2 density,strainInc,relSpinInc,tempOld,stretchOld,de fgradOld,
    3 fieldOld,stressOld,stateOld,enerInternOld,enerInel asOld,tempNew,
    4 stretchNew,defgradNew,fieldNew,stressNew,stateNew, enerInternNew,
    5 enerInelasNew)
    C
    include 'vaba_param.inc'
    C
    dimension props(nprops),density(nblock),coordMp(nblock,*),
    1 charLength(nblock),strainInc(nblock,ndir+nshr),
    2 relSpinInc(nblock,nshr),tempOld(nblock),
    3 stretchOld(nblock,ndir+nshr),fieldOld(nblock,nfiel dv),
    4 defgradOld(nblock,ndir+nshr+nshr),stressOld(nblock ,ndir+nshr),
    5 stateOld(nblock,nstatev),enerInternOld(nblock),
    6 enerInelasOld(nblock),tempNew(nblock),enerInelasNe w(nblock),
    7 stretchNew(nblock,ndir+nshr),defgradNew(nblock,ndi r+nshr+nshr),
    8 fieldNew(nblock,nfieldv),stressNew(nblock,ndir+nsh r),
    9 stateNew(nblock,nstatev),enerInternNew(nblock)
    c
    implicit none
    integer i,j,j1,k,k1,statusMP
    dimension StrainCurrent(nblock,ndir+nshr)
    real A,B,D,F,G,H,lambda,checkf
    real kk1,gk1,kk2,gk2,gv,E,nu,f1t
    real omega1,omega2,psi1,psi2,mu1,mu2
    real lambda1,lambda2,phi1,phi2,zeta1,zeta2,BM,mu
    real STRtrace,SIGtrace,DSTRtrace,DSIGtrace,S1
    C
    C
    character*80 cmname
    c
    C parameter(zero=0.,one=1.,two=2.,three=3.,four=4.,s ix=6.,nine=9.,half=.5,
    C 1 third=one/three,sixth=one/six,ninth=one/nine,fourth=one/four)
    c
    c
    c Model coefficients
    c
    kk1=props(1)
    gk1=props(2)
    kk2=props(3)
    gk2=props(4)
    gv=props(5)
    E=props(6)
    nu=props(7)
    c kv (infinite-volume conserved flow in dashpot)
    f1t=props(8)
    c
    omega1=0.0
    omega2=1.0+kk2/kk1
    psi1=1.0/(4.0*gv)
    psi2=1.0/(4.0*gk1)
    mu1=gk2/(2.0*gv)
    mu2=0.5*(1.0+gk2/gk1)
    lambda1=(-gk2/gv)/3.0
    lambda2=(1.0/3.0)*(kk2/kk1-gk2/gk1)
    phi1=-1.0/(6.0*gv)
    phi2=(1.0/(9.0*kk1))-(1.0/(6.0*gk1))
    zeta1=0.0
    zeta2=1.0/(3.0*kk1)
    BM=E/(3.0*(1.0-2.0*nu))
    mu=E/(2.0*(1.0+nu))
    c
    A=dt*lambda1+lambda2
    B=2.0*mu1*dt+2.0*mu2
    D=phi1*dt+phi2
    F=2.0*psi1*dt+2.0*psi2
    G=zeta1*dt+zeta2
    H=omega1*dt+omega2
    lambda=BM-(2.0*mu/3.0)
    C
    c
    if(totalTime.EQ.0) then
    do k = 1,nblock
    trace=strainInc(k,1)+strainInc(k,2)+strainInc(k,3)
    C update normal stresses
    stressNew(k,1)=stressOld(k,1)+(lambda*trace)
    1 + (2.0*mu*strainInc(k,1))
    stressNew(k,2)=stressOld(k,2)+(lambda*trace)
    1 + (2.0*mu*strainInc(k,2))
    stressNew(k,3)=stressOld(k,3)+(lambda*trace)
    1 + (2.0*mu*strainInc(k,3))
    C update shear stresses
    stressNew(k,4)=stressOld(k,4)
    1 + (2.0*mu*strainInc(k,4))
    stressNew(k,5)=stressOld(k,5)
    1 + (2.0*mu*strainInc(k,5))
    stressNew(k,6)=stressOld(k,6)
    1 + (2.0*mu*strainInc(k,6))
    end do
    else

    do 100 i = 1,nblock
    c
    c current normal strain tensor components from deformation gradient
    c
    StrainCurrent(i,1)=0.5*(defgradNew(i,1)*defgradNew (i,1))
    1 + (0.5*(defgradNew(i,7)*defgradNew(i,7)))
    1 + (0.5*((defgradNew(i,6)*defgradNew(i,6))-1))
    StrainCurrent(i,2)=0.5*(defgradNew(i,4)*defgradNew (i,4))
    1 + (0.5*(defgradNew(i,2)*defgradNew(i,2)))
    1 + (0.5*((defgradNew(i,8)*defgradNew(i,8))-1))
    StrainCurrent(i,3)=0.5*(defgradNew(i,9)*defgradNew (i,9))
    1 + (0.5*(defgradNew(i,5)*defgradNew(i,5)))
    1 + (0.5*((defgradNew(i,3)*defgradNew(i,3))-1))
    c
    c current shear strain tensor components from deformation gradient
    c
    StrainCurrent(i,4)=0.5*((defgradNew(i,1)*defgradNe w(i,4))
    1 + (defgradNew(i,7)*defgradNew(i,2))
    1 + (defgradNew(i,6)*defgradNew(i,8)))
    StrainCurrent(i,5)=0.5*((defgradNew(i,4)*defgradNe w(i,9))
    1 + (defgradNew(i,2)*defgradNew(i,5))
    1 + (defgradNew(i,8)*defgradNew(i,3)))
    StrainCurrent(i,6)=0.5*((defgradNew(i,1)*defgradNe w(i,9))
    1 + (defgradNew(i,7)*defgradNew(i,5))
    1 + (defgradNew(i,6)*defgradNew(i,3)))
    C
    statusMp=stateNew(i,1)
    c
    STRtrace=StrainCurrent(i,1)+StrainCurrent(i,2)
    1 + StrainCurrent(i,3)
    c
    c Trace of current stress tensor
    SIGtrace=stressOld(i,1)+stressOld(i,2)+stressOld(i ,3)
    c
    c Trace of current strain increment tensor
    DSTRtrace=strainInc(i,1)+strainInc(i,2)+strainInc( i,3)
    c
    c Evaluate the trace of the stress increment
    DSIGtrace=(H/G)*STRtrace+(dt/G)*((-zeta1*SIGtrace)
    1 + (omega1*STRtrace))
    c
    C Stress terms
    c
    c Update normal stresses
    c
    if(statusMP.EQ.1.0) then
    do j=1,3
    stressNew(i,j)= stressOld(i,j)+(A/F)*DSTRtrace
    1 + (B/F)*strainInc(i,j)-(D/F)*DSIGtrace
    2 + (dt/F)*((lamda1*STRtrace)+(2.0*mu1*StrainCurrent(i,j))
    3 - (phi1*SIGtrace)-(2.0*psi1*stressOld(i,j)))
    c Update shear stress
    stressNew(i,j+3)=stressOld(i,j+3)+(B/F)*strainInc(i,j+3)
    1 + ((2.0*dt)/F)*((mu1*StrainCurrent(i,j+3))
    2 - (psi1*stressOld(i,j+3)))
    end do
    else
    do j1=1,6
    stressNew(i,j1)=0.0
    end do
    endif

    c
    c Tensile failure criteria
    c
    S1=stressNew(i,3)
    checkf=S1/f1t
    if(checkf.GE.1) then
    stateNew(i,1)=0.0
    else
    stateNew(i,1)=1.0
    endif
    c
    100 continue
    endif
    c
    return
    end


    INPUT FILE

    *Heading
    ** Job name: oct15 Model name: Model-1
    *Preprint, echo=NO, model=NO, history=NO, contact=NO
    **
    ** PARTS
    **
    *Part, name=Part-1
    *Node
    1, 20., 20., 40.
    2, 20., -20., 40.
    3, 20., 20., 0.
    4, 20., -20., 0.
    5, -20., 20., 40.
    6, -20., -20., 40.
    7, -20., 20., 0.
    8, -20., -20., 0.
    *Element, type=C3D8R
    1, 5, 6, 8, 7, 1, 2, 4, 3
    *Nset, nset=_PickedSet5, internal, generate
    1, 8, 1
    *Elset, elset=_PickedSet5, internal
    1,
    ** Section: Section-1
    *Solid Section, elset=_PickedSet5, controls=EC-1, material=Material-1
    1.,
    *End Part
    **
    **
    ** ASSEMBLY
    **
    *Assembly, name=Assembly
    **
    *Instance, name=Part-1-1, part=Part-1
    *End Instance
    **
    *Nset, nset=_PickedSet13, internal, instance=Part-1-1
    3, 4, 7, 8
    *Elset, elset=__PickedSurf12_S3, internal, instance=Part-1-1
    1,
    *Surface, type=ELEMENT, name=_PickedSurf12, internal
    __PickedSurf12_S3, S3
    *End Assembly
    **
    ** ELEMENT CONTROLS
    **
    *Section Controls, name=EC-1, ELEMENT DELETION=YES
    1., 1., 1.
    **
    ** MATERIALS
    **
    *Material, name=Material-1
    *Density
    1.44e-09,
    *Depvar, delete=1
    1,
    *User Material, constants=8
    933., 1400.,22700.,34000., 1.,90000., 0.3, 2610.
    ** ----------------------------------------------------------------
    **
    ** STEP: applyload
    **
    *Step, name=applyload
    *Dynamic, Explicit
    , 0.0005
    *Bulk Viscosity
    0.06, 1.2
    **
    ** BOUNDARY CONDITIONS
    **
    ** Name: BC-2 Type: Displacement/Rotation
    *Boundary
    _PickedSet13, 1, 1
    _PickedSet13, 2, 2
    _PickedSet13, 3, 3
    _PickedSet13, 4, 4
    _PickedSet13, 5, 5
    _PickedSet13, 6, 6
    **
    ** LOADS
    **
    ** Name: Load-1 Type: Surface traction
    *Dsload
    _PickedSurf12, TRVEC, 1300., 0., 0., 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
    LE, PE, PEEQ, S, STATUS
    *Contact Output
    CSTRESS,
    **
    ** HISTORY OUTPUT: H-Output-1
    **
    *Output, history, variable=PRESELECT
    *End Step
    Attached Files

  • #2
    Since you are using an explicit simulation the resulting stress & strain values are often noisy. Are you sure you are running the Explicit anlysis correctly?

    -Jorgen
    Jorgen Bergstrom, Ph.D. PolymerFEM Administrator

    Comment

    Working...
    X