Notifications
Clear all

VUMAT for Hyperelastic Materials - Update to Bergstrom's Example

10 Posts
4 Users
0 Likes
3,277 Views
Posts: 21
Topic starter
(@twhohen)
Eminent Member
Joined: 5 years ago

Hello,

Several posts discuss VUMATs for hyperelastic materials. Bergstrom (2015) provides an example for the Neo-Hookean strain-energy function (also copied to several forums) that was written for older versions of Abaqus. An updated, complete subroutine is published here:

Hohenberger et al (2019). A constitutive model for both low and high strain nonlinearities in highly filled elastomers and implementation with user-defined material subroutines in Abaqus. Rubber Chemistry and Technology, 92 (4), 653-686.

It uses a modified (generalized) Yeoh strain-energy function and is confirmed to work in Abaqus v6.14 and R2018. I am sharing this here since many posts inquiring about hyperelastic VUMATs are now closed.

Regards,
Travis

9 Replies
Posts: 3982
(@jorgen)
Member
Joined: 4 years ago

Hi Travis, that I great. Are you able to also share the source code here?

-Jorgen

Reply
2 Replies
(@twhohen)
Joined: 5 years ago

Eminent Member
Posts: 21

@jorgen

I will check with the publisher. In the meantime, anybody with immediate interest in the code can contact the corresponding author, James Busfield. He is easy to find with a web search.

-th

Reply
(@sayed241)
Joined: 2 years ago

New Member
Posts: 2

@jorgen can you please share a vumat for composit materials 

 

Reply
Posts: 21
Topic starter
(@twhohen)
Eminent Member
Joined: 5 years ago

Hello,

I have confirmed that I can post the code. For a detailed explanation of the statements, refer to the article above which summarizes the finite elasticity theory that is necessary to build a hyperelastic VUMAT. A discussion of implementation of the theory in Abaqus is given in Appendix C. The article was built largely from the more detailed accounts in the books by Holzapfel (2000) and Bergström (2015).

Some pertinent updates to the Neo-Hookean subroutine that is commonly shared online are:

- updated names for Abaqus variables.
- a gen-Yeoh SEF is used. Set:
----- (K2,K3) = (0,0) for Neo-Hookean SEF.
----- (m,p,q) = (1,2,3) for Yeoh SEF.
----- (K2,K3) = (0,2) for Davies-De-Thomas SEF (1994, Rubber Chem. & Tech. 67 (4), 716).
- includes linear elastic initialization step required per Abaqus VUMAT subroutine template.
- includes logic to handle indeterminate stress condition that can arise with the gen-Yeoh SEF.

The code only handles 3D elements. A 2D plane-strain or axisymmetric (but not plane-stress) version can be written by removing references to the 5th and 6th elements in the stress, strain, and strain increment vectors. Alternatively, the logic in Bergström's Neo-Hookean example can be included to handle 2D and 3D cases.

Regards,
Travis

