下载文档原格式
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。
三角形常应变单元程序的编制与使用有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。
有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。
对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好的适应性，而矩形或四边形单元较三节点三角形有更高的计算精度。
Matlab语言是进行矩阵运算的强大工具，因此，用Matlab语言编写有限元中平面问题的程序有优越性。
本章将详细介绍如何利用Matlab语言编制三角形常应变单元的计算程序，程序流程图见图1。
有限元法中三节点三角形分析结构的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2）单元分析，建立单元刚度矩阵。
3）整体分析，建立总刚矩阵。
4）建立整体结构的等效节点荷载和总荷载矩阵5）边界条件处理。
6）解方程，求出节点位移。
7）求出各单元的单元应力。
8）计算结果整理。
计算结果整理包括位移和应力两个方面；位移计算结果一般不需要特别的处理，利用计算出的节点位移分量，就可画出结构任意方向的位移云图；而应力解的误差表现在单元内部不满足平衡方程，单元与单元边界处应力一般不连续，在边界上应力解一般与力的边界条件不相符合。
图1 程序流程图1.1 程序说明% 三角形常应变单元求解结构主程序●功能：运用有限元法中三角形常应变单元解平面问题的计算主程序。
●基本思想：单元结点按右手法则顺序编号。
●荷载类型：可计算结点荷载。
●说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。
1 程序准备format short e %设定输出类型clear all %清除所有已定义变量clc %清屏●说明：format short e －设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法；clear all －清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况；clc －清屏，使屏幕在本程序运行开始时2 全局变量定义global NNODE NPION NELEM NVFIX NFORCE COORD LNODS YOUNG POISS THICKglobal FORCE FIXEDglobal BMATX DMATX SMATX AREAglobal ASTIF ASLOD ASDISPglobal FP1●说明：NNODE—单元结点数，NPION—总结点数，NELEM—单元数，NVFIX—受约束边界点数，NFORCE—结点力数，COORD—结构结点坐标数组，LNODS —单元定义数组，YOUNG—弹性模量，POISS—泊松比，THICK—厚度FORCE —节点力数组(n,3) n:受力节点数目,(n,1):作用点,(n,2):x方向,(n,3):y 方向；FIXED—约束信息数组(n,3) n:受约束节点数目, (n,1):约束点(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0BMATX—单元应变矩阵(3*6)，DMATX—单元弹性矩阵(3*3)，SMATX—单元应力矩阵(3*6)，AREA—单元面积ASTIF—总体刚度矩阵，ASLOD—总体荷载向量，ASDISP—结点位移向量FP1—数据文件指针3 打开文件FP1=fopen('input.txt','rt'); %打开输入数据文件存放初始数据●说明：FP1=fopen('input.txt','rt'); －打开已存在的输入数据文件input.txt，且设置其为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1来执行FP2=fopen('output.txt','wt'); －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。
目录程序一：平面刚架静力分析程序（PF.FOR） (17)程序二：平面三节点有限元程序 (17)程序三：四节点矩形薄板单元程序 (24)程序一：平面刚架静力分析程序（PF.FOR）已知各杆截面均为矩形，柱截面宽0.4m,高0.5m, 梁截面宽0.4m,高0.4m，各杆E=3.65×104 MPa。
图2节点、单元编号如下图3，1.2.3…..为节点号，①②③……为单元号：图3总共有13个节点，13个单元。
计算源程序如下：! PF.FOR (A program for analysis of plane frame)! Version 6.3 2004! Main program reads the control data & organizes the whole ! calculation by calling subroutines.DIMENSION W(80000)CHARACTER IDFN*20,TITLE(5)*72READ(*,'(A12)')IDFNOPEN(3,FILE=IDFN,STATUS='OLD')READ(3,'(A72)')(TITLE(M),M=1,5)READ(3,*)E,NM,NJ,NS,NLCL1=1L2=L1+NML3=L2+NML4=L3+NML11=L4+NML12=L11+NJL21=L12+NJL22=L21+NSL31=L22+NSL41=L31+6*NMCALL IOMJS(TITLE,E,NM,NJ,NS,NLC,W(L1),W(L2),W(L3),&W(L4),W(L11),W(L12),W(L21),W(L22))CALL LCVCT (NM,W(L1),W(L2),W(L31),NJ,N)CALL LCDIA (NM,N,W(L31),W(L41),W(L41),W(L41),MAXBDW,NA) L51=L41+NL52=L51+36L53=L52+NA*2L54=L53L61=L54+N*2NW=L61+6*NM-1WRITE (*,1)NA,NW1 FORMAT(/40X,'( NA=',I6,' )',/40X,'( NW=',I6,' )')CALL FORMA (E,NM,NJ,N,NA,W(L1),W(L2),W(L3),W(L4),&W(L11),W(L12),W(L31),W(L51),W(L41),W(L52)) CALL AS(NS,N,NA,W(L21),W(L41),W(L52))CALL LDLT(N,NA,W(L41),W(L52),W(L53))DO 100 LC=1,NLCREAD (3,*)NLJL62=L61+NLJL63=L62+NLJL64=L63+NLJL71=L61L81=L71+6*NMCALL B0(LC,N,NLJ,W(L54))IF(NLJ.NE.0) CALL IOLJB(N,NLJ,W(L61),W(L62),&W(L63),W(L64),W(L54))READ(3,*)NLML82=L81+NLML83=L82+NLML84=L83+NLMCALL F0(NLM,NM,W(L71))IF(NLM.NE.0) CALL IOLMFB(NM,NJ,N,NLM,W(L81),&W(L82),W(L83),W(L84),W(L1),W(L2),W(L11),&W(L12),W(L31),W(L71),W(L54))CALL BS(NS,N,W(L21),W(L22),W(L54))CALL SLVEQ(N,NA,MAXBDW,W(L41),W(L52),W(L54))CALL OJD(NJ,N,W(L54))CALL COTF(E,NM,NJ,N,W(L1),W(L2),W(L3),W(L4),&W(L11),W(L12),W(L31),W(L54),W(L71)) NW=L84+NLM-1WRITE(*,1)NA,NW100 CONTINUEWRITE(*,'(/)')ENDSUBROUTINE IOMJS(TITLE,E,NM,NJ,NS,NLC,IST,IEN,&AR,RI,X,Y,IS,VS)! Read data of members, joints, supports & print them DIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),&X(NJ),Y(NJ),IS(NS),VS(NS)CHARACTER TITLE(5)*72WRITE(*,'(/)')WRITE(*,1)(TITLE(M),M=1,5)1 FORMAT(1X,A72)WRITE(*,2)E,NM,NJ,NS,NLC2 FORMAT(/13X,'The Input Data'&//10X,'The General Information'&//6X,'E',9X,'NM',5X,'NJ',5X,'NS',5X,'NLC'& /1X,1PE10.3,4I7)READ(3,*)(IST(M),IEN(M),AR(M),RI(M),M=1,NM)WRITE(*,3)3 FORMAT(/10X,'The Information of Members'&//1X,'member',2X,'start',2X,'end',9X,'A',15X,'I') WRITE(*,4)(M,IST(M),IEN(M),AR(M),RI(M),M=1,NM)4 FORMAT(1X,I4,I8,I6,1P2E16.6)READ(3,*)(X(M),Y(M),M=1,NJ)WRITE(*,5)5 FORMAT(/10X,'The Joint Coordinates'&//1X,'joint',11X,'X',17X,'Y')WRITE(*,6)(M,X(M),Y(M),M=1,NJ)6 FORMAT(1X,I4,2F18.6)READ (3,*)(IS(M),VS(M),M=1,NS)WRITE(*,7)7 FORMAT(/10X,'The Information of Supports'&//4X,'IS',9X,'VS')WRITE(*,8)(IS(M),VS(M),M=1,NS)8 FORMAT(1X,I5,F16.6)RETURNENDSUBROUTINE LCVCT(NM,IST,IEN,LV,NJ,N)! Determine location vector of elementDIMENSION IST(NM),IEN(NM),LV(6,NM)DO 100 M=1,NMI=IST(M)*3J=IEN(M)*3LV(1,M)=I-2LV(2,M)=I-1LV(3,M)=ILV(4,M)=J-2LV(5,M)=J-1LV(6,M)=J100 CONTINUEN=NJ*3RETURNENDSUBROUTINE LCDIA(NM,N,LV,MIN,IBDW,LD,MAXBDW,NA)! Determine location of diagonal elements of global stiffness ! matrix ADIMENSION LV(6,NM),MIN(N),IBDW(N),LD(N)DO 100 I=1,NMIN(I)=I100 CONTINUEDO 400 M=1,NMMINLV=LV(1,M)DO 200 I=2,6IF (LV(I,M).LT.MINLV) MINLV=LV(I,M)200 CONTINUEDO 300 I=1,6IF (MINLV.LT.MIN(LV(I,M))) MIN(LV(I,M))=MINLV300 CONTINUE400 CONTINUEMAXBDW=0DO 500 I=1,NIBDW(I)=I-MIN(I)+1IF(IBDW(I).GT.MAXBDW) MAXBDW=IBDW(I)500 CONTINUELD(1)=IBDW(1)DO 600 I=2,NLD(I)=LD(I-1)+IBDW(I)600 CONTINUENA=LD(N)RETURNENDSUBROUTINE RLCS(M,NM,NJ,IST,IEN,X,Y,RL,C,S)! Calculate length, cosine & sine of member DIMENSION IST(NM),IEN(NM),X(NJ),Y(NJ)I=IST(M)J=IEN(M)X1=X(J)-X(I)Y1=Y(J)-Y(I)RL=SQRT(X1*X1+Y1*Y1)C=X1/RLS=Y1/RLRETURNENDSUBROUTINE KEBAR (M,E,NM,NJ,IST,IEN,AR,RI,&X,Y,C,S,E1,E2,E3,E4)! Calculate element stiffness matrix along local axes DIMENSION IST(NM),IEN(NM),X(NJ),Y(NJ),AR(NM),RI(NM) CALL RLCS (M,NM,NJ,IST,IEN,X,Y,RL,C,S)E1=E*AR(M)/RLE2=12.0*E*RI(M)/(RL*RL*RL)E3=0.5*E2*RLE4=0.6666667*E3*RLRETURNENDSUBROUTINE KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE)! Calculate element stiffness matrix along global axes DIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),&X(NJ),Y(NJ),AE(6,6)CALL KEBAR (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4) A1=E1*C*C+E2*S*SA2=(E1-E2)*C*SA3=E1*S*S+E2*C*CA4=E3*SA5=E3*CA6=E4AE(1,1)=A1AE(2,1)=A2AE(2,2)=A3AE(3,1)=-A4AE(3,2)=A5AE(3,3)=A6AE(4,1)=-A1AE(4,2)=-A2AE(4,3)=A4AE(4,4)=A1AE(5,1)=-A2AE(5,2)=-A3AE(5,3)=-A5AE(5,4)=A2AE(5,5)=A3AE(6,1)=-A4AE(6,2)=A5AE(6,3)=0.5*A6AE(6,4)=A4AE(6,5)=-A5AE(6,6)=A6RETURNENDSUBROUTINE FORMA (E,NM,NJ,N,NA,IST,IEN,AR,RI,&X,Y,LV,AE,LD,A)! Form global stiffness matrix ADIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ),&LV(6,NM),AE(6,6),LD(N)DOUBLE PRECISION A(NA)DO 300 M=1,NMCALL KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE)DO 200 I=1,6DO 100 J=1,IIF (LV(I,M).GE.LV(J,M)) THENA(LD(LV(I,M))-LV(I,M)+LV(J,M))&=A(LD(LV(I,M))-LV(I,M)+LV(J,M))+AE(I,J)ELSEA(LD(LV(J,M))-LV(J,M)+LV(I,M))&=A(LD(LV(J,M))-LV(J,M)+LV(I,M))+AE(I,J)END IF100 CONTINUE200 CONTINUE300 CONTINUERETURNENDSUBROUTINE AS (NS,N,NA,IS,LD,A)! Introduce support conditions into global stiffness matrix A DIMENSION IS(NS),LD(N)DOUBLE PRECISION A(NA)DO 100 M=1,NSI=3*(IS(M)/10)-3+MOD(IS(M),10)A(LD(I))=1D22100 CONTINUERETURNENDSUBROUTINE LDLT (N,NA,LD,A,T)! Solve equations (1) - decomposition of matrix A DIMENSION LD(N)DOUBLE PRECISION A(NA),T(N),SUMDO 300 I=2,NLDI=LD(I)I1=I-LDI+LD(I-1)+1DO 200 J=I1,I-1LDJ=LD(J)J1=J-LDJ+LD(J-1)+1IF(I1.GT.J1) J1=I1SUM=0.0D0DO 100 K=J1,J-1SUM=SUM+T(K)*A(LDJ-J+K)100 CONTINUET(J)=A(LDI-I+J)-SUMA(LDI-I+J)=T(J)/A(LDJ)A(LDI)=A(LDI)-T(J)*A(LDI-I+J)200 CONTINUE300 CONTINUERETURNENDSUBROUTINE SLVEQ (N,NA,MAXBDW,LD,A,B)! Solve equations (2) - forward & back substitution DIMENSION LD(N)DOUBLE PRECISION A(NA),B(N)DO 200 I=2,NLDI=LD(I)I1=I-LDI+LD(I-1)+1DO 100 J=I1,I-1B(I)=B(I)-A(LDI-I+J)*B(J)100 CONTINUE200 CONTINUEDO 300 I=1,NB(I)=B(I)/A(LD(I))300 CONTINUEDO 500 I=N-1,1,-1IMIN=I+MAXBDWIF(IMIN.GT.N) IMIN=NDO 400 J=I+1,IMINLDJ=LD(J)J1=J-LDJ+LD(J-1)+1IF(I.GE.J1) B(I)=B(I)-A(LDJ-J+I)*B(J)400 CONTINUE500 CONTINUERETURNENDSUBROUTINE B0 (LC,N,NLJ,B)! Initialize joint load vector BDOUBLE PRECISION B(N)WRITE (*,1)LC,NLJ1 FORMAT(/15X,'Loading Case',I3&//10X,'The Loadings at Joints'&//17X,'NLJ=',I4)DO 100 I=1,NB(I)=0.0D0100 CONTINUERETURNENDSUBROUTINE IOLJB (N,NLJ,ILJ,PX,PY,PM,B)! Read data of loading at joint & form joint load vector B DIMENSION ILJ(NLJ),PX(NLJ),PY(NLJ),PM(NLJ)DOUBLE PRECISION B(N)READ (3,*)(ILJ(M),PX(M),PY(M),PM(M),M=1,NLJ)WRITE (*,1)1 FORMAT(/2X,'ILJ',11X,'PX',14X,'PY',15X,'PM')WRITE (*,2)(ILJ(M),PX(M),PY(M),PM(M),M=1,NLJ)2 FORMAT(1X,I4,2F16.4,F18.5)DO 100 M=1,NLJI=ILJ(M)*3B(I-2)=PX(M)B(I-1)=PY(M)B(I)=PM(M)100 CONTINUERETURNENDSUBROUTINE F0(NLM,NM,F)! Initialize terminal forces of membersDIMENSION F(6,NM)WRITE (*,1)NLM1 FORMAT(/10X,'The Loadings at Members'&//17X,'NLM=',I4)DO 200 J=1,NMDO 100 I=1,6F(I,J)=0.0100 CONTINUE200 CONTINUERETURNENDSUBROUTINE IOLMFB(NM,NJ,N,NLM,ILM,ITL,PV,DST,&IST,IEN,X,Y,LV,F,B)! Read data of loading at member & calculate fixed-end forces, ! add equivalent joint loads to vector BDIMENSION ILM(NLM),ITL(NLM),PV(NLM),DST(NLM),IST(NM),& IEN(NM),X(NJ),Y(NJ),LV(6,NM),F(6,NM)DOUBLE PRECISION B(N)READ(3,*)(ILM(M),ITL(M),PV(M),DST(M),M=1,NLM)WRITE (*,1)1 FORMAT(/2X,'ILM',2X,'ITL',11X,'PV',12X,'DST')WRITE(*,2)(ILM(M),ITL(M),PV(M),DST(M),M=1,NLM)2 FORMAT(1X,I4,I5,F16.4,F16.6)DO 100 M=1,NLML=ILM(M)CALL RLCS (L,NM,NJ,IST,IEN,X,Y,RL,C,S) D1=DST(M)D2=RL-D1IF (ITL(M).EQ.1.OR.ITL(M).EQ.3)THENP1=PV(M)*CP2=-PV(M)*SENDIFIF(ITL(M).EQ.2.OR.ITL(M).EQ.4)THENP1=PV(M)*SP2=PV(M)*CENDIFIF(ITL(M).EQ.1.OR.ITL(M).EQ.2)THENF1=-P1*D2/RLF4=-P1-F1F2=-P2*D2*D2*(RL+2.0*D1)/(RL*RL*RL)F5=-P2-F2F3=-P2*D1*D2*D2/(RL*RL)F6=P2*D1*D1*D2/(RL*RL)ENDIFIF(ITL(M).EQ.3.OR.ITL(M).EQ.4)THENG=P2*D1*D1/(12.0*RL*RL)F3=-G*((6.0*RL-8.0*D1)*RL+3.0*D1*D1)F6=G*D1*(4.0*RL-3.0*D1)F5=-6.0*G*D1*(2.0-D1/RL)F2=-P2*D1-F5F4=-P1*D1*D1/(2.0*RL)F1=-P1*D1-F4ENDIFF(1,L)=F(1,L)+F1F(2,L)=F(2,L)+F2F(3,L)=F(3,L)+F3F(4,L)=F(4,L)+F4F(5,L)=F(5,L)+F5F(6,L)=F(6,L)+F6B(LV(1,L))=B(LV(1,L))-F1*C+F2*SB(LV(2,L))=B(LV(2,L))-F1*S-F2*CB(LV(3,L))=B(LV(3,L))-F3B(LV(4,L))=B(LV(4,L))-F4*C+F5*SB(LV(5,L))=B(LV(5,L))-F4*S-F5*CB(LV(6,L))=B(LV(6,L))-F6100 CONTINUERETURNENDSUBROUTINE BS(NS,N,IS,VS,B)! Introduce support conditions into joint load vector B DIMENSION IS(NS),VS(NS)DOUBLE PRECISION B(N)DO 100 M=1,NSI=3*(IS(M)/10)-3+MOD(IS(M),10)B(I)=VS(M)*1D22100 CONTINUERETURNENDSUBROUTINE OJD(NJ,N,B)! Print joint displacementsDOUBLE PRECISION B(N)WRITE (*,1)1 FORMAT(/13X,'The Results of Calculation'&//10X,'The Joint Displacements'&//1X,'joint',8X,'u',15X,'v',14X,'phi') WRITE (*,2)(M,B(3*M-2),B(3*M-1),B(3*M),M=1,NJ)2 FORMAT(1X,I4,1P3E16.6)RETURNENDSUBROUTINE COTF(E,NM,NJ,N,IST,IEN,AR,RI,X,Y,LV,B,F)! Calculate & print terminal forces of members DIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ),&LV(6,NM),F(6,NM)DOUBLE PRECISION B(N),U1,U2,U3,U4,U5,U6WRITE (*,1)1 FORMAT(/10X,'The Terminal Forces'&//1X,'member',4X,'N(st)',6X,'Q(st)',7X,'M(st)',& 6X,'N(en)',6X,'Q(en)',7X,'M(en)')DO 100 M=1,NMCALL KEBAR(M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4) U1=B(LV(1,M))*C+B(LV(2,M))*SU2=-B(LV(1,M))*S+B(LV(2,M))*CU3=B(LV(3,M))U4=B(LV(4,M))*C+B(LV(5,M))*SU5=-B(LV(4,M))*S+B(LV(5,M))*CU6=B(LV(6,M))F(1,M)=F(1,M)+E1*(U1-U4)F(2,M)=F(2,M)+E2*(U2-U5)+E3*(U3+U6)F(3,M)=F(3,M)+E3*(U2-U5)+E4*(U3+0.5*U6)F(4,M)=F(4,M)+E1*(U4-U1)F(5,M)=F(5,M)+E2*(U5-U2)-E3*(U3+U6)F(6,M)=F(6,M)+E3*(U2-U5)+E4*(0.5*U3+U6)WRITE(*,2)M,F(1,M),F(2,M),F(3,M),F(4,M),F(5,M),F(6,M) 2 FORMAT(1X,I4,2(2F11.3,F12.3))100 CONTINUERETURNEND输入数据文件************************************************** ** 114811150上机试题1 ** **************************************************3.65E7 13 13 12 11 2 0.2 4.1667E-32 3 0.2 4.1667E-33 6 0.16 2.1333E-32 5 0.16 2.1333E-34 5 0.2 4.1667E-35 6 0.2 4.1667E-36 9 0.16 2.1333E-35 8 0.16 2.1333E-38 9 0.2 4.1667E-37 8 0.2 4.1667E-38 10 0.16 2.1333E-310 12 0.16 2.1333E-311 13 0.2 4.1667E-30 00 100 2010 010 1010 2020 020 1020 2025 1030 030 1030 1011 012 013 041 042 043 071 072 073 0111 0112 0113 032 150 0 03 0 0 -50 10 0 0 -25 63 2 -250 54 4 -10 107 4 -10 108 2 -200 59 3 10 1010 3 10 10输出文件************************************************* * * * 114811150上机试题1 * * * *************************************************The Input DataThe General InformationE NM NJ NS NLC3.650E+07 13 13 12 1The Information of Membersmember start end A I1 12 2.000000E-01 4.166700E-032 23 2.000000E-01 4.166700E-033 3 6 1.600000E-01 2.133300E-034 25 1.600000E-01 2.133300E-035 4 5 2.000000E-01 4.166700E-036 5 6 2.000000E-01 4.166700E-037 6 9 1.600000E-01 2.133300E-038 5 8 1.600000E-01 2.133300E-039 8 9 2.000000E-01 4.166700E-0310 7 8 2.000000E-01 4.166700E-0311 8 10 1.600000E-01 2.133300E-0312 10 12 1.600000E-01 2.133300E-0313 11 13 2.000000E-01 4.166700E-03The Joint Coordinatesjoint X Y1 .000000 .0000002 .000000 10.0000003 .000000 20.0000004 10.000000 .0000005 10.000000 10.0000006 10.000000 20.0000007 20.000000 .0000008 20.000000 10.0000009 20.000000 20.00000010 25.000000 10.00000011 30.000000 .00000012 30.000000 10.00000013 30.000000 10.000000The Information of SupportsIS VS11 .00000012 .00000013 .00000041 .00000042 .00000043 .00000071 .00000072 .00000073 .000000111 .000000112 .000000113 .000000( NA= 258 )( NW= 927 )Loading Case 1The Loadings at JointsNLJ= 3ILJ PX PY PM2 150.0000 .0000 .000003 .0000 .0000 -50.00000 10 .0000 .0000 -25.00000The Loadings at MembersNLM= 6ILM ITL PV DST3 2 -250.0000 5.0000004 4 -10.0000 10.0000007 4 -10.0000 10.0000008 2 -200.0000 5.0000009 3 10.0000 10.00000010 3 10.0000 10.000000The Results of CalculationThe Joint Displacementsjoint u v phi1 9.727862E-21 -7.663517E-21 -5.946892E-202 8.890687E-02 -1.049797E-04 -7.120787E-033 1.470285E-01 -2.323802E-04 -7.465712E-034 9.683457E-21 -3.536778E-20 -5.930753E-205 8.886287E-02 -4.844902E-04 -7.160652E-036 1.469822E-01 -7.581843E-04 5.116445E-047 1.058868E-20 -2.196870E-20 -6.233835E-208 8.890696E-02 -3.009411E-04 -6.177452E-039 1.470137E-01 -3.792984E-04 -1.977168E-0310 8.890696E-02 -3.520154E-02 -7.782787E-0311 0.000000E+00 0.000000E+00 0.000000E+0012 8.890696E-02 -7.411547E-02 -7.782787E-0313 0.000000E+00 0.000000E+00 0.000000E+00The Terminal Forcesmember N(st) Q(st) M(st) N(en) Q(en) M(en)1 76.635 97.279 594.689 -76.635 -97.279 378.0972 93.002 -27.030 -129.904 -93.002 27.030 -140.3963 27.030 93.002 90.396 -27.030 156.998 -410.3724 25.691 -16.367 -248.192 -25.691 116.367 -415.4805 353.678 96.835 593.075 -353.678 -96.835 375.2706 199.797 45.396 110.296 -199.797 -45.396 343.6647 -18.366 42.799 66.708 18.366 57.201 -138.7178 -25.747 37.514 -70.087 25.747 162.486 -554.7759 57.201 81.634 177.624 -57.201 18.366 138.71710 219.687 155.887 706.717 -219.687 -55.887 352.15111 .000 .000 25.000 .000 .000 -25.00012 .000 .000 .000 .000 .000 .00013 .000 .000 .000 .000 .000 .000( NA= 258 )( NW= 951 )Press any key to continue程序二：平面三节点有限元程序，如图1所示的悬臂梁，受均布荷载q=1N/mm2 作用。
平面三角形单元有限元程序设计P9 m 9 m一、题目如图1所示,一个厚度均匀的三角形薄板,在顶点作用沿板厚方向均匀分布的竖向载荷。
已知:P=150N/m,E=200GPa,=0、25,t=0、1m,忽略自重。
试计算薄板的位移及应力分布。
要求:1.编写有限元计算机程序,计算节点位移及单元应力。
(划分三角形单元,单元数不得少于30个);2.采用有限元软件分析该问题(有限元软件网格与程序设计网格必须一致),详细给出有限元软件每一步的操作过程,并将结果与程序计算结果进行对比(任选取三个点,对比位移值);3.提交程序编写过程的详细报告及计算机程序;4.所有同学参加答辩,并演示有限元计算程序。
有限元法中三节点三角形分析结构的步骤如下:1)整理原始数据,如材料性质、荷载条件、约束条件等,离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2)单元分析,建立单元刚度矩阵。
3)整体分析,建立总刚矩阵。
4)建立整体结构的等效节点荷载与总荷载矩阵5)边界条件处理。
6)解方程,求出节点位移。
7)求出各单元的单元应力。
8)计算结果整理。
一、程序设计网格划分如图,将薄板如图划分为6行,并建立坐标系,则刚度矩阵的集成建立与总刚度矩阵等维数的空矩阵,已变单元刚度矩阵的集成。
由单元分析已知节点、单元的排布规律,继而通过循环计算求得每个单元对应的节点序号。
通过循环逐个计算:(1)每个单元对应2种单元刚度矩阵中的哪一种; (2)该单元对应总刚度矩阵的那几行哪几列 (3)将该单元的单元刚度矩阵加入总刚度矩阵的对应行列循环又分为3层循环:(1)最外层:逐行计算(2)中间层:该行逐个计算(3)最里层:区分为第 奇/偶 数个计算XYPXYP节点编号单元编号单元刚度的集成:[][][][][][]''''''215656665656266256561661eZeeeZeZeeeekkkKkkkkkk+⋯++=⇓=⇒==⇒==⇒=⨯⨯⨯⨯⨯⨯M边界约束的处理:划0置1法适用:这种方法适用于边界节点位移分量为已知(含为0)的各种约束。
3结点三角形单元有限元程序（MATLAB语言）该程序包括以下6个部分：1．主程序tri_fem：用于数据的录入和其他程序的调用；2．总刚程序Kf：计算结构的总体刚度；3．各结点位移求解程序xf：求解各结点的位移；4．线性方程组求解程序Jordan：Gauss-Jordan法求解非约束结点的位移；5．应力应变程序ss：由各结点位移求解各单元内的三个结点的应力stress和应变strain；6．数据录入程序input：录入材料、几何尺寸、单元编号和结点编号、位移约束和已知载荷等。
以课本P25页例2.2为例，其input程序为function [E,v,t,EN,ecode,NN,node,RN,RC,PN,PC]=input()E=2.1e11; v=1/3; t=1; %杨氏模量Pa，泊松比，厚度EN=2; %单元数ecode=[3 1 2; %单元编号单元1 3-1-2;单元2 1-3-41 3 4];NN=4; %结点数node=[0 0; %各结点坐标2 0;2 1;0 1];RN=2; %被约束的位移数RC=[1 4]; %被约束的结点PN=2; %有载荷的结点数%PC（1）表示有载荷的结点，PC(2)表示各结点的力，PC（3）表示载荷方向，0为x方向，1为y方向PC=[2 3;-1/2 -1/2;1 1]; %结点2、3分别有y负方向上-1/2N的力作用在matlab环境下，输入则程序运行结果如下：该程序求解的结点位移结果和结点应力结果与课本给出的结果一致。
附录：%%-------------平面三角形单元有限元法---------------------%%function [x strain stress]=tri_fem()[E,v,t,EN,ecode,NN,node,RN,RC,PN,PC]=input; %调入已定材料、几何尺寸以及单元和结点编号及约束和载荷分布[n m]=size(ecode);if EN~=nerror('Wrong elementnumber EN or wrong elementcode ecode!');return;end[n m]=size(node);if NN~=nerror('Wrong nodenumber NN or wrong node-coordinate node!');return;ende=zeros(EN,6);A=zeros(EN,1); %面积for i=1:ENe(i,:)=[node(ecode(i,1),:),node(ecode(i,2),:),node(ecode(i,3),:)]; %各三角形单元的节点坐标D=[1,node(ecode(i,1),:)1,node(ecode(i,2),:)1,node(ecode(i,3),:)];A(i,1)=1/2*det(D);end%% 形成单元参数b=zeros(EN,3);c=zeros(EN,3); %各单元参数初始化for i=1:ENb(i,1)=e(i,4)-e(i,6); b(i,2)=e(i,6)-e(i,2); b(i,3)=e(i,2)-e(i,4);c(i,1)=-e(i,3)+e(i,5); c(i,2)=-e(i,5)+e(i,1); c(i,3)=-e(i,1)+e(i,3);end%% 求得总刚，并引入约束和载荷求得各结点位移K=Kf(E,v,t,EN,ecode,NN,A,b,c); %调用函数Kf，求得结构的总体刚度矩阵x=xf(NN,RN,RC,PN,PC,K); %调用函数xf,求得各结点位移%% 求解应力应变[strain stress]=ss(E,v,EN,ecode,A,b,c,x);%% 单元刚度矩阵与结构刚度矩阵function K=Kf(E,v,t,EN,ecode,NN,A,b,c)Ke=zeros(6,6); %单元的刚度矩阵，初始为6*6阶零矩阵K=zeros(NN*2,NN*2); %结构的总体刚度矩阵，初始为零矩阵for m=1:1:EN %m为单元号for i=1:1:3for j=1:1:3Ke(2*i-1,2*j-1)=b(m,i)*b(m,j)+(1-v)*c(m,i)*c(m,j)/2;Ke(2*i-1,2*j)=v*c(m,i)*b(m,j)+(1-v)*b(m,i)*c(m,j)/2;Ke(2*i,2*j-1)=v*b(m,i)*c(m,j)+(1-v)*c(m,i)*b(m,j)/2;Ke(2*i,2*j)=c(m,i)*c(m,j)+(1-v)*b(m,i)*b(m,j)/2;endendKe=E*t/(4*(1-v^2)*A(EN)).*Ke; %获得单元m的刚度矩阵Kb=mat2cell(Ke,ones(1,3)*2,ones(1,3)*2); %将单元矩阵Ke分为3*3块set1=ones(1,NN)*2;Ka=mat2cell(zeros(NN*2,NN*2),set1,set1); %将总刚K分为NN*NN块for i=1:1:3for j=1:1:3Ka(ecode(m,i),ecode(m,j))=Kb(i,j); %各单元刚度矩阵整体编号,并叠加endendK=K+cell2mat(Ka);end%分块矩阵K合成一个矩阵K%% 引入位移向量和右端项function x=xf(NN,RN,RC,PN,PC,K)x=ones(NN*2,1); %位移初始为0向量for i=1:RNx(RC(i)*2-1)=0;x(RC(i)*2)=0;end %被约束结点位移为0%%----------------引入已知结点载荷-------------%px=zeros(NN*2,1); %载荷初始为0向量for i=1:PNif PC(3,i)==1px(PC(1,i)*2)=PC(2,i);else if PC(3,i)==0px(PC(1,i)*2-1)=PC(2,i);endendend%%----------------引入已知结点载荷-------------%%%----------------求解非约束结点的位移X-------------%set1=ones(1,NN)*2;Ka=mat2cell(K,set1,set1);pxa=mat2cell(px,set1,1);AN=zeros(2*(NN-RN),2*(NN-RN));ANa=mat2cell(AN,ones(1,NN-RN)*2,ones(1,NN-RN)*2);bn=zeros(2*(NN-RN),1);bna=mat2cell(bn,ones(1,NN-RN)*2,1);BN=zeros(2*RN,2*(NN-RN));BNa=mat2cell(BN,ones(1,RN)*2,ones(1,NN-RN)*2);m=1;for i=1:1:NNif x(2*i)==1M(m)=i;m=m+1;endendfor i=1:1:NN-RNfor j=1:1:NN-RNANa(i,j)=Ka(M(i),M(j));bna(i,1)=pxa(M(i),1);endendfor i=1:RNfor j=1:NN-RNBNa(i,j)=Ka(RC(i),M(j));endendAN=cell2mat(ANa);bn=cell2mat(bna);BN=cell2mat(BNa);X=Jordan(AN,bn); %利用Gauss-Jordan法求解非约束结点的位移X %%----------------求解非约束结点的位移X-------------%%----------------由X可得所有结点位移x-------------%BN=BN*X;m=1; n=1;for i=1:1:NNif x(2*i)==1x(2*i-1)=X(m);x(2*i)=X(m+1);m=m+2;else if x(2*i)==0px(2*i-1)=BN(n);px(2*i)=BN(n+1);n=n+2;endendend%% 列主元Jordan消去法将系数矩阵化成对角矩阵求解方程组的数值解function x=Jordan(A,b)%开始计算，赋初值[n,m]=size(A);x=zeros(n,1);for k=1:n%选主元max1=0;for i=k:nif abs(A(i,k))>max1max1=abs(A(i,k));r=i;endend%交换两行if r>kfor j=k:nz=A(k,j); A(k,j)=A(r,j); A(r,j)=z;endz=b(k); b(k)=b(r); b(r)=z;end%消元计算b(k)=b(k)/A(k,k);for j=k+1:nA(k,j)=A(k,j)/A(k,k);endfor i=1:nif i~=kfor j=k+1:nA(i,j)=A(i,j)-A(i,k)*A(k,j);endb(i)=b(i)-A(i,k)*b(k);endendend%输出xfor i=1:nx(i)=b(i);end%% 求解应力应变function [strain stress]=ss(E,v,EN,ecode,A,b,c,x)strain=zeros(3,EN); %应变初始为0矩阵stress=zeros(3,EN); %应力初始为0矩阵D=E/(1-v^2)*[1 v 0;v 1 0;0 0 (1-v)/2];for m=1:1:ENB=[b(m,1) 0 b(m,2) 0 b(m,3) 0;0 c(m,1) 0 c(m,2) 0 c(m,3);c(m,1) b(m,1) c(m,2) b(m,2) c(m,3) b(m,3)];B=B/2/A(m,1); %应变矩阵S=D*B; %应力矩阵X=[x(2*ecode(m,1)-1),x(2*ecode(m,1)),x(2*ecode(m,2)-1),x(2*ecode(m,2)),x(2*ecode(m,3)-1),x (2*ecode(m,3))]';strain(:,m)=B*X;stress(:,m)=S*X;end。
平面三角形３结点有限元程序1、程序名：，2、程序功能本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种荷载的作用，在计算自重时y轴取垂直向上为正；能处理非零已知位移，如支座沉降的作用。
主要输出的内容包括：结点位移、单元应力、主应力、第一主应力与x轴的夹角以及约束结点的支座反力。
程序采用Visual Fortran 编制而成，输入数据全部采用自由格式。
3、程序流程及框图图１-1 程序流程图图１-2 程序框图其中，各子程序的功能如下：INPUT——输入结点坐标、单元信息和材料参数；MR——形成结点自由度序号矩阵；FORMMA——形成指标矩阵MA（N）并调用其他功能子程序，相当于主控程序；DIV——取出单元的3个结点号码和该单元的材料号并计算单元的b i，c i等；MGK——形成整体劲度矩阵并按一维存放在SK（NH）中；LOAD——形成整体结点荷载列阵F；OUTPUT——输出结点位移或结点荷载；TREAT——由于有非零已知位移，对K和F进行处理；DECOMP——整体劲度矩阵的分解运算；FOBA——前代、回代求出未知结点位移δ；ERFAC——计算约束结点的支座反力；KRS——计算单元劲度矩阵中的子块K rs。
4、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。
在运行程序之前，必须根据程序中输入要求建立一个存放原始数据的文件，这个文件的名字由少于8个字符或数字组成。
数据文件包括如下内容：⑴总控信息，共一条，9个数据NP，NE，NM，NR，NI，NL，NG，ND，NCNP——结点总数；NE——单元总数；NM——材料类型总数；NR——约束结点总数；NI ——问题类型标识，0为平面应力问题，1为平面应变问题；NL ——受荷载作用的结点的数目；NG ——考虑自重作用为1，不计自重为0；ND ——非零已知位移结点的数目；NC ——要计算支座约束反力的结点数目。
⑵ 材料信息，共NM 条，每条依次输入EO ，VO ，W ，tEO ——弹性模量（kN/m 2）；VO ——泊松比；W ——材料容重（kN/m 3）；t ——单元厚度（m ）。
