【精品】Fortran语言编写的有限元结构程序

Fortran语言编写弹性力学平面问题3节点三角形单元或4节点等参单元的有限元程序：c--------------------------------------------------------------------------c.....FEA2DP---A finite element analysis program for2D elastic problemscc Tangent matrix is stored with varioud band methodc This program is used to demonstrte the usage of vrious bandc Storage schem of symmetric and unsymmetric tangent matrixcc Wang shunjinc At chongqing vniversity(06/06/2013)c-------------------------------------------------------------------------program FEA2DPcc a(1)-a(n1-1):x(ndm,nummnp);a(n1)-a(n2-1):f(ndf,numnp)c a(n2)-a(n3-1):b(neq);a(n3)-a(n4-1):ad(neq)c a(n4)-a(n5-1):al(nad);a(n5)-a(n6-1):nu(nad)cc ia(1)-ia(n1-1):ix(nen1,numel);ia(n1)-ia(n2-1):id(ndf,numnp)c ia(n2)-ia(n3-1):jd((ndf*numnp);ia(n3)-ia(n4-1):idl(nen*numel*ndf)cimplicit real*8(a-h,o-z)dimension a(100000),ia(1000)character*80headcommon/cdata/numnp,numel,nummat,nen,neqcommon/sdata/ndf,ndm,nen1,nstcommon/iofile/ior,iowcnmaxm=100000imaxm=1000ior=1iow=2cc Open files for data input and outputcopen(ior,file='input.dat',form='formatted')open(iow,file='output.dat')cc.....Read titlecread(ior,'(a)')headwrite(iow,'(a)')headcc.....Read and print control informationcc numnp:number of nodesc numel:number of elementsc nummat:number of material typesc nload:number of loadsc ndm:number of coordinats of each nodec ndf:number of degrees of freedomc nen:number of nodes in each elementcread(ior,'(7i5)')numnp,numel,nummat,nload,ndm,ndf,nenwrite(iow,2000)numnp,numel,nummat,nload,ndm,ndf,nen cc.....Set poiters for allocation of data arrayscnen1=nen+4nst=nen*ndfnneq=ndf*numnpcn1=ndm*numnp+1n2=n1+ndf*numnp+1ci1=nen1*numel+1i2=i1+ndf*numnp+1i3=i2+ndf*numnp+1i4=i3+numel*nen*ndf+1cc.....Call mesh input subroutine to read all mesh dataccall pmesh(a(1),a(n1),ia(1),ia(i1),ndf,ndm,nen1,nload)cpute profileccall profil(ia(i2),ia(i3),ia(i1),ia(1),ndf,nen1,nad)cn3=n2+neq+1n4=n3+neq+1n5=n4+nad+1n6=n5+nad+1cc The lengthes of real and integer arrayscwrite(iow,2222)n6,i4cc The lengthes of array exceeds the limitationcif(n6>nmaxm.or.i4>imaxm)thenif(n6>nmaxm)write(iow,3333)n6,nmaxmif(i4>nmaxm)write(iow,4444)i4,imaxmstopend ifctute and aseemble element arraysccall assem(nad,ia(1),ia(i1),ia(i2),a(1),a(n2),a(n3),1a(n4),a(n5))cc Form load vectorccall pload(ia(i1),a(n1),a(n2),nneq,neq)cc.....Triangular decomposition of a matrix stored in profile formccall datri(ndf,numnp,ia(i2),neq,nad,.false.,a(n3),a(n5),a(n5))cc For unsymmtric tangent matirxc Call datri(ndf,numnp,ia(i2),neq,nad,.true.,a(n3),a(n4),a(n5))cc Solve equationsccall dasol(ndf,numnp,a(n2),ia(i2),neq,nad,aengy,a(n3),a(n5),a(n5)) cc For unsymmetric tangent matrixc Call dasol(ndf,numnp,a(n2),ia(i2),neq,nad,aengy,a(n3),a(n5),a(n5)) cc Output nodal displacementsccall prtdis(ia(i1),a(n2),ndf,numnp,neq)cc.....Close input and output files;destroy temporary disk filescclose(ior)close(iow)cc.....Input/output formatsc1000format(20a4)2000format(//x5x,'number of nodal points=',i6/15x,'number of elements=',i6/25x,'number of material sets=',i6/35x,'number of nodal loads=',i6/45x,'dimension of coordinate space=',i6/55x,'degree of freedoms/node=',i6/65x,'nodes per element(maximum)=',i6)2222format(//,10x,'the lengthe of real array is',i10,/,110x,'the lengthe of integer array is',i10)3333format(//,10x,'the lengthe of real array',i10,'exceed the',1'maximun value',i10)4444format(//,10x,'the lengthe of integer array',i10,'exceed the',1'maximun value',i10)cstopendccsubroutine pmesh(x,f,ix,id,ndf,ndm,nen1,nload)cc......Data input routine for mesh descriprioncimplicit real*8(a-h,o-z)dimension x(ndm,numnp),f(ndf,numnp),id(ndf,numnp),ix(nen1,numel)common/bdata/head(20)common/cdata/numnp,numel,nummat,nen,neqcommon/mater/ee,xnu,itypecommon/iofile/ior,iowcc.....Input constrain codes and nodal coordinate datacc id(k,j):constrain code of kth degree of freedom of node j,=0:free,=1:fixed c x(k,j):kth coordinate of node jcdo i=1,numnpread(ior,'(3i5,2f10.4)')j,(id(k,j),k=1,ndm),(x(k,j),k=1,ndm) end docwrite(iow,'(//17hnodal coordinates,/)')do i=1,numnpwrite(iow,'(3i5,2f10.4)')i,(id(k,i),k=1,ndm),(x(k,i),k=1,ndm) end docc.....element data inputcc ix(k,j):global node number of kth node in element jcdo i=1,numelread(ior,'(9i5)')j,(ix(k,j),k=1,nen)end docwrite(iow,'(//,18helement definition,/)')do i=1,numelwrite(iow,'(9i5)')j,(ix(k,j),k=1,nen)end docc.....Material data inputcc ee:young's modulus,xnu:poisson ratioc itype:type of problem,=1,:plane stress,=2:plane strain,=3:axi-symmetric cread(ior,'(2f10.4,i5)')ee,xnu,itypewrite(iow,'(//,19hmateial properties,/)')write(iow,'(2(e10.4,5x),i5)')ee,xnu,itypecc.....force/disp data inputcc f(k,j):concentrate load at node j in k directioncf=0.0d0do i=1,nloadread(ior,'(i5,2f10.4)')j,(f(k,j),k=1,ndf)end docwrite(iow,'(//,20happlied nodal forces,/)')do i=1,nloadwrite(iow,'(i5,2f10.4)')j,(f(k,j),k=1,ndf)end docreturncc format statementsc2000format('mesh1>',$)3000format(1x,'**warning**element connections necessary'1'to use block in macro program')4000format('**current problem valies**'/i6,'nodes,',1i5,'elmts,',i3,'matls,',i2,'dims,',i2,'dof/node,',2i3,'nodes/elmt')endccsubroutine assem(nad,ix,id,jd,x,b,ad,al,au)cc Call element subroutine and assemble global tangent matrixcimplicit real*8(a-h,o-z)dimension ilx(nen),xl(ndf,nen),ld(ndf,nen),s(nst,nst),p(nst)dimension ix(nen1,numel),id(ndf,numnp),jd(ndf*numnp)dimension x(ndm,numnp),b(neq),ad(neq),al(nad),au(nad)common/cdata/numnp,numel,nummat,nen,neqcommon/sdata/ndf,ndm,nen1,nstcnel=nencc elenment loopcdo320n=1,numels=0.0d0!element stiffness matrixp=0.0d0!nodal forcene=ndo310i=1,nenilx(i)=ix(i,ne)!current element definitiondo k=1,ndmxl(k,i)=x(k,ilx(i))!nodal coords in current elementend dokk=ilx(i)do k=1,ndfld(k,i)=id(k,kk)!equation numbersend do310continuecc Call element libccall elmt01(xl,ilx,s,p,ndf,ndm,nst)cc Asemmble tangent matrix and load vector if neededccall dasbly(ndf,nad,s,p,ld,jd,nst,b,ad,al,au)c320continue!end element loopcreturnendccsubroutine dasbly(ndf,nad,s,p,ld,jp,ns,b,ad,al,au)cc.....Assemble the symmetric or unsymmetric arrays for'dasol'cimplicit real*8(a-h,o-z)c logical alfl,aufl,bfldimension ad(neq),al(nad),au(nad)dimension ld(ns),jp(ndf*numnp),b(neq),s(ns,ns),p(ns)common/cdata/numnp,numel,nummat,nen,neqcommon/iofile/ior,iowcc alfl=true:for unsymmetric matirx assemblec alfl=false:for symmetric matirx assemblec s:element stiffness matrixc p:load or internal force vectorc ad:diagonal elementsc au:upper triangle elementsc al:lower triangle elementsc jp:pointer to last element in each row/column of al/au respectivec ld:equation numbers of each freedom degree in an element(get from id) cc.....Loop through the rows to perform the assemblycdo200i=1,nsii=ld(i)if(ii.gt.0)thenc if(aufl)then!assemble stiffness matrixcc.....Loop through the columns to perform the assemblycdo100j=1,nsif(ld(j).eq.ii)thenad(ii)=ad(ii)+s(i,j)elseif(ld(j).gt.ii)thenjc=ld(j)jj=ii+jp(jc)-jc+1au(jj)=au(jj)+s(i,j)c if(alfl)al(jj)=al(jj)+s(j,i)!unsymmetricendif100continueendifc if(bfl)b(ii)=b(ii)+p(i)!assemble nodal forcec endif200continuecreturnendccsubroutine dasol(ndf,numnp,b,jp,neq,nad,energy,ad,al,au)cc.....Solution of symmetric equations in profile formc.....Coeficient matrix must be decomposed into its triangularc.....Factor using datri beforce using dasol.cc jp:pointer to last element in each row/column of al/au respecive ccimplicit real*8(a-h,o-z)dimension ad(neq),al(nad),au(nad)dimension b(neq),jp(ndf*numnp)common/iofile/ior,iowdata zero/0.0d0/cc.....Find the first nonzero entry in the ring hand sidecdo is=1,neqif(b(is).ne.zero)go to200end dowrite(iow,2000)returnc200if(is.lt.neq)thencc.....Reduce the right hand sidecdo300j=is+1,neqjr=jp(j-1)jh=jp(j)-jrif(jh.gt.0)thenb(j)=b(j)-dot(al(jr+1),b(j-jh),jh)end if300continueend ifcc.....Multiply inverse of diagonal elementscenergy=zerodo400j=is,neqbd=b(j)b(j)=b(j)*ad(j)energy=energy+bd*b(j)400continuecc.....backsubstitutioncif(neq.gt.1)thendo500j=neq,2,-1jr=jp(j-1)jh=jp(j)-jrif(jh.gt.0)thencall saxpb(au(jr+1),b(j-jh),-b(j),jh,b(j-jh))end if500continueend ifcreturnc2000format('**dasol warning1**zero right-hand-side vector') endccsubroutine datest(au,jh,daval)cc.....test for rankcimplicit real*8(a-h,o-z)dimension au(jh)cdaval=0.0d0cdo j=1,jhdaval=daval+abs(au(j))end docreturnendccsubroutine datri(ndf,numnp,jp,neq,nad,flg,ad,al,au)cc.....Triangular decomposiontion of a matrix stored in profile form cimplicit real*8(a-h,o-z)logical flgdimension jp(ndf*numnp),ad(neq),al(nad),au(nad)common/iofile/ior,iowcc.....n.b.tol should be set to approximate half-word precision.cdata zero,one/0.0d0,1.0d0/,tol/0.5d-07/cc.....Set initial values for contditioning checkcdimx=zerodimn=zerocdo j=1,neqdimn=max(dimn,abs(ad(j)))end dodfig=zerocc.....Loop through the columns to perform the triangular decomposition cjd=1do200j=1,neqjr=jd+1jd=jp(j)jh=jd-jrif(jh.gt.0)thenis=j-jhie=j-1cc.....If diagonal is zeor compute a norm for singularity testcif(ad(j).eq.zero)call datest(au(jr),jh,daval)do100i=is,iejr=jr+1id=jp(i)ih=min(id-jp(i-1),i-is+1)if(ih.gt.0)thenjrh=jr-ihidh=id-ih+1au(jr)=au(jr)-dot(au(jrh),al(idh),ih)if(flg)al(jr)=al(jr)-dot(al(jrh),au(idh),ih)end if100continueend ifcc.....Reduce the diagonalcif(jh.ge.0)thendd=ad(j)jr=jd-jhjrh=j-jh-1call dredu(al(jr),au(jr),ad(jrh),jh+1,flg,ad(j))cc.....Check for possible errors and print warningscif(abs(ad(j)).lt.tol*abs(dd))write(iow,2000)jif(dd.lt.zero.and.ad(j).gt.zero)write(iow,2001)jif(dd.gt.zero.and.ad(j).lt.zero)write(iow,2001)jif(ad(j).eq.zero)write(iow,2002)jif(dd.eq.zero.and.jh.gt.0)thenif(abs(ad(j)).lt.tol*daval)write(iow,2003)jendifendifcc.....Stroe reciprocal of diagonal,compute condition checkscif(ad(j).ne.zero)thendimx=max(dimx,abs(ad(j)))dimn=min(dimn,abs(ad(j)))dfig=max(dfig,abs(dd/ad(j)))ad(j)=one/ad(j)end if200continuecc.....Print conditioning informationcdd=zeroif(dimn.ne.zero)dd=dimx/dimnifig=dlog10(dfig)+0.6write(iow,2004)dimx,dimn,dd,ifigcreturncc.....formatsc2000format('**datri warning1**loss of at least7digits in', 1'reducing diagonal of equation',i5)2001format('**datri warning2**sign of changed when', 1'reducing equation',i5)2002format('**datri warning3**reduced diagonal is zero zeri for', 1'equation',i5)2003format('**datri warning4**rank failure ffo zero unreduced', 1'diagonal in equation',i5)2004format(//'conditon check:d-max',e11.4,';d-min',e11.4, 1';ratio',e11.4/'maximim no.diagonal digits lost:',i3) 2005format('cond ck:dmax',1p1e9.2,';dmin',1p1e9.2,1';ratio',1p1e9.2)endccsubroutine dredu(al,au,ad,jh,flg,dj)cc.....Reduce diagonal element in triangular decompositioncimplicit real*8(a-h,o-z)logical flgdimension al(jh),au(jh),ad(jh)cdo j=1,jhud=au(j)*ad(j)dj=dj-al(j)*udau(j)=udend docc.....Finish computation of column of al for unsymmetric matricescif(flg)thendo j=1,jhal(j)=al(j)*ad(j)end doend ifcreturnendccsubroutine profil(jd,idl,id,ix,ndf,nen1,nad)cpute profile of global arrayscimplicit real*8(a-h,o-z)dimension jd(ndf*numnp),idl(numel*nen*ndf),id(ndf,numnp),1ix(nen1,numel)common/cdata/numnp,numel,nummat,nen,neqcommon/frdata/maxfcommon/iofile/ior,iowcc jd:column hight(address of diagonal elements)c id:boudary condition codes before this bubroutine's runningc id:equation numbers in global array(excluded restrained nodes)after running c idl:element strech orderc nad:total number of non-zero elements except diagonal elementsc in global tangent matrixcc.....Set up the equation numberscneq=0cdo10k=1,numnpdo10n=1,ndfj=id(n,k)if(j.eq.0)thenneq=neq+1id(n,k)=neqelseid(n,k)=0endif10continuecpute column heightsccall pconsi(jd,neq,0)cdo50n=1,numelmm=0nad=0do30i=1,nenii=iabs(ix(i,n))if(ii.gt.0)thendo20j=1,ndfjj=id(j,ii)if(jj.gt.0)thenif(mm.eq.0)mm=jjmm=min(mm,jj)nad=nad+1idl(nad)=jjendif20continueend if30continueif(nad.gt.0)thendo40i=1,nadii=idl(i)jj=jd(ii)jd(ii)=max(jj,ii-mm)40continueendif50continuecpute diagongal pointers for profilecnad=0jd(1)=0if(neq.gt.1)thendo60n=2,neqjd(n)=jd(n)+jd(n-1)60continuenad=jd(neq)end ifcc.....Set element search order to sequentialcdo70n=1,numelidl(n)=n70continuecc.....equation summarycmaxf=0mm=0if(neq.gt.0)mm=(nad+neq)/neqwrite(iow,2001)neq,numnp,mm,numel,nad,nummatcreturnc2001format(5x,'neq=',i5,5x,'numnp=',i5,5x,'mm=',i5,/5x, 1'numel=',i5,5x,'nad=',i5,5x,'nummat=',i5/) endcsubroutine saxpb(a,b,x,n,c)cc.....Vector times scalar added to second vectorcimplicit real*8(a-h,o-z)dimension a(n),b(n),c(n)cdo k=1,nc(k)=a(k)*x+b(k)end docreturnendcsubroutine pconsi(iv,nn,ic)cc.....Zero integer arraycdimension iv(nn)cdo n=1,nniv(n)=icend docreturnendcsubroutine elmt01(xl,ilx,s,p,ndf,ndm,nst)cc.....plane linear elastic element routinec ityp=1:plane stressc=2:plane strainc=3:axisymmetriccimplicit real*8(a-h,o-z)dimension xl(ndm,nen),ilx(nen),sigr(6)dimension d(18),s(nst,nst),p(nst),shp(3,9),sg(16),tg(16),wg(16)character wd(3)*12common/cdata/numnp,numel,nummat,nen,neqcommon/mater/ee,xnu,itypecommon/iofile/ior,iowdata wd/'plane stress','plane strain','axisymmetric'/cc xl(ndm,nen):coords of each node in current elementc ilx(nen):element definition of current elementc d(18):materials propertiesc s(nst,nst):element stiffness matrixc p(ns):nodal force and internal forcec shp(3,9):shape function and its derivativesc sg(16),tg(16),wg(16):weight coefficients of guass intergtation c l,k:integration pointscl=2k=2e=eenel=nencc d(14):thickness;d(11),d(12):body forcesc.....Set material patameter type and flagscityp=max(1,min(ityp,3))j=min(ityp,2)cd(1)=e*(1.+(1-j)*xnu)/(1.+xnu)/(1.-j*xnu)d(2)=xnu*d(1)/(1.+(1-j)*xnu)d(3)=e/2./(1.+xnu)d(13)=d(2)*(j-1)if((d(14).le.0.0d0).or.ityp.ge.2)d(14)=1.0d(15)=itypd(16)=ed(17)=xnud(18)=-xnu/el=min(4,max(1,l))k=min(4,max(1,k))d(5)=ld(6)=kc d(9)=t0c d(10)=e*alp/(1.-j*xnu)lint=0cwrite(iow,2000)wd(ityp),d(16),d(17),d(4),l,k,d(14),1d(11),d(12)cc.....stiffness/residual computationcl=kcc Compute cordinates and weights of integtation pointc`sg,tg:cootds;wg=wp*wqcif(l*l.ne.lint)call pguass(l,lint,sg,tg,wg)cpute integrals of shape functionscdo340l=1,lintcc Compute shape function and their derivatives to local and global coordinate systemccall shape(sg(l),tg(l),xl,shp,xsj,ndm,nen,ilx,.false.)cc Compute global coordinates of integration pointscxx=0.0yy=0.0do j=1,nenxx=xx+shp(3,j)*xl(1,j)yy=yy+shp(3,j)*xl(2,j)end doxsj=xsj*wg(l)*d(14)!xsj+|j|(sp,tq)*wp*wq*tcpute jacobian correction for plane stress and strain problemscif(ityp.le.2)thendv=xsjxsj=0.0zz=0.0c sigr4=-d(11)*dv!d(11)body forceelsecc For anisymmetric problemcdv=xsj*xx*3.1415926*2.zz=1./xxc sigr4=sigr(4)*xsj-d(11)*dvendifj1=1cc.....Loop over rowscdo330j=1,nelw11=shp(1,j)*dvw12=shp(2,j)*dvw22=shp(3,j)*xsjw22=shp(3,j)*dv*zzccpute the internal forces out of balancecc p(j1)=p(j1)-(shp(1,j)*sigr(1)+shp(2,j)*sigr(2))*dvc1-shp(3,j)*sigr4c p(j1+1)=p(j1+1)-(shp(1,j)*sigr(2)+shp(2,j)*sigr(3))*dvc1+d(12)*shp(3,j)*dv!d(12)body force cc.....Loop over columns(symmetry noted)c Compute stiffness matrixck1=j1a11=d(1)*w11+d(2)*w22a21=d(2)*w11+d(1)*w22a31=d(2)*(w11+w22)a41=d(3)*w12a12=d(2)*w12a32=d(1)*w12a42=d(3)*w11do320k=j,nelw11=shp(1,k)w12=shp(2,k)w22=shp(3,k)*zzs(j1,k1)=s(j1,k1)+w11*a11+w22*a21+w12*a41s(j1+1,k1)=s(j1+1,k1)+(w11+w22)*a12+w12*a42s(j1,k1+1)=s(j1,k1+1)+w12*a31+w11*a41s(j1+1,k1+1)=s(j1+1,k1+1)+w12*a32+w11*a42k1=k1+ndf320continuej1=ndf+j1330continue340continuecc.....Make stiffness symmetriccdo360j=1,nstdo360k=j,nsts(k,j)=s(j,k)360continuecreturncc.....Formats for input-outputc1000format(3f10.0,3i10)1001format(8f10.0)2000format(/5x,a12,'linear elastic element'//110x,'modulus',e18.5/10x,'poission ratio',f8.5/10x,'density',e18.5/ 210x,'guass ptr/dir',i3/10x,'stress pts',i6/10x,'thickness',e16.5/310x,'1-gravity',e16.5/10x,'2-gtavity',e16.5/10x,'alpha',e20.5/410x,'base temp',e16.5/)2001format(5x,'element stresses'//'elmt1-coord',2x,'11-stress',2x, 1'12-stress',2x,'22-stress',2x,'33-stress',3x,'1-coord',2x,3x,2'2-stress'/'matl2-coord',2x,'11-strain',2x,'12-strain'2x,3'22-strain',2x,'33-strain',6x,'angle'/39('-'))2002format(i4,0p1f9.3,1p6e11.3/i4,0p1f9.3,1p4e11.3,0p1f11.2/) 5000format('input:e,nu,rho,pts/stiff,pts/stre',1',type(1=stress,2=strain,3=axism)',/3x,'>',$)5001format('input:thickness,1-body force,1-body force,alpha,' 1,'temp-base'/3x,'>',$)endcsubroutine shape(ss,tt,xl,shp,xsj,ndm,nel,ilx,flg)cc.....Shape function routine for two dimension elementscimplicit real*8(a-h,o-z)logical flgdimension xl(ndm,nel),s(4),t(4),x(nel)dimension shp(3,nel),xs(2,2),sx(2,2)data s/-0.5d0,0.5d0,0.5d0,-0.5d0/,1t/-0.5d0,-0.5d0,0.5d0,0.5d0/cc.....Form4-node quatrilateral shape functionscc nel:nuber of nodes per elementcdo100i=1,4shp(3,i)=(0.5+s(i)*ss)*(0.5+t(i)*tt)shp(1,i)=s(i)*(0.5+t(i)*tt)shp(2,i)=t(i)*(0.5+s(i)*ss)100continuecc.....Form triangge bu adding their and fourth together for triangle element cif(nel.eq.3)thendo i=1,3shp(i,3)=shp(i,3)+shp(i,4)enddoendifcc.....Add quatratic terms if necessary for element with more than4nodes cif(nel.gt.4)call shap2(ss,tt,shp,ilx,nel)cc.....Construct jacobian and its inversecdo125i=1,2do125j=1,2xs(i,j)=0.0do120k=1,nelxs(i,j)=xs(i,j)+xl(i,k)*shp(j,k)120continue125continuecc xsj:determinate of jacob matrixcxsj=xs(1,1)*xs(2,2)-xs(1,2)*xs(2,1)cif(flg)returnc flg=false:form global derivativescif(xsj.le.0.0d0)xsj=1.0sx(1,1)=xs(2,2)/xsjsx(2,2)=xs(1,1)/xsjsx(1,2)=-xs(1,2)/xsjsx(2,1)=-xs(2,1)/xsjcc....Form global derivativescdo130i=1,neltp=shp(1,i)*sx(1,1)+shp(2,i)*sx(2,1)shp(2,i)=shp(1,i)*sx(1,2)+shp(2,i)*sx(2,2)shp(1,i)=tp130continuecreturnendcsubroutine shap2(s,t,shp,ilx,nel)cc....Add quadtatic function as necessarycimplicit real*8(a-h,o-z)dimension shp(3,9),ilx(nel)cs2=(1.-s*s)/2.t2=(1.-t*t)/2.do100i=5,9do100j=1,3shp(j,i)=0.0100continuecc.....Midsize nodes(serenipity)cif(ilx(5).eq.0)go to101shp(1,5)=-s*(1.-t)shp(2,5)=-s2shp(3,5)=s2*(1.-t)101if(nel.lt.6)go to107if(ilx(6).eq.0)go to102shp(1,6)=t2shp(2,6)=-t*(1.+s)shp(3,6)=t2*(1.+s)102if(nel.lt.7)go to107if(ilx(7).eq.0)go to103shp(1,7)=-s*(1.+t)shp(2,7)=s2shp(3,7)=s2*(1.+t)103if(nel.lt.8)go to107if(ilx(8).eq.0)go to104shp(1,8)=-t2shp(2,8)=-t*(1.-s)shp(3,8)=t2*(1.-s)cc.....Interior node(lagragian)c104if(nel.lt.9)go to107if(ilx(9).eq.0)go to107shp(1,9)=-4.*s*t2shp(2,9)=-4.*t*s2shp(3,9)=4.*s2*t2cc.....Correct edge nodes for interior node(lagrangian) cdo106j=1,3do105i=1,4105shp(j,i)=shp(j,i)-0.25*shp(j,9)do106i=5,8106if(ilx(i).ne.0)shp(j,i)=shp(j,i)-.5*shp(j,9)cc.....Correct corner nodes for presense of midsize nodes c107do108i=1,4k=mod(i+2,4)+5l=i+4do108j=1,3108shp(j,i)=shp(j,i)-0.5*(shp(j,k)+shp(j,l))returnendcsubroutine pguass(l,lint,r,z,w)cc.....Guass points and weights for two dimensionscimplicit real*8(a-h,o-z)dimension lr(9),lz(9),lw(9),r(16),z(16),w(16)c common/eldtat/dm,n,ma,mct,iel,neldata lr/-1,1,1,-1,0,1,0,-1,0/,lz/-1,-1,1,1,-1,0,1,0,0/data lw/4*25,4*40,64/cc lint:number of integration pointsc r,z:coordinates of integration pointsc w:wp*wq,product of the two weightsclint=l*lcc.....1x1integerationc1r(1)=0.z(1)=0.w(1)=4.creturncc.....2x2integerationc2g=1.0/sqrt(3.d0)do i=1,4r(i)=g*lr(i)z(i)=g*lz(i)w(i)=1.end docreturncc.....3x3integerationc3g=sqrt(0.60d0)h=1.0/81.0d0cdo i=1,9r(i)=g*lr(i)z(i)=g*lz(i)w(i)=h*lw(i)enddocreturncendcsubroutine pload(id,f,b,nneq,neq) cc.....Form load vector in compact formcimplicit real*8(a-h,o-z)dimension f(nneq),b(neq),id(nneq)common/iofile/ior,iowcb=0.0d0cj=id(n)if(j.gt.0)thenb(j)=f(n)endifenddocreturnendcsubroutine prtdis(id,b,ndf,numnp,neq)cc Print out nodal displacementscimplicit real*8(a-h,o-z)dimension id(ndf,numnp),b(neq),u(ndf,numnp)common/iofile/ior,iowcu=0.0d0do100i=1,numnpdo j=1,ndfn=id(j,i)if(n>0)u(j,i)=b(n)end do100continuecc Out nodal displacementscwrite(iow,'(//,19hnodal displacements,/)')do i=1,numnpwrite(iow,'(5x,i5,2x,3(e12.4,3x))')i,(u(k,i),k=1,ndf) end docreturnendcdouble precision function dot(a,b,n)implicit real*8(a-h,o-z)dimension a(n),b(n)cc.....Dot product functioncdot=0.0d0do10k=1,ndot=dot+a(k)*b(k)10continuereturn end。

八节点平面等参元Fortran源程序C 这是一个采用平面四边形八节点等参单元的平面有限元分析程序。

PROGRAM PLANEFEMIMPLICIT REAL*8(A-H,O-Z)CHARACTER*80 LINECHAR,NEWLINECHARCOMMON A(30000),L(4000)COMMON /SOL/NPOIN,NELEM,NTYPE,NMATSOPEN(5,FILE='FEMDATA',STATUS='OLD')OPEN(6,FILE='FEMOUT',STATUS='UNKNOWN')C 以下进行的是从数据文件FEMDATA中读入数据文件主标题信息。

READ(5,5000) LINECHARLOCATECHAR=INDEX(LINECHAR,'输入')IF(LOCATECHAR.NE.0) THENLINECHAR(LOCATECHAR:LOCATECHAR+3)='输出'ENDIFWRITE(6,5100) LINECHAR5000 FORMAT(A)5100 FORMAT(80('*')/A/80('*')/)C 以下进行的是先读入一行字符，然后从字符中读入网格单元总数NELEM。

READ(5,5000) LINECHARWRITE(6,5000) LINECHARLOCATECHAR=INDEX(LINECHAR,'=')LOCATECHAR=LOCATECHAR+1NEWLINECHAR=LINECHAR(LOCATECHAR:80)READ(NEWLINECHAR,5200) NELEM5200 FORMAT(I5)C 以下进行的是先读入一行字符，然后从字符中读入单元节点总数NPOIN。

READ(5,5000) LINECHARWRITE(6,5000) LINECHARLOCATECHAR=INDEX(LINECHAR,'=')LOCATECHAR=LOCATECHAR+1NEWLINECHAR=LINECHAR(LOCATECHAR:80)READ(NEWLINECHAR,5200) NPOINC 以下进行的是先读入一行字符，然后从字符中读入问题类型编号NTYPE。

近场动力学与有限单元法的混合模型与隐式求解格式郁杨天;章青;顾鑫【摘要】利用近场动力学方法(PD)在模拟不连续变形问题的独特优势和有限单元法(FEM)的计算效率,提出近场动力学与有限单元法混合建模的方法,并用于求解断裂力学问题.裂纹出现的区域采用改进的近场动力学微观弹脆性(PMB)模型进行离散,其他区域采用有限元离散,通过杆单元连接不同的子区域.在隐式求解体系下实现了两种方法的混合建模,该模型在求解静力学问题时无需引入阻尼项,有效提高了计算效率和计算精度.通过模拟计算简支梁的弹性变形和三点弯曲梁I型裂纹的扩展过程,与理论解吻合良好,验证了提出的混合模型和求解方法的准确性和有效性.%A hybrid model of peridynamics (PD) and finite element method (FEM) was proposed and applied to solve problems of fracture mechanics in order to combine the unique advantage of PD in solving discontinuities and the computational efficiency of FEM.The improved prototype microelastic brittle (PMB) model of peridynamics was utilized for the regions where material failure was expected.The region without failure was discretized by FEM.The truss element was introduced to bridge peridynamic subregions and finite element subregions.The hybrid model is based on the implicit schemes, and it need not consider a fictitious damping term in solving static problems.The computational efficiency and accuracy of the model were improved.The static elastic deformation of a simply supported beam and the propagation process of mode I fracture in a three points bend beam were simulated to verify the accuracy and utility of the presentedmodel.Results obtained by the model agreed well with the theoretical solutions.【期刊名称】《浙江大学学报（工学版）》【年(卷),期】2017(051)007【总页数】7页(P1324-1330)【关键词】近场动力学;有限单元法(FEM);混合模型;改进的近场动力学微观弹脆性模型;隐式求解【作者】郁杨天;章青;顾鑫【作者单位】河海大学工程力学系,江苏南京 210098;河海大学工程力学系,江苏南京 210098;河海大学工程力学系,江苏南京 210098【正文语种】中文【中图分类】TU375;TU311有限单元法已在工程中得到了广泛应用,但在求解裂纹扩展等不连续问题时,通常须借助断裂准则以判断开裂位置及扩展方向,在裂纹扩展后要不断进行网格重构,计算结果具有网格依赖性[1].设置黏结单元[2-3]虽然可以避免计算过程的复杂性,但必须预先判断裂纹扩展方向和所在区域.有限单元法在求解断裂破坏问题时所面临的局限性源自于连续性假设的理论基础与实际问题不连续的矛盾.近场动力学[4-6]作为一种新兴的基于非局部积分思想的数值方法,避免了传统的局部微分方程求解不连续问题时的奇异性,在破坏问题分析中具有独特的优势,取得了许多新的成果[7-14],已成为国际学术界研究的热点,在岩石破裂和土体破坏分析等方面得到了成功应用[15-19].近场动力学作为一种粒子类方法,计算效率较低,为了充分利用有限元与近场动力学各自的优势,国外一些学者对近场动力学与有限元的混合建模方法进行研究.Macek 等[20]将近场动力学中的物质点对用有限元中的杆单元来描述,采用有限元分析软件ABAQUS中的镶嵌单元(embedded element),实现了近场动力学与有限单元法的混合建模；为了避免镶嵌单元刚度过大,镶嵌单元的弹性模量和密度必须取很小的数值.Kilic等[21]通过设置重叠区域来实现近场动力学与有限元的耦合,在重叠区内,采用有限元结点位移插值的方法确定近场动力学物质点的位移,但无法通过物质点的位移求得有限元的结点位移.Liu等[22]通过引入界面单元实现了两种方法的耦合,给出两种耦合方案,但在界面单元处,由近场动力学物质点的受力情况计算界面单元结点力的过程较繁琐.此外,上述混合建模方法均建立在显式求解体系下,在求解静力学问题时必须引入阻尼项[21],而阻尼的大小直接影响计算的收敛速度[23],致使计算效率不高.本文建立近场动力学与有限单元法的混合模型,并应用于断裂力学问题.对可能出现破坏的区域,采用近场动力学方法离散,其他区域采用有限单元法离散.通过杆单元连接PD和FEM子域,在隐式求解体系下实现了两种方法的混合建模.该模型在求解静力学问题时无需引入阻尼项,提高了计算效率.此外,在PD子域,采用改进的微观弹脆性(PMB)模型,提高了计算精度.基于所建立的混合模型,模拟计算了简支梁的弹性变形和三点弯曲梁I型裂纹的扩展过程,验证了该方法的准确性和有效性.1.1 近场动力学方法如图1所示,假设在某一时刻t,空间域R中任一物质点x与其邻近一定范围内的其他物质点存在相互作用力f.根据牛顿第二定律,可得ρ (x,t)=f(u(x′,t)-u(x,t),x′-x)dVx′+b(x,t).式中：Hx为物质点x的近场范围,δ为近场范围尺寸,ρ为物质点的材料密度,u为物质点的位移,b为外荷载密度.物质点间的相互作用力函数f称为本构力函数.它包含了材料的本构信息,一般根据材料特性(如弹性、黏弹性、弹塑性、均匀各向同性、各向异性等)与长程作用力空间分布特征构造本构力函数[24-26].微观弹脆性(prototype microelastic brittle, PMB)材料[27]的本构力函数f为式中:ξ为相对位置,ξ=x′-x；η为相对位移,η=u′-u;c(ξ,δ)为微观模量函数；s为物质点对的相对伸长率,μ为一标量函数,其中,s0为临界伸长率,可以通过材料的抗拉强度ft和弹性模量E定义,当s超过s0时,点对间不再发生作用.在近场动力学理论中,定义标量函数φ(x,t)以反映物质点x的损伤：式中：0≤φ≤1,φ=0表示材料未损伤,φ=1表示该点完全损伤,不再与其他点发生相互作用.1.2 改进的PMB模型式(2)中的微观模量函数c(ξ,δ)可以表示为式中：g(ξ,δ)为核函数,反映物质点对间长程力的强度随两点间距离变化的规律. 在以往的PMB模型中,核函数通常取将物质点对的微观模量设为常数,忽略物质点对间距离对微观模量函数的影响,从而导致应用PMB模型计算弹性变形时存在较大的误差.在改进的PMB模型[23,28]中,取核函数的表达式为该核函数能够反映长程力的空间分布规律,计算精度高.根据近场动力学应变能密度与连续介质力学应变能密度相等的原则,可以导出改进的PMB模型中c(0,δ)的表达式：1.3 有限元方程有限元分析的支配方程为式中：K为整体劲度矩阵,u为整体位移列阵,F为整体等效结点荷载列阵,其中,为选择矩阵;k为单元劲度矩阵,其中B为应变转化矩阵,D为弹性矩阵；其中N 为形函数矩阵,r为体积力,为物体表面的面力.2.1 近场动力学隐式求解格式将物体均匀离散,取物质点间距为运动方程(1)的离散形式可以表示为式中：n为时间步长编号；Vj为j处物质点体积,在三维情况下在二维情况下对于静力问题,令则近场动力学平衡方程的离散形式为式中：n表示第n级荷载.将PMB模型的本构力函数式(2)写成矩阵形式：式中：KP为微观模量函数的矩阵形式.建立局部直角坐标系x′y′z′,并使x′轴与物质点对的长度方向平行,如图2所示.在局部坐标系下,物质点对i、j的作用力和位移分别为存在如下的关系：式中：对局部坐标系下的物质点对的微观模量矩阵Kp′进行坐标转换.令物质点对的长度方向(x′轴)与整体坐标x、y和z轴的余弦分别为l、m和n,如图3所示.由图3可知,通过坐标转换,可得整体坐标系下物质点对微观模量函数的矩阵形式：于是,式(14)可以写成可得与有限元方法求解静力问题形式相同的近场动力学求解方程.2.2 近场动力学与有限元的混合模型如图4所示,将所考虑的物体划分为近场动力学子域和有限元子域进行离散.在两种区域的交界面上,采用杆单元连接近场动力学中的物质点与有限元中的结点,实现近场动力学与有限单元法的混合建模.如图5所示,交界面上的有限元结点不仅与所在单元的其他结点发生作用,而且通过杆单元与其近场范围δ内的物质点相互作用.参照式(9)、(10),杆单元的劲度可以取为式中：ξ′为有限元结点与物质点的相对距离.在整体坐标系下,杆单元劲度矩阵的形式与式(20)相同.具体计算时,只需用杆单元劲度k和有限元结点与物质点的相对距离ξ′替换式(20)中的微观模量c和物质点对间距离‖ξ‖,再形成总体的整体劲度矩阵和荷载列阵进行求解.针对上述计算方法,采用Fortran语言编写了计算程序,主要的计算步骤如图6所示.3.1 简支梁的弹性变形考虑图7所示的简支梁,在梁的上表面中点受集中力P=10 kN的作用,材料的弹性模量为30 GPa,泊松比为0.3.在计算中,将梁划分为两个有限元子域和一个近场动力学子域,取物质点间距Δx=2.5 mm,近场范围δ=4Δx=10 mm,有限元网格采用10 mm×10 mm的四结点矩形单元.为了分析比较,对简支梁全部采用近场动力学模型和有限元模型进行计算.其中有限元采用ANSYS软件,单元类型为SOLID四边形单元,尺寸为10 mm×10 mm.如表1所示为采用不同方法得到的中性轴跨中挠度计算结果.计算机型号为Intel Core i3 CPU @2.53 GHz,内存为4.0 GB.表中，hm为跨中挠度，er为相对误差，tc为计算耗时.图8给出利用混合模型计算得到的简支梁中性轴竖向位移h与解析解、有限单元法、近场动力学法计算结果的对比.图中，L为中性轴位置.可见,采用混合模型方法计算得到的中性轴跨中挠度相对误差为0.31%,与有限元解的相对误差0.14%相当,并且计算耗时仅为近场动力学用时的一半,验证了提出的近场动力学与有限元混合模型的正确性和有效性. 在弹性问题的计算中,虽然近场动力学方法的计算精度和计算效率均低于有限元法,但近场动力学方法主要是用于有限元难以求解的材料和结构的破坏问题.采用提出的近场动力学与有限元混合模型,精度与有限元相当,计算效率得到了很大的提高,为近场动力学方法应用于工程结构破坏问题的计算分析提供了一种有效的途径.3.2 含初始裂缝三点弯曲梁的破坏分析考虑图9所示的三点弯曲梁,跨中位置设置一预制裂缝,初始裂缝长度a=80 mm,材料的弹性模量为30 GPa,泊松比为0.3,抗拉强度为1.65 MPa.计算中,将梁划分为两个有限元子域和一个近场动力学子域,取Δx=2.5 mm,δ=4Δx=10 mm,有限元单元采用大小为10 mm×10 mm的四结点矩形单元,加载方式为位移加载.三点弯曲梁I型裂纹张开口位移(crack-mouth-opening displacement, CMOD)的线弹性断裂力学(LEFM)的解答[30-31]为式中：图10给出计算得到的荷载-CMOD曲线,与线弹性断裂力学的结果有较好的一致性.图11给出不同时刻梁的破坏特征.从计算结果可以发现,荷载达到6.4 kN时结构出现损伤,认为裂纹开始稳定扩展；荷载继续增加,当荷载达到7.5 kN时,裂纹开始失稳扩展.随着裂纹的不断扩展,所需荷载逐渐减小.本文提出近场动力学与有限元混合建模的方法求解材料和结构的破坏问题.该方法在可能出现破坏的区域,采用改进的近场动力学PMB模型离散,其他区域采用有限元模型离散,通过杆单元来连接各子区域.在隐式求解体系下实现了近场动力学与有限元的混合建模和数值求解.在求解静力学问题时,无需引入阻尼项,提高了计算效率和计算精度.简支梁的弹性变形和三点弯曲梁I型裂纹的扩展过程的计算结果表明,采用提出的PD-FEM混合模型和静力求解方法能够模拟裂纹扩展问题,减少近场动力学的计算域,节约计算时间,为解决工程结构破坏问题提供了一种有效的分析方法. 近场动力学突破了传统连续介质力学方法求解破坏问题时的奇异性和困难,在各类材料和结构的静动力裂纹扩展、冲击破坏等许多力学问题中均得到成功应用.近场动力学处于理论体系完善和初步应用阶段,还有很大的研究空间,在基于状态的近场动力学理论和方法、复杂环境下本构力函数的构建、非均匀离散技术、高效的数值求解体系、软件研制和工程应用等方面亟需加大研究力度.【相关文献】[1] MUROTANI K, YAGAWA G, CHOI J B. Adaptivefinite elements using hierarchical mesh and its application to crack propagation analysis [J]. Computer Methods in Applied Mechanics and Engineering, 2013, 253: 1-14.[2] CAMACHO G T, ORTIZ M. Computational modeling of impact damage in brittle materials [J]. International Journal of Solids and Structures, 1996, 33(20):2899-2938. [3] XU X P, NEEDLEMAN A. Numerical simulations of fast crack growth in brittle solids [J]. Journal of the Mechanics and Physics of Solids, 1994, 42(9):1397-1434.[4] SILLING S A. Reformulation of elasticity theory for discontinuities and long-range forces [J]. Journal of the Mechanics and Physics of Solids, 2000, 48(1):175-209.[5] SILLING S A, EPTON M, WECKNER O, et al. Peridynamic states and constitutive modeling [J]. Journal of Elasticity, 2007, 88(2): 151-184.[6] SILLING S A. Linearized theory of peridynamic states [J]. Journal of Elasticity, 2010,99(1): 85-111.[7] SILLING S A, BOBARU F. Peridynamic modeling of membranes and fibers [J]. International Journal of Nonlinear Mechanics, 2005, 40(2): 395-409.[8] GERSTLE W, SAU N, SILLING S. Peridynamic modeling of concrete structures [J]. Nuclear Engineering and Design, 2007, 237(12): 1250-1258.[9] ASKARI E, BOBARU F, LEHOUCQ R B, et al. Peridynamics for multiscale materials modeling [C]∥Journal of Physics: Conference Series. [S. l.]: IOP Publishing, 2008, 125(1): 012078.[10] XU J, ASKARI A, WECKNER O, et al. Peridynamic analysis of impact damage in composite laminates [J]. Journal of Aerospace Engineering, 2008, 21(3):187-194.[11] 沈峰,章青,黄丹,等.冲击荷载作用下混凝土结构破坏过程的近场动力学模拟[J].工程力学,2012(增1): 12-15. SHEN Feng, ZHANG Qing, HUANG Dan, et al. Peridynamics modeling of failure process of concrete structure subjected to impact loading [J]. Engineering Mechanics, 2012(supple.1): 12-15.[12] 胡祎乐,余音,汪海.基于近场动力学理论的层压板损伤分析方法[J].力学学报,2013,45(4): 624-628. HU Yi-le, YU Yin, WANG Hai. Damage analysis method for laminates based onperidynamic theory [J]. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(4): 624-628.[13] 谷新保,周小平. 含圆孔拉伸板的近场动力学数值模拟[J]. 固体力学学报, 2015, 36(5): 376-383. GU Xin-bao, ZHOU Xiao-ping. The numerical simulation of tensile plate with circular hole using peridynamic theory [J]. Chinese Journal of Solid Mechanics, 2015, 36(5): 376-383. [14] 郁杨天,章青,顾鑫.含单边缺口混凝土梁破坏的近场动力学模拟[J].工程力学,2016,33(12):80-85. YV Yang-tian, ZHANG Qing, GU Xin. Impact failure simulation of a single-edged notched concrete beam based on peridynamics [J]. Engineering Mechanics, 2016, 33(12): 80-85. [15] LAI X, REN B, FAN H, et al. Peridynamics simulations of geomaterial fragmentation by impulse loads [J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2015, 39(12): 1304-1330.[16] LAI X, LIU L S, LIU Q W, et al. Slope stabilityanalysis by peridynamic theory[C]∥Applied Mechanics and Materials. [S. l.]: Trans Tech Publications, 2015, 744: 584-588.[17] ZHOU X P, GU X B, WANG Y T. Numerical simulations of propagation, bifurcation and coalescence of cracks in rocks [J]. International Journal of Rock Mechanics and Mining Sciences, 2015, 80: 241-254.[18] REN B, FAN H, BERGEL G L, et al. A peridynamics-SPH coupling approach to simulate soil fragmentation induced by shock waves [J]. Computational Mechanics, 2015, 55(2): 287-302.[19] FAN H, BERGEL G L, LI S. A hybrid peridynamics-SPH simulation of soil fragmentation by blast loads of buried explosive [J]. International Journal of Impact Engineering, 2016, 87: 14-27.[20] MACEK R W, SILLING S A. Peridynamics via finite element analysis [J]. Finite Elements in Analysis and Design, 2007, 43(15): 1169-1178.[21] KILIC B, MADENCI E. Coupling of peridynamic theory and the finite element method [J]. Journal of Mechanics of Materials and Structures, 2010, 5(5):707-733.[22] LIU W, HONG J W. A coupling approach of discretized peridynamics with finite element method [J]. Computer Methods in Applied Mechanics and Engineering, 2012, 245: 163-175.[23] HUANG D, LU G, QIAO P. An improved peridynamic approach for quasi-static elastic deformation and brittle fracture analysis [J]. International Journal of Mechanical Sciences, 2015, 94: 111-122.[24] AZIZI M A B. The peridynamic model of viscoelastic creep and recovery [J]. Multidiscipline Modeling in Materials and Structures, 2015, 11(4):579-597.[25] MADENCI E, OTERKUS S. Ordinary state-based peridynamics for plastic deformation according to von Mises yield criteria with isotropic hardening [J]. Journal of the Mechanics and Physics of Solids, 2015, 86: 192-219.[26] GHAJARI M, IANNUCCI L, CURTIS P. A peridynamic material model for the analysis of dynamic crack propagation in orthotropic media [J]. Computer Methods in Applied Mechanics and Engineering, 2014,276(7): 431-452.[27] SILLING S A, ASKARI E. A meshfree method based on the peridynamic model of solid mechanics [J]. Computers and Structures, 2005, 83(17): 1526-1535.[28] HUANG D, LU G, WANG C, et al. An extended peridynamic approach for deformation and fractureanalysis [J]. Engineering Fracture Mechanics, 2015, 141: 196-211.[29] 铁摩辛柯.弹性理论[M].北京：高等教育出版社,1990: 100-107.[30] JENQ Y S, SHAH S P. Mixed-mode fracture of concrete [J]. International Journal of Fracture, 1988,38(2): 123-142.[31] JOHN R, SHAH S P. Mixed-mode fracture of concrete subjected to impact loading [J]. Journal of Structural Engineering, 1990, 116(3): 585-602.。

程序框图：程序特点：问题类型：可用于计算结构力学的平面刚架问题单元类型：直接利用杆单元载荷类型：节点载荷及非节点载荷，其中非节点载荷包括均布荷载和垂直于杆件的集中荷载材料性质：所有杆单元几何性质相同，且由相同的均匀材料组成方程求解：结构刚度矩阵采用满阵存放，Gauss消元过程采用《数值分析》中的列主元素消去法输入文件：按先处理法的要求，由手工生成输入数据文件1.主要变量：ne: 单元个数nj: 结点个数n: 自由度e: 弹性模量（单位：KN/m2）a: 杆截面积zi: 惯性矩np: 结点荷载个数nf: 非结点荷载个数x(nj): 存放结点的x轴坐标y(nj): 存放结点的y轴坐标ij(ne,2): 存放单元结点编号，其中ij(nj,1)存放起始结点编号，ij(nj,2)存放终止结点编号jn(nj,3): 存放结点位移编号，以组成单元定位数组pj(np,3): 存放结点荷载信息，其中pj(np,1)存放结点荷载作用结点号，pj(np,2)存放荷载方向代码(1—x方向；2—y方向；3—转角)，pj(np,3)存放荷载大小pf(ne,4): 存放非结点荷载信息，其中pf(ne,1)存放荷载作用单元号，pf(ne,2)存放荷载代码(1—均布荷载，2—垂直集中荷载)，pf(ne,3)存放荷载大小，pf(ne,4)荷载作用距离（均布荷载，集中荷载均以单元起始结点为计算起始位置）。

2.子例行子程序哑元信息：第一部分：基本部分I. subroutine lsc(Length & Sin & Cos):输入哑元：m(单元号),nj,ne,x,y,ij输出哑元：bl(杆件长度),si(正弦值),co(余弦值)II. subroutine elv(Element Location Vector):输入哑元：m,ne,nj,ij,jn输出哑元：lv(单元定位数组)III. subroutine esm(Element Stiffness Matrix):输入哑元：e,a,zi,bl,si,co输出哑元：ek(整体坐标系下的单刚矩阵)IV. subroutine eff(Element Fixed-end Forces)输入哑元：i,pf,nf,bl输出哑元：fo(局部坐标系下单元固端力)第二部分：主程序直接调用部分I. subroutine tsm(Total Stiffness Matrix 计算总刚矩阵)输入哑元：ne,nj,n,e,x,y,ij,a,zi,jn输出哑元：tkII. subroutine jlp(Joint Load Vector 计算结点荷载)输入哑元：ne,nj,n,np,nf,x,y,ij,jn,pj,pf输出哑元：p(结点荷载列矩阵)III. subroutine gauss(带列主元素消去的高斯法)输入(输出)哑元：tk,p,n ;(注意，算出位移后，直接存储到结点荷载列矩阵)IV. subroutine mvn(Member-end forces of elements 计算各单元的杆端力)输入哑元：ne,nj,n,nf,e,x,y,ij,a,zi,jn,pf,p3.文件管理：源程序文件：pff.for程序需读入的数据文件：input.txt程序输出的数据文件：output4.数据文件格式：【输出文件格式】： 1. 第1部分： 每行数据依次为：结点号，结点x 方向位移，结点y 方向位移，结点转角位移 2. 第2部分：每行数据依次为：单元号，xi F ，yi F ，i M ，xj F ，yj F ，j M源程序：program PFF implicit nonereal tk(100,100),x(50),y(50),p(100),pj(50,3),pf(50,4) integer ij(50,2),jn(50,3) integer ne,nj,n,np,nf real e,a,ziopen(1,file="input.txt",status="old") open(2,file="output.txt",status="old")read(1,*) ne,nj,n,e,a,zi,np,nfcall input(ne,nj,x,y,ij,jn,np,nf,pj,pf)call tsm(ne,nj,n,e,x,y,ij,a,zi,jn,tk)call jlp(ne,nj,n,np,nf,x,y,ij,jn,pj,pf,p)call gauss(tk,p,n)call mvn(ne,nj,n,nf,e,x,y,ij,a,zi,jn,pf,p)endsubroutine input(ne,nj,x,y,ij,jn,np,nf,pj,pf)dimension x(nj),y(nj),ij(ne,2),jn(nj,3),pj(np,3),pf(nf,4) read(1,*)(x(i),y(i),i=1,nj)read(1,*)(ij(i,1),ij(i,2),i=1,ne)read(1,*)((jn(i,j),j=1,3),i=1,nj)if (np>0) read(1,*)((pj(i,j),j=1,3),i=1,np)if (nf>0) read(1,*)((pf(i,j),j=1,4),i=1,nf)endsubroutine tsm(ne,nj,n,e,x,y,ij,a,zi,jn,tk)dimension x(nj),y(nj),ij(ne,2),jn(nj,3),tk(n,n),ek(6,6),lv(6) do i=1,ndo j=1,ntk(i,j)=0enddoenddodo m=1,necall lsc(m,ne,nj,x,y,ij,bl,si,co)call esm(e,a,zi,bl,si,co,ek)call elv(m,ne,nj,ij,jn,lv)do l=1,6i=lv(l)if (i/=0) thendo k=1,6j=lv(k)if (j/=0) tk(i,j)=tk(i,j)+ek(l,k)enddoendifenddoenddoendsubroutine lsc(m,ne,nj,x,y,ij,bl,si,co) dimension x(nj),y(nj),ij(ne,2)i=ij(m,1)j=ij(m,2)dx=x(j)-x(i)dy=y(j)-y(i)bl=sqrt(dx*dx+dy*dy)si=dy/blco=dx/blendsubroutine esm(e,a,zi,bl,si,co,ek) dimension ek(6,6)c1=e*a/blc2=2.0*e*zi/blc3=3.0*c2/blc4=2.0*c3/bls1=c1*co*co+c4*si*sis2=(c1-c4)*si*cos3=c3*sis4=c1*si*si+c4*co*cos5=c3*cos6=c2ek(1,1)=s1ek(1,2)=s2ek(1,3)=-s3ek(1,4)=-s1ek(1,5)=-s2ek(1,6)=-s3ek(2,2)=s4ek(2,3)=s5ek(2,4)=-s2ek(2,5)=-s4ek(2,6)=s5ek(3,3)=2*s6ek(3,4)=s3ek(3,5)=-s5ek(3,6)=s6ek(4,4)=s1ek(4,5)=s2ek(4,6)=s3ek(5,5)=s4ek(5,6)=-s5ek(6,6)=2.0*s6do i=1,5do j=i+1,6ek(j,i)=ek(i,j)enddoenddoendsubroutine elv(m,ne,nj,ij,jn,lv)dimension ij(ne,2),jn(nj,3),lv(6)i=ij(m,1)j=ij(m,2)do k=1,3lv(k)=jn(i,k)lv(k+3)=jn(j,k)enddoendsubroutine jlp(ne,nj,n,np,nf,x,y,ij,jn,pj,pf,p)dimension x(nj),y(nj),ij(ne,2),jn(nj,3),pj(np,3),pf(nf,4),p(n),fo(6),pe(6),lv(6) do i=1,np(i)=0.0enddoif (np>0) thendo i=1,npj=int(pj(i,1))k=int(pj(i,2))l=jn(j,k)if (l/=0) p(l)=pj(i,3)enddoendifif(nf>0) thendo i=1,nfm=int(pf(i,1))call lsc(m,ne,nj,x,y,ij,bl,si,co)call eff(i,pf,nf,bl,fo)call elv(m,ne,nj,ij,jn,lv)pe(1)=-fo(1)*co+fo(2)*sipe(2)=-fo(1)*si-fo(2)*cope(3)=-fo(3)pe(4)=-fo(4)*co+fo(5)*sipe(5)=-fo(4)*si-fo(5)*cope(6)=-fo(6)do j=1,6l=lv(j)if (l/=0) p(l)=p(l)+pe(j) enddoenddoendifendsubroutine eff(i,pf,nf,bl,fo) dimension pf(nf,4),fo(6)no=int(pf(i,2))q=pf(i,3)c=pf(i,4)b=bl-cc1=c/blc2=c1*c1c3=c1*c2do j=1,6fo(j)=0.0enddogoto(10,20),no10 fo(2)=-q*c*(1.0-c2+c3/2.0)fo(3)=-q*c*c*(0.5-2.0*c1/3.0+0.25*c2) fo(5)=-q*c*c2*(1.0-0.5*c1)fo(6)=q*c*c*c1*(1.0/3.0-0.25*c1) return20 fo(2)=-q*b*b*(1.0+2.0*c1)/bl/blfo(3)=-q*c*b*b/bl/blfo(5)=-q*c2*(1.0+2.0*b/bl)fo(6)=q*c2*breturnendsubroutine gauss(e,d,n)dimension e(n,n),d(n),a(n,n+1)do i=1,ndo j=1,na(i,j)=e(i,j)enddoenddodo i=1,na(i,n+1)=d(i)enddodo k=1,n-1do i=k+1,nif (abs(a(i,k))>abs(a(k,k))) thendo j=1,n+1c=a(k,j)a(k,j)=a(i,j)a(i,j)=cenddoelseendifenddodo i=k+1,na(i,k)=a(i,k)/a(k,k)do j=k+1,n+1a(i,j)=a(i,j)-a(i,k)*a(k,j)enddoenddoenddoa(n,n+1)=a(n,n+1)/a(n,n)do i=n-1,1,-1do j=i+1,np=p+a(i,j)*a(j,n+1)enddoa(i,n+1)=(a(i,n+1)-p)/a(i,i)p=0enddodo i=1,nd(i)=a(i,n+1)enddoendsubroutine mvn(ne,nj,n,nf,e,x,y,ij,a,zi,jn,pf,p)dimension x(nj),y(nj),ij(ne,2),jn(nj,3),pf(nf,4),lv(6),fo(6),d(6),fd(6),f(6),ek(6,6),p(n) write(2,10)10 format(//2x,"结点位移"/5x,"结点号",9x,"u向位移",9x,"v向位移",9x,"角位移") do j=1,njdo i=1,3d(i)=0.0l=jn(j,i)if (l/=0) d(i)=p(l)enddowrite(2,20)j,d(1),d(2),d(3)20 format(2x,i6,4x,3e15.6)enddowrite(2,30)30 format(/2x,"单元杆端力及弯矩"/4x,"单元号",13x,"Fx",17x,"Fy",17x,"弯矩") do m=1,necall lsc(m,ne,nj,x,y,ij,bl,si,co)call esm(e,a,zi,bl,si,co,ek)call elv(m,ne,nj,ij,jn,lv)do i=1,6l=lv(i)d(i)=0.0if(l/=0) d(i)=p(l)enddodo i=1,6fd(i)=0.0do j=1,6fd(i)=fd(i)+ek(i,j)*d(j)enddoenddof(1)=fd(1)*co+fd(2)*sif(2)=-fd(1)*si+fd(2)*cof(3)=fd(3)f(4)=fd(4)*co+fd(5)*sif(5)=-fd(4)*si+fd(5)*cof(6)=fd(6)if (nf>0) thendo i=1,nfl=int(pf(i,1))if (m==l) thencall eff(i,pf,nf,bl,fo)do j=1,6f(j)=f(j)+fo(j)enddoendifenddoendifwrite(2,40)m,f40format(2x,i8,4x,"Ix=",f12.4,3x,"Iy=",f12.4,3x,"Mi=",f12.4/14x,"Jx=",f12.4,3x,"J y=",f12.4,3x,"Mj=",f12.4)enddoend【算例】：课题二：平面刚架有限元程序分析题目一：分析如图所示结构，其中5AB BC CD m ===， 3.5ED EF FG m ===，40GPa E =，20.02m A =，44410m I -=⨯。

fortran有限元程序课程设计一、课程目标知识目标：1. 掌握Fortran语言的基本语法和程序结构；2. 理解有限元方法的基本原理及其在工程问题中的应用；3. 学会使用Fortran编写有限元程序，解决简单的物理问题；4. 了解有限元程序的调试与优化方法。

技能目标：1. 能够运用Fortran语言编写简单的有限元程序；2. 能够对有限元程序进行调试和性能优化；3. 能够运用所学知识解决实际工程问题，具备一定的编程实践能力；4. 能够通过团队合作，共同完成较复杂的有限元程序编写。

情感态度价值观目标：1. 培养学生对编程和计算物理学的兴趣，激发学生的求知欲和探索精神；2. 培养学生严谨、细致、勤奋的学习态度，提高学生的问题解决能力；3. 培养学生的团队合作精神，提高沟通与协作能力；4. 增强学生的民族自豪感，认识我国在有限元领域的发展成果。

课程性质：本课程为高年级专业选修课，旨在使学生掌握Fortran有限元程序的编写和应用，提高学生的编程实践能力和解决实际问题的能力。

学生特点：学生已具备一定的数学、物理和编程基础，具有较强的逻辑思维能力和动手能力。

教学要求：结合课本内容，注重理论与实践相结合，强化编程实践，提高学生的实际操作能力。

同时，注重培养学生的团队合作精神，提高学生的综合素质。

通过本课程的学习，使学生能够独立编写和优化有限元程序，为后续学习和工作打下坚实基础。

二、教学内容1. Fortran语言基础：变量定义、数据类型、运算符、控制结构、数组、函数与子程序等；2. 有限元方法原理：有限元离散化、单元划分、形函数、刚度矩阵、载荷向量、边界条件处理等；3. 有限元程序编写：根据实际问题，运用Fortran语言编写有限元程序，包括前处理、核心计算和后处理；4. 程序调试与优化：调试技巧、性能分析、优化方法等；5. 实际工程案例：选取具有代表性的工程问题，运用所学的Fortran有限元程序解决。

有限元编程算例(Fortran)本程序通过Fortran语言编写，程序在Intel Parallel Studio XE 2013 with VS2013中成功运行，程序为《计算力学》（龙述尧等编）一书中的源程序，仅作研究学习使用，省去了敲写的麻烦。

源程序为：!Page149COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200) OPEN(5,FILE='DATAIN')!OPEN(6,FILE='DATAOUT',STATUS='NEW')CALL DATAIF(IND.EQ.0)GOTO 10EO=EO/(1.0-UN*UN)UN=UN/(1.0-UN)10 CALL TOTSTICALL LOADCALL SUPPORCALL SOLVEQCALL STRESSPAUSE!STOPENDSUBROUTINE DATACOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)READ(5,*)NJ,NE,NZ,NDD,NPJ,INDNJ2=NJ*2NPJ1=NPJ+1READ(5,*)EO,UN,GAMA,TEREAD(5,*)((JM(I,J),J=1,3),I=1,NE)READ(5,*)((CJZ(I,J),J = 1,2),I=1,NJ)!Page150READ(5,*)(NZC(I),I=1,NZ)READ(5,*)((PJ(I,J),J=1,2),I=1,NPJ1)WRITE(6,10)(I,(CJZ(I,J),J=1,2),I=1,NJ)10 FORMA T(4X,2HNO,6X,1HX,6X,1HY/(I6,2X,F7.2,F7.2))RETURNENDSUBROUTINE ELEST(MEO,IASK)COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200) IE=JM(MEO,1)JE=JM(MEO,2)ME=JM(MEO,3)CM=CJZ(JE,1)-CJZ(IE,1)BM=CJZ(IE,2)-CJZ(JE,2)CJ=CJZ(IE,1)-CJZ(ME,1)BJ=CJZ(ME,2)-CJZ(IE,2)AE=(BJ*CM-BM*CJ)/2.0IF(IASK.LE.1) GOTO 50DO 10 I=1,3DO 10 J=1,6B(I,J)=0.010 CONTINUEB(1,1)=-BJ-BMB(1,3)=BJB(1,5)=BMB(2,2)=-CJ-CMB(2,4)=CJB(2,6)=CMB(3,1)=B(2,2)B(3,2)=B(1,1)B(3,3)=B(2,4)B(3,4)=B(1,3)B(3,5)=B(2,6)!Page151B(3,6)=B(1,5)DO 20 I=1,3DO 20 J=1,6B(I,J)=B(I,J)/(2.0*AE) 20 CONTINUED(1,1)=EO/(1.0-UN*UN)D(1,2)=EO*UN/(1.0-UN*UN) D(2,1)=D(1,2)D(2,2)=D(1,1)D(1,3)=0.0D(2,3)=0.0D(3,1)=0.0D(3,2)=0.0D(3,3)=EO/(2.0*(1.0+UN))DO 30 I=1,3DO 30 J=1,6S(I,J)=0.0DO 30 K=1,3S(I,J)=S(I,J)+D(I,K)*B(K,J)30 CONTINUEIF(IASK.LE.2) GOTO 50DO 40 I=1,6DO 40 J=1,6EKE(I,J)=0.0DO 40 K=1,3!**********************************Exchange B And S***********************************************EKE(I,J)=EKE(I,J)+B(K,I)*S(K,J)*AE*TE40 CONTINUE50 CONTINUERETURNENDSUBROUTINE TOTSTICOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)!Page152DO 20 I=1,NJ2DO 20 J=1,NDDTKZ(I,J)=0.020 CONTINUE!*************Not Understanded*****************************DO 30 MEO=1,NECALL ELEST(MEO,3)DO 30 I=1,3DO 30 II=1,2LH=2*(I-1)+IILDH=2*(JM(MEO,I)-1)+IIDO 30 J=1,3DO 30 JJ=1,2L=2*(J-1)+JJLZ=2*(JM(MEO,J)-1)+JJLD=LZ-LDH+1IF(LD.LE.0) GOTO 30TKZ(LDH,LD)=TKZ(LDH,LD)+EKE(LH,L)30 CONTINUERETURNENDSUBROUTINE LOADCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)DO 10 I=1,NJ2P(I)=0.010 CONTINUEIF(NPJ.EQ.0) GOTO 30DO 20 I=1,NPJI1=I+1J=IFIX(PJ(I1,2))P(J)=PJ(I1,1)20 CONTINUE30 IF(GAMA.LE.0.0) GOTO 50!Page153DO 40 MEO=1,NECALL ELEST(MEO,1)PE=-GAMA*AE*TE/3.0IE=JM(MEO,1)JE=JM(MEO,2)ME=JM(MEO,3)P(2*IE)=P(2*IE)+PEP(2*JE)=P(2*JE)+PEP(2*ME)=P(2*ME)+PE40 CONTINUE50 CONTINUERETURNENDSUBROUTINE SUPPORCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)DO 60 I=1,NZMZ=NZC(I)TKZ(MZ,1)=1.0DO 10 J=2,NDDTKZ(MZ,J)=0.010 CONTINUEIF(MZ-NDD)20,20,3020 JO=MZGOTO 4030 JO=NDD40 DO 50 J = 2,JOJ1=MZ-JTKZ(J1+1,J)=0.050 CONTINUEP(MZ)=0.060 CONTINUERETURNEND!Page154SUBROUTINE SOLVEQCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)NJ1=NJ2-1DO 50 K=1,NJ1IF(NJ2-K-NDD+1)10,10,2010 IM=NJ2GOTO 3020 IM=K+NDD-130 K1=K+1DO 50 I=K1,IML=I-K+1C=TKZ(K,L)/TKZ(K,1)LD1=NDD-L+1DO 40 J=1,LD1M=J+I-KTKZ(I,J)=TKZ(I,J)-C*TKZ(K,M)40 CONTINUEP(I)=P(I)-C*P(K)50 CONTINUEP(NJ2)=P(NJ2)/TKZ(NJ2,1)DO 100 I1 = 1,NJ1I=NJ2-I1!************************************************************************下面一行可能出错IF(NDD-NJ2+I-1)60,60,7060 JO=NDDGOTO 8070 JO=NJ2-I+180 DO 90 J=2,JOLH=J+I-1P(I)=P(I)-TKZ(I,J)*P(LH)90 CONTINUEP(I)=P(I)/TKZ(I,1)100 CONTINUE!Page155WRITE(6,110)(I,P(2*I-1),P(2*I),I=1,NJ)!************************************************************************************ 110 FORMA T(2X,3HJD=,3X,2HU=,12X,2HV=/(I4,3X,F16.7,3X,F16.7))RETURNENDSUBROUTINE STRESSCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200) DIMENSION WY(6),YL(3)DO 60 MEO=1,NECALL ELEST(MEO,2)DO 10 I=1,3DO 10 J=1,2LH=2*(I-1)+JLDH=2*(JM(MEO,I)-1)+JWY(LH)=P(LDH)10 CONTINUEDO 20 I=1,3YL(I)=0.0DO 20 J=1,6YL(I)=YL(I)+S(I,J)*WY(J)20 CONTINUESIGX=YL(1)SIGY=YL(2)TOXY=YL(3)PYL=(SIGX+SIGY)/2.0SIG=(SIGX-SIGY)**2/4.0+TOXY*TOXYRYL=SQRT(SIG)SIG1=PYL+RYLSIG2=PYL-RYLIF(SIGY.EQ.SIG2) GOTO 30CETA1=TOXY/(SIGY-SIG2)CETA=90.0-57.29578*ATAN(CETA1)GOTO 40!Page15630 CETA=0.040 WRITE(6,50)MEO,SIGX,SIGY,TOXY,SIG1,SIG2,CETA50FORMA T(4X,2HE=,I3/2X,3HSX=,F11.3,3X,3HSY=,F11.3,3X,4HTAU=,F11.3/2X,3HS1=,F11.3,3X,3HS2=,F11. 3,3X,4HCET=,F11.3)!50FORMA T(4X,2HE=,I3/2X,3HSX=,Fll.3,3X,3HSY=,F11.3,3X,4HTAU=,F11.3/2X,3HSl=,Fll.3,3X,3HS2=,F11.3,3 X,4HCET=,F11.3)60 CONTINUERETURNEND输入文件为datain28,36,9,10,4,01,0.17,0,11,5,22,5,62,6,33,6,73,7,44,7,85,9,66,9,106,10,77,10,117,11,88,11,129,13,1010,13,1410,14,11 11,14,15 11,15,12 12,15,16 13,17,14 14,17,18 14,18,15 15,18,19 15,19,16 16,19,20 17,21,18 18,21,22 18,22,19 19,22,23 19,23,20 20,23,24 21,25,22 22,25,26 22,26,23 23,26,27 23,27,24 24,27,28 0,61,62,63,60,51,52,53,50,41,42,43,40,31,32,33,30,21,22,23,20,11,12,13,10,01,02,03,07,15,23,31,39,47,49,50,550,0-5E4,2-10E4,4-10E4,6-5E4,8输出结果为：DATAOUTNO X Y1 0.00 6.002 1.00 6.003 2.00 6.004 3.00 6.005 0.00 5.006 1.00 5.007 2.00 5.008 3.00 5.009 0.00 4.0010 1.00 4.0011 2.00 4.0012 3.00 4.0013 0.00 3.0014 1.00 3.0015 2.00 3.0016 3.00 3.0017 0.00 2.0018 1.00 2.0019 2.00 2.0020 3.00 2.0021 0.00 1.0022 1.00 1.0023 2.00 1.0024 3.00 1.0025 0.00 0.0026 1.00 0.0027 2.00 0.0028 3.00 0.00JD= U= V=1 -29766.873 -1173917.7502 -14003.185 -1174018.8753 -3753.270 -1179518.1254 0.000 -1181719.7505 -26382.471 -1072681.5006 -10746.993 -1073615.0007 -2064.593 -1082360.7508 0.000 -1085873.2509 -13536.995 -964010.12510 3372.794 -970055.12511 7268.415 -989269.12512 0.000 -998401.81213 7816.581 -835383.43814 27176.234 -861713.93815 22063.230 -905726.12516 0.000 -927165.18817 29514.479 -665602.87518 53419.637 -747340.43819 34876.832 -839806.81220 0.000 -881219.12521 29580.273 -416288.71922 52944.918 -632601.12523 17504.195 -803765.68824 0.000 -859481.93825 0.000 0.00026 -120102.820 -583505.37527 -76202.375 -787347.18828 0.000 -829170.812E= 1SX= -1489.530 SY=-101489.383 TAU= -1489.531 S1= -1467.348 S2=-101511.562 CET= 179.147 E= 2SX= -1475.844 SY=-100654.875 TAU= -1790.500 S1= -1443.531 S2=-100687.188 CET= 178.966 E= 3SX= -7021.670 SY=-101597.672 TAU= -3741.688 S1= -6873.875 S2=-101745.469 CET= 177.738 E= 4SX= -8067.500 SY= -98528.750 TAU= -4459.156 S1= -7848.227 S2= -98748.023 CET= 177.185 E= 5SX= -13143.328 SY= -99391.750 TAU= -1662.500 S1= -13111.293 S2= -99423.781 CET= 178.896 E= 6SX= -14652.781 SY= -98337.500 TAU= -1501.062S1= -14625.867 S2= -98364.414 CET= 178.973 E= 7SX= -2923.122 SY=-109168.297 TAU= -5888.469 S1= -2597.762 S2=-109493.656 CET= 176.837 E= 8SX= -716.078 SY=-103681.562 TAU= -8617.406 S1= 0.148 S2=-104397.789 CET= 175.249 E= 9SX= -9188.316 SY=-105121.867 TAU= -9771.594 S1= -8203.125 S2=-106107.062 CET= 174.243 E= 10SX= -12285.000 SY= -95180.250 TAU= -12199.594 S1= -10526.887 S2= -96938.359 CET= 171.799 E= 11SX= -14170.516 SY= -95500.750 TAU= -5489.531 S1= -13801.664 S2= -95869.602 CET= 176.156 E= 12SX= -22797.406 SY= -91347.000 TAU= -3902.844 S1= -22575.914 S2= -91568.492 CET= 176.752 E= 13SX= -5104.269 SY=-129494.438 TAU= -11708.750 S1= -4011.727 S2=-130586.977 CET= 174.669 E= 14SX= 969.672 SY=-108176.375 TAU= -21424.750 S1= 5024.582 S2=-112231.281 CET= 169.283 E= 15SX= -14954.572 SY=-110883.469 TAU= -18383.531 S1= -11552.273 S2=-114285.766 CET= 169.515 E= 16SX= -19890.141 SY= -86924.312 TAU= -25131.188 S1= -11514.844 S2= -95299.609 CET= 161.569 E= 17SX= -22109.688 SY= -87301.625 TAU= -10225.406 S1= -20543.453 S2= -88867.859 CET= 171.292 E= 18SX= -35190.453 SY= -77219.000 TAU= -9162.000 S1= -33280.023 S2= -79129.430 CET= 168.222 E= 19SX= -9785.850 SY=-171444.172 TAU= -20524.969 S1= -7220.594 S2=-174009.422 CET= 172.876 E= 20SX= 4594.438 SY=-113592.375 TAU= -46145.688 S1= 20477.398 S2=-129475.336 CET= 161.007 E= 21SX= -25287.307 SY=-118672.312 TAU= -30023.750 S1= -16467.512 S2=-127492.109 CET= 163.629 E= 22SX= -30634.422 SY= -71127.188 TAU= -44991.469 S1= -1543.715 S2=-100217.891 CET= 147.114 E= 23SX= -34259.609 SY= -71743.438 TAU= -14637.906 S1= -29220.699 S2= -76782.344 CET= 161.005 E= 24SX= -43958.047 SY= -53418.938 TAU= -17697.562 S1= -30369.627 S2= -67007.359 CET= 142.482 E= 25SX= -19028.160 SY=-252549.000 TAU= -34958.688 S1= -13907.055 S2=-257670.094 CET= 171.666 E= 26SX= 3973.812 SY=-114063.750 TAU= -92238.344 S1= 54459.047 S2=-164548.984 CET= 151.307 E= 27SX= -39180.809 SY=-121400.055 TAU= -39312.688 S1= -23409.074 S2=-137171.781 CET= 158.140 E= 28SX= -42804.766 SY= -43317.938 TAU= -65723.062 S1= 22662.211 S2=-108784.914 CET= 135.112 E= 29SX= -42224.094 SY= -43219.188 TAU= -10273.375 S1= -32436.225 S2= -53007.055 CET= 136.386 E= 30SX= -21830.422 SY= -25448.312 TAU= -23810.344 S1= 239.594 S2= -47518.328 CET= 137.172 E= 31SX= -48815.199 SY=-424587.344 TAU= -79800.078 S1= -32570.844 S2=-440831.688 CET= 168.494 E= 32SX=-132271.750 SY= -71582.000 TAU=-175409.250 S1= 76087.781 S2=-279941.531 CET= 130.093 E= 33SX= -45090.102 SY= -56761.105 TAU= 804.781 S1= -45034.867 S2= -56816.336 CET= 3.926 E= 34SX= 42332.711 SY= -9221.938 TAU= -47066.344 S1= 70218.328 S2= -37107.555 CET= 149.354 E= 35SX= -20899.344 SY= -19971.375 TAU= 16235.219 S1= -4193.512 S2= -36677.207 CET= 45.819E= 36SX= 73163.914 SY= -17873.250 TAU= -17873.344 S1= 76547.250 S2= -21256.586 CET= 169.281。