C *******************************************************************
C VUMAT (3D only) for gen-Yeoh Strain-Energy Function
C Compatibility confirmed for Abaqus v6.14 & R2018
C ---------------------------------------------------
C Authors: Travis Hohenberger & Richard Windslow
C Date: 2019-07-11
C E-mail: twhohen@gmail.com
C Source: Hohenberger et al (2019). Rubber Chem. & Tech., 92 (4), 653
C
*******************************************************************
C
C Strain-energy function:
C
C W = K1*(I1-3)^m + K2*(I1-3)^p + K3*(I1-3)^q + (1/D1)*(J-1)^2
C
C where K1, K2, K3, m, p, q are distortional fitting parameters and D1
C is a volumetric fitting parameter. I1 is the first invariant of the
C modified stretch tensor. J is the volumetric ratio.
C
C **********************************************************************
C
SUBROUTINE VUMAT(
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,
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),
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),
1 fieldNew(nblock,nfieldv),
2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
3 enerInternNew(nblock), enerInelasNew(nblock)
C
CHARACTER*80 cmname
C
C PARAMETERS
C ----------
REAL*8 oneThrd, half, twoThrd, one, two, three, thresh
PARAMETER(oneThrd=1.d0/3.d0, half=0.5d0, twoThrd=2.d0/3.d0,
$ one=1.d0, two=2.d0, three=3.d0, thresh=10.d0**-12.d0)
C
C LOCAL VARIABLES
C ---------------
REAL*8 g0 , k0 , twoG , lmda , trace , d1 ,
$ k1 , k2 , k3 , em , pe , qu ,
$ Bxx , Byy , Bzz , Bxy , Bxz , Byz ,
$ BbarXX , BbarYY , BbarZZ , BbarXY , BbarXZ , BbarYZ ,
$ dBbarXX , dBbarYY, dBbarZZ, dBbarXY, dBbarXZ, dBbarYZ ,
$ J , J23 , duDi1 , duDi3 , I1 , g1 ,
$ p0 , j1 , u1
C
k1 = props(1)
k2 = props(2)
k3 = props(3)
em = props(4)
pe = props(5)
qu = props(6)
d1 = props(7)
C
g0 = two * k1
k0 = two / d1
C
C ***************************************************************
C ------------ INITIALIZE MATERIAL AS LINEARLY ELASTIC -------------
C
***************************************************************
C
twoG = two * g0
lmda = k0 - twoG * oneThrd
C
IF (totalTime.EQ.0.0) THEN
C
DO k = 1,nblock
trace = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)
stressNew(k,1) = stressOld(k,1) + twoG*strainInc(k,1) +
$ lmda*trace
stressNew(k,2) = stressOld(k,2) + twoG*strainInc(k,2) +
$ lmda*trace
stressNew(k,3) = stressOld(k,3) + twoG*strainInc(k,3) +
$ lmda*trace
stressNew(k,4) = stressOld(k,4) + twoG*strainInc(k,4)
stressNew(k,5) = stressOld(k,5) + twoG*strainInc(k,5)
stressNew(k,6) = stressOld(k,6) + twoG*strainInc(k,6)
END DO
C
RETURN
C
END IF
C
C ***************************************************************
C ----------- START LOOP FOR MATERIAL POINT CALCULATIONS -----------
C
***************************************************************
C
DO k = 1,nblock
C
C CALCULATE LEFT CAUCHY-GREEN STRAIN TENSOR, B = U*U
C --------------------------------------------------
Bxx = stretchNew(k,1) * stretchNew(k,1) +
$ stretchNew(k,4) * stretchNew(k,4) +
$ stretchNew(k,6) * stretchNew(k,6)
Byy = stretchNew(k,2) * stretchNew(k,2) +
$ stretchNew(k,4) * stretchNew(k,4) +
$ stretchNew(k,5) * stretchNew(k,5)
Bzz = stretchNew(k,3) * stretchNew(k,3) +
$ stretchNew(k,5) * stretchNew(k,5) +
$ stretchNew(k,6) * stretchNew(k,6)
Bxy = stretchNew(k,1) * stretchNew(k,4) +
$ stretchNew(k,4) * stretchNew(k,2) +
$ stretchNew(k,6) * stretchNew(k,5)
Bxz = stretchNew(k,1) * stretchNew(k,6) +
$ stretchNew(k,4) * stretchNew(k,5) +
$ stretchNew(k,6) * stretchNew(k,3)
Byz = stretchNew(k,4) * stretchNew(k,6) +
$ stretchNew(k,2) * stretchNew(k,5) +
$ stretchNew(k,5) * stretchNew(k,3)
C
C CALCULATE J = |F| = |U|
C -----------------------
C
J = stretchNew(k,1) *
$ ( stretchNew(k,2) * stretchNew(k,3) -
$ stretchNew(k,5) * stretchNew(k,5) ) +
$ stretchNew(k,4) *
$ ( stretchNew(k,5) * stretchNew(k,6) -
$ stretchNew(k,3) * stretchNew(k,4) ) +
$ stretchNew(k,6) *
$ ( stretchNew(k,4) * stretchNew(k,5) -
$ stretchNew(k,2) * stretchNew(k,6) )
C
C CALCULATE MODIFIED STRAIN TENSOR, Bbar = J^(-2/3)*B
C ---------------------------------------------------
J23 = J**(-twoThrd)
C
BbarXX = J23 * Bxx
BbarYY = J23 * Byy
BbarZZ = J23 * Bzz
BbarXY = J23 * Bxy
BbarXZ = J23 * Bxz
BbarYZ = J23 * Byz
C
C FIRST INVARIANT OF Bbar = tr(Bbar)
C ----------------------------------
I1 = BbarXX + BbarYY + BbarZZ
C
C DEVIATORIC PART OF Bbar
C -----------------------
p0 = oneThrd * I1
C
dBbarXX = BbarXX - p0
dBbarYY = BbarYY - p0
dBbarZZ = BbarZZ - p0
dBbarXY = BbarXY
dBbarXZ = BbarXZ
dBbarYZ = BbarYZ
C
C DERIVATIVES OF STRAIN-ENERGY FUNCTION
C -------------------------------------
j1 = I1 - three
C
IF (j1.LT.thresh) THEN
duDi1 = zero
ELSE
duDi1 = em * k1 * j1**(em-one) +
$ pe * k2 * j1**(pe-one) +
$ qu * k3 * j1**(qu-one)
END IF
C
duDi3 = two/d1 * (J - one)
C
C COROTATIONAL CAUCHY (TRUE) STRESSES
C -----------------------------------
g1 = two/J * duDi1
C
stressNew(k,1) = g1 * dBbarXX + duDi3
stressNew(k,2) = g1 * dBbarYY + duDi3
stressNew(k,3) = g1 * dBbarZZ + duDi3
stressNew(k,4) = g1 * dBbarXY
stressNew(k,5) = g1 * dBbarYZ
stressNew(k,6) = g1 * dBbarXZ
C
C UPDATE SPECIFIC INTERNAL ENERGY
C -------------------------------
u1 = 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) +
$ two * ( (stressOld(k,4) + stressNew(k,4))*
$ strainInc(k,4) +
$ (stressOld(k,5) + stressNew(k,5))*
$ strainInc(k,5) +
$ (stressOld(k,6) + stressNew(k,6))*
$ strainInc(k,6) ) )
C
enerInternNew(k) = enerInternOld(k) + u1 / density(k)
C
END DO
C
RETURN
C
END SUBROUTINE VUMAT
Reply
2 Replies
(@fat32ntfs)
Joined: 6 years ago

