平面梁单元有限元fortran程序

合集下载

Fortran语言编写弹性力学平面问题3节点三角形单元或4节点等参单元的有限元程序

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源程序

八节点平面等参元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。

平面框架结构的有限元分析

三梁平面框架结构的有限元分析一、问题说明如图1所示的框架结构，其顶端受均布载荷作用，用有限元方法分析该结构的位移。

结构中各个部分的参数为：弹性模量E=300GPa，截面惯性矩I=6.5×105mm4，横截面积A=680mm2。

相应的有限元分析模型见图2，利用梁板壳分析程序完成该模型的力学分析。

图1框架结构图2有限元分析模型二．Fortran程序的输入数据（1）Facile.11 4 3 6 0 12 42 1 11 1 11 3 51 2 2 3 3 40 0 0 0 1000 01000 1000 0 1000 0 0（2）Facile.2111 211 1111 0 0 0 1 03E5 1.6E5680 6.5E5 6.5E5 6.5E50 0（3）Facile.312 41 02 03 04 05 06 0 19 0 20 0 21 0 22 023 0 24 08 -1200 12 -200000 14 -1200 18 200000输出的数据文件为：Facile7和Facile8，其中各节点位移结果在文件Facile8中。

三．计算结果各节点的位移计算结果见表1。

四．Ansys分析结果Ansys计算结果如下图所示，图3为节点x方向的位移云图，图4为节点y 方向的位移云图，图5为节点转角云图。

图3 节点x方向的位移图4 节点y方向的位移图5 节点转角各节点的位移值见表2。

五．结果对比通过对比表1和表2中的数据可以发现，Fortran程序与Ansys分析的结果十分接近。

Fortran平面钢架有限元分析

1 有限元分析软件的开发1.1 程序功能该程序为平面刚架静力分析程序，能针对平面刚架间问题进行有限元计算，计算杆端位移及杆端力大小。

程序从磁盘文件中读取单元编号、节点编号及坐标、材料属性、荷载、边界条件等信息；将杆端位移，杆端力等计算结果以磁盘文件的形式输出，采用等带宽二维数组存储整体刚度矩阵并使用高斯消去法进行求解。

1.2 程序结构及流程1.3 程序的输入与输出详细介绍输入输出数据的格式。

如：数据文件分几个部分，各有几行，分别包含哪些内容及其类型、先后次序，等等。

输入，共有九行。

第一行：7,13,5,1,2,2。

分别为，7个结点，13个自由度，5个单元，1个类型，2个结点荷载，2个非结点荷载。

第二行：1,2,3,0.0,0.0,0,0,,6.0,0.0。

分别为：一号结点的位移序号，x方向为1，y方向为2，转角为3，坐标为(0.0,0.0)，因为二号结点固结在地面，所以二号结点的位移序号，x方向为0，y方向为0，转角为0，坐标为(6.0,0.0)。

第三行：4,5,6,0.0,6.0,4,5,7,0.0,6.0。

分别为：三号结点的位移序号，x方向为4，y方向为5，转角为6，坐标为(0.0,6.0), 四号结点位移序号x方向和y相同，转角为7，坐标为(,0.0,6.0)。

第四行：8,9,10,6.0,6.0,0,0,11,0.0,12.0.五号结点位移序号，x方向为8，y方向为9，转角为10，坐标为（6.0,6.0）。

因为六号结点铰接在地面，所以六号结点的位移序号，x方向和y方向为0，转角为11，坐标为（0.0,12.0）。

第五行：12,0,13,6.0,12.0. 因为七号结点与地面用滑动支座固定，所以七号结点的位移序号，x方向为12，y方向为0，转角为13，坐标（6.0,12.0）.第六行：1,2,1,1,3,1,4,5,1,3,6,1,5,7,1,分别为，1号和2号结点组成的单元为1号类型。

1号和3号结点组成的单元为1号类型，4号和5号结点组成的为1号类型，3号和6号结点组成的单元为1号类型，5号和7号结点组成的单元为1号类型。

平面桁架杆单元fortran程序

