平面梁单元有限元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有限元程序解决。

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

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。