有限元程序设计的一些总结2012.11.051、program main因为我现在所看的有限元平面计算程序还未考虑施工过程，不需考虑分级加载的情况，所以相对来说主程序是比较简单的。

只需要控制整个程序的流程就可以了。

现在我所涉及到的主程序流程如下图所示：inputmmaaskfff solve stress xhs xhs1xhsoutput这样的话只需要熟悉各流程（子程序）的作用，就可以从宏观上把握最简单的有限元平面程序的框架。

在主程序中只需要调用各子程序就可以了。

2、subroutine inputInput子程序要实现的功能是将建立的有限元模型的基本信息输入到程序中去，以便进行信息处理。

因为我毕业设计做的就相当于一个有限元前处理过程，我对这个相对来说熟悉点。

一般平面应力应变问题，需要输入的信息包括：单元信息，节点信息，约束信息，材料信息，荷载信息（一般在fff子程序中输入）。

（1）单元信息：各个单元的编号、四个节点号（平面等参单元，注意编号顺序）、材料号。

（2）节点信息：节点号、节点x坐标、节点y坐标（平面问题）（3）约束信息：a.零位移约束。

一般1表示自由，0表示约束b.非零位移约束。

（我放在荷载信息中给出的，参看fff子程序）（4）材料信息：材料信息的输入与所采用的模型相关，如弹性模型、E-B模型等等（我理解的是考虑不同的本构关系）。