PROGRAM BAR3DDIMENSION MEA(100,4),JZ(50,2),AAI(10),AEU(10,2),R(2,6),RT(6,2) DIMENSION AKEP(2,2),AKE(6,6),P(300),X(100),Y(100),Z(100)DIMENSION LD(300),IS(6),A(50000),UPE(2),PPE(2),UE(6),RA(6,2)OPEN(1,FILE='FNAME1_t.txt')OPEN(2,FILE='FNAME2_t.txt')READ(1,*) NN,NE,NM,NA,NCwrite(*,*) NN,NE,NM,NA,NCNF=3ND=2NFD=NF*NDN=NN*NFDO 50 I=1,NN50 READ(1,*)K,X(I),Y(I),Z(I),(P(NF*(I-1)+J),J=1,NF)READ(1,*)((JZ(I,J),J=1,2),I=1,NC)DO 100 I=1,NE100 READ(1,*) IE,(MEA(I,J),J=1,4)READ(1,*) (AAI(I),I=1,NA)READ(1,*) ((AEU(I,J),J=1,2),I=1,NM)CLOSE(1)WRITE(2,460) NN,NE,NM,NA,NCWRITE(2,465) (I,X(I),Y(I),Z(I),P(3*I-2),P(3*I-1),P(3*I),I=1,NN)WRITE(2,470) ((JZ(I,J),J=1,2),I=1,NC)write(2,475) (I,(Mea(i,j),j=1,4),i=1,Ne)WRITE(2,480) (I,AAI(I),I=1,NA)WRITE(2,485) (I,(AEU(I,J),J=1,2),I=1,NM)400 FORMAT(5I5)405 FORMAT(I5,6F10.0)410 FORMAT(2I5)415 FORMAT(5I5)420 FORMAT(F10.0)425 FORMAT(2F10.0)460 FORMAT(/5X,'THE INPUT OF NN,NE,NM,NA,NC'//(5X,5I5)) 465 FORMAT(/5X,'THE INPUT OF X,Y,Z,P'//(5X,I5,6F10.2))470 FORMAT(/5X,'THE INPUT OF JZ'//(5X,I5,5X,I5))475 FORMAT(/5X,'THE INPUT OF MEA'//(5X,I5,5X,4I5))480 FORMAT(/5X,'THE INPUT OF AAI'//(5X,I5,5X,F20.6))485 FORMAT(/5X,'THE INPUT OF AEU'//(5X,I5,5X,2E15.6)) CALL FLD(NN,NE,MEA,NF,ND,N,NT,LD)DO 500 I=1,NT500 A(I)=0.0DO 600 IE=1,NECALL KEP(NN,NA,IE,MEA,AEU,AAI,X,Y,Z,AKEP)!CALL CR(NN,IE,MEA,X,Y,Z,R)CALL TRAN (2,6,R,RT)CALL DOT(6,2,2,RT,AKEP,RA)CALL DOT(6,2,6,RA,R,AKE)CALL FIS(IE,MEA,NF,ND,NFD,IS)DO 560 I=1,NFDDO 560 J=1,NFDIF (IS(I)-IS(J)) 560,520,520520 NI=IS(I)IJ= LD(NI)-NI+IS(J)A(ij)=A(ij)+AKE(I,J)560 CONTINUE600 CONTINUECALL FCC(NC,N,NT,NF,JZ,LD,A)CALL DECOM(N,NT,A,P,LD)WRITE(2,850) (I,P(3*I-2),P(3*I-1),P(3*I),I=1,NN)WRITE(2,860)DO 800 IE=1,NECALL KEP(NN,NA,IE,MEA,AEU,AAI,X,Y,Z,AKEP)CALL CR(NN,IE,MEA,X,Y,Z,R)CALL FIS(IE,MEA,NF,ND,NFD,IS)DO 750 I=1,NFDNI=IS(I)UE(I)=P(NI)750 CONTINUECALL DOT(2,NFD,1,R,UE,UPE)CALL DOT(2,2,1,AKEP,UPE,PPE)NA1=MEA(IE,3)AO=AAI(NA1)STRESS=PPE(2)/AOWRITE(2,880) IE,PPE(2),STRESS800 CONTINUE850 FORMAT(//25X,'NODAL DISPLACEMENT'// 5X,'NODE No.',12X,'U',12X,'V',12X,'W'/(7X,I5,3E15.6))860 FORMAT(//15X,'AXIAL FORCE AND STRESS OF ELEMENT'/11X,'ELEMENT No.',9X,'STRESS') 880 FORMAT(14X,I5,2E16.6)CLOSE (2)END!C BAR ELEMENT (3-D) COORDINATED TRANSFORMATIONSUBROUTINE CR(NN,IE,MEA,X,Y,Z,R)DIMENSION X(NN),Y(NN),Z(NN),MEA(100,4),R(2,6)CALL ZERO(2,6,R)I=MEA(IE,1)J=MEA(IE,2)SL=SQRT((X(J)-X(I))**2+(Y(J)-Y(I))**2+(Z(J)-Z(I))**2)R(1,1)=(X(J)-X(I))/SLR(1,2)=(Y(J)-Y(I))/SLR(1,3)=(Z(J)-Z(I))/SLR(2,4)=R(1,1)R(2,5)=R(1,2)R(2,6)=R(1,3)RETURNEND!C BAR ELENMENT LACAL STIFFNESSSUBROUTINE KEP(NN,NA,IE,MEA,AEU,AAI,X,Y,Z,AKEP)DIMENSION X(NN),Y(NN),Z(NN),MEA(100,4)DIMENSION AEU(10,2),AAI(NA),AKEP(2,2)NM1=MEA(IE,4)Eo=AEU(NM1,1)NA1=MEA(IE,3)AO=AAI(NA1)I=MEA(IE,1)J=MEA(IE,2)SL=SQRT((X(J)-X(I))**2+(Y(J)-Y(I))**2+(Z(J)-Z(I))**2)AKEP(1,1)=EO*AO/SLAKEP(2,2)=AKEP(1,1)AKEP(1,2)=-AKEP(1,1)AKEP(2,1)=AKEP(1,2)RETURNEND!C IS DimensionSUBROUTINE FIS(IE,MEA,NF,ND,NFD,IS)DIMENSION MEA(100,4),IS(NFD)DO 150 ID=1,NDDO 100 JF=1,NFI1=(ID-1)*NF+JFIS(I1)=(MEA(IE,ID)-1)*NF+JF100 CONTINUE150 CONTINUERETURNEND!C LD DimensionSUBROUTINE FLD(NN,NE,MEA,NF,ND,N,NT,LD)DIMENSION MEA(100,4),LD(N)LD(1)=1DO 300 K=1,NNNMIN=1000DO 200 I=1,NEDO 150 J=1,NDIF(MEA(I,J).NE.K) GOTO 150DO 100 L=1,NDIF (MEA(I,L).LT.NMIN) NMIN=MEA(I,L)100 CONTINUE150 CONTINUE200 CONTINUEDO 250 I=1,NFJ=(K-1)*NF+IIF(J.NE.1) LD(J)=LD(J-1)+(K-NMIN)*NF+I 250 CONTINUE300 CONTINUENT=LD(N)RETURNEND!C CONSTRAINED CONDITIONSUBROUTINE FCC(NC,N,NT,NF,JZ,LD,A)DIMENSION JZ(50,2),LD(N),A(NT)DO 350 K=1,NCI=JZ(K,1)J=JZ(K,2)NI=(I-1)*NF+JNJ=LD(NI)A(NJ)=1.0E25350 CONTINUERETURNEND! C LINEAR EQUATION BY LDLT METHOD SUBROUTINE DECOM(N,NT,A,B,LD)DIMENSION A(NT),B(N),LD(N)DO 20 I=1,NLDI=LD(I)IF (I.EQ.1) GOTO 10IO=LDI-IIM1=I-1MI=1-IO+LD(IM1)IF(MI.GT.IM1) GOTO 10DO 8 J=MI,IM1JO=LD(J)-JIJ=IO+JIF(J.EQ.1) GOTO 6JM1=J-1MJ=1-JO+LD(JM1)MIJ=MAX0(MI,MJ)IF(MIJ.GT.JM1) GOTO 6S=0.0DO 2 K=MIJ,JM1IK=IO+KJK=JO+KS=S+A(IK)*A(JK)2 CONTINUEA(IJ)=A(IJ)-S6 B(I)=B(I)-A(IJ)*B(J)8 CONTINUE10 ALDI=A(LDI)IF(I.EQ.1.OR.MI.GT.IM1) GOTO 13DO 12 J=MI,IM1IJ=IO+JLDJ=LD(J)S=A(IJ)A(IJ)=S/A(LDJ)ALDI=ALDI-S*A(IJ)12 CONTINUE13 A(LDI)=ALDIB(I)=B(I)/ALDI20 CONTINUEDO 40 LDI=2,NI=N-LDI+2IO=LD(I)-IMI=1-IO+LD(I-1)IM1=I-1IF(MI.GT.IM1) GOTO 40DO 30 J=MI,IM1IJ=IO+JB(J)=B(J)-A(IJ)*B(I)30 CONTINUE40 CONTINUERETURNENDSUBROUTINE ZERO(a,b,c)DIMENSION A(M,N)DO 10 I=1,MDO 10 J=1,NA(I,J)=0.010 CONTINUERETURNENDSUBROUTINE TRAN(M,N,A,B)DIMENSION A(M,N),B(N,M)DO 20 I=1,MDO 20 J=1,NB(J,I)=A(I,J)20 CONTINUERETURNENDSUBROUTINE DOT(M,N,L,A,B,C)DIMENSION A(M,N),B(N,L),C(M,L)DO 50 I=1,MDO 50 J=1,LC(I,J)=0.0DO 40 K=1,NC(I,J)=C(I,J)+A(I,K)*B(K,J)40 CONTINUE50 CONTINUERETURNEND8 13 1 1 111 0.0 0.0 0.0 0.0 0.0 0.02 2.0 0.0 0.0 0.0 0.0 0.03 2.0 1.0 0.0 0.0 -20000.0 0.04 4.0 2.0 0.0 0.0 -10000.0 0.05 4.0 0.0 0.0 0.0 0.0 0.06 6.0 0.0 0.0 0.0 0.0 0.07 6.0 1.0 0.0 0.0 0.0 0.08 8.0 0.0 0.0 0.0 0.0 0.01 11 21 32 33 34 35 36 37 38 28 31 12 1 12 13 1 13 2 3 1 14 25 1 15 3 5 1 16 3 4 1 17 5 4 1 18 5 6 1 19 5 7 1 110 4 7 1 111 6 7 1 112 6 8 1 113 7 8 1 1 5002.0E5 0.3。

Fortran语言-有限元程序分析-平面钢架

程序框图：程序特点：问题类型：可用于计算结构力学的平面刚架问题单元类型：直接利用杆单元载荷类型：节点载荷及非节点载荷，其中非节点载荷包括均布荷载和垂直于杆件的集中荷载材料性质：所有杆单元几何性质相同，且由相同的均匀材料组成方程求解：结构刚度矩阵采用满阵存放，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程序)飞箭系列