New Member
Posts: 2

@twhohen hi, i think you may have forgot to define zero after the line:

 IF (j1.LT.thresh) THEN
Reply
(@twhohen)
Joined: 5 years ago

Eminent Member
Posts: 21

@fat32ntfs,

Yes, it looks like zero is not defined. I never encountered problems validating this version of the code, so my compiler must have been assigning it a default value of 0. Best to define that zero parameter at the top of the code.

Thanks for catching that.

Regards,
Travis

Reply
Posts: 2
(@fat32ntfs)
New Member
Joined: 6 years ago

So it seems I cannot edit my post, so I apologize for the double post. I wanted to add the following general question:

In my VUMATs, I use matrix multiplications rather than computing the components individually, for example:

C	F^eE = U
U(1,1) = stretchNew(k,1)
U(2,2) = stretchNew(k,2)
U(3,3) = stretchNew(k,3)
U(1,2) = stretchNew(k,4)
U(2,3) = stretchNew(k,5)
U(1,3) = stretchNew(k,6)
U(2,1) = U(1,2)
U(3,1) = U(1,3)
U(3,2) = U(2,3)
C
C B^ee = F^eE * F^eE^{T} = U^2
B = zero
B = matmul(U,U)

Is there any reason to choose one over the other?

Reply
1 Reply
(@twhohen)
Joined: 5 years ago

