c...Plane strain elasto-viscoplasticity "main engine" c...for use with FELIPE version 3 - sample datafile is vpcyl4.dat program vplas c cv...This program performs elast-viscoplasticity analyses. cv...It is developed from the plasticity program plast.for cv...Added code for viscoplasticity is preceded by comment lines cv...which begin with cv (like these lines!) cv...VPLAS allows for non-associated flow: therefore a cv...frontal solution routine for nonsymmetric matrices is used. cv...plane strain conditions, with material thickness = 1.0 c CV...THE PROGRAM (AND FELSUB.FOR) SHOULD BE COMPILED IN DOUBLE PRECISION: CV... ftn77 /dreal vplas.for in Salford Fortran77 c parameter(melem=800,mpoin=2000,mnode=8, 1msets=4,mcomp=6,mtotv=4000,mfixv=150, 2mfron=150,mstif=22500,mbufa=100) c...these parameters fix the dimensions of the arrays: c...melem = max. no. of elements c...mpoin = max. no. of nodes c...mnode = max. no. of nodes in a single element c...msets = max. no. of mat. prop. sets (also used for loading sets) c...mcomp = max. no. of mat. props = 6 for (visco)plasticity c...mtotv = max no. of degrees of freedom (= 2*mpoin for elasticity) cv...mfixv = max no. of fixed degrees of freedom cv...mfron = max frontwidth in solver FRONT cv...mstif = mfron*mfron in non-symmetric case c c character char*1 character*80 title,infile,outfil,prtfil c dimension lnods(melem,mnode),lnodn(melem),lsidn(melem), 1 coord(mpoin,2),nkode(mpoin),aslod(mtotv),ndof(mpoin), 2 lpros(melem),props(msets,mcomp),asdis(mtotv), 4 surft(msets,2),lmtyp(melem) c...lnods(l,n) = global no. of n'th node of element l c...lnodn(l) = no. of nodes in element l ( = 8 for 8-noded quad.) c...lsidn = no. of sides in element l ( = 4 for 8-noded quad.) c...coord(np,1), coord(np,2) = x,y coords of node np c...nkode(np) = fixity code of node np c...aslod is global nodal load vector c...ndof(np) = no. of dof.s at node np (not used in this program) c...lpros(l) = mat. prop. set no. of element l c...props(m,i) = value of i'th mat. property in mat. prop. set m c...asdis is global nodal displacement vector c...surft(m,1),surft(m,2) = normal,tangential surface tractions in set m c...lmtyp(l) = material type of elt l (not used in this program) c real*8 pival,eqrhs(mbufa) 1 ,vecrv(mfron),gstif(mstif),gload(mfron),equat(mfron,mbufa) DIMENSION NPIVO(MBUFA),NAMEV(MBUFA),NACVA(MFRON),nfvar(mfixv), 1 spdis(mfixv),nvar(mpoin) cb...arrays above are used in the FRONT frontal solver c dimension strsg(3,4,melem),strsg1(3,4,melem),totld(mtotv), 1 accld(mtotv),accds(mtotv),aslod1(mtotv),resld(mtotv), 2 fincs(9),mohrg(4,melem),vivel(3,4,melem) cc..strsg holds stresses at Gauss points (Cartesian) cc..totld is the total load to be applied cc..accld is the accumulated load applied so far cc..aslod is the load being applied in the current increment cc..accds is the accumulated displacements cc..resld is the residual force vector cc..fincs is the set of load increments cc..mohrg is the yield code at Gauss pts: 0=elastic, 1=plastic cv...vivel stores the viscoplastic strain velocities ("epsilon dot") c c cc..tolerances for convergence tests: rtoler = 0.1 stoler = 1.e-8 cc..rtoler is pl. strain velocities tolerance (%) for the timesteps cc..stoler is tolerance for stresses in local (stress) iterations cc..stoler is used in subroutine stress c write(6,*)'Program VPLAS: elasto-viscoplasticity (Mohr-Coulomb).' write(6,*)'Handles 8-noded quadrilateral elements.' write(6,*)'Allows for non-associated flow.' write(6,*)'Dimensioned for',melem,' elements, and',mpoin,' nodes.' write(6,*) c c...read an existing .dat file c 1 continue c c...get the input file name in infile call askfil(infile,iflag) if (iflag.eq.1) then call yesno(char,'Do you want to exit?') if (char.eq.'y') stop goto 1 endif c write(6,*)'Opening file ',infile,'...' lenin = lenstr(infile) open(15,file=infile(1:lenin),status='old',err=999) write(6,*) lm4 = lenin - 4 prtfil = infile(1:lm4)//'.prt' write(6,*)'Creating file ',prtfil(1:lenin),' for results...' write(6,*) open(16,file=prtfil(1:lenin),status='unknown',err=999) outfil = infile(1:lm4)//'.out' write(6,*)'Creating file ',outfil(1:lenin),' for output...' write(6,*) open(17,file=outfil(1:lenin),status='unknown',err=999) c c...read the input mesh data, and echo it to the output files call input(15,melem,mpoin,mnode,msets,mcomp,nelem, 1 npoin,lnods,lnodn,lsidn,lmtyp,coord,nkode,lpros,props 2 ,ndofn,npros,title,ndof) call print(16,melem,mpoin,mnode,msets,mcomp,nelem,npoin, 1 lnods,lnodn,lsidn,lmtyp,coord,nkode,lpros,props,ndofn, 2 npros,title,ndof) call output(17,melem,mpoin,mnode,msets,mcomp,nelem,npoin, 1 lnods,lnodn,lsidn,lmtyp,coord,nkode,lpros,props,ndofn, 2 npros,title,ndof) call yesno(char,'Do you want the mesh details listed?') if (char.eq.'y') then call print(6,melem,mpoin,mnode,msets,mcomp,nelem,npoin, 1 lnods,lnodn,lsidn,lmtyp,coord,nkode,lpros,props,ndofn, 2 npros,title,ndof) endif c iflag = 0 do 3 l=1,nelem nnode = lnodn(l) nside = lsidn(l) nodsid = 10*nnode + nside if (nodsid.ne.84) then iflag = 1 write(6,*)'ERROR: Cannot handle element ',l write(6,*)'which has ',nside,' sides and ',nnode,' nodes' endif 3 continue if (iflag.eq.1) then call prompt stop endif c c...iprint tells if we have output any results yet iprint = 0 c cc..zero the accumulating variables, stresses, etc. cv...and vivel c do 10 i=1,mtotv accds(i) = 0.0 accld(i) = 0.0 resld(i) = 0.0 10 continue do 11 i=1,9 fincs(i) = 0.0 11 continue do 12 l=1,melem do 12 ij=1,4 do 12 i=1,3 strsg(i,ij,l) = 0.0 strsg1(i,ij,l) = 0.0 vivel(i,ij,l) = 0.0 12 continue do 13 l=1,melem do 13 ij=1,4 mohrg(ij,l) = 0 13 continue c cp...create the pointer array nvar. Needed because a node may cp...have 2 or 3 d.o.f.s; this info stored in ndof cp...also store info on fixed d.o.f.s in nfvar and spdis, cp...rather than in nkode ntotv = 0 nfixv = 0 do 5 np=1,npoin ndofnp = ndof(np) nvar(np) = ntotv kode = nkode(np) C...SORT OUT KODE AND STORE FIXED VARIABLE NO. IN NFVAR() KOD = KODE DO 4 M=1,NDOFNP KD = KOD/10 K = KOD - KD*10 KOD = KD IF(K.ne.0) then NFIXV=NFIXV+1 IF (NFIXV.GT.MFIXV) THEN WRITE(6,*)'No. of fixed dof.s exceeds parameter MFIXV' write(6,*)'Increase MFIXV and recompile.' call prompt stop endif M1 = NDOFNP + 1 - M NFVAR(NFIXV) = ntotv + M1 spdis(nfixv) = 0.0 endif 4 continue ntotv = ntotv + ndofnp 5 continue nsdis = nfixv cp...s/r LOAd will add spec displts, total no will be held in nsdis cp...Now local dof i of node np is global dof nvar(np)+i cp...and nfvar(i) is global dof of i'th fixed dof cp...and ntotv is total no. of degrees of freedom c cc..read the equivalent load vector, into totld, and then cc..read the load increment sizes call prompt write(6,*) write(6,*)'Starting subroutine LOAD...' call load(melem,mpoin,mnode,msets,mcomp,mtotv, 1 nelem,npoin,lnods,lnodn,coord, 2 surft,totld) c c c...read boundary data header (no boundary set data expected) read(15,*) read(15,*)nbrad,nbdir if (nbrad.gt.0.or.nbdir.gt.0) then write(6,*) 1 'ERROR: boundary set data present - cannot cope!' call prompt stop endif c...read data defining load increments read(15,*) read(15,*) nincs,(fincs(i),i=1,nincs) cv...read the timestepping control parameters read(15,*) read(15,*) dtint,ftime,allow,theta cv...theta is the time discretization parameter: =0 for explicit algorithm, cv... between 0.5 and 1.0 for stable implicit algorithm. cv...dtint is the initial timestep length cv...if no new g.pts. have yielded, next timestep is FTIME times previous one cv... (DNEXT = FTIME*DTIME) cv...max. timestep is ALLOW times initial one (ALLOW*DTINT) c c c...open a scratch file to write the element stiffness matrices on open(unit=18,status='scratch',form='unformatted') cv...open further scratch files for use in FRONT OPEN(19,STATUS='SCRATCH',FORM='UNFORMATTED') OPEN(20,STATUS='SCRATCH',FORM='UNFORMATTED') OPEN(21,STATUS='SCRATCH',FORM='UNFORMATTED') PIVAL = 1.d-8 c c c cc..start of load increment loop cc..============================ c cv...initialize the timestepping ttime = 0.0 dtime = dtint c irsol = 0 cv...irsol = 0 for first solution, in front or afront cv...irsol = 1 for re-soln with new matrix having same structure cv...irsol = 2 for re-solution with already-decomposed matrix isymm = 1 cv...isymm = 1 for symmetric matrix, isymm=0 for nonsymm matrix c do 200 inc=1,nincs finc = fincs(inc) write(6,*)' Load increment ',inc write(6,*) write(6,*)' Increment size: ',finc write(6,*) cc..increment load, adding residual from previous increment do 20 i=1,ntotv asdis(i) = 0.0 aslod(i) = finc*totld(i) accld(i) = accld(i) + aslod(i) aslod(i) = aslod(i) + resld(i) 20 continue cv...get norm of the loads (to use in testing if steady state reached) call norm2(mtotv,ntotv,aslod,fnorm) do 21 l=1,melem do 21 ij=1,4 do 21 i=1,3 strsg1(i,ij,l) = strsg(i,ij,l) 21 continue cc..strsg will stay fixed as the start-of-timestep stress cc..strsg1 will hold the end-of-timestep stress c kount = 0 cc..start of global timestepping cc..============================ 100 continue kount = kount + 1 cc..note: stiffness depends on current stress state call stiff(melem,mpoin,mnode,msets,mcomp, 1 nelem,npoin,lnods,lnodn,coord,lpros,props 2 ,strsg,mohrg,dtime,theta) c do 110 i=1,ntotv aslod1(i) = aslod(i) 110 continue c call solve(mpoin,mtotv,mhbw,npoin,nhbw,nkode, c 1 aslod1,asdis,gstif) c cb...the next subroutine replaces ASSMB and SOLVE in ELAST.FOR c call prompt if (isymm.eq.1) then write(6,*) write(6,*)'Starting subroutine FRONT...' CALL FRONT(MTOTV,MPOIN,MELEM,MFIXV,MNODE,MFRON,MSTIF, 1NPOIN,NELEM,NSDIS,LNODS,LNODN,NDOF,NVAR,NFVAR,ASLOD1,ASDIS,SPDIS, 2VECRV,GSTIF,GLOAD,EQUAT,NACVA, 3PIVAL,IRSOL,EQRHS,NPIVO,NAMEV,MBUFA) irsol = 1 cv...if time parameter theta=0.0, only elastic matrix is used if (theta.eq.0.0) irsol = 2 else cv...non-symmetric matrix eqn must be solved CALL AFRONT(MTOTV,MPOIN,MELEM,MFIXV,MNODE,MFRON,MSTIF, 1NPOIN,NELEM,NSDIS,LNODS,LNODN,NDOF,NVAR,NFVAR,ASLOD1,ASDIS,SPDIS, 2VECRV,GSTIF,GLOAD,EQUAT,NACVA, 3PIVAL,IRSOL,EQRHS,NPIVO,NAMEV,MBUFA) endif c c cc..now evaluate new stresses, from this displt increment, in strsg1 cv..In v.p. program, mohrg is not updated here write(6,*)'Evaluating new stresses...' call stress(melem,mpoin,mnode,msets,mcomp,mtotv, 1 nelem,npoin,lnods,lnodn,coord,lpros,props, 2 asdis,strsg,strsg1,mohrg,stoler 3 ,vivel,dtime,theta) c cv...update stresses and displacements do 160 i=1,ntotv accds(i) = accds(i) + asdis(i) 160 continue do 161 l=1,nelem do 161 ij=1,4 do 161 i=1,3 strsg(i,ij,l) = strsg1(i,ij,l) 161 continue c cv...calculate new plastic strain velocities in vivel, & update mohrg write(6,*)'Updating yield...' call strain(melem,msets,mcomp,nelem,lpros,props, 1 strsg,vivel,mohrg,dtime,theta,newfl) c cv...set next timestep - newfl tells if new g.pt. has yielded ttime = ttime + dtime dnext = ftime*dtime if (newfl.ne.0) dnext = dtint if (dnext.gt.allow*dtint) dnext = allow*dtint dtime = dnext c cv...also set irsol and isymm: cv...matrix becomes nonsymm if any yielding has started if (isymm.eq.1.and.newfl.ne.0.and.theta.ne.0.0) isymm=0 if (newfl.eq.0.and.irsol.eq.0) then irsol = 2 isymm = 1 else irsol = 1 endif write(6,*)'irsol, isymm:',irsol,isymm c cv...now form new r.h.s. for next timestep in aslod write(6,*)'Forming new right-hand side vector...' call newrhs(melem,mpoin,mnode,mtotv,msets,mcomp, 1 nelem,npoin,ntotv,lpros,props,lnods,lnodn,coord, 2 theta,dtime,nkode,strsg,aslod,vivel,mohrg) c cv...check if steady state reached write(6,*)'Checking for steady state...' call steady(melem,nelem,vivel,mohrg,kount,first,ratio) write(6,*)'Time = ',ttime,' pl. strain ratio = ',ratio if (ratio.gt.rtoler) goto 100 c cc..end of timestepping for this increment cc..====================================== 150 continue write(6,*) kount,' global timesteps' write(6,*) c cc..compute residual force vector call residu(melem,mpoin,mnode,mtotv,nelem,npoin,ntotv, 1 lnods,lnodn,coord,nkode,accld,strsg,resld) c if (inc.lt.nincs) then c...output results if required call yesno(char,'Do you want the results output?') if (char.eq.'y') then if (iprint.eq.1) then write(6,666) inc,ttime write(16,666) inc,ttime write(17,666) inc,ttime endif iprint = 1 call result(melem,mpoin,mnode,msets,mcomp,mtotv, 1 nelem,npoin,lnods,lnodn,coord,nkode,lpros,props, 2 accld,accds,strsg,mohrg) endif endif c 200 continue cc..end of load increment loop cc..========================== c if (iprint.eq.1) then write(6,666) nincs,ttime write(16,666) nincs,ttime write(17,666) nincs,ttime 666 format(' End of load inc ',i5,'. Time = ',e10.3) endif call result(melem,mpoin,mnode,msets,mcomp,mtotv, 1 nelem,npoin,lnods,lnodn,coord,nkode,lpros,props, 2 accld,accds,strsg,mohrg) c write(6,*) write(6,*)'VPLAS completed. Results are in file ',prtfil(1:lenin) write(6,*) write(6,*) call prompt stop 999 write(6,*) 1'Cannot open file with this name. Terminating program.' call prompt stop end