FEPG 中级教程
第一章
偏微分方程的“弱”形式 ——虚位移原理
§1.1 偏微分方程的弱解形式
1.1.1 问题的提出
工程或物理学中的许多问题，通常是以未知场函数应满足的偏微分方程和边界条件的形式提出来的，可以一般地表示为未知函数 u 应满足偏微分方程组
A1 (u ) A(u ) = A2 (u ) = 0 M
2.2.1 导数之间的变换………………………………………………………46
§2.3 数值积分………………………………………………………………49
2.3.1 高斯积分……………………………………………………………49 2.3.2 节点积分……………………………………………………………51
第三章有限元输入数据形式
§1.4 求解解梯度的最小二乘法………………………………………27
1.4.1 已知位移求应力………………………………………………………28 1.4.2 已知温度求热流密度…………………………………………………29 1.4.3 已知电势求电场强度…………………………………………………30
第二章分片多项式的形函数
（在Ω内） ,
(1.1.1)
Ω 域可以是体积域、面积域等，如图 1.1.1 所示。同时未知函数 u 还应满足边界条件
B1 (u ) B(u ) = B2 (u ) = 0 M Γ 是 Ω 域的边界。 y
B (u ) = 0
（在Γ上）
(1.1.2)
Ω 域
A(u ) = 0
4.2.1.4 源程序……………………………………………………………83 4.2.1.5 Fortran 源程序…………………………………………………88 4.2.2 BFT 元件程序……………………………………………………………89 4.2.2.1 功能………………………………………………………………89 4.2.2.2 命令行参数说明…………………………………………………90 4.2.2.3 参数及数组说明…………………………………………………90 4.2.2.4 源程序……………………………………………………………91 4.2.2.5 Fortran 源程序…………………………………………………94 4.2.3 E 元件程序………………………………………………………………96 4.2.3.1 功能……………………………………………………………96 4.2.3.2 命令行参数说明…………………………………………………97 4.2.3.3 参数及数组说明…………………………………………………97 4.2.3.4 源程序……………………………………………………………98 4.2.3.5 Fortran 源程序………………………………………………107 4.2.4 SOLV 求解器……………………………………………………………109 4.2.4.1 功能……………………………………………………………109 4.2.4.2 命令行参数说明………………………………………………109 4.2.4.3 源程序…………………………………………………………109 A.直接法求解…………………………………………………109 B.迭代法求解…………………………………………………119 4.2.4.4 Fortran 源程序………………………………………………129 4.2.5 U 元件程序……………………………………………………………131 4.2.5.1 功能……………………………………………………………131 4.2.5.2 命令行参数说明………………………………………………131 4.2.5.3 参数及数组说明………………………………………………132 4.2.5.4 源程序…………………………………………………………132 4.2.5.5.Fortran 源程序………………………………………………134

fortran有限元程序课程设计

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)本程序通过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。

