taishanbuzuo 发表于 2012-5-10 22:46 static/image/common/back.gif
VUMAT可以用全积分单元吗？不可以吧，可以计算，但是结果不对
楼主有碰到11楼的问题吗，我现在有和十一楼类似的问题~还请指教啊
你好，我最近也在看hashin判据，但是有一个疑问，就是关于特征长度的，自带的2维hashin判据中，是不是程序自己计算特征长度，不用人为定义。另外cohesive单元的失效和这个hashin的线性软化是不是很类似的，都定义了能量释放率，但是cohesive中定义了单元厚度，类似于单元特征长度。但是如果是自带计算特征长度的话，那么在能量释放率定义的时候，数值应该是有范围限制的。否则最大应力点对应的位移有可能大于总位移的。
请指教。
taishanbuzuo 发表于 2012-3-21 23:16
现在也没什么好办法调试 ，每次编号了就运行一下，看提示什么错误，根据错误再去修改程序 ...
您好，如何给程序编号呢，我是初学者，不会调试程序，您那能附上一张图片吗
您好，我这调用子程序的时候出现了这样的情况，貌似子程序没有运算，调用和没调用子程序算出的结果一样，怎么回事
taishanbuzuo 发表于 2012-5-10 22:46
VUMAT可以用全积分单元吗？不可以吧，可以计算，但是结果不对
“VUMAT不能用全积分单元”？这个怎么说？
我自己做了一个对比：三组模型，分别用杆单元、壳单元S4R以及实体单元C3D8R，用的是纯弹性材料。然后我把自己根据广义胡克定律编的VUMAT分别应用于这三组模型中，然后各自与ABAQUS自带的弹性材料模型作比较，发现计算结果中，用C3D8R的模型，VUMAT的与ABAQUS自带的弹性材料本构的结果完全一样，而其余两种单元的模型结果，却有差异，尤其是壳单元的。
不知道你有时间能否看一下我的帖子http://forum.simwe.com/thread-1118927-1-1.html
谢谢。
yespin 发表于 2014-12-31 13:31
“VUMAT不能用全积分单元”？这个怎么说？
我自己做了一个对比：三组模型，分别用杆单元、壳单元S4R以及 ...
完全积分单元可以实现对角质量阵吗
师兄，能加您qq吗？最近在做复材渐进损伤仿真，新手，很多不懂，请教您
大神，我在提示"THE PARAMETER COMPOSITE CAN NOT BE USED IN Abaqus/Explicit"，意思是explicit不支持复合材料体单元，你用vumat是怎么解决这个问题的？复合材料难道不是用体单元？
子程序不便上传，能否上传 一个inp文件？
e63liuling 发表于 2015-6-2 22:59
大神，我在提示"THE PARAMETER COMPOSITE CAN NOT BE USED IN Abaqus/Explicit"，意思是explicit不支持复合 ...
explict是支持体单元的，但是不支持abaqus自带的三维复合材料模型，需要自己编写vumat
楼主您好，目前正在做相关的学习研究，您可以具体分享下您的算例吗？或者做一个相对简单的模型分享个算例，可以吗？里面有很多小细节真的很难搞懂
问下楼主你之前发的一个unifiber.f的三维hashin子程序如何输入参数的，我不知道怎么使材料破坏
不发源代码，可以共享下编译好的for文件吗？？:victory:
sagitarhncs 发表于 2015-12-11 17:04
不发源代码，可以共享下编译好的for文件吗？？
想学习复合材料VUMAT的，可以参照以下帖子，有源代码：
http://forum.simwe.com/thread-1124382-1-1.html
前来学习
amani 发表于 2012-3-18 12:34
Would you please share your subroutine VUMAT?
既然你想让别人跟你分享，你自己又和网友们分享了啥，为你羞愧
本帖最后由 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
浙大 刘鹏飞
amani 发表于 2016-11-2 10:22
分享我最近2016年发表的复合材料层合板低速冲击ABAQUS-VUMAT有限元程序（HASHIN失效准则）
Liu PF, Liao B ...
老师您好!
最近细细的研读了您2017年发表的考虑塑性的复合材料层合板低速冲击损伤的文章，其中关于屈服判断相关的问题一直不太明白。是一开始冲击就产生屈服吗？开始冲击时等效塑性应变为零，这时屈服函数必定大于零。
本帖最后由 HXD900910 于 2017-3-10 21:02 编辑
amani 发表于 2016-11-2 10:22
分享我最近2016年发表的复合材料层合板低速冲击ABAQUS-VUMAT有限元程序（HASHIN失效准则）
Liu PF, Liao B ...
C:\Users\Administrator\Desktop
HXD900910 发表于 2017-3-10 20:59
老师您好!
最近细细的研读了您2017年发表的考虑塑性的复合材料层合板低速冲击损伤的文章，其中关于 ...
这个问题一直没弄懂，老师能解释解释吗，谢谢您