如果采用弹性模型，则一般输入信息为：材料号、容重、弹性模量、泊松比。

（5）荷载信息在fff 再谈。

注：真正input.txt 以什么样的格式编写，与各人的习惯有关，只需要把以上信息包含进去就可以了。

3、subroutine mma子程序mma 的作用是：对自由度编号、找出每行的半带宽、储存主元位置 因为整体刚度矩阵K 是一个高度稀疏、非零元素呈带状分布、对称的矩阵，在程序中是以一维变半带宽存储结构储存的。

什么叫一维变半带宽存储结构？首先，它是一维存储结构。

K 是一个nh*nh 的矩阵，将K 按一维储存就是用一个（nh*nh ）*1的向量去逐行存储K 中的元素。

有限元计算结构力学fortran程序计算结构力学程序计算结构力学编程大作业时间,2007年6月!!!***************************************************************** ***********!!! 关于程序的说明!!!***************************************************************** ***********!一、功能:! 1、可计算包括节点力，一般非节点力，支座沉降、温度荷载作用、制造误差的平! 面桁架、梁、刚架及其组合结构的节点位移与杆端力;! 2、可同时计算多种工况下的节点位移与杆端力。

!******************************************************************* **********!******************************************************************* ***********!! 二、变量说明:! NE——单元数;! N——结构中自由度数;! NJ——节点数;! NS——特殊节点数，包括支座节点、主从节点(1节点不做主节点)、连接桁架的铰节点(没有转角);! NAI——结构的单元截面类型数;! MT——单元截面类型号;! NL——荷载工况数;! H——截面高度;! E——弹性模量;! JC——单元定位向量数组;! X(NJ),Y(NJ)——节点的X，Y坐标值;! JE(NE,2)——单元两端节点码数组;! AI(NAI,2)——按单元类型顺序存放A与I，AI(I，1)—第I类单元的截面积，AI(I，2)—第I类单元的! 惯性矩;! MT(NE)——单元所属单元类型号;! JS(NS,4)——特殊节点信息，JS(I,1)—结点码;JS(I,2)，JS(I,3)，JS(I,4)—U，V，CETA约束信息，! 有约束为1，没有约束为0;从节点某位移同主节点时位移时，该位移约束信息填主节点码;!! PJ(NP,3)——节点荷载信息数组;PJ(I，1)—节点力所在节点号;PJ(I，2)—节点力作用坐标方向:! 坐标方向U，V，M分别为1，2，3; PJ(I，3)—节点力的大小(含正负号);U，V方向集中力时，! 与坐标轴正向同向为正，M按右手法则为正;本程序推导过程取y轴向下为正。