有限元编程算例(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。

有限元计算结构力学fortran程序

有限元计算结构力学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轴向下为正。

fortran算例程序

题目：《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 的逆矩阵。

Fortran 讲解-平面有限元源程序

Fortran 讲解FORTRAN是科学计算语言。

很多年不碰它了（87年用）。

你怎么不去FOR TRAN论坛？试着说说：SUBROUTINE KE(IO,NE,NWE,T,A1,A2,V,EK,BCA) *子例行程序KE(可以CALL调用)DIMENSION B(7),BCA(7,NE),EK(6,6) *定义1个一维数组和2个两维数组DO 10 I=1,7 *循环，10是下面的标号B(I)=BCA(I,IO) *给一维数组B(I)赋值。

但BCA函数我没见过，外部函数？10 CONTINUE *未完成7次循环继续A=A1/B(7)*T *表达式计算结果赋值给ADO 20 I=1,3 *下面是双重循环——外循环DO 20 J=I,3 *内循环I1=2*IJ1=2*JEK(I1-1,J1-1)=A*(B(I)*B(J)+A2*B(I+3)*B(J+3)) *2维数组通过表达式计算后赋值EK(I1-1,J1)=A*(V*B(I)*B(J+3)+A2*B(I+3)*B(J))EK(I1,J1-1)=A*(V*B(I+3)*B(J)+A2*B(I)*B(J+3))EK(I1,J1)=A*(B(I+3)*B(J+3)+A2*B(I)*B(J))20 CONTINUE *循环未完成继续DO 30 I=3,6 *有是一个双重循环DO 30 J=1,IEK(I,J)=EK(J,I)30 CONTINUE *循环未完成继续IF(NWE.EQ.0) GOTO 60 *如果NWE=0 转标号60处WRITE(6,40) IO *输出IO，6=显示器或打印机，40是表控格式，就是由40标号语句控制输出格式40 FORMAT(/1X,'EK NE='I5) *格式说明，似乎1X前多个/WRITE(6,50) EK *输出EK50 FORMAT(1X,6E11.4) *同上60 RETURN *返回操作系统END *程序结束EK(I1-1,J1-1)=A*(B(I)*B(J)+A2*B(I+3)*B(J+3))假如I=2,J=3I1=2*I=4,J1=2*J=6，这样上面的语句就相当于：EK(3,5)=A*(B(2)*B(3)+A2*B(5)*B(6))呵呵，就这意思。

梁单元有限元计算程序(matlab)

DOF(2)=3*i-1;
DOF(3)=3*i;
DOF(4)=3*j-2;
DOF(5)=3*j-1;
DOF(6)=3*j;
for n1=1:6
for n2=1:6
KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
end
end
z=KK;
%将该稀疏矩阵中元素放入整体刚度矩阵中的相应位置，
%使得整体刚度矩阵仍然为下三角阵形式的稀疏矩阵；
z=zபைடு நூலகம்triu(z,1)'-triu(z,1);
z=sparse(z);
[nm,nmn]=size(cod);
%输入各单元的角度
alpha=input('please input the angle (degree) of each element in order:');
%输入节点坐标，每一行代表该节点的坐标，按节点编号顺序排列，即用分号将节点分开
%用逗号将每个单元内的节点分开，按单元编号顺序排列
%如[1,2;2,3;1,3]表示三个杆中的节点编号分别为（1,2）、（2,3）、（1,3）
cod=input('please input the node of each element in order:');
%计算单元个数，nm为单元个数
%用逗号将每个节点的坐标分开，按单元编号顺序排列
%如[1,2;2,3;1,3]表示三个节点的横纵坐标分别为（1,2）、（2,3）、（1,3）
con=input('please input the coordinates (m) of each node in order:');

Fortran语言有限元程序分析报告平面钢架

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 PFFimplicit 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,nfreal 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) dimensionx(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*si s2=(c1-c4)*si*cos3=c3*sis4=c1*si*si+c4*co*co s5=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)dimensionx(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))) then do 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)dimensionx(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)do 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)enddoendifendifwrite(2,40)m,f40format(2x,i8,4x,"Ix=",f12.4,3x,"Iy=",f12.4,3x,"Mi=",f12.4/14x,"Jx=",f12.4,3x,"Jy=",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、平面深梁有限元前处理信息1、1有限元模型1、2源程序1、3输出结果及分析2、矩阵求逆子程序调试2、1 考题2、2 源程序（主程序+子程序） 2、3 子程序验证2、1、平面深梁有限元前处理信息1、1有限元模型长为4m被分为16份，宽为1m被分为4份1、2源程序DIMENSION 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)=(I-1)*4./16XY(M,2)=(J-1)/4.END DOEND DODO I=1,85WRITE(1,10) I,XY(I,1),XY(I,2)END DO10 FORMAT(1X,I4,2X,2F5.2)DO I=1,16DO J=1,4NS=(I-1)*8+JNX=(I-1)*8+4+JIJK(NS,1)=(I-1)*5+JIJK(NS,2)=I*5+J+1IJK(NS,3)=(I-1)*5+J+1IJK(NX,1)=(I-1)*5+JIJK(NX,2)=I*5+JIJK(NX,3)=I*5+J+1END DO END DO DO I=1,128WRITE(1,20) I,IJK(I,1),IJK(I,2),IJK(I,3)END DO20FORMAT(1X,I4,3X,3I4) PAUSEEND1、3输出结果及分析1、输出结果结点坐标：1 .00 .002 .00 .253 .00 .504 .00 .75 5 .00 1.006 .25 .007 .25 .258 .25 .509 .25 .75 10 .25 1.00 11 .50 .00 12 .50 .25 13 .50 .50 14 .50 .75 15 .50 1.00 16 .75 .00 17 .75 .25 18 .75 .50 19 .75 .75 20 .75 1.00 21 1.00 .0022 1.00 .25 23 1.00 .50 24 1.00 .75 25 1.00 1.00 26 1.25 .00 27 1.25 .25 28 1.25 .50 29 1.25 .75 30 1.25 1.00 31 1.50 .00 32 1.50 .25 33 1.50 .50 34 1.50 .75 35 1.50 1.00 36 1.75 .00 37 1.75 .25 38 1.75 .50 39 1.75 .75 40 1.75 1.00 41 2.00 .00 42 2.00 .2543 2.00 .5044 2.00 .7545 2.00 1.0046 2.25 .0047 2.25 .2548 2.25 .5049 2.25 .7550 2.25 1.0051 2.50 .0052 2.50 .2553 2.50 .5054 2.50 .7555 2.50 1.0056 2.75 .0057 2.75 .2558 2.75 .5059 2.75 .7560 2.75 1.0061 3.00 .0062 3.00 .2563 3.00 .5064 3.00 .75 65 3.00 1.0066 3.25 .0067 3.25 .2568 3.25 .5069 3.25 .7570 3.25 1.0071 3.50 .0072 3.50 .2573 3.50 .5074 3.50 .7575 3.50 1.0076 3.75 .0077 3.75 .2578 3.75 .5079 3.75 .7580 3.75 1.0081 4.00 .0082 4.00 .2583 4.00 .5084 4.00 .7585 4.00 1.00单元格：1 1 7 22 2 8 33 3 9 44 4 10 55 16 76 27 87 3 8 98 4 9 109 6 12 7 10 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 1522 12 17 1823 13 18 1924 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 38 56 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 4367 43 49 4468 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 6090 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 111 68 73 74 112 69 74 75 113 71 77 72 114 72 78 73 115 73 79 74 116 74 80 75 117 71 76 77 118 72 77 78 119 73 78 79 120 74 79 80 121 76 82 77 122 77 83 78 123 78 84 79 124 79 85 80 125 76 81 82 126 77 82 83 127 78 83 84 128 79 84 852、结果分析对于这次平面深梁有限元前处理信息编程，梁长4m被分为16等份，宽1m被分为4等份，用程序输出结点坐标与实际结果一致，对于每一个单元格的三个结点数，也与实际一致，充分说明可以用此程序解决平面深梁有限元前处理问题。