Eminent Member
Posts: 21

@fat32ntfs,

I have not personally tried matrix multiplication in my codes, but I have some that do (for example, search "Shawn Chester New Jersey Institute of Technology" to find his codes that use matrix multiplication). I am not sure which technique is faster, so if you test computation speeds, it would be nice to know the difference.

Regards,

Travis

Reply
Posts: 2
(@sayed241)
New Member
Joined: 2 years ago

subroutine vumat(
c Read only -
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 -
7 stressNew, stateNew, enerInternNew, enerInelasNew )
c
include 'vaba_param.inc'
c
c 3D Orthotropic Elasticity with Hashin 3d Failure criterion
c
c The state variables are stored as:
c state(*,1) = material point status
c state(*,2:7) = damping stresses
c
c User defined material properties are stored as
c * First line:
c props(1) --> Young's modulus in 1-direction, E1
c props(2) --> Young's modulus in 2-direction, E2
c props(3) --> Young's modulus in 3-direction, E3
c props(4) --> Poisson's ratio, nu12
c props(5) --> Poisson's ratio, nu13
c props(6) --> Poisson's ratio, nu23
c props(7) --> Shear modulus, G12
c props(8) --> Shear modulus, G13
c
c * Second line:
c props(9) --> Shear modulus, G23
c props(10) --> beta damping parameter
c props(11) --> "not used"
c props(12) --> "not used"
c props(13) --> "not used"
c props(14) --> "not used"
c props(15) --> "not used"
c props(16) --> "not used"
c
c * Third line:
c props(17) --> Ultimate tens stress in 1-direction, sigu1t
c props(18) --> Ultimate comp stress in 1-direction, sigu1c
c props(19) --> Ultimate tens stress in 2-direction, sigu2t
c props(20) --> Ultimate comp stress in 2-direction, sigu2c
c props(21) --> Ultimate tens stress in 2-direction, sigu3t
c props(22) --> Ultimate comp stress in 2-direction, sigu3c
c props(23) --> "not used"
c props(24) --> "not used"
c
c * Fourth line:
c props(25) --> Ultimate shear stress, sigu12
c props(26) --> Ultimate shear stress, sigu13
c props(27) --> Ultimate shear stress, sigu23
c props(28) --> "not used"
c props(29) --> "not used"
c props(30) --> "not used"
c props(31) --> "not used"
c props(32) --> "not used"
c

