c...Plane stress/strain elasticity "main engine" c...for use with FELIPE version 3 c... - sample datafile is elcyl4.dat or eldemo.dat c program elast c parameter(melem=100,mpoin=400,mnode=8, 1msets=4,mcomp=6,mtotv=800,mband=200) c 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 (= 2*mpoin for elasticity) c...mband = halfbandwidth+1 of global matrix gstif using band storage 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),gstif(mtotv,mband),lmtyp(melem), 5 nminvc(mtotv),astmp(mtotv),diagk(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...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...astmp is a scratch vector used in the equation solution substitutions c...surft(m,1),surft(m,2) = normal,tangential surface tractions in set m c...gstif holds assembled global stiffness matrix as a vector in skyline storage c... or as a matrix using band storage c...nminvc(i)=j means that G(i,j) is the first nonzero entry in row i of gstif c...lmtyp(l) = material type of elt l (not used in this program) c...diagk stores the diagonal of the stiffness matrix (used to find reactions) c c...ifull = 0 for demo version: restricted to 50 elts or 50 nodes ifull = 0 c write(6,*)'Program ELAST: solves 2D elasticity problems.' write(6,*)'Handles 8-noded quadrilateral elements.' write(6,*)'Dimensioned for',melem,' elements, and',mpoin,' nodes.' if (ifull.eq.0)write(6,*)'Demonstration datafile is eldemo.dat' 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 2 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 2 continue if (iflag.eq.1) then call prompt stop endif c if (ifull.eq.0.and.(nelem.gt.51.or.npoin.gt.51)) then write(6,*)'Demo version will not process meshes with' write(6,*)'51 elements or more than 51 nodes.' call prompt stop endif c write(6,*) 1'Is this a plane stress, plane strain or axisymmetric analysis?' 3 write(6,*) 1'Type 1 for plane stress, 2 for plane strain, or 3 for axisymm...' read(5,'(i1)',err=3) mstyp write(16,*) if (mstyp.eq.1) then write(6,*) 'Plane stress analysis' write(16,*)'Plane stress analysis' else if (mstyp.eq.2) then write(6,*) 'Plane strain analysis' write(16,*)'Plane strain analysis' else if (mstyp.eq.3) then write(6,*) 'Axisymmetric analysis' write(16,*)'Axisymmetric analysis' else write(6,*)'Invalid input!' goto 3 endif write(16,*) 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 call prompt write(6,*) write(6,*)'Starting subroutine STIFF...' call stiff(melem,mpoin,mnode,msets,mcomp,mstyp, 1 nelem,npoin,lnods,lnodn,coord,lpros,props) c call prompt write(6,*) write(6,*)'Starting subroutine LOAD...' call load(melem,mpoin,mnode,msets,mcomp,mtotv,mstyp, 1 nelem,npoin,lnods,lnodn,coord, 2 surft,aslod,lpros,props) c call prompt write(6,*) write(6,*)'Starting subroutine ASSMB...' call assmb(melem,mpoin,mnode,mtotv,mband, 1 nelem,npoin,nhbw,lnods,lnodn, 2 gstif,nminvc,diagk) c call prompt write(6,*) write(6,*)'Starting subroutine SOLVE...' call solve(mpoin,mtotv,mband,npoin,nkode,nhbw, 1 aslod,asdis,astmp,gstif,nminvc,diagk) c call prompt write(6,*) write(6,*)'Starting subroutine RESULT...' call result(melem,mpoin,mnode,msets,mcomp,mtotv,mstyp, 1 nelem,npoin,lnods,lnodn,coord,nkode,lpros,props, 2 aslod,asdis) c write(6,*) write(6,*)'ELAST completed. Results are in file ',prtfil(1:lenin) write(6,*) write(6,*) call prompt close(15) close(16) close(17) stop 999 write(6,*) 1'Cannot open file with this name. Terminating program.' call prompt stop end c c c subroutine stiff(melem,mpoin,mnode,msets,mcomp,mstyp, 1 nelem,npoin,lnods,lnodn,coord,lpros,props) c c...calculates the element stiffness matrices, and writes them c...to the scratch file. c dimension props(msets,mcomp),lnods(melem,mnode),encod(2,8), 1 estif(16,16),bmatx(4,16),dmatx(4,4),dbmat(4,16), 2 gpwts(4),lpros(melem),lnodn(melem),coord(mpoin,2), 3 shape(8),gploc(2,4),gpcod(2),deriv(2,8),cartd(2,8) c...encod holds the nodal coords for the nodes around the current element c...estif will be the stiffness matrix for the current element c...bmatx, dmatx are the B and D matrices, and dbmat is D*B c...gpwts,gploc are the Gauss point weights and locations (in xi,eta coords) c...gpcod will be the x,y coords of the current Gauss point c...shape holds the shape functions evaluated at the Gauss point c...deriv holds the derivatives (w.r.t. xi, eta) of the shape functions c...cartd holds the Cartesian derivatives (w.r.t x,y) of shape functions c rewind 18 c c...for axisymmetric analyses, D is a 4x4 matrix; otherwise it's 3x3 ndmatx = 3 if (mstyp.eq.3) ndmatx = 4 c...and we need the value of pi, for integrating. Form it using inverse tangent pi = 4.*atan(1.0) c...start loop over all elements do 80 l=1,nelem c c...extract Young's modulus & Poisson's ratio from props lset = lpros(l) ym = props(lset,1) pr = props(lset,2) c...also extract material thickness (default is 1.0) for plane stress/strain thick = props(lset,3) if (thick.eq.0.0) thick = 1.0 c c...put nodal coords in encod nnode = lnodn(l) do 15 n=1,nnode nn = lnods(l,n) encod(1,n) = coord(nn,1) encod(2,n) = coord(nn,2) 15 continue c...zero stiffness matrix do 20 j=1,16 do 20 i=1,16 estif(i,j) = 0.0 20 continue c c...get D matrix call dmat(mstyp,ym,pr,dmatx) c c...put local coords of 2x2 guass-pts in gploc call gauss(gploc,gpwts,ngaus) c c...start loop through gauss-pts for contribution to estif do 75 ij=1,ngaus c...put local coords of gauss-pt in s,t s = gploc(1,ij) t = gploc(2,ij) c...get shape fns (shape) & derivs (deriv) in local coords call sfr(s,t,shape,deriv) c...get xy-coords of gauss-pt into gpcod, using shape fn interpolation do 25 m=1,2 sum = 0.0 do 24 n=1,8 24 sum = sum + encod(m,n)*shape(n) gpcod(m) = sum 25 continue c...for axisymmetric analyses, we need the radius at the Gauss point if (mstyp.eq.3) radius = gpcod(1) c...get xy-derivs (cartd) & derterminant of Jacobian (djacb) call jacob(encod,shape,deriv,cartd,djacb) c...get B matrix call bmat(cartd,mstyp,radius,shape,bmatx) c...form matrix product D*B do 40 j=1,16 do 40 i=1,ndmatx sum=0.0 do 39 k=1,ndmatx 39 sum = sum + dmatx(i,k)*bmatx(k,j) dbmat(i,j) = sum 40 continue c...integrate B^T *D*B over volume of element darea = gpwts(ij)*djacb c...for axisymmetry, the 'thickness' is the circumference if (mstyp.eq.3) then dvolu = darea*2.*pi*radius c...otherwise, get volume by multiplying by thickness else dvolu = darea*thick endif c...B^T*D*B*dvolu = contribution to estif do 50 j=1,16 do 50 i=1,16 sum = 0.0 do 49 k=1,ndmatx 49 sum = sum + bmatx(k,i)*dbmat(k,j) estif(i,j) = estif(i,j) + sum*dvolu 50 continue 75 continue c...end of gauss loop c c...write estif to file write (18) estif write(6,*)'Stiffness matrix for elt ',l,' written to file.' 80 continue c...end of element loop c write(6,*)'Elt stiffness matrices written to file' endfile 18 c return end