Fortran语言编写的有限元结构程序

算例一计算简图及结果输出用平面刚架静力计算程序下图结构的内力。

各杆EA，EI相同。

已知：642EA=4.010KN,EI=1.610KN m⨯⨯∙计算简图如下：（1）输入原始数据控制参数3,5,8,7,1,2（NE,NJ,N,NW,NPJ,NPF）结点坐标集结点未知量编号0.0,0.0,0,0 0.0,4.0,1,2,3 0.0,4.0,1,2,4 4.0,4.0,5,6,7 4.0,0.0,0,0,8单元杆端结点编号及单元EA、EI 1,2,4.0E+06,1.6E+04 3,4,4.0E+06,1.6E+04 5,4,4.0E+06,1.6E+04结点荷载7.0,-15.非结点荷载1.0,2.0,2.0,-2.0,1.0,4.0,-25.0（2）输出结果NE= 3 NJ= 5 N= 8 NW= 7 NPJ= 1 NPF= 2 NODE X Y XX YY ZZ1 0.0000 0.0000 0 0 02 0.0000 4.0000 1 2 33 0.0000 4.0000 1 2 44 4.0000 4.000056 75 4.0000 0.0000 0 0 8ELEMENT NODE-I NODE-J EA EI1 12 0.400000E+07 0.160000E+052 3 4 0.400000E+07 0.160000E+053 54 0.400000E+07 0.160000E+05CODE PX-PY-PM7. -15.0000ELEMENT IND A Q1. 2. 2.0000 -18.00002. 1. 4.0000 -25.0000NODE U V CETA1 0.000000E+00 0.000000E+00 0.000000E+002 -0.221743E-02 -0.464619E-04 -0.139404E-023 -0.221743E-02 -0.464619E-04 0.357876E-024 -0.222472E-02 -0.535381E-04 -0.298554E-025 0.000000E+00 0.000000E+00 0.658499E-03ELEMENT N Q M1 N1= 46.4619 Q1= 10.7119 M1= -6.8477N2= -46.4619 Q2= 7.2881 M2= 0.00002 N1= 7.2881 Q1= 46.4619 M1= 0.0000N2= -7.2881 Q2= 53.5381 M2= 14.15233 N1= 53.5381 Q1= 7.2881 M1= 0.0000N2= -53.5381 Q2= -7.2881 M2= -29.1523算例二计算简图及结果输出用平面刚架静力计算程序下图结构的内力。

平面梁单元有限元fortran程序

