Here is an example of VUMAT using Hashin criteria for composites. Rumors said it was from the company. This is re-posted from nabble.com. I have not test it. And I also found the same VUMAT in the appendix of a PhD thesis (please refer to the foot of this post).
c http://abaqus-users.1086179.n5.nabble.com/Hashin-damage-theory-td12032.html
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)
*
f1tInv = zero
f2tInv = zero
f3tInv = zero
f1cInv = zero
f2cInv = zero
f3cInv = zero
f12Inv = zero
f23Inv = zero
f13Inv = zero
*
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
*
lFail = 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 .ge. one ) then
lDmg = 1
dmgFiberT(k) = one
end if
else if ( s11 .lt. zero ) then
* -- Compressive Fiber Mode
rfc = abs(s11) * f1cInv
if ( rfc .ge. 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.
* enomMin .lt. eMin .or.
* dmgFiberT(k) .eq. one ) then
statusMp(k) = zero
end if
*
nDmg = nDmk + lDmg
*
end if
*
end do
*
return
end
************************************************************
* betaDamping: Add beta damping *
************************************************************
subroutine betaDamping3d ( nblock,
* beta, dt, strainInc, sigOld, sigNew,
* statusMp, sigDampOld, sigDampNew )
*
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 sigOld(nblock,n_s33_Car),
* sigNew(nblock,n_s33_Car),
* strainInc(nblock,n_s33_Car),
* statusMp(nblock),
* sigDampOld(nblock,n_s33_Car),
* sigDampNew(nblock,n_s33_Car)
*
parameter ( zero = 0.d0, one = 1.d0, two=2.0d0,
* half = 0.5d0, third = 1.d0/3.d0 )
parameter ( asmall = 1.d-16 )
*
betaddt = beta / dt
*
do k =1 , nblock
sigDampNew(k,i_s33_Xx) = betaddt * statusMp(k) *
* ( sigNew(k,i_s33_Xx)
* - ( sigOld(k,i_s33_Xx) - sigDampOld(k,i_s33_Xx) ) )
sigDampNew(k,i_s33_Yy) = betaddt * statusMp(k) *
* ( sigNew(k,i_s33_Yy)
* - ( sigOld(k,i_s33_Yy) - sigDampOld(k,i_s33_Yy) ) )
sigDampNew(k,i_s33_Zz) = betaddt * statusMp(k) *
* ( sigNew(k,i_s33_Zz)
* - ( sigOld(k,i_s33_Zz) - sigDampOld(k,i_s33_Zz) ) )
sigDampNew(k,i_s33_Xy) = betaddt * statusMp(k) *
* ( sigNew(k,i_s33_Xy)
* - ( sigOld(k,i_s33_Xy) - sigDampOld(k,i_s33_Xy) ) )
sigDampNew(k,i_s33_Yz) = betaddt * statusMp(k) *
* ( sigNew(k,i_s33_Yz)
* - ( sigOld(k,i_s33_Yz) - sigDampOld(k,i_s33_Yz) ) )
sigDampNew(k,i_s33_Zx) = betaddt * statusMp(k) *
* ( sigNew(k,i_s33_Zx)
* - ( sigOld(k,i_s33_Zx) - sigDampOld(k,i_s33_Zx) ) )
*
sigNew(k,i_s33_Xx) = sigNew(k,i_s33_Xx)+sigDampNew(k,i_s33_Xx)
sigNew(k,i_s33_Yy) = sigNew(k,i_s33_Yy)+sigDampNew(k,i_s33_Yy)
sigNew(k,i_s33_Zz) = sigNew(k,i_s33_Zz)+sigDampNew(k,i_s33_Zz)
sigNew(k,i_s33_Xy) = sigNew(k,i_s33_Xy)+sigDampNew(k,i_s33_Xy)
sigNew(k,i_s33_Yz) = sigNew(k,i_s33_Yz)+sigDampNew(k,i_s33_Yz)
sigNew(k,i_s33_Zx) = sigNew(k,i_s33_Zx)+sigDampNew(k,i_s33_Zx)
*
end do
*
return
end
************************************************************
* EnergyInternal3d: Compute internal energy for 3d case *
************************************************************
subroutine EnergyInternal3d(nblock, sigOld, sigNew ,
* strainInc, curDensity, enerInternOld, enerInternNew)
*
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 )
*
parameter( two = 2.d0, half = .5d0 )
*
dimension sigOld (nblock,n_s33_Car), sigNew (nblock,n_s33_Car),
* strainInc (nblock,n_s33_Car), curDensity (nblock),
* enerInternOld(nblock), enerInternNew(nblock)
*
do k = 1, nblock
stressPower = half * (
* ( sigOld(k,i_s33_Xx) + sigNew(k,i_s33_Xx) )
* * ( strainInc(k,i_s33_Xx) )
* + ( sigOld(k,i_s33_Yy) + sigNew(k,i_s33_Yy) )
* * ( strainInc(k,i_s33_Yy))
* + ( sigOld(k,i_s33_Zz) + sigNew(k,i_s33_Zz) )
* * ( strainInc(k,i_s33_Zz))
* + two * ( sigOld(k,i_s33_Xy) + sigNew(k,i_s33_Xy) )
* * strainInc(k,i_s33_Xy)
* + two * ( sigOld(k,i_s33_Yz) + sigNew(k,i_s33_Yz) )
* * strainInc(k,i_s33_Yz)
* + two * ( sigOld(k,i_s33_Zx) + sigNew(k,i_s33_Zx) )
* * strainInc(k,i_s33_Zx) )
*
enerInternNew(k) = enerInternOld(k) + stressPower/curDensity(k)
end do
*
return
end
************************************************************
* CopyR: Copy from one array to another *
************************************************************
subroutine CopyR(nCopy, from, to )
*
include 'vaba_param.inc'
*
dimension from(nCopy), to(nCopy)
*
do k = 1, nCopy
to(k) = from(k)
end do
*
return
end
*********************************************************************
*****
* eig33Anal: Compute eigen values of a 3x3 symmetric matrix analytically *
*********************************************************************
*****
subroutine eig33Anal( nblock, sMat, eigVal )
*
include 'vaba_param.inc'
*
parameter(i_s33_Xx=1,i_s33_Yy=2,i_s33_Zz=3 )
parameter(i_s33_Xy=4,i_s33_Yz=5,i_s33_Zx=6 )
parameter(i_s33_Yx=i_s33_Xy )
parameter(i_s33_Zy=i_s33_Yz )
parameter(i_s33_Xz=i_s33_Zx,n_s33_Car=6 )
*
parameter(i_v3d_X=1,i_v3d_Y=2,i_v3d_Z=3 )
parameter(n_v3d_Car=3 )
*
parameter ( zero = 0.d0, one = 1.d0, two = 2.d0,
* three = 3.d0, half = 0.5d0, third = one / three,
* pi23 = 2.094395102393195d0,
* fuzz = 1.d-8,
* preciz = fuzz * 1.d4 )
*
dimension eigVal(nblock,n_v3d_Car), sMat(nblock,n_s33_Car)
*
do k = 1, nblock
sh = third*(sMat(k,i_s33_Xx)+sMat(k,i_s33_Yy)+sMat(k,i_s33_Zz))
s11 = sMat(k,i_s33_Xx) - sh
s22 = sMat(k,i_s33_Yy) - sh
s33 = sMat(k,i_s33_Zz) - sh
s12 = sMat(k,i_s33_Xy)
s13 = sMat(k,i_s33_Xz)
s23 = sMat(k,i_s33_Yz)
*
fac = max(abs(s11), abs(s22), abs(s33))
facs = max(abs(s12), abs(s13), abs(s23))
if( facs .lt. (preciz*fac) ) then
eigVal(k,i_v3d_X) = sMat(k,i_s33_Xx)
eigVal(k,i_v3d_Y) = sMat(k,i_s33_Yy)
eigVal(k,i_v3d_Z) = sMat(k,i_s33_Zz)
else
q = third*((s12**2+s13**2+s23**2)+half*(s11**2+s22**2+s33**2))
fac = two * sqrt(q)
if( fac .gt. fuzz ) then
ofac = two/fac
else
ofac = zero
end if
s11 = ofac*s11
s22 = ofac*s22
s33 = ofac*s33
s12 = ofac*s12
s13 = ofac*s13
s23 = ofac*s23
r = s12*s13*s23
* + half*(s11*s22*s33-s11*s23**2-s22*s13**2-s33*s12**2)
if( r .ge. one-fuzz ) then
cos1 = -half
cos2 = -half
cos3 = one
else if( r .le. fuzz-one ) then
cos1 = -one
cos2 = half
cos3 = half
else
ang = third * acos(r)
cos1 = cos(ang)
cos2 = cos(ang+pi23)
cos3 =-cos1-cos2
end if
eigVal(k,i_v3d_X) = sh + fac*cos1
eigVal(k,i_v3d_Y) = sh + fac*cos2
eigVal(k,i_v3d_Z) = sh + fac*cos3
end if
end do
*
return
end
The mentioned dissertation is “ADVANCED MESOMECHANICAL MODELING OF TRIAXIALLY BRAIDED
COMPOSITES FOR DYNAMIC IMPACT ANALYSIS WITH FAILURE” by Zifeng Nie from The Graduate Faculty of The University of Akron. This thesis is public accessible. Ask Google 😀©
本文发表于水景一页。永久链接：<https://cnzhx.net/fe/2015/05/09/an-example-of-vumat-using-hashin-criteria-for-composites/>。转载请保留此信息及相应链接。