Announcement

Collapse
No announcement yet.

VUMAT EXAMPLES for Hyperelastic materials

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

  • VUMAT EXAMPLES for Hyperelastic materials

    Hello,

    I am new with VUMAT. Recently I have finished reading VUMAT example for Neo-Hookean material from the "Mechanics of Solid Polymers" textbook. Could anyone share more VUMAT examples of Hyperelastic materials for learning purposes? Thank you very much.

    Regards,

    Willsen

  • #2
    Hello,

    I have written a hyperelastic VUMAT that can handle strain-energy functions of the form:

    W = K1(I1-3)^m + K2(I1-3)^p + K3(I1-3)^q

    which covers the Neo-Hookean, Yeoh, and (almost) Davies-De-Thomas (1994) models. I think it will be published in an appendix of a related work soon. If you want to take a look before then, let me know and I may be allowed to share it with you before it hits the press.

    Regards,
    Travis

    Comment


    • #3
      I forgot to add the compressibility term which is required for VUMAT. The form is:

      W = K1(I1-3)^m + K2(I1-3)^p + K3(I1-3)^q + 1/D1(J-1)^2

      Comment


      • #4
        Hi,

        Thank you very much for your kind offer. I am quite confident with VUMAT now so all good!
        Congratulation with your publication!

        Regards,

        Willsen

        Comment


        • #5
          Hi,

          I am new to VUMAT too and appreciate if anyone share me a VUMAT subroutine and its CAE file so I can go through it carefully.

          Thanks so much,
          Somaye

          Comment


          • #6
            Hi Somaye,

            Recalling my experience on how to learn VUMAT, I did a lot of in depth study on the basic of Continuum Mechanics. Through a proper understanding of deformation gradient, different strain and stress measures in large deformation, it becomes really straight forward to understand the different physical quantities available in VUMAT. Jorgen's VUMAT for Neo Hookean material was super valuable to me.
            My suggestion will be to understand his VUMAT in depth and apply it in a single element test. Do an identical test using ABAQUS' in-built Neo-Hookean material model and compare the results.
            Build upon that. Try to write a more complicated VUMAT to describe Mooney-Rivlin material and compare it with the equivalent mateiral model from ABAQUS, and so on.


            Attached is an example of VUMAT for Mooney Rivlin model:

            ! This is a VUMAT subroutine for Mooney-Rivlin


            subroutine vumat(nblock, ndi, nshr, nstatev, nfieldv, nprops, &
            lanneal, stepTime, totTime, dt, cmname, coordMp, &
            charLen, props, density, strainInc, relSpinInc, temp0, &
            U0, F0, field0, stress0, state0, &
            intEne0, inelaEn0, temp1, U1, &
            F1, field1, stress1, state1, intEne1, inelaEn1)
            implicit none

            integer :: nblock, ndi, nshr, nstatev, nfieldv, nprops, lanneal
            real :: stepTime, totTime, dt
            character*80 :: cmname
            real :: coordMp(nblock,*)
            real :: charLen, props(nprops), density(nblock), &
            strainInc(nblock,ndi+nshr), relSpinInc(nblock,nshr), &
            temp0(nblock), U0(nblock,ndi+nshr), &
            F0(nblock,ndi+nshr+nshr), field0(nblock,nfieldv), &
            stress0(nblock,ndi+nshr), state0(nblock,nstatev), &
            intEne0(nblock), inelaEn0(nblock), temp1(nblock), &
            U1(nblock,ndi+nshr), F1(nblock,ndi+nshr+nshr), &
            field1(nblock,nfieldv), stress1(nblock,ndi+nshr), &
            state1(nblock,nstatev), intEne1(nblock), inelaEn1(nblock)

            ! local variables
            real :: J, U(3,3), U2dev(3,3), U4dev(3,3), C10, C01, D1 ,traceU2dev, traceU4dev , constant1, constant2
            integer :: i
            C10 = props(1)
            C01 = props(2)
            D1 = props(3)

            ! loop through all blocks
            do i = 1, nblock
            ! setup F (upper diagonal part)
            U(1,1) = U1(i,1)
            U(2,2) = U1(i,2)
            U(3,3) = U1(i,3)
            U(1,2) = U1(i,4)
            U(2,1) = U(1,2)
            if (nshr .eq. 1) then
            U(2,3) = 0.0
            U(1,3) = 0.0
            else
            U(2,3) = U1(i,5)
            U(1,3) = U1(i,6)
            end if
            U(3,2)=U(2,3)
            U(3,1)=U(1,3)
            ! calculate J
            J = (U(1,1)*(U(2,2)*U(3,3)-U(2,3)**2)) + (U(1,2)*(U(2,3)*U(1,3)-U(1,2)*U(3,3))) + (U(1,3)*(U(1,2)*U(2,3)-U(2,2)*U(1,3)))


            ! Define the square of the deviatoric stretch tensor
            U2dev=(J**(-2.0/3.0))*matmul(U,U)

            ! Define the forth power of the deviatoric stretch tensor
            U4dev=matmul(U2dev,U2dev)

            ! Define some constants to simplify the stress expression
            traceU2dev=(U2dev(1,1)+U2dev(2,2)+U2dev(3,3))
            traceU4dev=(U4dev(1,1)+U4dev(2,2)+U4dev(3,3))
            constant1=2.0*C10/J
            constant2=2.0*C01/J


            !The corotational stress
            stress1(i,1) = (constant1+constant2*traceU2dev)*(U2dev(1,1)-1.0/3.0*traceU2dev) - constant2*(U4dev(1,1)-1.0/3.0*traceU4dev) + 2.0/D1*(J-1.0)
            stress1(i,2) = (constant1+constant2*traceU2dev)*(U2dev(2,2)-1.0/3.0*traceU2dev) - constant2*(U4dev(2,2)-1.0/3.0*traceU4dev) + 2.0/D1*(J-1.0)
            stress1(i,3) = (constant1+constant2*traceU2dev)*(U2dev(3,3)-1.0/3.0*traceU2dev) - constant2*(U4dev(3,3)-1.0/3.0*traceU4dev) + 2.0/D1*(J-1.0)
            stress1(i,4) = (constant1+constant2*traceU2dev)*U2dev(1,2) - constant2*U4dev(1,2)
            if (nshr .eq. 3) then
            stress1(i,5) = (constant1+constant2*traceU2dev)*U2dev(2,3) - constant2*U4dev(2,3)
            stress1(i,6) = (constant1+constant2*traceU2dev)*U2dev(3,1) - constant2*U4dev(3,1)
            end if
            end do
            end subroutine vumat

            Comment


            • #7
              Thank you so much Willsen!
              Can the Jorgen's VUMAT for Neo Hookean material or yours be run as .f or .for file in ABAQUS? It seems those formats are kind of different from the umat file I see here :
              https://abaqus-docs.mit.edu/2017/Eng...ces/umathrt2.f

              Comment


              • #8
                Hi There,
                I am having a problem with inputting my Vumat into Abaqus explicit.
                I have already implemented many subroutines into Abaqus and got results. However, my recent vumat gives this error :"problem during compilation"
                my vumat is just simple neo-hookean model.
                I appreciate if someone helps me with that. the code is written here:

                subroutine vumat(
                C Read only (unmodifiable)variables -
                1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
                2 stepTime, totalTime, dt, cmname, coordMp, charLength,
                3 props, density, strainInc, relSpinInc,
                4 tempOld, stretchOld, defgradOld, fieldOld,
                5 stressOld, stateOld, enerInternOld, enerInelasOld,
                6 tempNew, stretchNew, defgradNew, fieldNew,
                C Write only (modifiable) variables -
                7 stressNew, stateNew, enerInternNew, 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),
                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),
                8 defgradNew(nblock,ndir+nshr+nshr),
                9 fieldNew(nblock,nfieldv),
                1 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
                2 enerInternNew(nblock), enerInelasNew(nblock),
                C
                character*80 cmname
                C
                parameter( zero = 0.D0, one = 1.D0, two = 2.D0, three = 3.D0,
                1 third = one/three, half = 0.5D0, twoThirds = two/three,
                2 threeHalfs = 1.5D0 )
                C
                C elastic constants
                C
                c10 = props(1)
                d1 = props(2)
                shrMod = two*c10
                blkMod = two/d1
                C
                C
                do k = 1,nblock
                C
                C calculate left cauchy-green tensor in terms of stretch tensor
                C and jacobian
                bxx =stretchNew(k,1) * stretchNew(k,1)
                1 + stretchNew(k,4) * stretchNew(k,4)
                byy =stretchNew(k,4) * stretchNew(k,4)
                1 + stretchNew(k,2) * stretchNew(k,2)
                bzz =stretchNew(k,3) * stretchNew(k,3)
                bxy =stretchNew(k,1) * stretchNew(k,4)
                1 + stretchNew(k,4) * stretchNew(k,2)
                bxz = zero
                byz = zero
                detu=stretchNew(k,3)*(stretchNew(k,1)*stretchNew(k ,2)
                1 - stretchNew(k,4)*stretchNew(k,4))
                C
                if ( nshr .gt. 1 ) then
                bxx =bxx + stretchNew(k,6) * stretchNew(k,6)
                byy =byy + stretchNew(k,5) * stretchNew(k,5)
                bzz =bzz + stretchNew(k,6) * stretchNew(k,6)
                1 + stretchNew(k,5) * stretchNew(k,5)
                bxy =bxy + stretchNew(k,6) * stretchNew(k,6)
                1 + stretchNew(k,4) * stretchNew(k,2)
                bxz = stretchNew(k,1) * stretchNew(k,6)
                1 + stretchNew(k,4) * stretchNew(k,5)
                + stretchNew(k,6) * stretchNew(k,3)
                byz = stretchNew(k,4) * stretchNew(k,6)
                1 + stretchNew(k,2) * stretchNew(k,5)
                2 + stretchNew(k,5) * stretchNew(k,3)
                detu = detu + stretchNew(k,6) *
                1 ( stretchNew(k,4) * stretchNew(k,5)
                2 - stretchNew(k 6) * stretchNew(k,2) )
                3 - stretchNew(k,5) * ( stretchNew(k,1) * stretchNew(k,5)
                4 - stretchNew(k,6) * stretchNew(k,4) )
                end if
                xpow = exp ( - log(detu) * two_thirds )
                bxx = bxx * xpow
                byy = byy * xpow
                bzz = bzz * xpow
                bxy = bxy * xpow
                bxz = bxz * xpow
                byz = byz * xpow
                C BI1 ( first invariant of BIJ )
                C
                bi1 = bxx + byy + bzz
                C
                C BDIJ ( deviatoric part of BIJ )
                C
                bdxx = bxx - third * bi1
                bdyy = byy - third * bi1
                bdyy = byy - third * bi1
                bdzz = bzz - third * bi1
                bdxy = bxy
                bdxz = bxz
                bdyz = byz
                C
                duDi1 = half * shrMod
                duDi3 = blkMod * ( detu - one )
                C
                C Calculate Cauchy (true) stresses
                C
                detuInv = one / detu
                factor = two * duDi1 * detuInv
                C
                stressNew(k,1) = factor * bdxx + duDi3
                stressNew(k,2) = factor * bdyy + duDi3
                stressNew(k 3) = factor * bdzz + duDi3
                stressNew(k,3) = factor * bdzz + duDi3
                stressNew(k,4) = factor * bdxy
                if ( nshr .gt. 1 ) then
                stressNew(k,5) = factor * bdyz
                stressNew(k,6) = factor * bdxz
                end if
                C
                end do
                C
                return
                end




                Thanks!

                Comment

                Working...
                X