c...thermoelasticity "main engine" c...for use with FELIPE version 3 - sample file is thermqd2.dat c program therm cp...program for steady-state solution of 2D thermoelasticity problems. cp...uses 4-noded or 8-noded quad elts where corner nodes have (u,v,p) cp...degrees of freedom (p=temperature). Example of a cp...coupled problem, and of additional d.o.f.s per node. Also uses cp...a Gaussian elimination solver for the nonsymmetric stiffness cp...matrix (stored as a band). Theory: solve cp cp / T \ / \ / \ cp | K Q | | u | | f | cp | | | | = | | cp | 0 -H | | p | | r | cp \ / \ / \ / cp...see manual for definition of terms. p is the temperature. 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,mbw=300) 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...mbw = max. bandwidth of global stiffness matrix cp...We have to use bandwidth and store whole band, because K is not cp...symmetric! c c character*1 char,dbdcmp(8) character*80 title,infile,outfil,prtfil c dimension lnods(melem,mnode),lnodn(melem),lsidn(melem), 1 coord(mpoin,2),nkode(mpoin),aslod(mtotv), 2 lpros(melem),props(msets,mcomp),asdis(mtotv), 4 surft(msets,2),gstif(mtotv,mbw),lmtyp(melem) 5 ,ndof(mpoin),nvar(mpoin),dbdcrd(8),dbdval(8) 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 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 vector c... 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 c...dbdcmp,dbdval,dbdcrd are component, coord, value for Dirichlet c... boundary value sets. c write(6,*)'Program THERM: 2D thermoelasticity problems.' write(6,*)'Handles 4-noded or 8-noded quadrilateral elements,' write(6,*)'with (u,v,t) d.o.f.s at all nodes.' 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 2 l=1,nelem nnode = lnodn(l) nside = lsidn(l) nodsid = 10*nnode + nside if (nodsid.ne.84.and.nodsid.ne.44) 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 cp...create the pointer array nvar. Needed when there are variable cp...dof's per node. Not strictly needed in THERM, but will be needed cp...in POREL. The no of dof's in node n is stored in ndof(n) nvar(1) = 0 nvtot = 0 do 5 np=1,npoin ndofnp = ndof(np) nvar(np) = nvtot nvtot = nvtot + ndofnp 5 continue cp...Now local dof i of node np is global dof nvar(np)+i cp...and nvtot is total no. of degrees of freedom c write(6,*)'Is this a plane stress or plane strain analysis?' 3 write(6,*)'Type 1 for plane stress, or 2 for plane strain...' 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 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,lsidn,coord,lpros,props 2 ,lmtyp,ndof,nvar) cp...previous line added for this and following routines 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,lmtyp,ndof,nvar,lpros,props) c c c...read data for Dirichlet boundaries, if present 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 c...read Dirichlet boundary values, if present, from end of .dat file call bndry(ndbdry,dbdcmp,dbdcrd,dbdval) endif c call prompt write(6,*) write(6,*)'Starting subroutine ASSMB...' call assmb(melem,mpoin,mnode,mtotv,mbw, 1 nelem,npoin,nbw,lnods,lnodn, 2 gstif,lmtyp,ndof,nvar) c call prompt write(6,*) write(6,*)'Starting subroutine SOLVE...' call solve(mpoin,mtotv,mbw,npoin,nbw,nkode, 1 coord,aslod,asdis,gstif,ndof,nvar, 2 ndbdry,dbdcmp,dbdcrd,dbdval) c call prompt write(6,*) write(6,*)'Starting subroutine RESULT...' call result(melem,mpoin,mnode,msets,mcomp,mtotv,mstyp, 1 nelem,npoin,lnods,lnodn,lsidn,coord,nkode,lpros,props, 2 aslod,asdis,lmtyp,ndof,nvar) c write(6,*) write(6,*)'THERM 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