subroutine vumat( c Read only - 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' 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), charLength(nblock), 1 coordMp(nblock,*), 2 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) * 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 = 17, * i_pro_sigu1c = 18, * i_pro_sigu2t = 19, * i_pro_sigu2c = 20, * i_pro_sigu3t = 21, * i_pro_sigu3c = 22, * i_pro_sigu12 = 25, * i_pro_sigu13 = 26, * i_pro_sigu23 = 27 ) * Temporary arrays dimension eigen(maxblk*3) parameter ( Gft = 133.d0, Gfc = 40.d0 ) parameter ( Gmt = 0.6d0, Gmc = 2.1d0 ) * * 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) * 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 = E1 * ( xnu21 + xnu31*xnu23 ) * gg C13 = E1 * ( xnu31 + xnu21*xnu32 ) * gg C23 = E2 * ( xnu32 + xnu12*xnu31 ) * gg * 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, G12, G23, G13, * 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, G12, G23, G13, * stateNew(1,i_svd_strain), * stressNew ) * * Failure evaluation * call copyr ( nblock, * stateOld(1,i_svd_DmgFiberT), stateNew(1,i_svd_DmgFiberT) ) call copyr ( nblock, * stateOld(1,i_svd_DmgFiberC), stateNew(1,i_svd_DmgFiberC) ) 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) ) nDmg = 0 call eig33Anal ( nblock, stretchNew, eigen ) call Hashin3d ( nblock, nDmg, charLength, * f1t, f2t, f3t, f1c, f2c, f3c, f12, f23, f13, * E1, E2, G12, G13, G23, C11, C22, * stateNew(1,i_svd_DmgFiberT), * stateNew(1,i_svd_DmgFiberC), * stateNew(1,i_svd_DmgMatrixT), * stateNew(1,i_svd_DmgMatrixC), * stateNew(1,i_svd_statusMp), * stateNew(1,i_svd_strain),stressNew,eigen ) * -- 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, G12, G23, G13, * 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 * ************************************************************ subroutine OrthoEla3dExp ( nblock, * dmgFiberT, dmgFiberC, dmgMatrixT, dmgMatrixC, * C11, C22, C33, C12, C23, C13, G12, G23, G13, * strain, stress ) * include 'vaba_param.inc' * Orthotropic elasticity, 3D case - * parameter( zero = 0.d0, one = 1.d0, two = 2.d0, half = 0.5d0) parameter ( smt = 0.9d0, smc = 0.7d0 ) * 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), * dmgFiberT(nblock), dmgFiberC(nblock), * dmgMatrixT(nblock), dmgMatrixC(nblock), * stress(nblock,n_s33_Car) * -- shear fraction in matrix tension and compression mode * do k = 1, nblock * -- Compute damaged stiffness 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 ) * G12 dG23 = ( one - df ) * * ( one - smt*dmt ) * ( one - smc*dmc ) * G23 dG13 = ( one - df ) * * ( one - smt*dmt ) * ( one - smc*dmc ) * G13 * -- 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) 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 * * 缺少除了失效区域的hashin系数 ************************************************************ subroutine Hashin3d ( nblock, nDmg, charLength, * f1t, f2t, f3t, f1c, f2c, f3c, f12, f23, f13, * E1, E2, G12, G13, G23, C11, C22, * dmgFiberT, dmgFiberC, dmgMatrixT, dmgMatrixC, * statusMp, strain, stress, eigen ) * include 'vaba_param.inc' parameter( * zero = 0.d0, one = 1.d0, * half = 0.5d0, two= 2.d0, three =3.d0, alpha = 0.9d0) * parameter( cri = .99d0, thickness = 0.52d0 ) parameter( smt = 0.d0, smc = 0.d0, sft = 0.d0, sfc = 0.d0) parameter ( Gft = 133.d0, Gfc = 40.d0 ) parameter ( Gmt = 0.6d0, Gmc = 2.1d0 ) * 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 dmgFiberT(nblock), dmgFiberC(nblock), * dmgMatrixT(nblock), dmgMatrixC(nblock), * strain(nblock,n_s33_Car), * stress(nblock,n_s33_Car), * eigen(nblock,n_v3d_Car), * statusMp(nblock), charLength(nblock) * Real :: e_11P, e_11N, e_22P, e_22N Real :: e_eq1t, e_eq1c, e_eq2t, e_eq2c Real :: e_ft0, e_fft, e_fc0, e_ffc, e_mt0, e_fmt, e_mc0, e_fmc * f1tInv = zero f2tInv = zero f3tInv = zero f1cInv = zero f2cInv = zero f3cInv = zero f12Inv = zero f23Inv = zero f13Inv = zero * if ( f1t .gt. zero ) f1tInv = one / f1t if ( f2t .gt. zero ) f2tInv = one / f2t if ( f3t .gt. zero ) f3tInv = one / f3t if ( f1c .gt. zero ) f1cInv = one / f1c if ( f2c .gt. zero ) f2cInv = one / f2c if ( f3c .gt. zero ) f3cInv = one / f3c if ( f12 .gt. zero ) f12Inv = one / f12 if ( f23 .gt. zero ) f23Inv = one / f23 if ( f13 .gt. zero ) f13Inv = one / f13 * e1t0 = f1t / E1 e1c0 = f1c / E1 e2t0 = f2t / E2 e2c0 = f2c / E2 * do k = 1, nblock if ( statusMp(k) .eq. one ) then * lDmg = 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) * s11 = stress(k,i_s33_Xx) s22 = stress(k,i_s33_Yy) s33 = stress(k,i_s33_Zz) s12 = stress(k,i_s33_Xy) s23 = stress(k,i_s33_Yz) s13 = stress(k,i_s33_Zx) * e_11P = (abs(strain(k,i_s33_Xx)) + strain(k,i_s33_Xx))/two e_11N = (abs(strain(k,i_s33_Xx)) - strain(k,i_s33_Xx))/two e_22P = (abs(strain(k,i_s33_Yy)) + strain(k,i_s33_Yy))/two e_22N = (abs(strain(k,i_s33_Yy)) - strain(k,i_s33_Yy))/two * e_eq1t = ((e_11P)**two + (e12)**two + (e13)**two)**half e_eq1c = e_11N e_eq2t = ((e_22P)**two + (e12)**two)**half e_eq2c = ((e_22N)**two + (e12)**two)**half * charLength(k) = (charLength(k)**three / thickness)**half e_ft0 = f1t / C11 e_fft = two*Gft/charLength(k)/f1t e_fc0 = f1c / C11 e_ffc = two*Gfc/charLength(k)/f1c e_mt0 = f2t / C22 e_fmt = two*Gmt/charLength(k)/f2t e_mc0 = f2c / C22 e_fmc = two*Gmc/charLength(k)/f2c * * Evaluate Fiber modes if ( s11 .gt. zero ) then * -- Tensile Fiber Mode rft = (e11/e1t0)**two * + alpha * (e12 * f12Inv * G12)**two * + alpha * (e13 * f12Inv * G12)**two if ( rft .ge. one ) then lDmg = 1 dmgFiberT(k) = e_fft * (e_eq1t-e_ft0) / * e_eq1t / (e_fft - e_ft0) if (dmgFiberT(k) .lt.zero) then dmgFiberT(k) = zero end if if( dmgFiberT(k) .gt.cri )then dmgFiberT(k) = cri end if end if else if ( s11 .lt. zero ) then * -- Compressive Fiber Mode rfc = (e11 /e1c0)**two if ( rfc .ge. one ) then lDmg = 1 dmgFiberC(k) = e_ffc * (e_eq1c - e_fc0)/ * e_eq1c / (e_ffc - e_fc0) if (dmgFiberC(k).lt.zero) then dmgFiberC(k) = zero end if if( dmgFiberC(k).gt.cri )then dmgFiberC(k) = cri end if end if end if * * Evaluate Matrix Modes if ( ( s22 + s33 ) .gt. zero ) then * -- Tensile Matrix mode rmt = ((e22 * E2 + e33 * E2)*(f2tInv))**two * + (e12 * G12 * f12Inv)**two * + (e13 * G13 * f13Inv)**two * - (e22 * E2 * e33 * E2)*(f23Inv)**two * + (e23 * G23 * f23Inv)**two if ( rmt .ge. one ) then lDmg = 1 dmgMatrixT(k) = e_fmt * (e_eq2t - e_mt0)/ * e_eq2t/ (e_fmt - e_mt0) if (dmgMatrixT(k).lt.zero) then dmgMatrixT(k) = zero end if if(dmgMatrixT(k).gt.cri)then dmgMatrixT(k) = cri end if end if else if ( ( s22 + s33 ) .lt. zero ) then * -- Compressive Matrix Mode rmc = (e22*E2+e33*E2)*((E2 * e2c0 * half * f23Inv)**two * - one)*f2cInv * + ((e22*E2+e33*E2)*half*f23Inv)**two * + (e23 * G23 * f23Inv)**two * - (e22 * E2 * e33 * E2)*(f23Inv)**two * + (e12 * G12 * f12Inv)**two * + (e13 * G13 * f13Inv)**two if ( rmc .ge. one ) then lDmg = 1 dmgMatrixC(k) = e_fmc * (e_eq2c - e_mc0)/ * e_eq2c / (e_fmc - e_mc0) if (dmgMatrixC(k).lt.zero) then dmgMatrixC(k) = zero end if if( dmgMatrixC(k).gt.cri )then dmgMatrixC(k) = cri end if end if end if * * eigMax=max(eigen(k,i_v3d_X),eigen(k,i_v3d_Y),eigen(k,i_v3d_Z)) eigMin=min(eigen(k,i_v3d_X),eigen(k,i_v3d_Y),eigen(k,i_v3d_Z)) enomMax = eigMax - one enomMin = eigMin - one * * if ( enomMax .gt. eMax .or. * enomMin .lt. eMin .or. dmgFiberT(k) .eq. one) then statusMp(k) = zero end if nDmg = nDmg + lDmg * 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' * parameter(one = 1.d0) dimension from(nCopy), to(nCopy) * do k = 1, nCopy to(k) = from(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