| 2 | 1/1 | 返回列表 |
| 查看: 1371 | 回復(fù): 1 | ||
| 【懸賞金幣】回答本帖問題,作者zhuhao3802將贈送您 10 個金幣 | ||
zhuhao3802鐵蟲 (初入文壇)
|
[求助]
根據(jù)simwe論壇最大應(yīng)力準(zhǔn)則改寫VUMAT,計算結(jié)果不對,大佬們幫忙看看
|
|
|
************************************************************ * 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 * * 應(yīng)變更新 call strainUpdate ( nblock,strainInc, * stateOld(1,i_svd_strain), stateNew(1,i_svd_strain) ) * * 應(yīng)力更新 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 ************************************************************ * (失效判據(jù)) ************************************************************ 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 |
鐵蟲 (初入文壇)
| 2 | 1/1 | 返回列表 |
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 環(huán)境調(diào)劑 +8 | chenhanheng 2026-03-02 | 8/400 |
|
|---|---|---|---|---|
|
[考研] 0856材料工程,初試313調(diào)劑 +7 | 賣個關(guān)子吧 2026-03-03 | 7/350 |
|
|
[考研] 理學(xué),工學(xué),農(nóng)學(xué)調(diào)劑,少走彎路,這里歡迎您! +8 | likeihood 2026-03-02 | 11/550 |
|
|
[考研] 085700資環(huán)求調(diào)劑,初始279,六級已過,英語能力強 +3 | 085700資環(huán)調(diào)劑 2026-03-03 | 4/200 |
|
|
[考研] 266材料化工求調(diào)劑 +3 | 哇塞王帥 2026-03-03 | 3/150 |
|
|
[考研] 一志愿中科大能動297求調(diào)劑,本科川大 +3 | 邵11 2026-03-03 | 3/150 |
|
|
[考研] 材料學(xué)碩318求調(diào)劑 +15 | February_Feb 2026-03-01 | 17/850 |
|
|
[考研] 環(huán)境工程專碩307求調(diào)劑 +3 | ccc! 2026-03-03 | 3/150 |
|
|
[考研] 307求調(diào)劑 +6 | wyyyqx 2026-03-01 | 6/300 |
|
|
[考研] 調(diào)劑材料學(xué)碩 +4 | 詞凝Y 2026-03-02 | 4/200 |
|
|
[考研] 0854復(fù)試調(diào)劑 276 +5 | wmm9 2026-03-01 | 7/350 |
|
|
[考研] 275求調(diào)劑 +7 | 明遠(yuǎn)求學(xué) 2026-03-01 | 7/350 |
|
|
[考研] 289求調(diào)劑 +8 | yang婷 2026-03-02 | 9/450 |
|
|
[考研] 303求調(diào)劑 +5 | 今夏不夏 2026-03-01 | 5/250 |
|
|
[考研] 264求調(diào)劑 +4 | 巴拉巴拉根556 2026-02-28 | 4/200 |
|
|
[考研] 材料學(xué)調(diào)劑 +10 | 提神豆沙包 2026-02-28 | 12/600 |
|
|
[基金申請]
|
Doma 2026-03-01 | 7/350 |
|
|
[考研] 290求調(diào)劑 +9 | 材料專碩調(diào)劑; 2026-02-28 | 11/550 |
|
|
[考研] 295復(fù)試調(diào)劑 +3 | 簡木ChuFront 2026-03-01 | 3/150 |
|
|
[考研] 311求調(diào)劑 +9 | 南迦720 2026-02-28 | 10/500 |
|