作业2：有限单元法的编程实现2.1有限元概述有限元分析的基本概念是将求解域离散为若干子域，并通过它们边界上的节点相互联结成为组合体，对每一单元假定一个合适的近似解，然后推导求解这个域总的满足条件，从而得到问题的解。

这个解只是近似解，因为实际问题被较简单的问题所代替。

由于大多数实际问题难以得到准确解，而有限元不仅计算精度高，而且能适应各种复杂形状，因而成为行之有效的工程分析手段。

对于不同物理性质和数学模型的问题，有限元求解法的基本步骤是相同的，只是具体公式推导和运算求解不同。

有限元求解问题的基本步骤通常为：2.1.1确定求解域及求解域的离散化，即子域的剖分根据实际问题近似确定求解域的物理性质和几何区域。

将求解域近似为具有不同有限大小和形状且彼此相连的有限个单元组成的离散域，习惯上称为有限元网络划分，求解域的离散化是有限元法的核心技术之一。

选择单元的形式，确定单元数，节点数及单元及节点的编号。

2.1.2单元分析一个具体的物理问题通常可以用一组包含问题状态变量边界条件的微分方程式表示，为适合有限元求解，通常将微分方程化为等价的泛函形式。

作为弹性力学微分方程的等效积分的形式，虚位移原理与虚应力原理分别是平衡方程与力的边界条件和几何方程与位移边界条件的等效积分形式。