dimension props(nprops), density(nblock),
1 coordMp(nblock,*),
2 charLength(*), strainInc(nblock,ndir+nshr),
3 relSpinInc(nblock,nshr), tempOld(nblock),
4 stretchOld(nblock,ndir+nshr), defgradOld(nblock,ndir+nshr+nshr),
5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(*),
8 stretchNew(nblock,ndir+nshr), defgradNew(nblock,ndir+nshr+nshr),
9 fieldNew(nblock,nfieldv), stressNew(nblock,ndir+nshr),
1 stateNew(nblock,nstatev),
2 enerInternNew(nblock), enerInelasNew(nblock)
*
character*80 cmname
*
parameter( zero = 0.d0, one = 1.d0, two = 2.d0, half = .5d0 )
*
parameter(
* i_svd_DmgFiberT = 1,
* i_svd_DmgFiberC = 2,
* i_svd_DmgMatrixT = 3,
* i_svd_DmgMatrixC = 4,
* i_svd_statusMp = 5,
* i_svd_dampStress = 6,
c * i_svd_dampStressXx = 6,
c * i_svd_dampStressYy = 7,
c * i_svd_dampStressZz = 8,
c * i_svd_dampStressXy = 9,
c * i_svd_dampStressYz = 10,
c * i_svd_dampStressZx = 11,
* i_svd_Strain = 12,
c * i_svd_StrainXx = 12,
c * i_svd_StrainYy = 13,
c * i_svd_StrainZz = 14,
c * i_svd_StrainXy = 15,
c * i_svd_StrainYz = 16,
c * i_svd_StrainZx = 17,
* n_svd_required = 17 )
*
parameter(
* i_s33_Xx = 1,
* i_s33_Yy = 2,
* i_s33_Zz = 3,
* i_s33_Xy = 4,
* i_s33_Yz = 5,
* i_s33_Zx = 6 )
*
* Structure of property array
parameter (
* i_pro_E1 = 1,
* i_pro_E2 = 2,
* i_pro_E3 = 3,
* i_pro_nu12 = 4,
* i_pro_nu13 = 5,
* i_pro_nu23 = 6,
* i_pro_G12 = 7,
* i_pro_G13 = 8,
* i_pro_G23 = 9,
*
* i_pro_beta = 10,
*
* i_pro_sigu1t = 17,
* i_pro_sigu1c = 18,
* i_pro_sigu2t = 19,
* i_pro_sigu2c = 20,
* i_pro_sigu3t = 21,
* i_pro_sigu3c = 22,
* i_pro_sigu12 = 25,
* i_pro_sigu13 = 26,
* i_pro_sigu23 = 27 )
* Temporary arrays
dimension eigen(maxblk*3)
*
* Read material properties
*
E1 = props(i_pro_E1)
E2 = props(i_pro_E2)
E3 = props(i_pro_E3)
xnu12 = props(i_pro_nu12)
xnu13 = props(i_pro_nu13)
xnu23 = props(i_pro_nu23)
G12 = props(i_pro_G12)
G13 = props(i_pro_G13)
G23 = props(i_pro_G23)
*
xnu21 = xnu12 * E2 / E1
xnu31 = xnu13 * E3 / E1
xnu32 = xnu23 * E3 / E2
*
*
* Compute terms of stiffness matrix
gg = one / ( one - xnu12*xnu21 - xnu23*xnu32 - xnu31*xnu13
* - two*xnu21*xnu32*xnu13 )
C11 = E1 * ( one - xnu23*xnu32 ) * gg
C22 = E2 * ( one - xnu13*xnu31 ) * gg
C33 = E3 * ( one - xnu12*xnu21 ) * gg
C12 = E1 * ( xnu21 + xnu31*xnu23 ) * gg
C13 = E1 * ( xnu31 + xnu21*xnu32 ) * gg
C23 = E2 * ( xnu32 + xnu12*xnu31 ) * gg
*
f1t = props(i_pro_sigu1t)
f1c = props(i_pro_sigu1c)
f2t = props(i_pro_sigu2t)
f2c = props(i_pro_sigu2c)
f3t = props(i_pro_sigu3t)
f3c = props(i_pro_sigu3c)
f12 = props(i_pro_sigu12)
f13 = props(i_pro_sigu13)
f23 = props(i_pro_sigu23)
*
beta = props(i_pro_beta)
*
* Assume purely elastic material at the beginning of the analysis
*
if ( totalTime .eq. zero ) then
if (nstatev .lt. n_svd_Required) then
call xplb_abqerr(-2,'Subroutine VUMAT requires the '//
* 'specification of %I state variables. Check the '//
* 'definition of *DEPVAR in the input file.',
* n_svd_Required,zero,' ')
call xplb_exit
end if
call OrthoEla3dExp ( nblock,
* stateOld(1,i_svd_DmgFiberT),
* stateOld(1,i_svd_DmgFiberC),
* stateOld(1,i_svd_DmgMatrixT),
* stateOld(1,i_svd_DmgMatrixC),
* C11, C22, C33, C12, C23, C13, G12, G23, G13,
* strainInc,
* stressNew )
return
end if
*
* Update total elastic strain
call strainUpdate ( nblock, strainInc,
* stateOld(1,i_svd_strain), stateNew(1,i_svd_strain) )
*
* Stress update
call OrthoEla3dExp ( nblock,
* stateOld(1,i_svd_DmgFiberT),
* stateOld(1,i_svd_DmgFiberC),
* stateOld(1,i_svd_DmgMatrixT),
* stateOld(1,i_svd_DmgMatrixC),
* C11, C22, C33, C12, C23, C13, G12, G23, G13,
* stateNew(1,i_svd_strain),
* stressNew )
*
* Failure evaluation
*
call copyr ( nblock,
* stateOld(1,i_svd_DmgFiberT), stateNew(1,i_svd_DmgFiberT) )
call copyr ( nblock,
* stateOld(1,i_svd_DmgFiberC), stateNew(1,i_svd_DmgFiberC) )
call copyr ( nblock,
* stateOld(1,i_svd_DmgMatrixT), stateNew(1,i_svd_DmgMatrixT) )
call copyr ( nblock,
* stateOld(1,i_svd_DmgMatrixC), stateNew(1,i_svd_DmgMatrixC) )
nDmg = 0
call eig33Anal ( nblock, stretchNew, eigen )
call Hashin3d ( nblock, nDmg,
* f1t, f2t, f3t, f1c, f2c, f3c, f12, f23, f13,
* stateNew(1,i_svd_DmgFiberT),
* stateNew(1,i_svd_DmgFiberC),
* stateNew(1,i_svd_DmgMatrixT),
* stateNew(1,i_svd_DmgMatrixC),
* stateNew(1,i_svd_statusMp),
* stressNew, eigen )
* -- Recompute stresses if new Damage is occurring
if ( nDmg .gt. 0 ) then
call OrthoEla3dExp ( nblock,
* stateNew(1,i_svd_DmgFiberT),
* stateNew(1,i_svd_DmgFiberC),
* stateNew(1,i_svd_DmgMatrixT),
* stateNew(1,i_svd_DmgMatrixC),
* C11, C22, C33, C12, C23, C13, G12, G23, G13,
* stateNew(1,i_svd_strain),
* stressNew )
end if
*
* Beta damping
if ( beta .gt. zero ) then
call betaDamping3d ( nblock,
* beta, dt, strainInc,
* stressOld, stressNew,
* stateNew(1,i_svd_statusMp),
* stateOld(1,i_svd_dampStress),
* stateNew(1,i_svd_dampStress) )
end if
*
* Integrate the internal specific energy (per unit mass)
*
call EnergyInternal3d ( nblock, stressOld, stressNew,
* strainInc, density, enerInternOld, enerInternNew )
*
return
end

