c************************************************************************* c************************************************************************* c** ** c** WARNING WARNING WARNING WARNING WARNING WARNING WARNING ** c** ** c** This is a 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 MSU. ** c** ** c** Mark F. Horstemeyer ** c** Chair for Advanced Manufacture and Design ** c** Center for Advanced Vehicular System ** c** Mississippi State University ** c** ** c** Phone: (662) 325-7308 ** c** Fax: (662) 325-7223 ** c** email: mfhorst@me.msstate.edu ** c** ** c** WARNING WARNING WARNING WARNING WARNING WARNING WARNING ** c** ** 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************************************************************************* C************************************************************************* c c References: c c 1. D. J. Bammann, 1990, Modeling temperature and strain rate dependent c large deformation of metals., Appl. Mech. Rev., 43, S312-19. c 2 D. J. Bammann, M. L. Chiesa, M. F. Horstemeyer, L. I. Weingarten, 1993, c Failure in ductile materials using finite element methods, In: c Structural Crashworthiness and Failure, N, Jones and T. Wierzbicki (eds.), c Elsevier Applied Science, 1-54. c 2. M.F. Horstemeyer, M.M. Matalanis, A.M. Sieber, M.L. Botos, 2000, c Micromechanical finite element calculations of temperature and void c configuration effects on void growth and coalescence, Int. J. c of Plasticity, 16, 979-1015 c 3. M.F. Horstemeyer, J. Lathrop, A. M. Gokhale, M. Dighe, 2000, c Modeling stress state dependent damage evolution in a cast Al-Si c -Mg aluminum alloy, Theoretical and Appl. Fracture Mechanics, 33, 31-47 c 3. M. F. Horstemeyer, A. M. Gokhale, 1999, A void-crack nucleation model c for ductile metals, Int. J. of Solids adn Structures, 36, 5029-5055 c 4. D. J. Bammann, M. L. Chiesa, M. F. Horstemeyer, L. I. Weingarten, 1993, c Failure in ductile materials using finite element methods, In: c Structural Crashworthiness and Failure, N, Jones and T. Wierzbicki (eds.), c Elsevier Applied Science, 1-54. c c C************************************************************************* C************************************************************************* C User subroutine VUMAT subroutine vumat ( C Read only - * nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, * stepTime, totalTime, dt, cmname, coordMp, charLength, * props, density, strainInc, relSpinInc, * tempOld, stretchOld, defgradOld, fieldOld, * stressOld, stateOld, enerInternOld, enerInelasOld, * tempNew, stretchNew, defgradNew, fieldNew, C Write only - * stressNew, stateNew, enerInternNew, enerInelasNew ) C include 'vaba_param.inc' C dimension coordMp(nblock,*), charLength(nblock), props(nprops), 1 density(nblock), strainInc(nblock,ndir+nshr), 2 relSpinInc(nblock,nshr), tempOld(nblock), 3 stretchOld(nblock,ndir+nshr), 4 defgradOld(nblock,ndir+nshr+nshr), 5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr), 6 stateOld(nblock,nstatev), enerInternOld(nblock), 7 enerInelasOld(nblock), tempNew(nblock), 8 stretchNew(nblock,ndir+nshr), 9 defgradNew(nblock,ndir+nshr+nshr), 1 fieldNew(nblock,nfieldv), 2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev), 3 enerInternNew(nblock), enerInelasNew(nblock) dimension xi(6) C character*80 cmname 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 ) data con1 /.81649658092773/ data pi /3.1415927/ 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 c c****************************************************************** c ntens = ndir+nshr c c------- BCJ parameters c c set material constants 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) c htcp = props(29) c cc21 = props(45) cc22 = props(46) cc23 = props(47) cc24 = props(48) cc25 = props(49) cc26 = props(50) c 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 c-------------------------------------------------------------------- c If stepTime equals to zero, assume the material pure elastic c and use initial elastic modulus c if ( TotalTime .eq. zero ) then do k = 1, nblock c c--------- initialize state variables if (props(54).eq.one)then randx = rand(0)/(1-2*rand(0)) vvf = props(44)*randx endif stateOld(k,8) = props(28) stateOld(k,10) = pi*props(31)**two stateOld(k,13) = props(35)*exp(-props(51)/stateOld(k,8)) stateOld(k,18) = props(44) if (props(47) .eq. one) then stateOld(k,18) = vvf else stateOld(k,18) = props(44) endif cf = (cd1+cd2*stateOld(k,13)*stateOld(k,10)) * * exp(props(52)*stateOld(k,8))*(dcs0/dcs)**zz stateOld(k,14) = cf*(stateOld(k,13)*stateOld(k,10) * +stateOld(k,18)) c--------- temperature is set first by *initial condition command in abaqus c or by props(28) stateOld(k,8) = props(28) if(tempOld(k) .eq. zero) then if(props(28) .eq. zero) then write(*,*)' error - temperature is zero' stop else stateOld(k,8) = props(28) endif else stateOld(k,8) = tempOld(k) endif if(props(28).eq.zero) stateOld(k,8) = tempOld(k) end do end if if ( StepTime .eq. zero ) then do k = 1, nblock c--------- g = shear modulus twog = 2*g blk = bulk modulus g = props(1) blk = props(3) c---------- Trial stress blkg = blk-twothds*g trace = strainInc(k,1) + strainInc(k,2) + strainInc(k,3) blkgt = blkg*trace do i=1,ntens stressNew(k,i) = stressOld(k,i) + twog * strainInc(k,i) end do stressNew(k,1) = stressNew(k,1) * + twog * strainInc(k,1) + blkgt stressNew(k,2) = stressNew(k,2) * + twog * strainInc(k,2) + blkgt stressNew(k,3) = stressNew(k,3) * + twog * strainInc(k,3) + blkgt end do return end if c-------------------------------------------------------------------- c do 100 k = 1, nblock c c--------- temperature if(props(28).eq.zero) stateOld(k,8) = tempOld(k) theta = stateOld(k,8) c--------- state variables do i=1,nstatev stateNew(k,i) = stateOld(k,i) end do c-------- damage c stateOld(k,14) = zero dam1 = one-stateOld(k,14) dam2 = one-min(one,dt*stateOld(k,16)/dam1) phi1 = one-stateOld(k,18) c c--------- g = shear modulus twog = 2*g blk = bulk modulus if (props(5).eq.zero) then tratio = zero g = props(1) else tratio = theta/props(5) tratio = min(tratio,0.9999) g = props(1)*(one-tratio*exp(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 davg = third*(strainInc(k,1) + strainInc(k,2) + strainInc(k,3)) pold = third*(stressOld(k,1) + stressOld(k,2) + stressOld(k,3)) p = pold*dam2 + dam1*blk * davg * three c c---- check for melt c if(theta.gt.props(5)) then do i=1,ntens stressNew(k,i) = zero stateNew(k,i) = zero end do p = min(zero,p) stressNew(k,1) = p stressNew(k,2) = p stressNew(k,3) = p stateNew(k,7) = zero stateNew(k,9) = zero go to 200 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 = plastic hardeing moduli c rs1,rs2 = static recovery functions c rd1,rd2 = dynamic recovery functions c c deviatoric stress ds11 = stressOld(k,1)-pold ds22 = stressOld(k,2)-pold ds33 = stressOld(k,3)-pold ds12 = stressOld(k,4) ds23 = stressOld(k,5) ds13 = stressOld(k,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*exp(-cc2/theta) ytheta = cc3*exp(cc4/theta)*adj ftheta = cc5*exp(-cc6/theta) if(dj2.eq.0)then rd1=cc7*(1-ca*(4./27.))*exp(-cc8/theta) h1=(cc9-cc10*theta)*(1+ca*(4./27.)) rs1 = cc11 * exp(-cc12/theta) rd2 = cc13*(1-ca*(4./27.))*exp(-cc14/theta) h2 = (cc15-cc16*theta)*(1+ca*(4./27.)) rs2 = cc17 * exp(-cc18/theta) rd3 = cc21*(1-ca*(4./27.))*exp(-cc22/theta) h3 = (cc23-cc24*theta)*(1+ca*(4./27.)) rs3 = cc25 * exp(-cc26/theta) else rd1 = cc7*(1-ca*(4./27.-dj3**2./dj2**3.) * -cb*dj3/dj2**1.5)*exp(-cc8/theta) h1 = (cc9-cc10*theta)*(1+ca*(4./27.-dj3**2./dj2**3.) * +cb*dj3/dj2**1.5) rs1 = cc11 * exp(-cc12/theta) rd2 = cc13*(1-ca*(4./27.-dj3**2./dj2**3.) * -cb*dj3/dj2**1.5)*exp(-cc14/theta) h2 = (cc15-cc16*theta)*(1+ca*(4./27.-dj3**2./dj2**3.) * +cb*dj3/dj2**1.5) rs2 = cc17 * exp(-cc18/theta) rd3 = cc21*(1-ca*(4./27.-dj3**2./dj2**3.) * -cb*dj3/dj2**1.5)*exp(-cc22/theta) h3 = (cc23-cc24*theta)*(1+ca*(4./27.-dj3**2./dj2**3.) * +cb*dj3/dj2**1.5) rs3 = cc25 * exp(-cc26/theta) endif c c---- compute effective strain rate c if(dt .ne. zero)then dum = zero do i=4,ntens dum = dum + strainInc(k,i)**two end do ddd = sqrt(strainInc(k,1)**two+strainInc(k,2)**two * + strainInc(k,3)**two + half*dum) * con1 / dt else ddd = zero endif c c---- calculate trial alpha, kappa and yield radius c alphaxx = stateOld(k,1)+stateOld(k,20) alphayy = stateOld(k,2)+stateOld(k,21) alphazz = stateOld(k,3)+stateOld(k,22) alphaxy = stateOld(k,4)+stateOld(k,23) alphayz = stateOld(k,5)+stateOld(k,24) alphazx = stateOld(k,6)+stateOld(k,25) alpm = con1 * sqrt(alphaxx**two + alphayy**two + alphazz**two * + two*(alphaxy**two + alphayz**two + alphazx**two)) sto = dt*rs1*alpm*(dcs0/dcs)**zz sto2 = dt*rs2*(dcs0/dcs)**zz if(iparam1.eq.0)then sto = dt * (rs1+rd1*ddd+rs3+rd3*ddd)*alpm*(dcs0/dcs)**zz sto2 = dt * (rs2+rd2*ddd)*(dcs0/dcs)**zz endif do i=1,ntens stateNew(k,i) = stateOld(k,i) * (one - sto) stateNew(k,i+19) = stateOld(k,i+19) * (one - sto) end do if(iparam3.eq.0)then trialk = stateOld(k,7) else * trialk = (-one+sqrt(one+two*sto2*(stateOld(k,7) * + half*h2*ddd*dt)))/ max(1.e-30,sto2) endif stateNew(k,7) = stateOld(k,7) - sto2 * trialk * trialk ak = (vtheta * log((ddd+sqrt(ddd**two+ftheta**two))/ftheta) * + ytheta + stateNew(k,7))*dam1 c c---- calculate trial elastic deviatoric stresses c do i=1,3 stressNew(k,i) = dam2*(stressOld(k,i)-pold) * + dam1*twog*(strainInc(k,i)-davg) end do do i=4,ntens stressNew(k,i) = dam2*stressOld(k,i) + dam1*g*strainInc(k,i) end do c c---- compute xi (deviatoric stress - 2/3 alpha) c do i=1,ntens xi(i) = stressNew(k,i) - twothds * stateNew(k,i) end do c c---- compute (magnitude of xi) squared c dum = zero do i=4,ntens dum = dum + xi(i)**two end do 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) then c c---- elastic process begins here do i=1,3 stressNew(k,i) = stressNew(k,i) + p end do go to 200 end if c c---- plasticity process begins here c ximag = sqrt(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+h2*dam1)) endif dgam2 = dgam / ximag c dsig = dam1*twog * dgam2 do i=1,ntens stressNew(k,i) = stressNew(k,i) - dsig * xi(i) end do stressNew(k,1) = stressNew(k,1) + p stressNew(k,2) = stressNew(k,2) + p stressNew(k,3) = stressNew(k,3) + p c stateNew(k,7) = stateNew(k,7) + dgam * con1 * h2*(dcs0/dcs)**zz stateNew(k,7) = max(zero,stateNew(k,7)) c dalph = (h1+h3) * dgam2 do i=1,ntens stateNew(k,i) = stateNew(k,i) + dalph * xi(i)*(dcs0/dcs)**zz end do c c---- update plastic strain (for output purposes only) c stateNew(k,9) = stateOld(k,9) + dgam * con1 c c---- update temperature for adiabatic problems c dum = zero do i=4,ntens dum = dum + stressNew(k,i)*xi(i) end do stateNew(k,8) = stateOld(k,8) * + htcp*dgam2*( stressNew(k,1)*xi(1) * + stressNew(k,2)*xi(2) + stressNew(k,3)*xi(3) * + two*dum) c c---- update damage c epsdot = dgam*con1/dt sigeqv = (stressNew(k,1)-stressNew(k,2))**two * + (stressNew(k,2)-stressNew(k,3))**two * + (stressNew(k,3)-stressNew(k,1))**two * + 6.0d0*(stressNew(k,4)**2+stressNew(k,5)**2 * +stressNew(k,6)**2) sigeqv = sqrt(half*sigeqv) sige = max(1.d-15,sigeqv) c c Cocks-Ashby large pore growth term cacon = abs(vtheta/ytheta) if(cacon.lt.props(54)) cacon=props(54) dterm=two*(two*cacon-1)/(two*cacon+1) arg = min(15.,p*dterm/sige) beta = sinh(max(zero,arg) ) c90 = one + cacon psi = min(15.,beta*dt*epsdot*c90) tmp = max(zero,(one+(phi1**c90-one)*exp(psi))) stateNew(k,18) = min((one-tmp**(one/c90)),.99) c Cocks-Ashby void growth rate vy = vtheta/ytheta stateNew(k,19) = beta*epsdot*(one/(one-stateNew(k,18))**vy * -(one-stateNew(k,18))) c McClintock form of void growth sqtth = three**half abc = sqtth/(two*(one-props(30))) * * sinh(sqtth*half*(one-props(30))*(two*p/sige*third)) if(props(31).gt.zero) then vrad = props(31)*exp(stateNew(k,9)*abc/con1) else vrad = zero end if stateNew(k,10) = pi*vrad**two stateNew(k,11) = three*stateNew(k,10)*abc*epsdot c Nucleation of voids c deviatoric stress ds11 = stressNew(k,1)-p ds22 = stressNew(k,2)-p ds33 = stressNew(k,3)-p ds12 = stressNew(k,4) ds23 = stressNew(k,5) ds13 = stressNew(k,6) c invariants of stress di1 = three*p 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) if(dj2.le.zero) then r1 = zero r2 = zero r3 = zero else r1 = (4.d0/27.d0-dj3**2/dj2**3) r2 = dj3/dj2**(op5) r3 = di1/(dj2**half) endif r3 = abs(r3) zzz = (props(32)*r1+props(33)*r2+props(34)*r3) zzz = abs(zzz) zzzz = (props(37)**half/(props(36)*props(38)**third))*zzz stateNew(k,17) = stateOld(k,13) stateNew(k,13) = props(35)*exp(stateNew(k,9)*zzzz/con1) * * exp(-props(51)/stateNew(k,8)) c added for nonmonotonic path sequences, statev(17) is old nucleation if(stateNew(k,13).lt.stateNew(k,17)) then stateNew(k,13) = abs(stateNew(k,17)+stateNew(k,13)) endif c Coalescence factor cf = (cd1+cd2*stateNew(k,13)*stateNew(k,10))* * * exp(props(52)*stateNew(k,8))*(dcs0/dcs)**zz c Damage damage=cf*(stateNew(k,13)*stateNew(k,10)+stateNew(k,18)) if(damage.gt.0.6) damage = .99 stateNew(k,14) = min(damage,0.99) c stateNew(k,14) = zero c Nucleation Rate epsdot = abs(epsdot) stateNew(k,15)=zzzz*stateNew(k,13)*epsdot c Damage Rate zsecond = cf*(stateNew(k,15)*stateNew(k,10) * +stateNew(k,13)*stateNew(k,11) * +stateNew(k,19)) zthird = (stateNew(k,13)*stateNew(k,10)+stateNew(k,18)) * * cd2*(dcs0/dcs)**zz*exp(props(52)*stateNew(k,8)) * * (stateNew(k,15)*stateNew(k,10) * +stateNew(k,13)*stateNew(k,11)) stateNew(k,16)=zsecond+zthird c c Triaxiality stateNew(k,12) = p/sige c 200 continue c stateNew(k,7) = max(stateNew(k,7),zero) c 100 continue return end c subroutine vdisp( c Read only variables - 1 nblock, nDof, nCoord, kstep, kinc, 2 stepTime, totalTime, dtNext, dt, 3 cbname, jBCType, jDof, jNodeUid, amp, 4 coordNp, u, v, a, rf, rmass, rotaryI, c Write only variable - 5 rval ) C-----------------------------------------------------------------* C-----------------------------------------------------------------* C----- CALCULATION OF IMPOSED DISPLACEMENT AND VELOCITY --------* C-------------- TO OBTAIN CONSTANT STRAIN RATE -----------------* C-----------------------------------------------------------------* C-----------------------------------------------------------------* C--------- Variables to be defined : C C rval(nDof,*): Values of the prescribed variable for degrees of freedom 1–6 C (translation and rotation) at the nodes. C--------------------------------------------------* C--------------------------------------------------* C--------- Variables passed in for information : C C nblock : Number of nodal points to be processed in this call to VDISP. C nDof : Number of degrees of freedom (equals 6). C nCoord : Number of coordinate components (equals 3). C kstep : Step number. C kinc : Increment number. C stepTime : Value of time since the step began. C totalTime : Value of total time. The time at the beginning of the step is given by totalTime-stepTime. C dtNext : Next time increment size. C dt : Current time increment size. C cbname : User-specified name corresponding to the associated boundary condition. C jBCType : Indicator for type of prescribed variable: 0 for displacement, C 1 for velocity, and C 2 for acceleration. C jDof(nDof) : Indicator for prescribed degrees of freedom. C The values given by rval(j,k) are prescribed only if jDof(j) equals 1. C jNodeUid(*) : Node numbers. C amp(nblock) : Amplitude values corresponding to the associated amplitude functions. C coordMp(*,*) : Nodal point coordinates. C u(nDof,*) : Nodal point displacements at stepTime. C v(nDof, *) : Nodal point velocities at zero step time during initialization C or nodal point mean velocities at stepTime-dt/2 during time incrementation. C a(nDof, *) : Nodal point accelerations at stepTime before the boundary condition is prescribed. C rf(nDof, *) : Nodal point reaction at stepTime-dt. C rmass(*) : Nodal point masses. C rotaryI(3,3,*): Nodal point rotary inertia. C-----------------------------------------------------------------* C-----------------------------------------------------------------* C include 'vaba_param.inc' c character*80 cbname dimension jDof(nDof), jNodeUid(nblock), 1 amp(nblock), coordNp(nCoord,nblock), 2 u(nDof,nblock), v(nDof,nblock), a(nDof,nblock), 3 rf(nDof,nblock), rmass(nblock), rotaryI(3,3,nblock), 4 rval(nDof,nblock) data zero,half,one,two/0.0d0,0.5d0,1.0d0,2.0d0/ C-----------------------------------------------------------------* C-----------------------------------------------------------------* C--------- Values of the prescribed variable for DOF ------------* C-----------------------------------------------------------------* C xl0 = one ! Initial Element Length rate = 5.d-01 ! Strain Rate xret = xl0*exp(rate*totalTime) do 100 k = 1, nblock do 100 j = 1, nDof if( jDof(j) .gt. 0 ) then if(jBCType.eq.0) rval(j, k) = xret-xl0 ! Displacement if(jBCType.eq.1) rval(j, k) = rate*xret ! Velocity if(jBCType.eq.2) rval(j, k) = rate*rate*xret ! Acceleration end if 100 continue C C-----------------------------------------------------------------* C-----------------------------------------------------------------* return end C