c************************************************************************* c c USCAR-USAMP CRADA property model developed at Sandia. c The BAMMANN - CHIESA - JOHNSON Plasticity Model c was modified to include dendrite cell size effects and c damage nucleation, growth, and coalescence were c added by Mark F. Horstemeyer c c************************************************************************* c************************************************************************* c** ** c** WARNING WARNING WARNING WARNING WARNING WARNING WARNING ** c** ** c** this is proprietary software ** c** ** c** recipients of this routine are not allowed to distribute ** c** their source to other facilities. rather, requests for ** c** the source should be directed to the author at sandia. ** c** ** c** Mark F. Horstemeyer ** c** Sandia National Labs ** c** Mail Stop 9405 ** c** PO Box 969 ** c** Livermore, CA 94551-0969 ** c** ** c** email: mfhorst@sandia.gov ** c** ** c** WARNING WARNING WARNING WARNING WARNING WARNING WARNING ** c** ** c************************************************************************* c************************************************************************* C C Copyright C Sandia Corporation 1998 C All Rights Reserved C C Notice: The Government is granted for itself and others acting on C its behalf a paid-up, nonexclusive, irrevocable worldwide license in C this data to reproduce, prepare derivative works, and perform publicly C and display publicly. Beginning five (5) years after permission to C assert copyright was obtained, the Government is granted for itself C and others acting on its behalf a paid-up, nonexclusive, irrevocable C worldwide license in this data to reproduce, prepare derivative works, C distribute copies to the public, perform publicly and display publicly, C and to permit others to do so. NEITHER THE UNITED STATES NOR THE C UNITED STATES DEPARTMENT OF ENERGY, NOR ANY OF THEIR EMPLOYEES, MAKES C ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR C RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR USEFULNESS OF ANY C INFORMATION, APPARATUS, PRODUCT, OR PROCESS DISCLOSED, OR REPRESENTS C THAT ITS USE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS. C c************************************************************************* c c Sandia strain rate and temperature dependent plasticity model c was originally implemented by Mike Chiesa. c This subroutine can be used for axisymmetric, plane strain and c 3d solid elements. c Porosity was added by Mark Horstemeyer 1/8/97. c Void Nucleation was added by Mark Horstemeyer 9/28/97. c Void coalescence was added by Mark Horstemeyer 3/1/98. c Dendrite cell size was added by Mark Horstemeyer 3/15/98. c c************************************************************************* c c w a r n i n g c c this is proprietary software c c recipients of this routine are not allowed to distribute c their source to other facilities. rather, request for c the source should be directed to the author at snll. c phone: 510-294-2816 c c c*********************************************************************** c c l e g a l n o t i c e c c this computer code material was prepared as an account c of work sponsored by the united states government. c neither the united states nor the united states c department of energy, nor any of their employees, c nor any of their contractors, subcontractors, or their c employees, makes any warranty, express or implied, or c assumes any legal liability or responsibility for the c accuracy, completeness or usefulness of any information, c apparatus, produce or process disclosed, or represents c that its use would not infringe privately-owned rights c C************************************************************************* subroutine umat(stress,statev,ddsdde,sse,spd,scd, 1 rpl,ddsddt,drplde,drpldt, 2 stran,dstran,time,dtime,temp,dtemp,predef,dpred,cmname, 3 ndi,nshr,ntens,nstatv,props,nprops,coords,drot,pnewdt, 4 celent,dfgrd0,dfgrd1,noel,npt,layer,kspt,kstep,kinc) INCLUDE 'ABA_PARAM.INC' c *************************************************** character*8 cmname dimension stress(ntens),statev(nstatv), 1 ddsdde(ntens,ntens), 2 ddsddt(ntens),drplde(ntens), 3 stran(ntens),dstran(ntens),time(2),predef(1),dpred(1), 4 props(nprops),coords(3),drot(3,3),t(6),de(6), 5 dfgrd0(3,3),dfgrd1(3,3),xi(6) c parameter ( zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0, * third = one / three, half = 0.5d0, twothds = two / three, * op5 = 1.5d0 ) parameter ( four = 4.d0, twtsev = 27.d0, f427 = four / twtsev ) c data con1 /.81649658092773d0/ data pi /3.1415927d0/ data iparam1,iparam2,iparam3 / 0,0,0 / c c*********************************************************************** c c * ntens = number of non-zero stress components (4 for 2d, 6 for 3d) c * nshr = number of non-zero shear components (1 for 2d, 3 for 3d) c * ndi = number of non-zero normal stresses (always 3) c * nstatv = number of state variables (25) c * nprops = number of material parameters (50) c c************************************************************************* c c * statev(1) = alpha-xx c * statev(2) = alpha-yy c * statev(3) = alpha-zz c * statev(4) = alpha-xy c * statev(5) = alpha-yz c * statev(6) = alpha-zx c * statev(7) = kappa c * statev(8) = temperature c * statev(9) = effective plastic strain c * statev(10) = McClintock void growth(second phase pores) c * statev(11) = rate of change of M porosity c * statev(12) = stress triaxiality c * statev(13) = nucleation c * statev(14) = damage c * statev(15) = nucleation rate c * statev(16) = damage rate c * statev(17) = nucleation from previous time step c * statev(18) = Cocks-Ashby void growth(large pores) c * statev(19) = rate of change of CA porosity c * statev(20) = alpha-xx long range c * statev(21) = alpha-yy long range c * statev(22) = alpha-zz long range c * statev(23) = alpha-xy long range c * statev(24) = alpha-yz long range c * statev(25) = alpha-zx long range c c************************************************************************* c c * props(1) thru props(5) are constants for Johnson/Bammann c * formulas for shear and bulk moduli c * props(1) = mu zero , props(2) = a c * props(3) = K zero , props(4) = b c * props(5) = T melt , props(6) = C1 c * props(7) = C2 , props(8) = C3 c * props(9) = C4 , props(10)= C5 c * props(11)= C6 , props(12)= C7 c * props(13)= C8 , props(14)= C9 c * props(15)= C10 , props(16)= C11 c * props(17)= C12 , props(18)= C13 c * props(19)= C14 , props(20)= C15 c * props(21)= C16 , props(22)= C17 c * props(23)= C18 , props(24)= C19 c * props(25)= C20 , props(26)= CA c * props(27)= CB c * props(28)= initial temperature c * props(29)= heat generation coefficient c * props(30)= McClintock damage constant, n c * props(31)= initial void radius c * props(32)= torsional constant a in nucleation model c * props(33)= tension/comp constant b in nucleation model c * props(34)= triaxiality constant c in nucleation model c * props(35)= coefficient constant in nucleation model c * props(36)= fracture toughness, related to nucleation model c * props(37)= ave size of particles, related to nucleation model c * props(38)= particles vol fraction, related to nucleation model c * props(39)= coalescence factor, D=nucleation*void volume*coal. c * props(40)= coalescence factor, D=nucleation*void volume*coal. c * props(41)= reference grain size or dendrite cell size,dcs0 c * props(42)= grain size or dendrite cell size of material,dcs c * props(43)= grain size or dendrite cell size exponent,zz c * props(44)= initial void volume fraction for CA void growth c * props(45)= C21 , props(46)= C22 c * props(47)= C23 , props(48)= C24 c * props(49)= C25 , props(50)= C26 c * props(51)= nucleation temperature dependence c * props(52)= coalescence temperature dependence c * props(53)= flag to use vvfr4.dat file, 0=no, 1=yes c * props(54)= Cacon c * props(55)= elastic modulus - porosity exponent c****************************************************************** c c * iparam1 = 0 for linear return (uses total strain in recovery) c * = 1 for quadratic return (uses plastic strain in recovery) c * (only linear return implemented in this version) c c * iparam2 = 0 for Simo tangent stiffness matrix c * = 1 for Lathrop tangent stiffness matrix c * (only Simo stiffness implemented in this version) c c * iparam3 = 0 for trial kappa = kappa(n) c * = 1 for trial kappa = kappa(n + 1/2) c c************************************************************************ c---- initialize void volume fraction c this initization reads void volume fraction from the file c vvfr4.dat which is a text file in which each line contains c two entries: c element number, void volume fraction c the first line must be element 1, the next element 2 etc. c this file may be using the program mkvvf c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if (props(53).eq.one)then if(time(2).eq.zero)then write(6,'(a,i4)') "(initdmg.f) init umat, el = ",noel open(11,file='/cavs/compdesign/users/doglesby/vvfr4.dat', * status='OLD') do i=1,noel read(11,*)idum, vvf c vvf = zero end do close(11) props(44) = vvf end if else go to 13 end if c print *, "after initialization" c---- initialize state variables c temperature is set first by *initial condition command in abaqus c or by props(28) 13 if (props(54).eq.one)then if(time(2).eq.zero)then iseed=3 do i=1,noel iseed = 2045*iseed + 1 iseed = iseed -((iseed/1048576)*1048576) randx = real(iseed + 1)/1048577.0 vvf = (.02)*randx end do props(44) = vvf end if else go to 15 end if c 15 if(time(2).eq.zero)then statev(1) = zero statev(2) = zero statev(3) = zero statev(4) = zero statev(5) = zero statev(6) = zero statev(7) = zero statev(8) = props(28) statev(9) = zero statev(10) = (four/three)*pi*props(31)**three statev(10) = pi*props(31)**two statev(11) = zero statev(12) = zero statev(13) = props(35)*exp(-props(51)/statev(8)) statev(15) = zero statev(16) = zero statev(17) = zero if(props(47) .eq. one) then statev(18) = vvf else statev(18) = props(44) endif cf = (props(39)+props(40)*statev(13)*statev(10))* * (props(41)/props(42))**props(43)*exp(props(52)*statev(8)) statev(14) = cf*(statev(13)*statev(10)+statev(18)) statev(19) = zero statev(20) = zero statev(21) = zero statev(22) = zero statev(23) = zero statev(24) = zero statev(25) = zero if(temp.eq.zero)then if(props(28).eq.zero)then write(*,*)' error - temperature is zero' stop else statev(8) = props(28) endif else statev(8) = temp endif endif c c****************************************************************** c cc1 = props(6) cc2 = props(7) cc3 = props(8) cc4 = props(9) cc5 = props(10) cc6 = props(11) cc7 = props(12) cc8 = props(13) cc9 = props(14) cc10 = props(15) cc11 = props(16) cc12 = props(17) cc13 = props(18) cc14 = props(19) cc15 = props(20) cc16 = props(21) cc17 = props(22) cc18 = props(23) cc19 = props(24) cc20 = props(25) cc21 = props(45) cc22 = props(46) cc23 = props(47) cc24 = props(48) cc25 = props(49) cc26 = props(50) ca = props(26) cb = props(27) cd1 = props(39) cd2 = props(40) dcs0 = props(41) dcs = props(42) zz = props(43) zeta = props(55) c htcp = props(29) c c---- damage c dam1=one-statev(14) dam2 = one-min(one,dtime*statev(16)/dam1) phi1 = one-statev(18) c c---- g = shear modulus twog = 2*g blk = bulk modulus c if(props(28).eq.zero) statev(8) = temp theta = statev(8) if (props(5).eq.0) then tratio = zero g = props(1) else tratio = theta/props(5) tratio = min(tratio,0.9999d0) g = props(1)*(one-tratio*dexp(props(2)*(one-one/tratio))) endif twoga = two * g blka = props(3) - props(4)*tratio twog = (twoga * dam1**zeta) blk = (blka * dam1**zeta) c c---- calculate pressure c davg = third * (dstran(1) + dstran(2) + dstran(3)) pold = third * (stress(1) + stress(2) + stress(3)) p = pold*dam2 + dam1*blk * davg * three c c---- check for melt c if(theta.gt.props(5)) then do 35 i=1,ntens stress(i) = zero 35 statev(i) = zero p = min(zero,p) stress(1) = p stress(2) = p stress(3) = p statev(7) = zero statev(9) = zero go to 315 endif c c---- compute function evaluations c theta = temperature c ytheta = static yield stress c vtheta,ftheta = functions to define rate dependence of yield c h1,h2,h3 = plastic hardeing moduli c rs1,rs2,h3 = static recovery functions c rd1,rd2,h3 = dynamic recovery functions c c deviatoric stress ds11 = stress(1)-pold ds22 = stress(2)-pold ds33 = stress(3)-pold ds12 = stress(4) ds23 = stress(5) ds13 = stress(6) c invariants of stress dj2 = half*(ds11**two+ds22**two+ds33**two+ * two*(ds12**two+ds23**two+ds13**two)) dj3 = ds11*(ds22*ds33-ds23*ds23)-ds12*(ds12*ds33-ds23*ds13) * +ds13*(ds12*ds23-ds22*ds13) c c stress state dependent material constants adj = half*(one+tanh(cc19*(cc20-theta))) if(cc19.eq.zero)adj=one vtheta = cc1*dexp(-cc2/theta) ytheta = cc3*dexp(cc4/theta)*adj ftheta = cc5*dexp(-cc6/theta) if(dj2.eq.0)then rd1 = cc7*(one-ca*f427)*dexp(-cc8/theta) h1 = (cc9-cc10*theta)*(one+ca*f427) rs1 = cc11 * dexp(-cc12/theta) rd2 = cc13*(one-ca*f427)*dexp(-cc14/theta) h2 = (cc15-cc16*theta)*(one+ca*f427) rs2 = cc17 * dexp(-cc18/theta) rd3 = cc21*(one-ca*f427)*dexp(-cc22/theta) h3 = (cc23-cc24*theta)*(one+ca*f427) rs3 = cc25 * dexp(-cc26/theta) else rd1 = cc7*(one-ca*(f427-dj3**two/dj2**three) * -cb*dj3/dj2**op5)*dexp(-cc8/theta) h1 = (cc9-cc10*theta)*(one+ca*(f427-dj3**two/dj2**three) * +cb*dj3/dj2**op5) rs1 = cc11 * dexp(-cc12/theta) rd2 = cc13*(one-ca*(f427-dj3**two/dj2**three) * -cb*dj3/dj2**op5)*dexp(-cc14/theta) h2 = (cc15-cc16*theta)*(one+ca*(f427-dj3**two/dj2**three) * +cb*dj3/dj2**op5) rs2 = cc17 * dexp(-cc18/theta) rd3 = cc21*(one-ca*(f427-dj3**two/dj2**three) * -cb*dj3/dj2**op5)*dexp(-cc22/theta) h3 = (cc23-cc24*theta)*(one+ca*(f427-dj3**two/dj2**three) * +cb*dj3/dj2**op5) rs3 = cc25 * dexp(-cc26/theta) endif c c---- update alpha using abaqus rotation matrix c if(rs1.ne.zero .or. rd1.ne.zero .or. h1.ne.zero)then if(ntens.eq.4)then term1 = drot(1,1)*statev(1) + drot(1,2)*statev(4) term2 = drot(1,1)*statev(4) + drot(1,2)*statev(2) term3 = drot(2,1)*statev(1) + drot(2,2)*statev(4) term4 = drot(2,1)*statev(4) + drot(2,2)*statev(2) term5 = statev(1) + statev(2) statev(1) = drot(1,1)*term1+drot(1,2)*term2-term5*drot(1,3)**two statev(2) = drot(2,1)*term3+drot(2,2)*term4-term5*drot(2,3)**two statev(3) = - (statev(1) + statev(2)) statev(4) = drot(2,1)*term1+drot(2,2)*term2 - $ term5*drot(1,3)*drot(2,3) else term1=drot(1,1)*statev(1)+drot(1,2)*statev(4)+drot(1,3)*statev(6) term2=drot(1,1)*statev(4)+drot(1,2)*statev(2)+drot(1,3)*statev(5) term3=drot(1,1)*statev(6)+drot(1,2)*statev(5)+drot(1,3)*statev(3) term4=drot(2,1)*statev(1)+drot(2,2)*statev(4)+drot(2,3)*statev(6) term5=drot(2,1)*statev(4)+drot(2,2)*statev(2)+drot(2,3)*statev(5) term6=drot(2,1)*statev(6)+drot(2,2)*statev(5)+drot(2,3)*statev(3) term7=drot(3,1)*statev(1)+drot(3,2)*statev(4)+drot(3,3)*statev(6) term8=drot(3,1)*statev(4)+drot(3,2)*statev(2)+drot(3,3)*statev(5) term9=drot(3,1)*statev(6)+drot(3,2)*statev(5)+drot(3,3)*statev(3) statev(1)=term1*drot(1,1)+term2*drot(1,2)+term3*drot(1,3) statev(2)=term4*drot(2,1)+term5*drot(2,2)+term6*drot(2,3) statev(3)=term7*drot(3,1)+term8*drot(3,2)+term9*drot(3,3) statev(4)=term1*drot(2,1)+term2*drot(2,2)+term3*drot(2,3) statev(5)=term4*drot(3,1)+term5*drot(3,2)+term6*drot(3,3) statev(6)=term1*drot(3,1)+term2*drot(3,2)+term3*drot(3,3) endif endif c---- update long range alpha using abaqus rotation matrix c if(rs3.ne.zero .or. rd3.ne.zero .or. h3.ne.zero)then if(ntens.eq.4)then term1 = drot(1,1)*statev(20) + drot(1,2)*statev(23) term2 = drot(1,1)*statev(23) + drot(1,2)*statev(21) term3 = drot(2,1)*statev(20) + drot(2,2)*statev(23) term4 = drot(2,1)*statev(23) + drot(2,2)*statev(21) term5 = statev(20) + statev(21) statev(20) = drot(1,1)*term1+drot(1,2)*term2-term5*drot(1,3)**two statev(21) = drot(2,1)*term3+drot(2,2)*term4-term5*drot(2,3)**two statev(22) = - (statev(20) + statev(21)) statev(23) = drot(2,1)*term1+drot(2,2)*term2 - $ term5*drot(1,3)*drot(2,3) else term1=drot(1,1)*statev(20)+drot(1,2)*statev(23)+drot(1,3)* * statev(25) term2=drot(1,1)*statev(23)+drot(1,2)*statev(21)+drot(1,3)* * statev(24) term3=drot(1,1)*statev(25)+drot(1,2)*statev(24)+drot(1,3)* * statev(22) term4=drot(2,1)*statev(20)+drot(2,2)*statev(23)+drot(2,3)* * statev(25) term5=drot(2,1)*statev(23)+drot(2,2)*statev(21)+drot(2,3)* * statev(24) term6=drot(2,1)*statev(25)+drot(2,2)*statev(24)+drot(2,3)* * statev(22) term7=drot(3,1)*statev(20)+drot(3,2)*statev(23)+drot(3,3)* * statev(25) term8=drot(3,1)*statev(23)+drot(3,2)*statev(21)+drot(3,3)* * statev(24) term9=drot(3,1)*statev(25)+drot(3,2)*statev(24)+drot(3,3)* * statev(22) statev(20)=term1*drot(1,1)+term2*drot(1,2)+term3*drot(1,3) statev(21)=term4*drot(2,1)+term5*drot(2,2)+term6*drot(2,3) statev(22)=term7*drot(3,1)+term8*drot(3,2)+term9*drot(3,3) statev(23)=term1*drot(2,1)+term2*drot(2,2)+term3*drot(2,3) statev(24)=term4*drot(3,1)+term5*drot(3,2)+term6*drot(3,3) statev(25)=term1*drot(3,1)+term2*drot(3,2)+term3*drot(3,3) endif endif c c---- compute effective strain rate c if(dtime.ne.zero)then dum = zero do 40 i=4,ntens 40 dum = dum + dstran(i)**two ddd = dsqrt(dstran(1)**two+dstran(2)**two+dstran(3)**two + $ half*dum) * con1 / dtime else ddd = zero endif c c---- calculate trial alpha, kappa and yield radius c alphaxx = statev(1)+statev(20) alphayy = statev(2)+statev(21) alphazz = statev(3)+statev(22) alphaxy = statev(4)+statev(23) alphayz = statev(5)+statev(24) alphazx = statev(6)+statev(25) alpm = con1 * dsqrt(alphaxx**two + alphayy**two + alphazz**two+ $ two*(alphaxy**two + alphayz**two + alphazx**two)) sto = dtime*rs1*alpm*(dcs0/dcs)**zz sto2 = dtime*rs2*(dcs0/dcs)**zz if(iparam1.eq.0)then sto = dtime * (rs1+rd1*ddd+rs3+rd3*ddd)*alpm*(dcs0/dcs)**zz sto2 = dtime * (rs2+rd2*ddd)*(dcs0/dcs)**zz endif do 50 i=1,ntens 50 statev(i) = statev(i) * (one - sto) do 51 i=20,ntens+20-1 51 statev(i) = statev(i) * (one - sto) if(iparam3.eq.0)then trialk = statev(7) else trialk = (-one+dsqrt(one+two*sto2*(statev(7) $ +half*h2*ddd*dtime)))/ $ max(1.d-30,sto2) endif statev(7) = statev(7) - sto2 * trialk * trialk ak = (vtheta * log((ddd+dsqrt(ddd**two+ftheta**two))/ftheta) + $ ytheta + statev(7))*dam1 c c---- calculate trial elastic deviatoric stresses c do 60 i=1,3 60 stress(i) = dam2*(stress(i)-pold)+dam1*twog*(dstran(i)-davg) do 70 i=4,ntens 70 stress(i) = dam2*stress(i) + dam1*g*dstran(i) c c---- compute xi (deviatoric stress - 2/3 alpha) c do 80 i=1,ntens 80 xi(i) = stress(i) - twothds * statev(i) c c---- compute (magnitude of xi) squared c dum = zero do 85 i=4,ntens 85 dum = dum + xi(i)**two ximag2 = xi(1)**two + xi(2)**two + xi(3)**two + two*dum c c---- check for plasticity c ak2 = ximag2 - twothds * ak * abs(ak) if(ak2 .le. zero .or. ddd.eq.zero) go to 300 c c---- plasticity process begins here c ximag = dsqrt(ximag2) c c---- return trial stresses to yield surface, add pressure term c and update state variables c if(iparam1.eq.0)then dgam = (ximag-con1*ak)/(dam1*twog+twothds*(dcs0/dcs)**zz * *(h1+h3+h2*dam1)) endif dgam2 = dgam / ximag c dsig = dam1*twog * dgam2 do 90 i=1,ntens 90 stress(i) = stress(i) - dsig * xi(i) stress(1) = stress(1) + p stress(2) = stress(2) + p stress(3) = stress(3) + p c statev(7) = statev(7) + dgam * con1 * h2*(dcs0/dcs)**zz statev(7) = max(zero,statev(7)) c dalph = (h1+h3) * dgam2 do 100 i=1,ntens 100 statev(i) = statev(i) + dalph * xi(i)*(dcs0/dcs)**zz c c---- update plastic strain (for output purposes only) c statev(9) = statev(9) + dgam * con1 c c---- update temperature for adiabatic problems c dum = zero do 110 i=4,ntens 110 dum = dum + stress(i)*xi(i) statev(8) = statev(8) + htcp*dgam2*( stress(1)*xi(1) + $ stress(2)*xi(2) + stress(3)*xi(3) + two*dum) c c---- update damage c epsdot=dgam*con1/dtime sige=max(1.d-15,dsqrt(half*((stress(1)-stress(2))**two+ * (stress(2)-stress(3))**two+(stress(3)-stress(1))**two+ * 6.d0*(stress(4)**two+stress(5)**two+stress(6)**two)))) c c Cocks-Ashby large pore growth term cacon = abs(vtheta/ytheta) if(cacon.lt.18.0d0) then cacon=18.0d0 endif dterm=two*(two*cacon-one)/(two*cacon+one) arg = min(15.d0,p*dterm/sige) beta = dsinh(max(zero,arg) ) c90 = one + cacon psi = min(15.d0,beta*dtime*epsdot*c90) tmp = max(zero,(one+(phi1**c90-one)*dexp(psi))) statev(18) = min((one-tmp**(one/c90)),.99d0) c Cocks-Ashby void growth rate statev(19) = beta*epsdot*(one/(one-statev(18)) * **(vtheta/ytheta) * -(one-statev(18))) c McClintock form of void growth sqtth = three**half abc=sqtth/(two*(one-props(30)))*dsinh(sqtth* * half*(one-props(30))*(2*p/sige+third)) vrad=props(31)*dexp(statev(9)*abc/con1) statev(10)=pi*vrad**two statev(11)=three*statev(10)*abc*epsdot c Nucleation of voids c deviatoric stress ds11=stress(1)-p ds22=stress(2)-p ds33=stress(3)-p ds12=stress(4) ds23=stress(5) ds13=stress(6) c invariants of stress di1=three*p dj2=half*(ds11**two+ds22**two+ds33**two+ * 2*(ds12**two+ds23**two+ds13**two)) dj3=ds11*(ds22*ds33-ds23**two)-ds22*(ds11*ds33-ds13**two) * +ds33*(ds22*ds11-ds12**two) if(dj2.le.zero)then r1=zero r2=zero r3=zero else r1=(f427-dj3**two/dj2**three) r2=dj3/dj2**op5 r3=di1/(dj2**half) endif r3=abs(r3) c if(r3.lt.0)then c r3=-r3 c endif zzz=(props(32)*r1+props(33)*r2+props(34)*r3) zzz=abs(zzz) zzzz=(props(37)**half/(props(36)*props(38)**third))*zzz statev(17)=statev(13) statev(13)=props(35)*dexp(statev(9)*zzzz/con1)* * dexp(-props(51)/statev(8)) c added for nonmonotonic path sequences, statev(17) is old nucleation if(statev(13).lt.statev(17)) then statev(13) = abs(statev(17)+statev(13)) endif c Coalescence factor cf=(cd1+cd2*statev(13)*statev(10))* * dexp(props(52)*statev(8)) * *(dcs0/dcs)**zz c Damage damage=cf*(statev(13)*statev(10)+statev(18)) if(damage .gt. 0.6d0) damage = .99d0 statev(14)=min(damage,0.99d0) c Nucleation Rate epsdot=abs(epsdot) statev(15)=zzzz*statev(13)*epsdot c Damage Rate zsecond=cf*(statev(15)*statev(10)+statev(13)*statev(11)+ * statev(19)) zthird=(statev(13)*statev(10)+statev(18))* * cd2*(dcs0/dcs)**zz*dexp(props(52)*statev(8))* * (statev(15)*statev(10)+statev(13)*statev(11)) statev(16)=zsecond+zthird c c Triaxiality statev(12)=p/sige c c---- form elastic-plastic constitutive matrix c if bbb=1 then conventional stiffness is calculated c otherwise simo stiffness is calculated c if(iparam2.eq.0)then alpm2 = con1 * dsqrt(alphaxx**two + alphayy**two + alphazz**two $ + two*(alphaxy**two+alphayz**two+alphazx**two)) bbb = (con1*(ak+h2*con1*dgam) + alpm2 - alpm) / ximag bbb = min(one,bbb) slope = h1+h3+h2-(rd1+rd3)*alpm2*alpm2-rd2*statev(7)*statev(7) slope = max((1.d-5*(h1+h2+h3)),slope) gamma = one/(one+slope/three/g) - (one-bbb) r = dam1*(twog * gamma / ximag**two) z1 = dam1*(blk + twothds * twog * bbb) z2 = dam1*(blk - third * twog * bbb) z3 = dam1*(g * bbb) do 140 i=1,ntens do 140 j=1,ntens 140 ddsdde(i,j) = - r * xi(i) * xi(j) ddsdde(1,1) = z1 + ddsdde(1,1) ddsdde(2,2) = z1 + ddsdde(2,2) ddsdde(3,3) = z1 + ddsdde(3,3) ddsdde(1,2) = z2 + ddsdde(1,2) ddsdde(1,3) = z2 + ddsdde(1,3) ddsdde(2,1) = z2 + ddsdde(2,1) ddsdde(3,1) = z2 + ddsdde(3,1) ddsdde(2,3) = z2 + ddsdde(2,3) ddsdde(3,2) = z2 + ddsdde(3,2) do 130 i=4,ntens 130 ddsdde(i,i) = z3 + ddsdde(i,i) endif 2000 format(' debug',5e13.5/4e13.5) go to 400 c c---- elastic process begins here c 300 continue do 310 i=1,3 310 stress(i) = stress(i) + p dgam=zero c c---- form elastic stiffness matrix c 315 continue do 320 i=1,ntens do 320 j=1,ntens 320 ddsdde(i,j) = zero z1 = dam1*(blk + twothds * twog) z2 = dam1*(blk - third * twog) ddsdde(1,1) = z1 ddsdde(2,2) = z1 ddsdde(3,3) = z1 ddsdde(1,2) = z2 ddsdde(1,3) = z2 ddsdde(2,1) = z2 ddsdde(3,1) = z2 ddsdde(2,3) = z2 ddsdde(3,2) = z2 do 330 i=4,ntens 330 ddsdde(i,i) = dam1*g c c---- clean up c 400 continue statev(7) = max(statev(7),zero) c 101 format(1x,6(d12.6,2x)) 102 format(1x,a24,6(d12.6,2x)) return end C SUBROUTINE DISP(U,KSTEP,KINC,TIME,NODE,NOEL,JDOF,COORDS) C-----------------------------------------------------------------* C-----------------------------------------------------------------* C----- CALCULATION OF IMPOSED DISPLACEMENT AND VELOCITY --------* C-------------- TO OBTAIN CONSTANT STRAIN RATE -----------------* C-----------------------------------------------------------------* C-----------------------------------------------------------------* C ----------- C Input : C ----------- C KSTEP : Step number C KINC : Increment number C TIME(1) : Step time C TIME(2) : Total time C NODE : Node number C JDOF : Degree of freedom C ----------- C Output : C ----------- C U(1) : Total value of the prescribed value at this point. C The variable may be displacement, rotation, pore C pressure, temperature, etc., depending on the degree C of freedom constrained. C U(2) : dU(1)/dTIME C U(3) : d2U(1)/(dTIME)2 C =============================================================== C INCLUDE 'ABA_PARAM.INC' C DIMENSION U(3),TIME(2),COORDS(3) DATA ZERO,ONE,TWO/0.D0,1.D0,2.D0/ C----------------------------------------------------------------- C----------------------------------------------------------------- C--------- Computation of displacement C----------------------------------------------------------------- C XL0 = ONE RATE = 1.D-01 U(1) = XL0*DEXP(RATE*TIME(2))-XL0 U(2) = RATE*XL0*DEXP(RATE*TIME(2)) U(3) = RATE*U(2) C C====================================================================== RETURN END C