************************************************************
* OrthoEla3dExp: Orthotropic elasticity - 3d *
************************************************************
subroutine OrthoEla3dExp ( nblock,
* dmgFiberT, dmgFiberC, dmgMatrixT, dmgMatrixC,
* C11, C22, C33, C12, C23, C13, G12, G23, G13,
* strain, stress )
*
include 'vaba_param.inc'

* Orthotropic elasticity, 3D case -
*
parameter( zero = 0.d0, one = 1.d0, two = 2.d0)
parameter(
* i_s33_Xx = 1,
* i_s33_Yy = 2,
* i_s33_Zz = 3,
* i_s33_Xy = 4,
* i_s33_Yz = 5,
* i_s33_Zx = 6,
* n_s33_Car = 6 )
*
dimension strain(nblock,n_s33_Car),
* dmgFiberT(nblock), dmgFiberC(nblock),
* dmgMatrixT(nblock), dmgMatrixC(nblock),
* stress(nblock,n_s33_Car)
* -- shear fraction in matrix tension and compression mode
parameter ( smt = 0.9d0, smc = 0.5d0 )
*
do k = 1, nblock
* -- Compute damaged stiffness
dft = dmgFiberT(k)
dfc = dmgFiberC(k)
dmt = dmgMatrixT(k)
dmc = dmgMatrixC(k)
df = one - ( one - dft ) * ( one - dfc )
*
dC11 = ( one - df ) * C11
dC22 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C22
dC33 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C33
dC12 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C12
dC23 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C23
dC13 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C13
dG12 = ( one - df )
* * ( one - smt*dmt ) * ( one - smc*dmc ) * G12
dG23 = ( one - df )
* * ( one - smt*dmt ) * ( one - smc*dmc ) * G23
dG13 = ( one - df )
* * ( one - smt*dmt ) * ( one - smc*dmc ) * G13
* -- Stress update
stress(k,i_s33_Xx) = dC11 * strain(k,i_s33_Xx)
* + dC12 * strain(k,i_s33_Yy)
* + dC13 * strain(k,i_s33_Zz)
stress(k,i_s33_Yy) = dC12 * strain(k,i_s33_Xx)
* + dC22 * strain(k,i_s33_Yy)
* + dC23 * strain(k,i_s33_Zz)
stress(k,i_s33_Zz) = dC13 * strain(k,i_s33_Xx)
* + dC23 * strain(k,i_s33_Yy)
* + dC33 * strain(k,i_s33_Zz)
stress(k,i_s33_Xy) = two * dG12 * strain(k,i_s33_Xy)
stress(k,i_s33_Yz) = two * dG23 * strain(k,i_s33_Yz)
stress(k,i_s33_Zx) = two * dG13 * strain(k,i_s33_Zx)
end do
*
return
end

