根据论坛最大应力准则改写VUMAT，计算结果不对，大佬们帮忙看看
************************************************************* define the material of matrix *
************************************************************
subroutine VUMAT(
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'
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)
c
character*80 cmname
parameter( zero = 0.d0, one = 1.d0, two = 2.d0, half = .5d0 )
*
parameter(
* i_svd_DmgMatrixT = 1,
* i_svd_DmgMatrixC = 2,
* i_svd_DmgMatrixS = 3,
* i_svd_statusMp = 4,
* i_svd_dampStress= 5,
c * i_svd_dampStressXx = 5,
c * i_svd_dampStressYy = 6,
c * i_svd_dampStressZz = 7,
c * i_svd_dampStressXy = 8,
c * i_svd_dampStressYz = 9,
c * i_svd_dampStressZx = 10,
* i_svd_Strain = 11,
c * i_svd_StrainXx = 11,
c * i_svd_StrainYy = 12,
c * i_svd_StrainZz = 13,
c * i_svd_StrainXy = 14,
c * i_svd_StrainYz = 15,
c * i_svd_StrainZx = 16,
* n_svd_required = 16 )
*
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 )
*
*the metarial of matrix is ragard as isotropic
*
parameter (
* i_pro_E = 1,
* i_pro_nu = 2,
* i_pro_sigu1t = 3,
* i_pro_sigu1c = 4,
* i_pro_sigus = 5,
*
* i_pro_beta= 6)
*
* * i_pro_rate = 7,
* * i_pro_c = 8 )
* Temporary arrays
*
* Read material properties
*
E = props(i_pro_E)
xnu = props(i_pro_nu)
*
beta = props(i_pro_beta)
*
G=half*E/(one+xnu)
gg = one / ( one - xnu*xnu - xnu*xnu - xnu*xnu
* - two*xnu*xnu*xnu )
C11= E * ( one - xnu*xnu ) * gg
C22= E * ( one - xnu*xnu ) * gg
C33= E * ( one - xnu*xnu ) * gg
C12= E * ( xnu + xnu*xnu ) * gg
C13= E * ( xnu + xnu*xnu ) * gg
C23= E * ( xnu + xnu*xnu ) * gg
*
f1t = props(i_pro_sigu1t)
f1c = props(i_pro_sigu1c)
f1s = props(i_pro_sigus)
*
* 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 isoelastic ( nblock,
* stateOld(1,i_svd_DmgMatrixT),
* stateOld(1,i_svd_DmgMatrixC),stateOld(1,i_svd_DmgMatrixS),
* strainInc, stressNew,
* C11, C22, C33, C12, C23, C13, G)
return
end if
*
*应变更新
call strainUpdate ( nblock,strainInc,
* stateOld(1,i_svd_strain), stateNew(1,i_svd_strain) )
*
* 应力更新
call isoelastic ( nblock,
* stateOld(1,i_svd_DmgMatrixT),
* stateOld(1,i_svd_DmgMatrixC),
* stateOld(1,i_svd_DmgMatrixS),
* stateNew(1,i_svd_strain), stressNew,
* C11, C22, C33, C12, C23, C13, G)
* Failure evaluation
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) )
call copyr ( nblock,
* stateOld(1,i_svd_DmgMatrixS), stateNew(1,i_svd_DmgMatrixS) )
nDmg = 0
call Failure3d( nblock, nDmg,
* f1t, f1c, f1s,
* stateOld(1,i_svd_DmgMatrixT),
* stateOld(1,i_svd_DmgMatrixC),
* stateOld(1,i_svd_DmgMatrixS),
* stateOld(1,i_svd_statusMp),
* stressNew)
* -- Recompute stresses if new Damage is occurring
if ( nDmg .gt. 0 ) then
call isoelastic ( nblock,
* stateNew(1,i_svd_DmgMatrixT),
* stateNew(1,i_svd_DmgMatrixC),
* stateNew(1,i_svd_DmgMatrixS),
* stateNew(1,i_svd_strain), stressNew,
* C11, C22, C33, C12, C23, C13, G)
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 isoelastic ( nblock, strain, stress,
* DmgMatrixC, DmgMatrixS,DmgMatrixT,
* C11, C22, C33, C12, C23, C13, G)
*
include 'vaba_param.inc'
DOUBLE PRECISION ::DAMAGE,MaxDam
*isotropic elasticity, 3D case - （刚度退化计算）
*
parameter( zero = 0.d0, one = 1.d0, two = 2.d0,ten=10.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 )
*
dimensionstrain(nblock,n_s33_Car), stress(nblock,n_s33_Car),
* DmgMatrixT(nblock), DmgMatrixC(nblock), DmgMatrixS(nblock)
*
DAMAGE = ZERO
MaxDam = ZERO
do k = 1, nblock
* -- Stress update
MaxDam=max(DmgMatrixT(k),DmgMatrixT(k),DmgMatrixT(k))
if (MaxDam.eq.one)then
DAMAGE=one/ten
else
DAMAGE=ZERO
endif
*
dC11 = ( one - DAMAGE ) * C11
dC22 = ( one - DAMAGE ) * C22
dC33 = ( one - DAMAGE ) * C33
dC12 = ( one - DAMAGE ) * C12
dC23 = ( one - DAMAGE ) * C23
dC13 = ( one - DAMAGE ) * C13
dG = ( one - DAMAGE ) * G
dG = ( one - DAMAGE ) * G
dG = ( one - DAMAGE ) * G
* -- 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 * dG * strain(k,i_s33_Xy)
stress(k,i_s33_Yz) = two * dG * strain(k,i_s33_Yz)
stress(k,i_s33_Zx) = two * dG * 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
************************************************************
* （失效判据）
************************************************************
subroutine Failure3d ( nblock, nDmg,
* f1t, f1c, f1s, DmgMatrixT, DmgMatrixC,DmgMatrixS,
* statusMp, stress)
*
include 'vaba_param.inc'
PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,NINE=9.0D0,FOUR=4.D0)
PARAMETER (PI=3.141592653589,TOL=0.0001)
DOUBLE PRECISION :: I1,I2,I3,A,B,C,E1,E2,Fai,SP1,SP2,SP3,SSP1,
1 SSP2,SSP3,S1,S2,ST
parameter ( eMax = 1.d0, eMin = -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 )
*
parameter(i_v3d_X=1,i_v3d_Y=2,i_v3d_Z=3 )
parameter(n_v3d_Car=3 )
*
dimensionDmgMatrixS(nblock),
* DmgMatrixT(nblock), DmgMatrixC(nblock),
* stress(nblock,n_s33_Car),
* statusMp(nblock)
c
c
do k = 1, nblock
if ( statusMp(k) .eq. one ) then
I1 = stress(k,1)+stress(k,2)+stress(k,3)
I2 = stress(k,1)*stress(k,2)
1 +stress(k,2)*stress(k,3)
2 +stress(k,1)*stress(k,3)
3 -stress(k,4)**TWO-stress(k,5)**TWO
4 -stress(k,6)**TWO
I3 = stress(k,1)*stress(k,2)*stress(k,3)
1 -stress(k,1)*stress(k,5)**TWO
2 -stress(k,2)*stress(k,6)**TWO
3 -stress(k,3)*stress(k,4)**TWO
4 +TWO*stress(k,4)*stress(k,5)*stress(k,6)
C
C calculate principal stress
A = I1/THREE
E1 = I1**TWO-THREE*I2
IF(E1.GT.TOL)THEN
B=SQRT(E1)
ELSE
C WRITE(*,*) 'Warning: B is less than 0.01,set to 0'
B=ZERO
ENDIF
C
IF(B.GT.0)THEN
E2=(TWO*I1**THREE-NINE*I1*I2+NINE*THREE*I3)/
1 (B**(THREE/TWO)*TWO)
IF(E2.GT.ONE)THEN
C WRITE(*,*) 'Warning:cos(Fai)>1,E2= ', E2, 'adjust to 1'
E2=ONE
ELSEIF(E2.LT.-ONE)THEN
C WRITE(*,*) 'Warning:cos(Fai)<-1,E2= ', E2, 'adjust to -1'
E2=-ONE
ENDIF
Fai=ACOS(E2)/THREE
SP1=A+TWO*B*COS(Fai)/THREE
SP2=A+TWO*B*COS(Fai+TWO*PI/THREE)/THREE
SP3=A+TWO*B*COS(Fai+FOUR*PI/THREE)/THREE
ELSE
SP1=A
SP2=A
SP3=A
ENDIF
S1=SP1+SP2+SP3
S2=SP1*SP2+SP2*SP3+SP3*SP1
ST=SQRT(ABS(TWO*S1**TWO/NINE-TWO*S2/THREE))
C
Cmake sure SP1>=SP2>=SP3
IF((SP1.GT.SP2).and.(SP2.GT.SP3))THEN
SSP1=SP1
SSP2=SP2
SSP3=SP3
ELSEIF ((SP1.GT.SP3).and.(SP3.GT.SP2))THEN
SSP1=SP1
SSP2=SP3
SSP3=SP2
ELSEIF ((SP2.GT.SP1).and.(SP1.GT.SP3))THEN
SSP1=SP2
SSP2=SP1
SSP3=SP3
ELSEIF ((SP2.GT.SP3).and.(SP3.GT.SP1))THEN
SSP1=SP2
SSP2=SP3
SSP3=SP1
ELSEIF ((SP3.GT.SP1).and.(SP1.GT.SP2))THEN
SSP1=SP3
SSP2=SP1
SSP3=SP2
ELSEIF ((SP3.GT.SP2).and.(SP2.GT.SP1))THEN
SSP1=SP3
SSP2=SP2
SSP3=SP1
ENDIF
C write(*,*),'A=',SSP1,SSP2,SSP3
c
IF (SSP3.GT.ZERO.AND.SSP1.GT.f1t) THEN
matrixdmg = one
DmgMatrixT(k) = one
END IF
C
IF (SSP1.LT.ZERO.AND.ABS(SSP3).GT.f1c) THEN
matrixdmg = one
DmgMatrixC(k) = one
END IF
c
IF (ABS(ST).GT.f1s) THEN
matrixdmg = one
DmgMatrixS(k) = one
END IF
*
nDmg =matrixdmg
*
if (DmgMatrixT(k).eq.one.or.
* DmgMatrixT(k).eq.one.or.
* DmgMatrixT(k).eq.one)then
statusMp(k) = zero
end if
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
************************************************************
* 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
有大佬帮忙看下么
页:
[1]