| 查看: 1344 | 回复: 1 | ||
| 【悬赏金币】回答本帖问题，作者zhuhao3802将赠送您 10 个金币 | ||
[求助]
根据simwe论坛最大应力准则改写VUMAT，计算结果不对，大佬们帮忙看看
|
||
|
************************************************************ * define the material of matrix * ************************************************************ 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, c Write only - 7 stressNew, stateNew, enerInternNew, enerInelasNew ) c include 'vaba_param.inc' 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) 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 :AMAGE,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 ) * dimension strain(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 ) * dimension DmgMatrixS(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 C make 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 ) * 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 |
回复此楼» 猜你喜欢
《岩石力学与工程学报》修改后增刊复审，申请正刊的几率大吗？
已经有11人回复
-
Sustainable Engineering Materials首届青年编委招募启事
已经有17人回复
-
建筑环境与结构工程论文润色/翻译怎么收费?
已经有200人回复
-
混凝土材料类SCI期刊求推荐
已经有7人回复
-
哈工大申请考核制招收博士生9月28日截止
已经有20人回复
-
哈工大26年度博士申请考核第二轮
已经有6人回复
-
advance in structural engineering 期刊Awaiting EIC Assignmeng阶段十天了
已经有3人回复
-
ansys有限元分析求助
已经有2人回复
2楼2020-10-21 19:23:25