复合材料VUMAT之三维hashin准则
本帖最后由 taishanbuzuo 于 2015-12-14 18:37 编辑VUMAT终于调试成功，现将计算结果贴一下：
复合材料的失效准则采用的三维HASHIN准则，损伤演化段位最简单的线性演化。与试验结果吻合的非常不错。
再次提醒大家，原有的场变量子程序有很大的缺点，虽然编程序比较简单，但是计算结果不好。当单元达到损伤准则以后，不是一下子就彻底失去承载能力，而是需要一段下降段。考虑下降段与不考虑下降段结果相差很大。
子程序不便上传，先传几幅图，用以交流。
一张图片是前四层基体拉伸损伤情况
一张为前四层纤维拉伸损伤情况
另一张为某单元应力应变曲线。
想学习复合材料VUMAT的，可以参照以下帖子，有源代码：http://forum.simwe.com/thread-1124382-1-1.html
本帖最后由 amani 于 2016-11-2 10:24 编辑
971619083 发表于 2016-10-31 11:02
既然你想让别人跟你分享，你自己又和网友们分享了啥，为你羞愧
分享我最近2016年发表的复合材料层合板低速冲击ABAQUS-VUMAT有限元程序（HASHIN失效准则）
Liu PF, Liao BB, Jia LY, Peng XQ. Finite element analysis of dynamic progressive failure of carbon fiber composite laminates under low velocity impact. Composite Structures, 2016, 149: 408-422.
subroutine vumat(
c Read only -
1nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
2stepTime, totalTime, dt, cmname, coordMp, charLength,
3props, density, strainInc, relSpinInc,
4tempOld, stretchOld, defgradOld, fieldOld,
5stressOld, stateOld, enerInternOld, enerInelasOld,
6tempNew, stretchNew, defgradNew, fieldNew,
c Write only -
7stressNew, 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),
1coordMp(nblock,*),
2charLength(*), strainInc(nblock,ndir+nshr),
3relSpinInc(nblock,nshr), tempOld(nblock),
4stretchOld(nblock,ndir+nshr), defgradOld(nblock,ndir+nshr+nshr),
5fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6stateOld(nblock,nstatev), enerInternOld(nblock),
7enerInelasOld(nblock), tempNew(*),
8stretchNew(nblock,ndir+nshr), defgradNew(nblock,ndir+nshr+nshr),
9fieldNew(nblock,nfieldv), stressNew(nblock,ndir+nshr),
1stateNew(nblock,nstatev),
2enerInternNew(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 = 11,
* i_pro_sigu1c = 12,
* i_pro_sigu2t = 13,
* i_pro_sigu2c = 14,
* i_pro_sigu3t = 15,
* i_pro_sigu3c = 16,
* i_pro_sigu12 = 17,
* i_pro_sigu13 = 18,
* i_pro_sigu23 = 19,
* i_pro_G1t = 20,
* i_pro_G1c = 21,
* i_pro_G2t = 22,
* i_pro_G2c = 23,
* i_pro_Gs = 24)
* 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)
G1t=props(i_pro_G1t)
G1c=props(i_pro_G1c)
G2t=props(i_pro_G2t)
G2c=props(i_pro_G2c)
Gs=props(i_pro_Gs)
*
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= E2 * ( xnu12-xnu32*xnu13 ) * gg
C13= E3 * ( xnu13 - xnu12*xnu23 ) * gg
C23= E3 * ( xnu23 - xnu21*xnu13 ) * gg
C21= E1 * ( xnu21-xnu31*xnu23 ) * gg
C31= E1 * ( xnu31 - xnu21*xnu32 ) * gg
C32= E2 * ( xnu32 - xnu12*xnu31 ) * gg
C44=G12
C55=G23
C66=G13
*
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, C21,C32,C31,C44,C55,C66,
* 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, C21,C32,C31,C44,C55,C66,
* stateNew(1,i_svd_strain),
* stressNew )
*
nDmg = 0
call eig33Anal ( nblock, stretchNew, eigen )
call Hashin3d( nblock, nDmg,
* f1t, f2t, f3t, f1c, f2c, f3c, f12, f23, f13,
* C11, C22, C33, C12, C23, C13, C21,C32,C31,C44,C55,C66,
* G1t,G1c,G2t,G2c,Gs,
* stateNew(1,i_svd_statusMp),
* stateNew(1,i_svd_strain), eigen,stressNew,charLength,
* stateOld(1,i_svd_DmgFiberT),
* stateOld(1,i_svd_DmgFiberC),
* stateOld(1,i_svd_DmgMatrixT),
* stateOld(1,i_svd_DmgMatrixC),
* stateNew(1,i_svd_DmgFiberT),
* stateNew(1,i_svd_DmgFiberC),
* stateNew(1,i_svd_DmgMatrixT),
* stateNew(1,i_svd_DmgMatrixC))
* -- 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, C21,C32,C31,C44,C55,C66,
* 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, C21,C32,C31,C44,C55,C66,
* strain, stress )
*
include 'vaba_param.inc'
*Orthotropic elasticity, 3D case -
*
parameter( zero = 0.d0, one = 1.d0, two = 2.d0)
parameter ( smt = 0.9d0, smc = 0.5d0 )
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 )
*
dimensionstrain(nblock,n_s33_Car),
* dmgFiberT(nblock), dmgFiberC(nblock),
* dmgMatrixT(nblock), dmgMatrixC(nblock),
* stress(nblock,n_s33_Car),stresstr(nblock,n_s33_Car)
* -- shear fraction in matrix tension and compression mode
*
do k = 1, nblock
* -- Compute damaged stiffness
c WRITE(*,'(1X,''dmgFiberT(k)'',E12.5)')dmgFiberT(k)
c WRITE(*,'(1X,''dmgFiberC(k)'',E12.5)')dmgFiberC(k)
c WRITE(*,'(1X,''dmgMatrixT(k)'',E12.5)')dmgMatrixT(k)
c WRITE(*,'(1X,''dmgMatrixC(k)'',E12.5)')dmgMatrixC(k)
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 ) * C44
dG23 = ( one - df )
* * ( one - smt*dmt ) * ( one - smc*dmc ) * C55
dG13 = ( one - df )
* * ( one - smt*dmt ) * ( one - smc*dmc ) * C66
* -- 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)
c WRITE(*,'(1X,''stress(k,i_s33_Xx)'',E12.5)')stress(k,i_s33_Xx)
c WRITE(*,'(1X,''stress(k,i_s33_Yy)'',E12.5)')stress(k,i_s33_Yy)
c WRITE(*,'(1X,''stress(k,i_s33_Zz)'',E12.5)')stress(k,i_s33_Zz)
c WRITE(*,'(1X,''stress(k,i_s33_Xy)'',E12.5)')stress(k,i_s33_Xy)
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 compress *
************************************************************
subroutine Hashin3d( nblock,nDmg,
* f1t,f2t, f3t,f1c, f2c, f3c, f12, f23, f13,
* C11, C22, C33, C12, C23, C13, C21,C32,C31,C44,C55,C66,
* G1t,G1c,G2t,G2c,Gs,
* statusMp, strain, eigen,stressNew,charLength,
* dmgftold,dmgfcold,dmgmtold,dmgmcold,
* dmgftnew,dmgfcnew,dmgmtnew,dmgmcnew)
*
include 'vaba_param.inc'
parameter( zero = 0.d0, one = 1.d0,two=2.d0, half = 0.5d0,
* three = 3.d0,pi=3.14,fz=0.99 )
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 dmgftold(nblock),dmgfcold(nblock),
* dmgmtold(nblock), dmgmcold(nblock),
* dmgftnew(nblock),dmgfcnew(nblock),
* dmgmtnew(nblock), dmgmcnew(nblock),
* dft(nblock),dfc(nblock),
* dmt(nblock), dmc(nblock),
* stressNew(nblock,n_s33_Car),
* strain(nblock,n_s33_Car),statusMp(nblock),
* eigen(nblock,n_v3d_Car),charLength(nblock)
*
if ( f1t .gt. zero ) E01t = f1t/C11
if ( f2t .gt. zero ) E02t = f2t/C22
if ( f3t .gt. zero ) E03t = f3t/C33
if ( f1c .gt. zero ) E01c = f1c/C11
if ( f2c .gt. zero ) E02c = f2c/C22
if ( f3c .gt. zero ) E03c = f3c/C33
if ( f12 .gt. zero ) E012 = f12/(2*C44)
if ( f23 .gt. zero ) E023 = f23/(2*C55)
if ( f13 .gt. zero ) E013 = f13/(2*C66)
c WRITE(*,'(1X,''f1t'',E12.5)')f1t
c WRITE(*,'(1X,''E1'',E12.5)')E1
c WRITE(*,'(1X,''f1c'',E12.5)')f1c
c WRITE(*,'(1X,''E2'',E12.5)')E2
c WRITE(*,'(1X,''f2t'',E12.5)')f2t
c WRITE(*,'(1X,''E3'',E12.5)')E3
c WRITE(*,'(1X,''E01t'',E12.5)')E01t
c WRITE(*,'(1X,''E02t'',E12.5)')E02t
c WRITE(*,'(1X,''E01c'',E12.5)')E01c
c WRITE(*,'(1X,''E02c'',E12.5)')E02c
*
lDmg = 0
do k = 1, nblock
dft(k) = zero
dfc(k) = zero
dmt(k) = zero
dmc(k) = zero
dmgftnew(k) = zero
dmgfcnew(k) = zero
dmgmtnew(k) = zero
dmgmcnew(k) = zero
END DO
c WRITE(*,'(1X,''statusMp(k)'',E12.5)')statusMp(k)
c do k = 1, nblock
c statusMp(k) = one
c end do
do k = 1, nblock
c if ( statusMp(k) .eq. one ) then
*
c lFail = 0
*
E11 = strain(k,i_s33_Xx)
E22 = strain(k,i_s33_Yy)
E33 = strain(k,i_s33_Zz)
E12 = strain(k,i_s33_Xy)
E23 = strain(k,i_s33_Yz)
E13 = strain(k,i_s33_Zx)
C WRITE(*,'(1X,''E11'',E12.5)')E11
C WRITE(*,'(1X,''E22'',E12.5)')E22
*
*
* Evaluate Fiber modes
if ( E11 .gt. zero ) then
* -- Tensile Fiber Mode
rft =(E11/E01t)**2
if ( rft .ge. one ) then
lDmg1 = 1
Ef1t=2*G1t/(f1t*charLength(k))
dft(k)=Ef1t*(E11-E01t)/(E11*(Ef1t-E01t))
if (dft(k).ge. one) then
dft(k)=1
end if
if (dft(k).le. zero) then
dft(k)=0
end if
c WRITE(*,'(1X,''Ef1t'',E12.5)')Ef1t
c WRITE(*,'(1X,''dmgFiberT(k)'',E12.5)')dmgFiberT(k)
end if
c WRITE(*,'(1X,''Ef1t'',E12.5)')Ef1t
c WRITE(*,'(1X,''dmgFiberT(k)'',E12.5)')dmgFiberT(k)
c WRITE(*,'(1X,''lDmg1'',E12.5)')lDmg1
else if ( E11 .lt. zero ) then
* -- Compressive Fiber Mode
rfc = (E11/E01C)**2
if ( rfc .ge. one ) then
lDmg2 = 1
Ef1c=2*G1c/(f1c*charLength(k))
dfc(k)=Ef1c*(E11-E01c)/(E11*(Ef1c-E01c))
if (dfc(k).ge. one) then
dfc(k)=1
end if
if (dfc(k).le. zero) then
dfc(k)=0
end if
end if
c WRITE(*,'(1X,''Ef1c'',E12.5)')Ef1c
c WRITE(*,'(1X,''dmgFiberC(k)'',E12.5)')dmgFiberC(k)
c WRITE(*,'(1X,''lDmg2'',E12.5)')lDmg2
end if
*
* Evaluate Matrix Modes
if ( ( E22 + E33 ) .gt. zero ) then
* -- Tensile Matrix mode
rmt = (E22+E33)**2/(E02t**2)
* - (E22*E33)/(E023**2)
* + ( E12/E012)**2+(E13/E013)**2+(E23/E023)**2
if ( rmt .ge. one ) then
lDmg3 = 1
Ef2t=2*G2t/(f2t*charLength(k))
dmt(k)=Ef2t*(E22-E02t)/(E22*(Ef2t-E02t))
if (dmt(k).ge. one) then
dmt(k)=1
end if
if (dmt(k).le. zero) then
dmt(k)=0
end if
end if
else if ( ( E22+E33 ) .lt. zero ) then
* -- Compressive Matrix Mode
rmc = ( E22+E33 )**2/(E02c*E03c)
* + (E22+E33 )*(E02C/(2*E012)-1)/E02C
* - (E22*E33)/(E023**2)
* + ( E12/E012)**2+(E13/E013)**2+(E23/E023)**2
if ( rmc .ge. one ) then
lDmg4 = 1
Ef2c=2*G2c/(f2c*charLength(k))
dmc(k)=Ef2c*(E22-E02c)/(E22*(Ef2c-E02c))
if (dmc(k).ge. one) then
dmc(k)=1
end if
if (dmc(k).le. zero) then
dmc(k)=0
end if
end if
end if
dmgftnew(k)=MAX(dmgftold(k),dft(k))
dmgfcnew(k)=MAX(dmgfcold(k),dfc(k))
dmgmtnew(k)=MAX(dmgmtold(k),dmt(k))
dmgmcnew(k)=MAX(dmgmcold(k),dmc(k))
C WRITE(*,'(1X,''rft'',E12.5)')rft
C WRITE(*,'(1X,''rfc'',E12.5)')rfc
*
c eigMax=max(eigen(k,i_v3d_X),eigen(k,i_v3d_Y),eigen(k,i_v3d_Z))
c eigMin=min(eigen(k,i_v3d_X),eigen(k,i_v3d_Y),eigen(k,i_v3d_Z))
c enomMax = eigMax - one
c enomMin = eigMin - one
*
if ( dmgftnew(k).gt.fz.or.
* dmgfcnew(k).gt.fz) then
statusMp(k) = zero
else
statusMp(k) = one
end if
c WRITE(*,'(1X,''statusMp(k)'',E12.5)')statusMp(k)
*
nDmg = nDmg + lDmg1+ lDmg2+ lDmg3+ lDmg4
*
c 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 )
*
dimensionsigOld(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
**************************************************************************
* 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
浙大 刘鹏飞
yespin 发表于 2014-12-31 13:31
“VUMAT不能用全积分单元”？这个怎么说？
我自己做了一个对比：三组模型，分别用杆单元、壳单元S4R以及 ...
完全积分单元可以实现对角质量阵吗 shannonlixi 发表于 2012-6-26 22:58 static/image/common/back.gif
你好，我有个问题，abaqus自带的hashin准则计算原理也是你的图中的曲线，为什么还要自己编写vumat程序？ ...
hashin准则是二维的损伤准则，不适合用于厚板，也不适于用于分析考虑层间效应的问题。 你好，我有个问题，abaqus自带的hashin准则计算原理也是你的图中的曲线，为什么还要自己编写vumat程序？ Would you please share your subroutine VUMAT? 什么是“原有的场变量子程序”，是帮助里面的吗？ 从图上看效果还是不错的，可以说得具体一些吗？期待向你学习 给予支持和鼓舞 divinelove 发表于 2012-3-19 10:33 static/image/common/back.gif
什么是“原有的场变量子程序”，是帮助里面的吗？
以前的场变量子程序都是瞬间衰减的，将模量瞬间衰减到0或者一个很小的值，
与实际情况还是有区别的。
UMAT和VUMAT可以根据实际情况编写下降段。 本帖最后由 divinelove 于 2012-3-21 15:51 编辑
taishanbuzuo 发表于 2012-3-21 12:15 http://forum.simwe.com/static/image/common/back.gif
以前的场变量子程序都是瞬间衰减的，将模量瞬间衰减到0或者一个很小的值，
与实际情况还是有区别的。
UMA ...
我下一步打算编一个vumat出来，带有下降的，今天刚刚把那个自己带umat的改成自己需要的了，但是不知道该怎么调试，你是怎么调的？ divinelove 发表于 2012-3-21 15:50 static/image/common/back.gif
我下一步打算编一个vumat出来，带有下降的，今天刚刚把那个自己带umat的改成自己需要的了，但是不知道该 ...
现在也没什么好办法调试 ，每次编号了就运行一下，看提示什么错误，根据错误再去修改程序 请问，在vumat中怎样查询到层号和单元号呢？ szk8617 发表于 2012-4-4 16:56 static/image/common/back.gif
请问，在vumat中怎样查询到层号和单元号呢？
这个 不知道，需要用到这个吗？ taishanbuzuo 发表于 2012-4-4 21:12 static/image/common/back.gif
这个 不知道，需要用到这个吗？
恭喜了，我发过一个帖子http://forum.simwe.com/forum.php?mod=viewthread&tid=1029132&highlight=vumat
不知道是否遇到过相似问题？我试用了full integration 和 reduced integration, 冲击过程中单元都出现了不规则扭曲，这可能是shear locking 或 hourglass作用，而且所有内部单元均定义了接触，实在找不出原因，现在还在调试vumat，我不清楚是不是只要定义了damage evolution， 在冲击（bending）载荷下就会出现单元的不规则扭曲呢？ 恩，我的进展还没达到楼上的水平，停顿在查询abaqus调用vumat时，所处的层号和单元号上了。
因为我的刚度退化模式是基于单层中基体裂纹密度的。在获取单元编号时，从国外论坛上查得：
nblock :Number of material points to be processed in this call to VUMAT.
The material points are the integration points (Gauss points).You
can't control or change NBLOCK -- ABAQUS selects a value based, I
assume, on achieving a balance between memory/storage requirements and
efficiency.You can't get element numbers in VUMAT, but you do get the
material name (CMNAME).If you really need the element number, you can
use a field or state variable to store the element number.
可以通过定义场变量来存储获得，有谁知道这怎样定义吗？
请问楼主用的是什么evolution law？用到characteristic length没？
如果用到，是自己定义的还是用的vumat里的charlength？ 楼主和大神们好，我最近在写复合材料的umat. 不知道能否加下QQ讨论讨论。
QQ: 51380358 shiyu_vincent 发表于 2012-4-6 06:05 static/image/common/back.gif
恭喜了，我发过一个帖子http://forum.simwe.com/forum.php?mod=viewthread&tid=1029132&highlight=vumat
...
VUMAT可以用全积分单元吗？不可以吧，可以计算，但是结果不对 新手请教楼主：定义层合板属性时因为要定义密度，depvar里面设置的两个变量数是不是要比User Material里面设置的变量数多一个啊？ taishanbuzuo 发表于 2012-6-26 23:17 static/image/common/back.gif
hashin准则是二维的损伤准则，不适合用于厚板，也不适于用于分析考虑层间效应的问题。 ...
哦，明白了，自带的只能用平面应力单元做，呵呵，师兄是西工大的，以后有问题还要请教啊，呵呵 taishanbuzuo 发表于 2012-6-26 23:17 static/image/common/back.gif
hashin准则是二维的损伤准则，不适合用于厚板，也不适于用于分析考虑层间效应的问题。 ...
你好，我一直不知道自带的hashin准则计算的时候，刚度退化的系数是怎么来的，请问可以指教一下么？谢谢啊