在导出它们的的过程中都未涉及到物理方程所以它们适用于线弹性、非线性线弹性及弹塑性的问题。

现有的有限元计算多采用以位移为未知量的形式，可用虚位移原理来描述其平衡方程，其矩阵形式为：(u )u 0T T T V S f dV TdS σδεσδδ--=⎰⎰ (2.1) 我们需要对单元构造一个适合的近似解，即推导有限单元的格式，其中包括选择合理的单元坐标系，建立单元试函数，以某种方法给出单元各状态变量的离散关系，从而形成单元矩阵（结构力学中称刚度阵或柔度阵）。

这里以平面问题为例，单元位移与节点位移的关系表示为：[,,...]e i e e i i i j j a u u N a N N a Na ⎧⎫⎪⎪====⎨⎬⎪⎪⎩⎭∑ (2.2) 应变与节点位移的关系为：[][]e e e e i j i j Lu LNa L N N a B B a Ba ε===== (2.3)应力与节点位移的关系为：e e D DBa Sa σε=== (2.4)单元刚度矩阵的表达式为：T ee K B DBtdxdy Ω=⎰ (2.5) 单元等效节点载荷矩阵表示为T T e ee e ef s P P P N ftdxdy N TtdS ΩΩ=+=+⎰⎰ (2.6) 2.1.3 整体分析将单元总装形成离散域的总矩阵方程（联合方程组），反映对近似求解域的离散域的要求，即单元函数的连续性要满足一定的连续条件。