program beam2d!character*64 fname1,fname2dimension mea(100,4),jz(50,2),aai(10,2),aeu(10,2),qq(100) dimension r(6,6),rt(6,6),ra(6,6),pop(6),po(6)dimension akep(6,6),ake(6,6),p(300),x(100),y(100) dimension ld(300),is(6),a(50000),upe(6),ppe(6),ue(6)open(1,file='fname1_b.txt')open(2,file='fname2_b.txt')read(1,*) nn,ne,nm,na,ncwrite(*,*) NN,NE,NM,NA,NCnf=3nd=2nfd=nf*ndn=nn*nfdo 50 i=1,nn50 read(1,*) k,x(i),y(i),(p(nf*(i-1)+j),j=1,nf)read(1,*) ((jz(i,j),j=1,2),i=1,nc)do 100 i=1,ne100 read(1,*) ie,(mea(i,j),j=1,4),qq(i)read(1,*) ((aai(i,j),j=1,2),i=1,na)read(1,*) ((aeu(i,j),j=1,2),i=1,nm)close(1)write(2,460) nn,ne,nm,na,ncwrite(2,465) (i,x(i),y(i),p(3*i-2),p(3*i-1),p(3*i),i=1,nn)write(2,470) ((jz(i,j),j=1,2),i=1,nc)write(2,475) (i,(mea(i,j),j=1,4),i=1,ne)write(2,480) ((aai(i,j),j=1,2),i=1,na)write(2,485) ((aeu(i,j),j=1,2),i=1,nm)460 format(/5x,'the input of nn,ne,nm,na,nc'//(5x,5i5)) 465 format(/5x,'the input of x,y,p'//(5x,i5,5f10.2))470 format(/5x,'the input of jz'//(5x,i5,5x,i5))475 format(/5x,'the input of mea,qq'//(5x,i5,5x,4i5,f10.2)) 480 format(/5x,'the input of aai'//(5x,i5,5x,2e15.6))485 format(/5x,'the input of aeu'//(5x,i5,5x,2e15.6))!c structure stiffness statementcall fld(nn,ne,mea,nf,nd,n,nt,ld)do 500 i=1,nt500 a(i)=0do 600 ie=1,necall kep(nn,ie,mea,aeu,aai,x,y,akep)call cr(nn,ie,mea,x,y,r)call tran(6,6,r,rt)call dot(6,6,6,rt,akep,ra)call dot(6,6,6,ra,r,ake)call fis(ie,mea,nf,nd,nfd,is)do 560 i=1,nfddo 560 j=1,nfdif(is(i)-is(j)) 560,520,520520 ni=is(i)ij=ld(ni)-ni+is(j)a(ij)=a(ij)+ake(i,j)560 continue600 continue! c equivalent node forcedo 700 ie=1,necall fixf(nn,ne,ie,mea,x,y,qq,pop)call cr(nn,ie,mea,x,y,r)call tran(6,6,r,rt)call dot(6,6,1,rt,pop,po)call fis(ie,mea,nf,nd,nfd,is)do 650 i=1,nfdni=is(i)p(ni)=p(ni)-po(i)650 continue700 continuecall fcc(nc,n,nt,nf,jz,ld,a)call decom(n,nt,a,p,ld)write(2,850) (i,p(3*i-2),p(3*i-1),p(3*i),i=1,nn) write(2,860)do 800 ie=1,necall kep(nn,ie,mea,aeu,aai,x,y,akep)call cr(nn,ie,mea,x,y,r)call fis(ie,mea,nf,nd,nfd,is)do 750 i=1,nfdni=is(i)ue(i)=p(ni)750 continuecall dot(6,nfd,1,r,ue,upe)call dot(6,6,1,akep,upe,ppe)call fixf(nn,ne,ie,mea,x,y,qq,pop)do 760 i=1,nfdppe(i)=ppe(i)+pop(i)760 continuewrite(2,880) ie,ppe(1),ppe(2),ppe(3)write(2,900) ppe(4),ppe(5),ppe(6)800 continue850 format(//25x,'nodal displacement'//5x,'node no.',12x,'u',12x,'v',12x,'theta'/(7x,i5,3e15.6)) 860 format(//15x,'axial forces of element',/11x,'element no.','i-j',10x,'nx',12x,'ny',12x,'mz') 880 format(4x,i6,2x,'i',3e15.6)900 format(12x,'j',3e15.6)close(2)end!c is dimensionSUBROUTINE FIS(IE,MEA,NF,ND,NFD,IS)DIMENSION MEA(100,4),IS(NFD)DO 150 ID=1,NDDO 100 JF=1,NFI1=(ID-1)*NF+JFIS(I1)=(MEA(IE,ID)-1)*NF+JF100 CONTINUE150 CONTINUERETURNEND!c ld dimensionSUBROUTINE FLD(NN,NE,MEA,NF,ND,N,NT,LD)DIMENSION MEA(100,4),LD(N)LD(1)=1DO 300 K=1,NNNMIN=1000DO 200 I=1,NEDO 150 J=1,NDIF(MEA(I,J).NE.K) GOTO 150DO 100 L=1,NDIF (MEA(I,L).LT.NMIN) NMIN=MEA(I,L)100 CONTINUE150 CONTINUE200 CONTINUEDO 250 I=1,NFJ=(K-1)*NF+IIF(J.NE.1) LD(J)=LD(J-1)+(K-NMIN)*NF+I 250 CONTINUE300 CONTINUENT=LD(N)RETURNEND! c constrained conditionSUBROUTINE FCC(NC,N,NT,NF,JZ,LD,A)DIMENSION JZ(50,2),LD(N),A(NT)DO 350 K=1,NCI=JZ(K,1)J=JZ(K,2)NI=(I-1)*NF+JNJ=LD(NI)A(NJ)=1.0E25350 CONTINUERETURNEND! c linear equation by ldlt method SUBROUTINE DECOM(N,NT,A,B,LD)DIMENSION A(NT),B(N),LD(N)DO 20 I=1,NLDI=LD(I)IF (I.EQ.1) GOTO 10IO=LDI-IIM1=I-1MI=1-IO+LD(IM1)IF(MI.GT.IM1) GOTO 10DO 8 J=MI,IM1JO=LD(J)-JIJ=IO+JIF(J.EQ.1) GOTO 6JM1=J-1MJ=1-JO+LD(JM1)MIJ=MAX0(MI,MJ)IF(MIJ.GT.JM1) GOTO 6S=0.0DO 2 K=MIJ,JM1IK=IO+KJK=JO+KS=S+A(IK)*A(JK)2 CONTINUEA(IJ)=A(IJ)-S6 B(I)=B(I)-A(IJ)*B(J)8 CONTINUE10 ALDI=A(LDI)IF(I.EQ.1.OR.MI.GT.IM1) GOTO 13DO 12 J=MI,IM1IJ=IO+JLDJ=LD(J)S=A(IJ)A(IJ)=S/A(LDJ)ALDI=ALDI-S*A(IJ)12 CONTINUE13 A(LDI)=ALDIB(I)=B(I)/ALDI20 CONTINUEDO 40 LDI=2,NI=N-LDI+2IO=LD(I)-IMI=1-IO+LD(I-1)IM1=I-1IF(MI.GT.IM1) GOTO 40DO 30 J=MI,IM1IJ=IO+JB(J)=B(J)-A(IJ)*B(I)30 CONTINUE40 CONTINUERETURNENDSUBROUTINE zero(M,N,A)DIMENSION A(M,N)DO 10 I=1,MDO 10 J=1,NA(I,J)=0.010 CONTINUERETURNENDsubroutine tran(m,n,a,b)dimension a(m,n),b(n,m)do 20 i=1,mdo 20 j=1,nb(j,i)=a(i,j)20 continuereturnendsubroutine dot(m,n,l,a,b,c)! dimension a(m,n),b(n,l),c(m,l)do 50 i=1,mdo 50 j=1,lc(i,j)=0.0do 40 k=1,nc(i,j)=c(i,j)+a(i,k)*b(k,j)40 continue50 continuereturnend! c beam element(2-d) coordinate transformation subroutine cr(nn,ie,mea,x,y,r)dimension x(nn),y(nn),mea(100,4),r(6,6)call zero(6,6,r)i=mea(ie,1)j=mea(ie,2)sl=sqrt((x(j)-x(i))**2+(y(j)-y(i))**2)r(1,1)=(x(j)-x(i))/slr(1,2)=(y(j)-y(i))/slr(2,1)=-r(1,2)r(2,2)=r(1,1)r(3,3)=1.0r(4,4)=r(1,1)r(4,5)=r(1,2)r(5,4)=r(2,1)r(5,5)=r(2,2)r(6,6)=1.0returnend! c beam element(2-d)subroutine kep(nn,ie,mea,aeu,aai,x,y,akep)dimension x(nn),y(nn),mea(100,4),aeu(10,2),aai(10,2),akep(6,6) real izcall zero(6,6,akep)nm1=mea(ie,4)eo=aeu(nm1,1)na1=mea(ie,3)ao=aai(na1,1)iz=aai(na1,2)i=mea(ie,1)j=mea(ie,2)sl=sqrt((x(j)-x(i))**2+(y(j)-y(i))**2)akep(1,1)=eo*ao/slakep(4,4)=akep(1,1)akep(1,4)=-akep(1,1)akep(4,1)=akep(1,4)akep(2,2)=12.0*eo*iz/sl**3akep(2,5)=-akep(2,2)akep(5,2)=akep(2,5)akep(5,5)=akep(2,2)akep(2,3)=6.0*eo*iz/sl**2akep(3,2)=akep(2,3)akep(2,6)=akep(2,3)akep(6,2)=akep(2,6)akep(3,5)=-akep(2,3)akep(5,3)=akep(3,5)akep(5,6)=akep(5,3)akep(6,5)=akep(5,6)akep(3,3)=4.0*eo*iz/slakep(6,6)=akep(3,3)akep(3,6)=2.0*eo*iz/slakep(6,3)=akep(3,6)returnend! c fixed and forcesubroutine fixf(nn,ne,ie,mea,x,y,qq,pop)dimension x(nn),y(nn),qq(ne),mea(100,4),pop(6)i=mea(ie,1)j=mea(ie,2)sl=sqrt((x(j)-x(i))**2+(y(j)-y(i))**2)pop(1)=0.0pop(4)=0.0pop(2)=-qq(ie)*sl/2.0 pop(5)=pop(2)pop(3)=-qq(ie)*sl*sl/12.0 pop(6)=-pop(3)returnend4 3 1 1 51 0 0 0 0 02 3 0 0 -200 03 6 0 0 0 -504 12 0 0 0 01 11 21 33 24 21 12 1 1 02 23 1 1 03 34 1 1 -25.08.611E-3 2.17e-42.0E8 0.3。

平面梁单元MATLAB有限元程序

平面梁单元MATLAB有限元程序function []=beam2n()clear allclose allclearclose%------------------------ BEAM2 ---------------------------disp('========================================'); disp(' PROGRAM BEAM2 '); disp(' Beam Bending Analysis '); disp(' T.R.Chandrupatla and A.D.Belegundu '); disp('========================================');InputData;Bandwidth;Stiffness;ModifyForBC;BandSolver;ReactionCalc;Output;%------------------------ function InputData ---------------------------function [] = InputData();global NN NE NM NDIM NEN NDNglobal ND NL NCH NPR NMPC NBWglobal X NOC F AREA MAT SMI Sglobal PM NU U MPC BT REACTglobal CNSTglobal TITLE FILE1 FILE2global LINP LOUTglobal NQdisp(blanks(1));FILE1 = input('Input Data File Name ','s'); LINP = fopen(FILE1,'r');FILE2 = input('Output Data File Name ','s'); LOUT = fopen(FILE2,'w');DUMMY = fgets(LINP);TITLE = fgets(LINP);DUMMY = fgets(LINP);TMP = str2num(fgets(LINP));[NN, NE, NM, NDIM, NEN, NDN] =deal(TMP(1),TMP(2),TMP(3),TMP(4),TMP(5),TMP(6));NQ = NDN * NN;DUMMY = fgets(LINP);TMP = str2num(fgets(LINP));[ND, NL, NMPC]= deal(TMP(1),TMP(2),TMP(3));NPR=1; % E%----- Coordinates -----DUMMY = fgets(LINP);for I=1:NNTMP = str2num(fgets(LINP));[N, X(N,:)]=deal(TMP(1),TMP(2:1+NDIM)); end%----- Connectivity -----DUMMY = fgets(LINP);for I=1:NETMP = str2num(fgets(LINP));[N,NOC(N,:), MAT(N), SMI(N)] = ...deal(TMP(1),TMP(2:1+NEN), TMP(2+NEN), TMP(3+NEN));end%----- Specified Displacements ----- DUMMY = fgets(LINP); for I=1:NDTMP = str2num(fgets(LINP));[NU(I,:),U(I,:)] = deal(TMP(1), TMP(2)); end%----- Component Loads -----DUMMY = fgets(LINP);F = zeros(NQ,1);for I=1:NLTMP = str2num(fgets(LINP));[N,F(N)]=deal(TMP(1),TMP(2)); end%----- Material Properties ----- DUMMY = fgets(LINP);for I=1:NMTMP = str2num(fgets(LINP));[N, PM(N,:)] = deal(TMP(1), TMP(2:NPR+1));end%----- Multi-point Constraints B1*Qi+B2*Qj=B0if NMPC > 0DUMMY = fgets(LINP);for I=1:NMPCTMP = str2num(fgets(LINP));[BT(I,1), MPC(I,1), BT(I,2), MPC(I,2), BT(I,3)] = ...deal(TMP(1),TMP(2),TMP(3),TMP(4),TMP(5));endendfclose(LINP);%------------------------ function Bandwidth ---------------------------function []=Bandwidth();global NN NE NM NDIM NEN NDNglobal ND NL NCH NPR NMPC NBW global X NOC F AREA MAT SMI S global PM NU U MPC BT STRESS REACT global CNSTglobal TITLE FILE1 FILE2global LINP LOUT%----- Bandwidth Evaluation ----- NBW = 0;for N=1:NENABS = NDN*(abs(NOC(N, 1) - NOC(N, 2)) + 1);if (NBW < NABS)NBW = NABS;endendfor I=1:NMPCNABS = abs(MPC(I, 1) - MPC(I, 2)) + 1;if (NBW < NABS)NBW = NABS;endenddisp(sprintf('Bandwidth = %d', NBW));%------------------------ function Stiffness ---------------------------function []=Stiffness();global NN NE NM NDIM NEN NDNglobal ND NL NCH NPR NMPC NBW global X NOC F AREA MAT SMI S global PM NU U MPC BT REACTglobal CNSTglobal NQ%----- Global Stiffness Matrix S = zeros(NQ,NBW);for N = 1:NEdisp(sprintf('Forming Stiffness Matrix of Element %d', N));%-------- Element Stiffness and Temperature Load -----N1 = NOC(N, 1);N2 = NOC(N, 2);M = MAT(N);EL = abs(X(N1) - X(N2));EIL = PM(M, 1) * SMI(N) / EL^3;SE(1, 1) = 12 * EIL;SE(1, 2) = EIL * 6 * EL;SE(1, 3) = -12 * EIL;SE(1, 4) = EIL * 6 * EL;SE(2, 1) = SE(1, 2);SE(2, 2) = EIL * 4 * EL * EL;SE(2, 3) = -EIL * 6 * EL;SE(2, 4) = EIL * 2 * EL * EL;SE(3, 1) = SE(1, 3);SE(3, 2) = SE(2, 3);SE(3, 3) = EIL * 12;SE(3, 4) = -EIL * 6 * EL;SE(4, 1) = SE(1, 4);SE(4, 2) = SE(2, 4);SE(4, 3) = SE(3, 4);SE(4, 4) = EIL * 4 * EL * EL;disp('.... Placing in Global Locations'); for II = 1:NENNRT = NDN * (NOC(N, II) - 1);for IT = 1:NDNNR = NRT + IT;I = NDN * (II - 1) + IT;for JJ = 1:NENNCT = NDN * (NOC(N, JJ) - 1);for JT = 1:NDNJ = NDN * (JJ - 1) + JT;NC = NCT + JT - NR + 1;if (NC > 0)S(NR, NC) = S(NR, NC) + SE(I, J);endendendendendend%------------------------ function ModifyForBC ---------------------------function []=ModifyForBC();global NN NE NM NDIM NEN NDNglobal ND NL NCH NPR NMPC NBWglobal X NOC F AREA MAT SMI Sglobal PM NU U MPC BT REACTglobal CNSTglobal NQ%----- Decide Penalty Parameter CNST ----- CNST = 0;for I = 1:NQif CNST < S(I, 1); CNST = S(I, 1); end endCNST = CNST * 10000;%----- Modify for Boundary Conditions ----- % --- Displacement BC ---for I = 1:NDN = NU(I);S(N, 1) = S(N, 1) + CNST;F(N) = F(N) + CNST * U(I);end%--- Multi-point Constraints ---for I = 1:NMPCI1 = MPC(I, 1);I2 = MPC(I, 2);S(I1, 1) = S(I1, 1) + CNST * BT(I, 1) * BT(I, 1);S(I2, 1) = S(I2, 1) + CNST * BT(I, 2) * BT(I, 2);IR = I1;if IR > I2; IR = I2; endIC = abs(I2 - I1) + 1;S(IR, IC) = S(IR, IC) + CNST * BT(I, 1) * BT(I, 2);F(I1) = F(I1) + CNST * BT(I, 1) * BT(I, 3);F(I2) = F(I2) + CNST * BT(I, 2) * BT(I, 3); end%------------------------ function BandSolver ---------------------------function []=BandSolver();global NN NE NM NDIM NEN NDNglobal ND NL NCH NPR NMPC NBWglobal X NOC F AREA MAT SMI Sglobal PM NU U MPC BT REACTglobal CNSTglobal NQ%----- Equation Solving using Band Solver ----- disp('Solving using Band Solver(bansol.m)');[F] = bansol(NQ,NBW,S,F);%------------------------ function ReactionCalc ---------------------------function []=ReactionCalc();global NN NE NM NDIM NEN NDNglobal ND NL NCH NPR NMPC NBWglobal X NOC F AREA MAT SMI Sglobal PM NU U MPC BT REACTglobal CNSTfor I = 1:NDN = NU(I);REACT(I) = CNST * (U(I) - F(N));end%------------------------ function Output ---------------------------function []=Output();global NN NE NM NDIM NEN NDNglobal ND NL NCH NPR NMPC NBWglobal X NOC F AREA MAT SMI Sglobal PM NU U MPC BT REACTglobal CNSTglobal TITLE FILE1 FILE2global LINP LOUTdisp(sprintf('Output for Input Data from file %s\n',FILE1)); fprintf(LOUT,'Output for Input Data from file %s\n',FILE1);disp(TITLE);fprintf(LOUT,'%s\n',TITLE);disp(' Node# X-Displ Rotation');fprintf(LOUT,' Node# X-Displ Rotation\n'); I=[1:NN]';% print a matrixdisp(sprintf(' %4d %15.4E %15.4E\n',[I,F(2*I-1),F(2*I)]')); fprintf(LOUT,' %4d %15.4E %15.4E\n',[I,F(2*I-1),F(2*I)]');%----- Reaction Calculation -----disp(sprintf(' DOF# Reaction'));fprintf(LOUT,' DOF# Reaction\n');for I = 1:NDN = NU(I);R = CNST * (U(I) - F(N));disp(sprintf(' %4d %15.4E',N,REACT(I)));fprintf(LOUT,' %4d %15.4E\n',N,REACT(I)); endfclose(LOUT);disp(sprintf('The Results are available in the text file %s', FILE2));。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。