有限元编程算例(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平面钢架有限元分析

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语言-有限元程序分析-平面钢架

程序框图：程序特点：问题类型：可用于计算结构力学的平面刚架问题单元类型：直接利用杆单元载荷类型：节点载荷及非节点载荷，其中非节点载荷包括均布荷载和垂直于杆件的集中荷载材料性质：所有杆单元几何性质相同，且由相同的均匀材料组成方程求解：结构刚度矩阵采用满阵存放，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.ALLOCATABLE::COORD(:,:),PROPS(:,:,:) !声明两个可变大小的数组，COORD(:,:)是二维数组，PROPS(:,:,:)是三维数组。

2.Fortran程序行首为C代表改行为注释，不会被编译3.全局变量（common），不同的程序之间，也就是在不同的函数之间或者是主程序跟函数之间，除了可以通过传递参数的方法来共享内存，还可以通过“全局变量”来让不同程序中声明出来的变量使用相同的内存位置。

4.Dimensional维的，viscoplastic塑性的，elastic有弹力的，finite有限的，element元素，program程序。

5.THREE DIMENSIONAL ELASTIC-VISCOPLASTIC FINITE ELEMENT PROGRAM三维弹塑性有限元程序6.Module可以用来封装程序模块，通常是用来把程序中，具备相关功能的函数及变量封装在一起。

程序在开始定义了一个module模块，在模块中定义了MXKKK=，MXGSJ=1000，MXGSJ=1000三个常量（PARAMETER表示常量），并且每个常量都赋了值。

在module模块中定义了NELEM,NPOIN,NPROP,MXDFN,NSTEP,IDEVP,IDDP,LTYPE以及NFIX1,NPL,NVL,NSL,NHL,NTL,IDCVG,NTOTV,NKK以及DTIME,TOLER,SCALE,DSCALE这些全局变量（common表示全局变量），定义了ICM(3,8),CGAUS(2),VSHAP(8,8),DERIV(3,8,8)以及POSGP(3),COPG(3),EJ(3,3),EJACI(3,3),R(8,8)这些维数与大小都确定的全局数组变量，定义了COORD(:,:),PROPS(:,:,:)以及STRSG(:,:,:),DJ(:,:),CARTD(:,:,:,:)以及TRANJ(:,:,:,:),DJRMX(:,:,:)以及DREMX(:,:,:),DJEMX(:,:,:,:)以及CREMX(:,:,:),CJEMX(:,:,:,:)以及MELEM(:,:),MPROP(:),ISSOR(:,:),NNDEX(:)以及MPFIX(:,:),MPSJ(:),MMATP(:),MPIV(:)以及TSTIF(:)以及ADISP(:),TDISP(:),ALOAD(:)以及PSNBR(:,:,:),PSNBJ(:,:)以及PSTNR(:,:,:),PSTNJ(:,:)以及STRSP(:,:),STRSJ(:,:)这些维数确定但是大小不确定的可变大小的数组，ALLOCATABLE表示可变大小的数组变量。

读有限元Fortran程序笔记

读有限元程序笔记1.ALLOCATABLE::COORD(:,:),PROPS(:,:,:) !声明两个可变大小的数组，COORD(:,:)是二维数组，PROPS(:,:,:)是三维数组。

4.Dimensional维的，viscoplastic塑性的，elastic有弹力的，finite有限的，element元素，program程序。

程序在开始定义了一个module模块，在模块中定义了MXKKK=50000000，MXGSJ=1000，MXGSJ=1000三个常量（PARAMETER表示常量），并且每个常量都赋了值。

Fortran---平面四边形八结点等参元有限元程序

% -------------------------------------------------------------------------% 这是一个平面四边形八结点等参元有限元主程序，本程序具有以下功能：% 1,本程序能计算平面应力问题。

% 2，自动画出变形图和荷载图。

% 3，输出节点位移，节点力,节点应力，和节点主应力。

% 时间：2004年12月%--------------------------------------------------------------------------function N8plane% 变量说明列表:% NfreeNum 节点自由度数(本程序为2);EleNodeNum 单元节点数(本程序为8) % ---------------------------单变量---------------------------------------% FP1 数据文件指针(输入);NPOIN 节点总数;NELEM 单元总数;NVFIX 约束自由度数;% E 弹性模量;POISSON 泊松比;h 单元厚度;method 问题类型(1.单元应力 2.单元应变)% ----------------------------数组----------------------------------------% LNODS 单元定义;COORD 坐标;FPOIN 点荷载定义;k 单刚矩阵;K 总刚矩阵% Oxy 单元中点应力;O12单元中点主应力% ----------------------------矢量----------------------------------------% FIXED 约束信息;NFORCE 荷载信息;ASLOD 等效节点荷载;Ecoord 单元坐标% DISP 总体位移;EDisp 单元位移% ------------------------------------------------------------------------% 主调函数:% N8ek 生成单刚N8Assemble 组装总刚;N8stress 求解应力;% ------------------------------------------------------------------------%%%clear%-----输入数据存储的文件名------------file_in = input('数据输入的文件: ', 's' ) ;file_out = input( '数据输出的文件: ', 's') ;%--------------------------------------------------------------------------% 打开文件读取数据%--------------------------------------------------------------------------FP1=fopen(file_in,'r');NfreeNum=2;EleNodeNum=8;NPOIN=fscanf(FP1,'%d',1); NELEM=fscanf(FP1,'%d',1); NVFIX=fscanf(FP1,'%d',1);E=fscanf(FP1,'%f',1);POISSON=fscanf(FP1,'%f',1);h=fscanf(FP1,'%f',1);method=fscanf(FP1,'%d' ,1);%-------读取结构信息--------LNODS=fscanf(FP1,'%f',[8,NELEM])'; %单元定义COORD=fscanf(FP1,'%f',[2,NPOIN])'; %坐标FIXED=fscanf(FP1,'%f',NVFIX)'; %约束信息NFORCE=fscanf(FP1,'%d',2)'; %荷载信息if NFORCE(1)~=0FPOIN=fscanf(FP1,'%f',[3,NFORCE(1)])'; %节点力end%--------------定义非结点力-------------------------------------------------if NFORCE(2)~=0NF=NFORCE(1);np_number=fscanf( FP1, '%d', 1 ) ; % 读入均布侧压的单元边数PBorder=zeros(np_number,6);for i=1:np_numberPBorder(i,1)=fscanf( FP1, '%d', 1 ) ; % 单元编号PBorder(i,2)=fscanf( FP1, '%d', 1 ) ; % 侧压边左结点号PBorder(i,3)=fscanf( FP1, '%d', 1 ) ; % 侧压边中结点号PBorder(i,4)=fscanf( FP1, '%d', 1 ) ; % 侧压边右结点号PBorder(i,5)=fscanf( FP1, '%f', 1 ) ; % 侧压力的大小PBorder(i,6)=fscanf( FP1, '%d', 1 ) ; % 侧压力的方向：1--x方向，2--y方向%--------------------------------------------------------------------------% 以下语句将均布荷载转化为集中力%--------------------------------------------------------------------------if PBorder(i,6)==1 % 判断均布力的方向1--x方向，2--y方向xi=COORD(PBorder(i,2),1); % 提取侧压边两端点坐标yi=COORD(PBorder(i,2),2);xj=COORD(PBorder(i,4),1);yj=COORD(PBorder(i,4),2);Blength=sqrt((xj-xi)^2+(yj-yi)^2); % 边界长度FPOIN(NF+1,1)=PBorder(i,2); % 给侧压力下端点分配力，并输入到节点力矩阵中FPOIN(NF+1,2)=PBorder(i,5)*Blength*h/6;FPOIN(NF+1,3)=0;FPOIN(NF+2,1)=PBorder(i,3); % 给侧压力中端点分配力，并输入到节点力矩阵中FPOIN(NF+2,2)=PBorder(i,5)*Blength*h*2/3;FPOIN(NF+2,3)=0;FPOIN(NF+3,1)=PBorder(i,4); % 给侧压力上端点分配力，并输入到节点力矩阵中FPOIN(NF+3,2)=PBorder(i,5)*Blength*h/6;FPOIN(NF+3,3)=0;NF=NF+3;elsexi=COORD(PBorder(i,2),1); % 提取侧压边两端点坐标，为计算边界长度作准备yi=COORD(PBorder(i,2),2);xj=COORD(PBorder(i,4),1);yj=COORD(PBorder(i,4),2);Blength=sqrt((xj-xi)^2+(yj-yi)^2); % 边界长度FPOIN(NF+1,1)=PBorder(i,2); % 给侧压力左端点分配力，并输入到节点力矩阵中FPOIN(NF+1,2)=0;FPOIN(NF+1,3)=PBorder(i,5)*Blength*h/6;FPOIN(NF+2,1)=PBorder(i,3); % 给侧压力中端点分配力，并输入到节点力矩阵中FPOIN(NF+2,2)=0;FPOIN(NF+2,3)=PBorder(i,5)*Blength*h*2/3;FPOIN(NF+3,1)=PBorder(i,3); % 给侧压力右端点分配力，并输入到节点力矩阵中FPOIN(NF+3,2)=0;FPOIN(NF+3,3)=PBorder(i,5)*Blength*h/6;NF=NF+3; % NF为转化后集中力的个数endendend%--------------------------------------------------------------------------% 均布荷载转化完毕%--------------------------------------------------------------------------fclose(FP1);K=zeros(2*NPOIN,2*NPOIN);%-----生成单刚,组装总刚-----for i=1:NELEMfor j=1:EleNodeNumEcoord((j-1)*NfreeNum+1:j*NfreeNum)=COORD(LNODS(i,j),1:NfreeNum); % 提取i单元各节点的坐标endk=N8ek(E,POISSON,h,Ecoord,method); % 高斯积分求解单刚K=N8Assemble(K,k,LNODS,NfreeNum,EleNodeNum,i); % 叠加总刚if k~=k'error('单刚矩阵不对称，请检查') % 判断刚度对称性endendif K~=K'error('总刚矩阵不对称，请检查')end%-------生成荷载向量---------ASLOD(1:2*NPOIN)=0;ASLOD1(1:2*NPOIN)=0;for i=1:NF % NF为转化后集中力的个数ASLOD1((FPOIN(i,1)*2-1):FPOIN(i,1)*2)=FPOIN(i,2:3);ASLOD=ASLOD+ASLOD1; % 实现相同节点力的叠加ASLOD1(1:2*NPOIN)=0;%NPOIN(i,1)为作用点,NPOIN(i,2:3)为x,y方向的节点力end%----------------------------------------------K1=K;%-----------加约束-----------for j=1:NVFIXa=FIXED(j);K1(1:NfreeNum*NPOIN,a)=0; K1(a,1:NfreeNum*NPOIN)=0; K1(a,a)=1;ASLOD(a)=0;endDISP=K1\ASLOD'; % 求出节点位移%-----------------Fnode=K*DISP; % 求出节点力%------求单元节点应力--------EDisp=zeros(16,1);sigma=zeros(NPOIN,3);nodetimes=zeros(NPOIN); % 统计节点求解的次数,因不同单元可能共用某个节点for i=1:NELEMfor j=1:EleNodeNumEDisp(j*2-1:j*2)=DISP(2*LNODS(i,j)-1:2*LNODS(i,j));endfor j=1:EleNodeNumEcoord((j-1)*NfreeNum+1:j*NfreeNum)=COORD(LNODS(i,j),1:NfreeNum);endEnodesigma=N8stress(E,POISSON,Ecoord,EDisp,method); % 求解i单元的节点应力for j=1:EleNodeNumnum=LNODS(i,j);nodetimes(num)=nodetimes(num)+1;if nodetimes(num)==0sigma(num,1:3)=Enodesigma(j,1:3);elsesigma(num,1:3)=sigma(num,1:3)+Enodesigma(j,1:3); % 公共节点处应力叠加endendendfor i=1:NPOINsigma(i,1:3)=sigma(i,1:3)/nodetimes(i); % i节点应力mainsigma(i,1:3)=mainstress2d(sigma(i,1:3)); % i节点主应力end%-------------------------------------------------------------------------%%%--------------------------------------------------------------------------% 程序输出部分%--------------------------------------------------------------------------%FP2 = fopen( file_out, 'w' ); % 打开输出文件fprintf( FP2,'----------------------------------------------------------------------------------------\n');fprintf( FP2,'这是一个平面四边形八结点等参元有限元主程序，本程序具有以下功能：\n'); fprintf( FP2,' 1,本程序能计算平面应力问题。

有限元方法(有FORTRAN程序)飞箭系列

§4.2 元件程序的结构………………………………………………………81
4.2.1 START 元件程序…………………………………………………………81 4.2.1.1 功能………………………………………………………………82 4.2.1.2 命令行参数说明…………………………………………………82 4.2.1.3 参数及数组说明…………………………………………………83
2.2.1 导数之间的变换………………………………………………………46
§2.3 数值积分………………………………………………………………49
2.3.1 高斯积分……………………………………………………………49 2.3.2 节点积分……………………………………………………………51
第三章有限元输入数据形式
Γ 0
图 1.1.1
x
要求解的未知函数 u 可以是标量函数场(例如温度)，也可以是几个变量组成的向量函数场(例如位移、应变、应力等)。 A, B 是表示对于独立变量(例如空间坐标、时间坐标等)的微分算子。偏微分方程数应和未知场函数的数目相对应，因此，上述偏微分方程可以是单个的方程，也可以是一组方程。所以在式(1.1.1) 和(1.1.2)中采用了矩阵形式。对于工程或物理学中遇到的偏微分方程一般是没有理论解的，即未知函数没有解析表达式。而工程上又需要了解这些未知函数，所以我们一般用数值的方法来求解这些偏微分方程。有限元方法就是一种数值求解偏微分方程的方法，它实
§1.4 求解解梯度的最小二乘法………………………………………27
1.4.1 已知位移求应力………………………………………………………28 1.4.2 已知温度求热流密度…………………………………………………29 1.4.3 已知电势求电场强度…………………………………………………30

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

(完整word版)FORTRAN经典入门程序20例

对于FORTRA的初学者。

这些例子可作为小练习1.例题：计算工资问题。

每小时工资为RATE如果工作超过40小时，加班呢部分工资是正常时间工资的1.5倍。

C Payroll with overtimeprogram payrollreal rate, hours, payread (*,*) rate, hoursif (hours>40) thenpay=40*rate+(hours-40)*1.5*rateelsepay=hours*rateEND IFprint *, "rate=" , rateprint *, "hours=" , hoursprint *, "pay=" ,payend2•学生成绩问题。

大于80为A级。

大于60小于80为B级。

小于60为C级。

IF的嵌套。

注意空格可以看清楚else if ,e nd if,pri nt 的内容•PROGRAGRADESTUDENTREA0GRADEIF (GRADE .GE. 80) THENPRINT*,GRADE, "GRADE=>A"ELSEIF (GRADELT.60) THENPRINT*,GRADE"GRADE=>C"ELSEPRINT*,GRADE"GRADE=>B"END IFEND IFEND3. 三个数按从小到大排序。

PROGRA M AXMINREALA,B,C,TREA0A,B,CIF (A.GT.B) THENT=AA=BB=TELSEEND IFIF (B.GT.C) THENT=BB=CC=TELSEEND IFIF (A.GT.B) THENT=AA=BB=TEND IFPRINT*,A,B,CEND4. 运用EISE IF语句。

重做例子2PROGRAM2READ*,*) GRADEIF (GRADE .GE. 80.0) THENPRINT*, GRADE, "=>A"ELSE IF(GRADE .GE. 70.0) THENPRINT*, GRADE, "=>B"ELSE IF(GARDE .GE. 60.0) THENPRINT*, GRADE, "=>C"ELSEPRINT*, GARDE, "=>D"END IFEND3x 6,x 05. 计算y 2x 2x 8,x 0PROGRAEQUATIONREAD*,*) XIF (X .GE. 0.0) Y=3*X+6IF (X .LT. 0.0) Y=-X**2+2*X-8PRINT*, "X=" ,X, "Y=" ,YEND6. CONTINUED句。

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语言有限元程序分析报告平面钢架

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算例程序

题目：《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源代码

作业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语言编写的有限元结构程序

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

各杆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-03 ELEMENT 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 讲解-平面有限元源程序

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))呵呵，就这意思。