题目：《FORTRAN语言》算例程序专业：班级：学号：姓名：时间：目录1.平面深梁有限元前处理信息生成1.1深梁有限元模型(3)1.2源程序(3)1.3输出结果和分析(4)2.矩阵求逆子程序调试2.1子程序源程序(9)2.2考题(11)2.2.1主程序(11)2.3考题及验证(11)1.平面深梁有限元前处理信息生成1.1深梁有限元模型16117176 81 每个小方格长和宽都为0.25。

1.2源程序c programdimension xy(85,2),ijk(128,3)open(1,file="xyijk.dat",status="unknown")do i=1,17do j=1,5m=(i-1)*5+jxy(m,1)=0.25*(i-1)xy(m,2)=0.25*(j-1)end doend dodo i=1,85write(1,*) i, (xy(i,j),j=1,2)end dodo i=1,16do j=1,4ns=8*(i-1)+jnx=ns+4ijk(ns,1)=(i-1)*5+jijk(ns,2)=ijk(ns,1)+6ijk(ns,3)=ijk(ns,1)+1ijk(nx,1)=ijk(ns,1)ijk(nx,2)=ijk(ns,1)+5ijk(nx,3)=ijk(ns,1)+6end doend dodo i=1,128write(1,*) i, (ijk(i,j),j=1,3)end doclose(1)end1.3输出结果及分析1 0.000000E+00 0.000000E+002 0.000000E+00 2.500000E-013 0.000000E+00 5.000000E-014 0.000000E+00 7.500000E-015 0.000000E+00 1.0000006 2.500000E-01 0.000000E+007 2.500000E-01 2.500000E-018 2.500000E-01 5.000000E-019 2.500000E-01 7.500000E-0110 2.500000E-01 1.00000011 5.000000E-01 0.000000E+0012 5.000000E-01 2.500000E-0113 5.000000E-01 5.000000E-0114 5.000000E-01 7.500000E-0115 5.000000E-01 1.00000016 7.500000E-01 0.000000E+0017 7.500000E-01 2.500000E-0118 7.500000E-01 5.000000E-0119 7.500000E-01 7.500000E-0121 1.000000 0.000000E+0022 1.000000 2.500000E-0123 1.000000 5.000000E-0124 1.000000 7.500000E-0125 1.000000 1.00000026 1.250000 0.000000E+0027 1.250000 2.500000E-0128 1.250000 5.000000E-0129 1.250000 7.500000E-0130 1.250000 1.00000031 1.500000 0.000000E+0032 1.500000 2.500000E-0133 1.500000 5.000000E-0134 1.500000 7.500000E-0135 1.500000 1.00000036 1.750000 0.000000E+0037 1.750000 2.500000E-0138 1.750000 5.000000E-0139 1.750000 7.500000E-0140 1.750000 1.00000041 2.000000 0.000000E+0042 2.000000 2.500000E-0143 2.000000 5.000000E-0144 2.000000 7.500000E-0145 2.000000 1.00000046 2.250000 0.000000E+0047 2.250000 2.500000E-0148 2.250000 5.000000E-0149 2.250000 7.500000E-0150 2.250000 1.00000051 2.500000 0.000000E+0052 2.500000 2.500000E-0153 2.500000 5.000000E-0154 2.500000 7.500000E-0155 2.500000 1.00000056 2.750000 0.000000E+0057 2.750000 2.500000E-0158 2.750000 5.000000E-0159 2.750000 7.500000E-0160 2.750000 1.00000061 3.000000 0.000000E+0062 3.000000 2.500000E-0163 3.000000 5.000000E-0165 3.000000 1.00000066 3.250000 0.000000E+0067 3.250000 2.500000E-0168 3.250000 5.000000E-0169 3.250000 7.500000E-0170 3.250000 1.00000071 3.500000 0.000000E+0072 3.500000 2.500000E-0173 3.500000 5.000000E-0174 3.500000 7.500000E-0175 3.500000 1.00000076 3.750000 0.000000E+0077 3.750000 2.500000E-0178 3.750000 5.000000E-0179 3.750000 7.500000E-0180 3.750000 1.00000081 4.000000 0.000000E+0082 4.000000 2.500000E-0183 4.000000 5.000000E-0184 4.000000 7.500000E-0185 4.000000 1.0000001 1 7 22 2 8 33 3 9 44 4 10 55 16 76 27 87 3 8 98 4 9 109 6 12 710 7 13 811 8 14 912 9 15 1013 6 11 1214 7 12 1315 8 13 1416 9 14 1517 11 17 1218 12 18 1319 13 19 1420 14 20 1521 11 16 1722 12 17 1824 14 19 2025 16 22 1726 17 23 1827 18 24 1928 19 25 2029 16 21 2230 17 22 2331 18 23 2432 19 24 2533 21 27 2234 22 28 2335 23 29 2436 24 30 2537 21 26 2738 22 27 2839 23 28 2940 24 29 3041 26 32 2742 27 33 2843 28 34 2944 29 35 3045 26 31 3246 27 32 3347 28 33 3448 29 34 3549 31 37 3250 32 38 3351 33 39 3452 34 40 3553 31 36 3754 32 37 3855 33 38 3956 34 39 4057 36 42 3758 37 43 3859 38 44 3960 39 45 4061 36 41 4262 37 42 4363 38 43 4464 39 44 4565 41 47 4266 42 48 4368 44 50 4569 41 46 4770 42 47 4871 43 48 4972 44 49 5073 46 52 4774 47 53 4875 48 54 4976 49 55 5077 46 51 5278 47 52 5379 48 53 5480 49 54 5581 51 57 5282 52 58 5383 53 59 5484 54 60 5585 51 56 5786 52 57 5887 53 58 5988 54 59 6089 56 62 5790 57 63 5891 58 64 5992 59 65 6093 56 61 6294 57 62 6395 58 63 6496 59 64 6597 61 67 6298 62 68 6399 63 69 64 100 64 70 65 101 61 66 67 102 62 67 68 103 63 68 69 104 64 69 70 105 66 72 67 106 67 73 68 107 68 74 69 108 69 75 70 109 66 71 72 110 67 72 73111 68 73 74112 69 74 75113 71 77 72114 72 78 73115 73 79 74116 74 80 75117 71 76 77118 72 77 78119 73 78 79120 74 79 80121 76 82 77122 77 83 78123 78 84 79124 79 85 80125 76 81 82126 77 82 83127 78 83 84128 79 84 852.矩阵求逆子程序调试2.1子程序源程序SUBROUTINE gjcp(A,N)integer Nreal A(N,N) !存放矩阵AINTEGER IP(N) !记录主列号REAL P !工作单元，放主元INTEGER I0,R !工作单元，放主列号EPS=0.001write(*,*)'gjcp ok0000000' DO K=1,NP=0I0=KIP(K)=KDO I=K,NIF(ABS(A(I,K)).GT.ABS(P))THENP=A(I,K)I0=IIP(K)=IENDIFENDDOIF(ABS(P).LE.EPS)THENWRITE(*,*)'DET=0'stop !后来加的GOTO 10ENDIFIF(I0.NE.K)THENDO J=1,NS=A(K,J)A(K,J)=A(I0,J)A(I0,J)=SENDDOENDIFA(K,K)=1./PDO I=1,NIF(I.NE.K)THENA(I,K)=-A(I,K)*A(K,K)DO J=1,NIF(J.NE.K)THENA(I,J)=A(I,J)+A(I,K)*A(K,J)ENDIFENDDOENDIFENDDODO J=1,NIF(J.NE.K)THENA(K,J)=A(K,K)*A(K,J)ENDIFENDDOENDDOwrite(*,*)'gjcp ok01'DO K=N-1,1,-1R=IP(K)IF(R.NE.K)THENDO I=1,NS=A(I,R)A(I,R)=A(I,K)A(I,K)=SENDDOENDIFENDDO10 write(*,*)'the end'END2.2考题1 1 2求矩阵-1 2 0 的逆矩阵。