************************************************************
* strainUpdate: Update total strain *
************************************************************
subroutine strainUpdate ( nblock,
* strainInc, strainOld, strainNew )
*
include 'vaba_param.inc'
*
parameter(
* i_s33_Xx = 1,
* i_s33_Yy = 2,
* i_s33_Zz = 3,
* i_s33_Xy = 4,
* i_s33_Yz = 5,
* i_s33_Zx = 6,
* n_s33_Car = 6 )
*
dimension strainInc(nblock,n_s33_Car),
* strainOld(nblock,n_s33_Car),
* strainNew(nblock,n_s33_Car)
*
do k = 1, nblock
strainNew(k,i_s33_Xx)= strainOld(k,i_s33_Xx)
* + strainInc(k,i_s33_Xx)
strainNew(k,i_s33_Yy)= strainOld(k,i_s33_Yy)
* + strainInc(k,i_s33_Yy)
strainNew(k,i_s33_Zz)= strainOld(k,i_s33_Zz)
* + strainInc(k,i_s33_Zz)
strainNew(k,i_s33_Xy)= strainOld(k,i_s33_Xy)
* + strainInc(k,i_s33_Xy)
strainNew(k,i_s33_Yz)= strainOld(k,i_s33_Yz)
* + strainInc(k,i_s33_Yz)
strainNew(k,i_s33_Zx)= strainOld(k,i_s33_Zx)
* + strainInc(k,i_s33_Zx)
end do
*
return
end

************************************************************
* Hashin3d w/ Modified Puck: Evaluate Hashin 3d failure *
* criterion for fiber, Puck for matrix *
************************************************************
subroutine Hashin3d ( nblock, nDmg,
* f1t, f2t, f3t, f1c, f2c, f3c, f12, f23, f13,
* dmgFiberT, dmgFiberC, dmgMatrixT, dmgMatrixC,
* statusMp, stress, eigen )
*
include 'vaba_param.inc'

