c...poroelasticity "main engine" c...for use with FELIPE version 3 - sample file is conslex2.dat c program consl cp...program for poroelasticity (consolidation) problems. cp...uses 8-noded quad elts where corner nodes have (u,v,p) cp...degrees of freedom (p=pore pressure). Example of a cp...coupled problem. Theory: solve cp _ __ __ _ cp (X + Y) u = th f + f - (th X - Y) u cp n+1 n n+1 n cp __ cp th = (1.0-theta)/theta cp cp...Note: to make sign of pore pressures compatible with the cp...compression-positive sign convention for u,f, -Q in cp...standard theory is replaced by Q cp...changes to elast.for are prefaced by cp... comment lines c parameter(melem=200,mpoin=400,mnode=8, 1msets=4,mcomp=6,mtotv=1200,mstif=100000) 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 (e.g. Young's mod, Poisson ratio) c...mtotv = max no. of degrees of freedom (corner nodes have 3 dofs) c...mstif = max. length of vector gstif storing global matrix using c skyline storage c c c character*1 char character*80 infile,outfil,prtfil,title c dimension lnods(melem,mnode),lnodn(melem),lsidn(melem), 1 coord(mpoin,2),nkode(mpoin),aslod(mtotv),asrhs(mtotv), 2 lpros(melem),props(msets,mcomp),asdis(mtotv), 4 surft(msets,2),gstif(mstif),lmtyp(melem) 5 ,ndof(mpoin),nvar(mpoin),nnout(5) 6 ,nminvc(mtotv),lpoint(mtotv),astmp(mtotv) 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...asrhs is right-hand-side vector at each timestep c...ndof(np) = no. of dof.s at node np cp..nvar(np) is a pointer array for global dof of node np 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/pressure vector in current timestep c...surft(m,1),surft(m,2) = normal,tangential surface tractions in set m c...gstif holds assembled global stiffness matrix in compact form c...lmtyp(l) = material type of elt l (1 for structure, 2 for soil) c...nminvc(i)=j means that G(i,j) is the first nonzero entry in row i of gstif c...lpoint(i)=ii means that G(i,i) is stored at gstif(ii) in skyline storage c write(6,*)'Program CONSL: soil consolidation problems.' write(6,*)'Handles 8-noded quadrilateral elements,' write(6,*)'with (u,v,t) d.o.f.s at corner nodes, and (u,v) at' write(6,*)'midside nodes. Time-dependent analysis. Writes' write(6,*)'output for undrained and drained solutions.' 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 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 nvar(1) = 0 ntotv = 0 do 5 np=1,npoin ndofnp = ndof(np) nvar(np) = ntotv ntotv = ntotv + ndofnp 5 continue cp...Now local dof i of node np is global dof nvar(np)+i cp...and ntotv is total no. of degrees of freedom c c...analyse the "mask" of the global stiffness matrix, i.e. the range c...of non-zero values in each row. This information is stored in c...nminvc: row i starts in column nminvc(i). This is used to construct c...lpoint: g(i,i) will be located at gstif(lpoint(i)) c...this is a generalisation of the s/r in PLADV, to allow for variable c...d.o.f.s per node c call mask(melem,mpoin,mnode,mtotv,mstif,ntotv, 1 nelem,npoin,nhbw,lnods,lnodn,nminvc,lpoint,ndof,nvar) c...doing this now, avoids repeating it within assmb at each iteration! c write(6,*)'Type the numbers of 5 nodes ' write(6,*)'for which you would like the displacements at each' write(6,*)'timestep to be output:' read(5,*) (nnout(i),i=1,5) c c...open a scratch file to write the element stiffness matrices on open(unit=18,status='scratch',form='unformatted') write(6,*)'Scratch file opened ...' write(6,*) c cp...STIFF forms the X and Y matrices on file 18 for each element call prompt write(6,*) write(6,*)'Starting subroutine STIFF...' call stiff(melem,mpoin,mnode,msets,mcomp, 1 nelem,npoin,lnods,lnodn,lsidn,coord,lpros,props 2 ,lmtyp,ndof,nvar) c 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,aslod,lmtyp,ndof,nvar) c c...read data for boundaries (should be blank) read(15,*) read(15,'(2i5)')nbrad,nbdir c if (nbrad.gt.0) then write(6,*) 'ERROR: nbrad > 0.' write(6,*)'This program cannot handle radiating boundaries.' call prompt stop endif c if (nbdir.gt.0) then write(6,*) 'ERROR: nbdir > 0.' write(6,*)'This program cannot handle Dirichlet boundaries.' cp...Permeable boundaries (where the excess pore pressure p is fixed cp...at zero) are handled in the same way as fixed displacement dof's, cp...via the fixity code only. call prompt stop endif c c...read data defining load increments (not expecting any) read(15,*) read(15,'(i5)') nincs if (nincs.gt.0) then write(6,*) 1 'Error: cannot handle load increments - data ignored.' endif c cv...read the timestepping control parameters read(15,*) read(15,*) dtint,ftime,allow,theta,ttend,dtone cv...theta is the time discretization parameter: =0 for explicit algorithm, cv... between 0.5 and 1.0 for stable implicit algorithm. cv...dtone is the small timestep used for the initial solution cv...dtint is the initial timestep length in timestepping 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) cp...stop when time reaches TTEND c tht = (1.0-theta)/theta c cv...initialize the timestepping ttime = 0.0 dtime = dtone c c...form r.h.s. vector from initial load aslod. We take f_0 = 0 and c...f_1 = aslod. So r.h.s. is tht*f_0 + f_1 = f_1 c...There is no contibution from v_0, which is zero. do 20 i=1,ntotv asrhs(i) = aslod(i) 20 continue c write(6,*)'Undrained solution' write(6,*)'------------------' write(6,*) write(6,*)'Undrained timestep = ',dtone write(6,*) write(6,*)'Starting subroutine ASSMB...' call assmb(melem,mpoin,mnode,mtotv,mstif, 1 nelem,npoin,lnods,lnodn, 2 gstif,lmtyp,ndof,nvar,theta,dtime,nminvc,lpoint) c...add 10e12 to diagonal elements for fixed dof (including cp...pore pressure dof's on permeable boundaries). write(6,*)'Starting subroutine FIXDOF...' call fixdof(mtotv,mstif,ntotv,mpoin,npoin,gstif,nkode,coord, 1 asrhs,ndof,nvar,lpoint) c... T c...solution of equations by L D L decomposition write(6,*) write(6,*)'Starting matrix decomposition in LDLDEC...' call ldldec(mtotv,mstif,ntotv,gstif,lpoint) c...solve write(6,*) write(6,*)'Starting equation substitution in LDLSUB...' call ldlsub(mtotv,mstif,ntotv,gstif,asrhs,asdis,astmp, 1 lpoint) c ttime = ttime + dtime write(6,*) write(6,*)'Time = ',ttime write(6,*) c call result(melem,mpoin,mnode,msets,mcomp,mtotv, 1 nelem,npoin,lnods,lnodn,lsidn,coord,nkode,lpros,props, 2 aslod,asdis,lmtyp,ndof,nvar) c kount = 0 dtime = dtint write(6,*)'Starting timestepping, with timestep = ',dtime write(6,*) c write(6,*)'Forming new right-hand side vector...' cv...now form new r.h.s. for next timestep in asrhs c...now f_1 = aslod and f_2 = aslod, so rhs = tht*f_1 + f_2 do 30 i=1,ntotv asrhs(i) = tht*aslod(i) + aslod(i) 30 continue call newrhs(melem,mpoin,mnode,mtotv, 1 nelem,npoin,lnods,lnodn,nkode, 2 asrhs,asdis,lmtyp,ndof,nvar,theta,dtime) c cc..start of global timestepping cc..============================ 100 continue kount = kount + 1 c write(6,*)'Assembling...' call assmb(melem,mpoin,mnode,mtotv,mstif, 1 nelem,npoin,lnods,lnodn, 2 gstif,lmtyp,ndof,nvar,theta,dtime,nminvc,lpoint) write(6,*)'Applying fixities...' call fixdof(mtotv,mstif,ntotv,mpoin,npoin,gstif,nkode, 1 coord,asrhs,ndof,nvar,lpoint) write(6,*)'Decomposing matrix...' call ldldec(mtotv,mstif,ntotv,gstif,lpoint) c cp...perform 5 timesteps using the timestep dtime do 180 kk=1,5 c cp...solve - no need to re-assemble or decompose, as if dtime is cp...unchanged then stiffness matrix will be unchanged write(6,*)'Performing substitutions...' call ldlsub(mtotv,mstif,ntotv,gstif,asrhs,asdis,astmp, 1 lpoint) c ttime = ttime + dtime write(6,*) write(6,*)'Time = ',ttime write(6,*) c do 170 i=1,5 cp...output results at selected nodes nn = nnout(i) ndofnn = ndof(nn) nvarnn = nvar(nn) write(6,*)'Displts/pressure at node ',nn dispx = asdis(nvarnn+1) dispy = asdis(nvarnn+2) if (ndofnn.eq.3) then pore = asdis(nvarnn+3) write(6,610) dispx,dispy,pore else write(6,610) dispx,dispy endif 610 format(3(5x,e12.5)) 170 continue call prompt c c...form new rhs vector (if we will not be changing the timestep) if (kk.lt.5) then write(6,*)'Forming new right-hand side vector...' c...both f_n and f_n+1 are = aslod in subsequent timesteps do 175 i=1,ntotv asrhs(i) = tht*aslod(i) + aslod(i) 175 continue call newrhs(melem,mpoin,mnode,mtotv, 1 nelem,npoin,lnods,lnodn,nkode, 2 asrhs,asdis,lmtyp,ndof,nvar,theta,dtime) endif 180 continue c call yesno(char,'Do you want to stop?') if (char.eq.'y') goto 200 c c...output results to file write(6,665) ttime write(16,665) ttime write(17,665) ttime 665 format('Consolidation progressing: time = ',e10.3) c call result(melem,mpoin,mnode,msets,mcomp,mtotv, 1 nelem,npoin,lnods,lnodn,lsidn,coord,nkode,lpros,props, 2 aslod,asdis,lmtyp,ndof,nvar) c cv...set next timestep dnext = ftime*dtime if (dnext.gt.allow*dtint) dnext = allow*dtint dtime = dnext write(6,*) write(6,*)'Changing timestep length to ',dtime write(6,*) c cv...now form new r.h.s. for next timestep in asrhs write(6,*)'Forming new right-hand side vector...' do 185 i=1,ntotv asrhs(i) = tht*aslod(i) + aslod(i) 185 continue call newrhs(melem,mpoin,mnode,mtotv, 1 nelem,npoin,lnods,lnodn,nkode, 2 asrhs,asdis,lmtyp,ndof,nvar,theta,dtime) c if (ttime.lt.ttend) goto 100 c cc..end of timestepping cc..=================== c 200 write(6,*)'END OF TIMESTEPPING' write(6,*) write(6,*) kount,' global timesteps' write(6,*) c c write(6,666) ttime write(16,666) ttime write(17,666) ttime 666 format('Drained solution. Time = ',e10.3) c call result(melem,mpoin,mnode,msets,mcomp,mtotv, 1 nelem,npoin,lnods,lnodn,lsidn,coord,nkode,lpros,props, 2 aslod,asdis,lmtyp,ndof,nvar) c write(6,*) write(6,*)'POREL 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