parameter( zero = 0.d0, one = 1.d0, half = 0.5d0, three =3.d0 )
parameter(
* i_s33_Xx = 1,
* i_s33_Yy = 2,
* i_s33_Zz = 3,
* i_s33_Xy = 4,
* i_s33_Yz = 5,
* i_s33_Zx = 6,
* n_s33_Car = 6 )
*
parameter(i_v3d_X=1,i_v3d_Y=2,i_v3d_Z=3 )
parameter(n_v3d_Car=3 )
*
parameter ( eMax = 1.00d0, eMin = -0.8d0 )
*
dimension dmgFiberT(nblock), dmgFiberC(nblock),
* dmgMatrixT(nblock), dmgMatrixC(nblock),
* stress(nblock,n_s33_Car),
* eigen(nblock,n_v3d_Car),
* statusMp(nblock)
*

*
if ( f1t .gt. zero ) f1tInv = one / f1t
if ( f2t .gt. zero ) f2tInv = one / f2t
if ( f3t .gt. zero ) f3tInv = one / f3t
if ( f1c .gt. zero ) f1cInv = one / f1c
if ( f2c .gt. zero ) f2cInv = one / f2c
if ( f3c .gt. zero ) f3cInv = one / f3c
if ( f12 .gt. zero ) f12Inv = one / f12
if ( f23 .gt. zero ) f23Inv = one / f23
if ( f13 .gt. zero ) f13Inv = one / f13
*
do k = 1, nblock
if ( statusMp(k) .eq. one ) then
*
ldmg = 0
*
s11 = stress(k,i_s33_Xx)
s22 = stress(k,i_s33_Yy)
s33 = stress(k,i_s33_Zz)
s12 = stress(k,i_s33_Xy)
s23 = stress(k,i_s33_Yz)
s13 = stress(k,i_s33_Zx)
*
* Evaluate Fiber modes
if ( s11 .gt. zero ) then
* -- Tensile Fiber Mode
rft = (s11*f1tInv )**2 + (s12*f12Inv )**2 + (s13*f13Inv )**2
if ( rft .lt. one ) then
lDmg = 1
dmgFiberT(k) = one
end if
else if ( s11 .lt. zero ) then
* -- Compressive Fiber Mode
rfc = abs(s11) * f1cInv
if ( rfc .lt. one ) then
lDmg = 1
dmgFiberC(k) = one
end if
end if
*
* Evaluate Matrix Modes
if ( ( s22 + s33 ) .gt. zero ) then
* -- Tensile Matrix mode
rmt = ( s11 * half * f1tInv )**2
* + ( s22**2 * abs(f2tInv * f2cInv) )
* + ( s12 * f12Inv )**2
* + ( s22 * (f2tInv + f2cInv) )
if ( rmt .ge. one ) then
lDmg = 1
dmgMatrixT(k) = one
end if
else if ( ( s22 + s33 ) .lt. zero ) then
* -- Compressive Matrix Mode
rmc = ( s11 * half * f1tInv )**2
* + ( s22**2 * abs(f2tInv * f2cInv) )
* + ( s12 * f12Inv )**2
* + ( s22 * (f2tInv + f2cInv) )
if ( rmc .ge. one ) then
lDmg = 1
dmgMatrixC(k) = one
end if
end if
*
eigMax=max(eigen(k,i_v3d_X),eigen(k,i_v3d_Y),eigen(k,i_v3d_Z))
eigMin=min(eigen(k,i_v3d_X),eigen(k,i_v3d_Y),eigen(k,i_v3d_Z))
enomMax = eigMax - one
enomMin = eigMin - one
*
if ( enomMax .gt. eMax .or.
1 enomMin .lt. eMin .or.
2 dmgFiberT(k) .eq. one .or.
3 dmgFiberC(k) .eq. one .or.
4 dmgMatrixT(k) .eq. one .or.
5 dmgMatrixC(k) .eq. one ) then
statusMp(k) = zero
end if
*
nDmg = nDmg + lDmg
*
end if
*
end do
*
return
end

AM USING THIS VUMAT ??? HOW CAN I ADDD DAMAGE EVOLUTION IN THIS VUMAT PLEASE

Reply
Share: