平面三角形单元有限元程序的设计说明
有限元实验2-三角形单元的形函数性质
实验三：三角形单元的形函数性质验证
一、 实验目的
1、加深对平面三角形单元有限元分析过程的理解；
2、掌握平面三角形单元形函数矩阵的求解过程和性质。
二、 实验要求
1、明确形函数矩阵的含义和求法；
2、根据题目要求，给出具体的计算过程；
3、编制相应的matlab 计算程序并调试运行。
三、 实验内容
用有限元法求图示平面三角形单元的形函数。
已知节点i,j,m 在xoy 平面中。
用所编程序验证形函数特性：
1、在单元任一点上，三个形函数之和等于1；
2、形函数N i 在i 点的函数值为1，在j 点及m 点的函数值为零；
3、三边上任一点的形函数与第三个顶点的坐标无关。
四、 实验提示
1、()() 21,y c x b a y x N i i i i ++∆=
，其中下标i , j , m 轮换; 2、()()m i j m i j i m m j j i y x y x y x y x y x y x ++++=∆2
1 -21; 3、j m i m m i j m m j i x x c y y b y x y x a -=-=-=;;，其中下标i , j , m 轮换。
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。
平面三角形单元有限元程序设计
.. 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行，并建立坐标系，那么X Y P X YP刚度矩阵的集成建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。
由单元分析节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。
通过循环逐个计算：〔1〕每个单元对应2种单元刚度矩阵中的哪一种； 〔2〕该单元对应总刚度矩阵的那几行哪几列〔3〕将该单元的单元刚度矩阵参加总刚度矩阵的对应行列循环又分为3层循环：〔1〕最外层：逐行计算〔2〕中间层：该行逐个计算〔3〕最里层：区分为第奇/偶数个计算单元刚度的集成：[][][][][][]''''''215656665656266256561661eZeeeZeZeeeekkkKkkkkkk+⋯++=⇓=⇒==⇒==⇒=⨯⨯⨯⨯⨯⨯边界约束的处理：划0置1法适用：这种方法适用于边界节点位移分量为(含为0)的各种约束。
有限元程序设计--第五章 线性三角形单元
16
单元应变和应力矩阵
由于同一单元中的D、B矩阵都是常数矩阵，所以S矩阵也是常 数矩阵。也就是说，三角形三节点单元内的应力分量也是常 量。 当然，相邻单元的E， ， A和bi、ci(i，j， m)一般不完全相同， 因而具有不同的应力，这就造成在相邻单元的公共边上存在 着应力突变现象。但是随着网格的细分，这种突变将会迅速 减小。
b2 B2 0 c2
b1 , c1 , b2 , c2 , b3 , c3 与x、y无关，都是常量，因此B矩
阵也是常量。单元中任一点的应变分量是 B矩阵与单元节点位
移的乘积，因而也都是常量。因此，这种单元被称为常应变单
元。
03:25
15
单元应变和应力矩阵
平面应力： x σ ( x, y ) y
03:25
12
形函数的性质
当单元的位移函数满足完备性要求时，称单元是完备的（通
常较容易满足）。当单元的位移函数满足协调性要求时，称 单元是协调的。
当势能泛函中位移函数的导数是2阶时，要求位移函数在单元
的交界面上具有C1或更高的连续性，这时构造单元的插值函 数往往比较困难。在某些情况下，可以放松对协调性的要求 ，只要单元能够通过分片试验 (Patch test)，有限元分析的解 答仍然可以收敛于正确的解。这种单元称为非协调单元。
3
平面三角形单元
显然，三角形三个节点的的位移可由下列方程给出，在各节点上 的水平位移方程为： u1=1+2 x1 +3 y1 u2=1+2 x2 +3 y2 u3=1+2 x3 +3 y3
1 1 x1 解出 2 1 x2 1 x 3 3 y1 y2 y3
基于Matlab语言的按平面三角形单元划分的结构有限元程序设计
基于Matlab说话的按平面三角形单元划分的构造有限元程序设计专业：建筑与土木匠程班级：建工研12-2姓名：韩志强学号：471220580基于Matlab说话的按平面三角形单元划分构造有限元程序设计一、有限单元发及Matlab说话概述1. 有限单元法跟着现代工业.临盆技巧的成长,不竭请求设计高质量.高程度的大型.庞杂和周详的机械及工程构造.为此目标,人们必须预先经由过程有用的盘算手腕,确实的猜测即将诞生的机械和工程构造,在将来工作时所产生的应力.应变和位移是以,须要追求一种简略而又准确的数值剖析办法.有限单元法恰是顺应这种请求而产生和成长起来的一种十分有用的数值盘算办法.有限元法把一个庞杂的构造分化成相对简略的“单元”,各单元之间经由过程结点互相衔接.单元内的物理量由单元结点上的物理量按必定的假设内插得到,如许就把一个庞杂构造从无穷多个自由度简化为有限个单元构成的构造.我们只要剖析每个单元的力学特征,然后按照有限元法的规矩把这些单元“拼装”成整体,就可以或许得到整体构造的力学特征.有限单元法根本步调如下：(1)构造离散：构造离散就是树立构造的有限元模子,又称为网格划分或单元划分,即将构造离散为由有限个单元构成的有限元模子.在该步调中,须要依据构造的几何特征.载荷情形等肯定单元体内随意率性一点的位移插值函数.(2)单元剖析：依据弹性力学的几何方程以及物理方程肯定单元的刚度矩阵.(3)整体剖析:把各个单元按本来的构造从新衔接起来,并在单元刚度矩阵的基本上肯定构造的总刚度矩阵,形成如下式所示的整体有限元线性方程：式中阵.(4)载荷移置:依据静力等效道理,将载荷移置到响应的节点上,形成节点载荷矩阵.(5)鸿沟前提处理:对式①所示的有限元线性方程进行鸿沟前提处理.(6)求解线性方程:求解式①所示的有限元线性方程,得到节点的位移.在该步调中,如有限元模子的节点越多,则线性方程的数目就越多,随之有限元剖析的盘算量也将越大.(7)求解单元应力及应变依据求出的节点位移求解单元的应力和应变.(8)成果处理与显示进入有限元剖析的后处理部分,对盘算出来的成果进行加工处理,并以各类情势将盘算成果显示出.在用有限元法进行构造剖析时,将会碰到大量的数值盘算,因而在实用上是必定要借助于盘算机和有限元程序,才干完成这些庞杂而沉重的数值盘算工作.而Matlab是当今国际科学界最具影响力和活气的软件.它来源于矩阵运算,并已经成长成一种高度集成的盘算机说话.它供给了壮大的科学盘算,灵巧的程序设计流程,高质量的图形可视化与界面设计,便捷的与其他程序和说话接口的功效.Matlab在列国高校与研讨单位起侧重大的感化.“工欲善其事,必先利其器”.假如有一种十分有用的对象能解决在教授教养与研讨中碰到的问题,那么Matlab说话恰是如许的一种对象.它可以将运用者从繁琐.无谓的底层编程中解放出来,把有限的珍贵时光更多地花在解决问题中,如许无疑会进步工作效力. 今朝,Matlab已经成为国际上最风行的科学与工程盘算的软件对象,如今的Matlab已经不但仅是一个“矩阵试验室”了,它已经成为了一种具有普遍运用远景的全新的盘算机高等编程说话了,有人称它为“第四代”盘算机说话,它在国表里高校和研讨部分正扮演侧重要的脚色.Matlab说话的功效也越来越壮大,不竭顺应新的请求提出新的解决办法.可以预感,在科学运算.主动掌握与科学画图范畴Matlab说话将长期保持其举世无双的地位.为此,本例采取Matlab说话编程,以运用其快捷壮大的矩阵数值盘算功效.二、问题描写一矩形薄板,一边固定,推却150kN分散荷载,构造简图及按平面三角形单元划分的有限元模子图如下所示.薄板厚在本例中,所取构造模子及数据重要用于程序设计理论剖析,与工程现实无关.三、参数输入：单元个数NELEM=4节点个数NNODE=6受束缚鸿沟点数NVFIX=2节点荷载个数NFORCE=1弹性模量YOUNG=2e8泊松比POISS=0.2厚度THICK=0.002单元节点编码数组节点坐标数组节点力数组FORCE =[4 0 -150]束缚信息数组以上数值数据为程序运行前输入的初始数据,存为“”文本格局,同时必须放在Matlab工作目次下,路径不合错误程序不克不及主动读取指定初始文件,运行出错.初始数据文本文件输入格局如下图：四、Matlab说话程序源代码：1.程序中变量解释NNODE 单元节点数NPION 总结点数NELEM 单元数NVFIX 受束缚鸿沟点数FIXED 束缚信息数组NFORCE 节点力数FORCE 节点力数组COORD 构造节点坐标数组LNODS 单元界说数组YOUNG 弹性模量POISS 泊松比THICK 厚度B 单元应变矩阵(3*6) D 单元弹性矩阵(3*3) S 单元应力矩阵(3*6) A单元面积ESTIF 单元刚度矩阵ASTIF 总体刚度矩阵ASLOD 总体荷载向量ASDISP 节点位移向量ELEDISP 单元节点位移向量STRESS 单元应力FP1 数据文件指针2.Matlab说话程序代码%****************************************************************************** %初始化及数据挪用format short e %设定输出类型clear %消除内存变量FP1=fopen('','rt'); %打开输入数据文件,读入掌握数据NELEM=fscanf(FP1,'%d',1), %单元个数（单元编码总数）NPION=fscanf(FP1,'%d',1), %结点个数（结点编码总数）NVFIX=fscanf(FP1,'%d',1)%受束缚鸿沟点数NFORCE=fscanf(FP1,'%d',1), %结点荷载个数YOUNG=fscanf(FP1,'%e',1), %弹性模量POISS=fscanf(FP1,'%f',1), %泊松比THICK=fscanf(FP1,'%f',1) %厚度LNODS=fscanf(FP1,'%d',[3,NELEM])' %单元界说数组（单元结点号）%响应为单元结点号（编码）.按逆时针次序输入COORD=fscanf(FP1,'%f',[2,NPION])' %结点坐标数组%坐标：x,y坐标（共NPOIN组）FORCE=fscanf(FP1,'%f',[3,NFORCE])' %结点力数组（受力结点编号,x偏向,y偏向）FIXED=fscanf(FP1,'%d',[3,NVFIX])' %束缚信息（束缚点,x束缚,y束缚）%有束缚为1,无束缚为0%***************************************************************************** %生成单元刚度矩阵并构成总体刚度矩阵ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0%***************************************************************************** for i=1:NELEM%生成弹性矩阵DD= [1 POISS 0;POISS 1 0;0 0 (1-POISS)/2]*YOUNG/(1-POISS^2)%***************************************************************************** %盘算当前单元的面积AA=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)])/2%***************************************************************************** %生成应变矩阵Bfor j=0:2b(j+1)= COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1);endB=[b(1) 0 b(2) 0 b(3) 0;...0 c(1) 0 c(2) 0 c(3);...c(1) b(1) c(2) b(2) c(3) b(3)]/(2*A);B1( :,:,i)=B;%***************************************************************************** %求应力矩阵S=D*BS=D*B;ESTIF=B'*S*THICK*A; %求解单元刚度矩阵a=LNODS(i,:); %暂时向量,用来记载当前单元的节点编号for j=1:3for k=1:3ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)…=ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2); %依据节点编号对应关系将单元刚度分块叠加到总刚%度矩阵中endendend%***************************************************************************** %将束缚信息参加总体刚度矩阵（对角元素改一法）for i=1:NVFIXif FIXED(i,2)==1ASTIF(:,(FIXED(i,1)*2-1))=0; %一列为零ASTIF((FIXED(i,1)*2-1),:)=0; %一行动零ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1;%对角元素为1endif FIXED(i,3)==1ASTIF( :,FIXED(i,1)*2)=0; %一列为零ASTIF(FIXED(i,1)*2,:)=0; %一行动零ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; %对角元素为1endend%***************************************************************************** %生成荷载向量ASLOD(1:2*NPION)=0; %总体荷载向量置零for i=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end%***************************************************************************** %求解内力ASDISP=ASTIF\ASLOD' %盘算节点位移向量ELEDISP(1:6)=0;%当前单元节点位移向量for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);%掏出当前单元的节点位移向量endiSTRESS=D*B1(:, :, i)*ELEDISP' %求内力endfclose(FP1); %封闭数据文件五、程序运行成果NELEM =4NPION =6NVFIX =2NFORCE =1YOUNG =200000000POISS =THICK =LNODS =1 2 62 3 42 4 52 5 6COORD =0 01 02 02 11 10 1FORCE =4 0 -150FIXED =1 1 16 1 1D =2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 00 0 8.3333e+007A =D =2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 00 0 8.3333e+007A =D =2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 00 0 8.3333e+007A =D =2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007A =ASDISP =（解释：以上12个ASDISP输出成果数据暗示节点位移向量,即各节点分离在X,Y偏向上产生的位移.）i =1STRESS =-1.6793e+005-3.3586e+004-1.3207e+005i =2STRESS =-5.5503e+004-5.5503e+004-5.5503e+004i =3STRESS =5.5503e+004-4.9526e+004-9.4497e+004i =4STRESS =1.6793e+005-2.7040e+004-1.7932e+004（解释：以上四组STRESS输出成果数据暗示单元应力SX,SY,SXY,i为单元号.）六、ANSYS建模比较剖析运用ANSYS有限元剖析软件,完整按照前面Matlab程序设计的各项前提参数以及单元划分方法树立ANSYS模子,其求解成果与以上程序盘算成果比较,验证程序是否准确.ANSYS模子节点单元如下图所示：ANSYS模子求解变形图如下所示：ANSYS求解节点位移成果列表显示如下：（单位：m）PRINT U NODAL SOLUTION PER NODE***** POST1 NODAL DEGREE OF FREEDOM LISTING *****THE FOLLOWING DEGREE OF FREEDOM RESULTS ARE IN THE GLOBAL COORDINATESYSTEMNODE UX UY UZ USUM1 0.0000 0.0000 0.0000 0.00006 0.0000 0.0000 0.0000 0.0000MAXIMUM ABSOLUTE VALUESNODE 4 4 0 4ANSYS求解单元应力成果列表显示如下：（单位：KN/m2）PRINT S ELEMENT SOLUTION PER ELEMENT***** POST1 ELEMENT NODAL STRESS LISTING *****THE FOLLOWING X,Y,Z VALUES ARE IN GLOBAL COORDINATESELEMENT= 1 PLANE182NODE SX SY SZ SXY SYZ1 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.00002 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 6 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 ELEMENT= 2 PLANE182NODE SX SY SZ SXY SYZ2 -55503. -55503. 0.0000 -55503. 0.00003 -55503. -55503. 0.0000 -55503. 0.00004 -55503. -55503. 0.0000 -55503. 0.0000 ELEMENT= 3 PLANE182NODE SX SY SZ SXY SYZ2 55503. -49526. 0.0000 -94497. 0.00004 55503. -49526. 0.0000 -94497. 0.00005 55503. -49526. 0.0000 -94497. 0.0000 ELEMENT= 4 PLANE182NODE SX SY SZ SXY SYZ2 0.16793E+06 -27040. 0.0000 -17932. 0.00005 0.16793E+06 -27040. 0.0000 -17932. 0.00006 0.16793E+06 -27040. 0.0000 -17932. 0.0000结论经由过程比较Matlab说话设计程序运行成果和ANSYS建模剖析成果可知,两种方法算出的成果完整一致,说程序设计准确.所以本程序实用于按三角形单元划分的平面构造有限元剖析.format short eclearFP1=fopen('LinearTriangleElement of Thin plate.txt','rt');NELEM=fscanf(FP1,'%d',1)NPION=fscanf(FP1,'%d',1)NVFIX=fscanf(FP1,'%d',1)NFORCE=fscanf(FP1,'%d',1)YOUNG=fscanf(FP1,'%e',1)POISS=fscanf(FP1,'%f',1)THICK=fscanf(FP1,'%f',1)LNODS=fscanf(FP1,'%d',[3,NELEM])COORD=fscanf(FP1,'%f',[2,NPION])FORCE=fscanf(FP1,'%f',[3,NFORCE])FIXED=fscanf(FP1,'%d',[3,NVFIX])ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0for i=1:NELEM%生成弹性矩阵DD= YOUNG/(1-POISS^2)*[1 POISS 0;POISS 1 0;0 0 (1-POISS)/2]%盘算当前单元的面积AA=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)])/2%生成应变矩阵Bfor j=0:2b(j+1)= COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1); endB=[b(1) 0 b(2) 0 b(3) 0;...0 c(1) 0 c(2) 0 c(3);...c(1) b(1) c(2) b(2) c(3) b(3)]/(2*A);B1( :,:,i)=B;%求应力矩阵S=D*BS=D*B;ESTIF=B'*S*THICK*A;a=LNODS(i,:);for j=1:3for k=1:3ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)…=ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2);endendend%将束缚信息参加总体刚度矩阵for i=1:NVFIXif FIXED(i,2)==1ASTIF(:,(FIXED(i,1)*2-1))=0;ASTIF((FIXED(i,1)*2-1),:)=0;ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1;endif FIXED(i,3)==1ASTIF( :,FIXED(i,1)*2)=0;ASTIF(FIXED(i,1)*2,:)=0;ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1;endend%生成荷载向量ASLOD(1:2*NPION)=0;for i=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3)end%求解内力ASDISP=ASTIF\ASLOD'ELEDISP(1:6)=0;for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);endiSTRESS=D*B1(:, :, i)*ELEDISP'endfclose(FP1);format short eclearFP1=fopen('LinearTriangleElement of Thin plate.txt','rt'); NELEM=fscanf(FP1,'%d',1)NPION=fscanf(FP1,'%d',1)NVFIX=fscanf(FP1,'%d',1)NFORCE=fscanf(FP1,'%d',1)YOUNG=fscanf(FP1,'%e',1)POISS=fscanf(FP1,'%f',1)THICK=fscanf(FP1,'%f',1)LNODS=fscanf(FP1,'%d',[NELEM,3])COORD=fscanf(FP1,'%f',[NPION,2])FORCE=fscanf(FP1,'%f',[NFORCE,3])FIXED=fscanf(FP1,'%d',[NVFIX,3])ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0 for i=1:NELEM%生成弹性矩阵DD= YOUNG/(1-POISS^2)*[1 POISS 0;POISS 1 0;0 0 (1-POISS)/2]%盘算当前单元的面积AA=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)])/2%生成应变矩阵Bfor j=0:2b(j+1)= COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1); endB=[b(1) 0 b(2) 0 b(3) 0;...0 c(1) 0 c(2) 0 c(3);...c(1) b(1) c(2) b(2) c(3) b(3)]/(2*A);B1( :,:,i)=B;%求应力矩阵S=D*BS=D*B;ESTIF=B'*S*THICK*A;a=LNODS(i,:);for j=1:3for k=1:3ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)…=ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2);endendend%将束缚信息参加总体刚度矩阵for i=1:NVFIXif FIXED(i,2)==1ASTIF(:,(FIXED(i,1)*2-1))=0;ASTIF((FIXED(i,1)*2-1),:)=0;ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1;endif FIXED(i,3)==1ASTIF( :,FIXED(i,1)*2)=0;ASTIF(FIXED(i,1)*2,:)=0;ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1;endend%生成荷载向量ASLOD(1:2*NPION)=0;for i=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3)end%求解内力ASDISP=ASTIF\ASLOD'ELEDISP(1:6)=0;for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2); endiSTRESS=D*B1(:, :, i)*ELEDISP'endfclose(FP1);。
平面三角形单元有限元程序设计
平面三角形单元有限元程序设计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)的各种约束。
三角形单元有限元程序设计
三角形单元有限元程序设计一、引言
⑴文档背景
⑵文档目的
⑶读者对象
⑷名词解释
二、程序结构设计
⑴程序流程图
⑵程序组成模块描述
⑶主要数据结构设计
⑷代码逻辑设计
三、数据预处理
⑴数据输入格式与读取
⑵数据清洗与去噪
⑶数据预处理方法
⑷数据分割与划分
四、网格
⑴网格算法介绍
⑵网格质量评估与改善
⑶网格的实现方法
五、有限元方法
⑴有限元离散的原理
⑵有限元单元的选择
⑶有限元离散的网格节点选取
⑷有限元插值函数的计算六、模型求解
⑴线性方程组的求解方法
⑵模型参数的设置与调整
⑶迭代算法的选择与实现七、模型评估与验证
⑴结果评估指标的选择
⑵模型验证方法
⑶结果可视化与分析
八、性能优化
⑴程序性能分析与评估
⑵程序高效化的方法与技巧
⑶并行计算与优化
九、实例与案例
⑴实例介绍与问题描述
⑵实例数据处理过程
⑶实例模型求解与结果分析十、总结与展望
⑴工作总结
⑵程序改进与优化展望
⑶研究方向与发展趋势
十一、附件
●附件1、程序源代码
●附件2、输入数据样例文件
●附件3、网格结果数据文件附录：法律名词及注释
●法律名词1、注释1
●法律名词2、注释2
●法律名词3、注释3
请注意，附件的具体文件名和数据内容将根据你的实际情况进行调整。
法律名词及注释部分需要根据实际需要进行扩充和修改。
请合理对文档结构和章节内容进行调整，以符合你的具体要求。
2012硕研(平面问题三角形单元有限元程序设计)
0、
主程序段—1、2
REAL KS(200,100) （定义数组总刚度矩阵）
OPEN(5,FILE=‘FEMT3.IN’)(打开原始数据文件) DIMENSION LND(50,3),X(100),Y(100),JR(20,3),PJ(20,3),P(200)
OPEN(6,FILE=‘FEMT3.OUT’,STATUS=‘NEW’)（输出文件）
READ(5,*)((PJ(I,J),J=1,3),I=1,NPJ)（荷载信息。依 次为：荷载所在结点号，x、y方向荷载大小）
100
NW=0（半带宽） DO 110 IE=1,NE DO 110 I=1,3 DO 110 J=1,3 IW=IABS(LND(IE,I)-LND(IE,J)) IF(NW.LT.IW) THEN NW=IW ENDIF 110 CONTINUE NW=(NW+1)*2 IF(IPS.NE.0) THEN E=E/(1.0-PR*PR) PR=PR/(1.0-PR) ENDIF END
• 有关带宽，储存方式请自修《有限元法概论》 中§5-12；《结构矩阵分析原理》第五章。
6．子程序6（形成总刚度矩阵[KS]）
SUBROUTINE STIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS) DIMENSION LND(NE,3),X(NJ),Y(NJ) REAL KS(NJ2,NW),KE(6,6) DO 5 I=1,NJ2 DO 5 J=1,NW 5 KS(I,J)=0. DO 10 IE=1,NE CALL STE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)
半带宽计算示例：
半带宽d的一般公式 半带宽d=(相邻结点码的最大差值+1)×结点位移未知量数 三角形单元：半带宽d=(相邻结点码的最大差值+1)×2 图中相邻结点码的最大差值是4,故d=(4+1)×2=10
第六章三角形单元的有限元法程序设计
从第j列的主对角线元素起到该列上方第一个非零 元素为止，所含元素的个数称为第j列的列高，记为hj ； 如果把第j列上方第1个非零元素的行号记为mj，则第j 列的列高为 hj = j - mj + 1 其实，hj就是第j行的左带宽，因而必有 UBW= max(hj)
j=1,2, …,N
利用节点位移信息数组 ID （去约束后节点位移自 由度编码），可容易地确定刚度矩阵 [K] 任何一列的 列高。
{ } [S ]{ }
例1：对角受压的正方形薄板，载 荷沿厚度均匀分布，为2N/m。由 于对称性，取1/4部分作为计算对
2N/m y x
象，试用有限元程序进行计算。
2N/m 2m 2m
h 1.0m, E 1.0, 0.0
例2：简支梁，梁高3m,跨度18m，厚度1m，承 受均布荷载10N/m2。已知 E 2 1010 N / m2 , 0.167
y为挤压应力，属次要应 力。材料力学不考虑 , FEM与弹性力学误差较大。
6.2 提高计算精度的方法
（1） 计算结果的整理 计算结果包括位移和应力两个方面。在位移方 面，一般无须进行整理工作。应力结果则需要 整理。通常认为计算出的应力是三角形单元形 心处的应力。而相邻单元之间的应力存在突变， 甚至正、负符号都不相同。为了由计算结果推 算出结构内某一点的接接实际的应力，必须通 过某种平均计算。通常可采用两单元平均法或 绕结点平均法。
平均法整理单元应力
两单元平均法：把两个相邻单元中的常应力加 以平均，用来表示公共边界中点处的应力。 绕结点平均法：把环绕某一结点的各单元常应 力加以平均，用以表示该结点的应力。在内结 点效果较好，而在边界结点可能很差，一般改 为应由内结点的应力外推计算出来。 （2）网格的细分 通过网格的细分，使每个单元的面积缩小，那 么尽管每个单元是应变、常应力单元，仍可较 好地反映结构中的应力变化，使得到的解答收 敛于问题的精确解。
有限元方法第五章平面三角形单元
这样，位移模式 (e) 和 (f) 就可以写为
u Ni ui N j u j N mum v Nivi N jv j Nmvm
(5-11)
也可写成矩阵形式
f
u v
Ni I
NjI
NmI e N e
(5-12)
式中 I是二阶单位矩阵；Ni 、Nj 、Nm 是坐标的函数， 它们反映了单元的位移状态，所以一般称之为形状函数，简 称形函数。矩阵 [N] 叫做形函数矩阵。三节点三角形单元的 形函数是坐标的线性函数。单元中任一条直线发生位移后仍 为一条直线，即只要两单元在公共节点处保持位移相等。则 公共边线变形后仍为密合。
式中 1、2、…6是待定常数。因三角形单元共有六个
自由度，且位移函数u、v在三个节点处的数值应该等于 这些点处的位移分量的数值。假设节点i、j、m的坐标分 别为（xi , yi ）、（xj , yj ）、（xm , ym ），代入 (b) 式， 得：
ui 1 2 xi 3 yi
, v j 4 5xi 6 yi
二、经典解与有限元解的区别：
微分 经 典 解 法 —— （解析法）
数目增到∞ 大小趋于 0
建立一个描述连续体 性质的偏微分方程
有限单元 离散化 集合
总体分析解
有限元法——连续体——单元——代替原连续体
（近似法）
（单元分析）
线性方程组
三、有限元法算题的基本步骤
1. 力学模型的选取 （平面问题，平面应变问题，平面应力问题，轴对称问题， 空间问题，板，梁，杆或组合体等，对称或反对称等）
0
(b)
Ni xm
,
ym
1 2
ai
bi xm
ci ym
0
平面三角形3节点有限元程序
平面三角形３结点有限元程序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 ）。
有限元方法-第五章--平面三角形单元
D
E
1 2
1
0
对 1 0
称
1
(i)
2
所以，[S]的子矩阵可记为
Si DBi
E
2 1 2
bi
1
bi
2
ci
ci
1
ci
2
bi
（ i
,
j
,
m轮换） (5-19)
对于平面应变问题，只要将 (i) 式中的E换成E/1-2 ， 换成 /1-，即得到其弹性矩阵
D
1
E1 1 2
1
1
起来，便可近似地表示整个区域的真实位移函数。这种 化繁为简、联合局部逼近整体的思想，正是有限单元法 的绝妙之处。
基于上述思想，我们可以选择一个单元位移模式，
单元内各点的位移可按此位移模式由单元节点位移通过
插值而获得。线性函数是一种最简单的单元位移模式，
故设
u 1 2x 3y
v 4 5x 6y
(b)
0
(b)
Ni xm
,
ym
1 2
ai
bi xm
ci ym
0
(c)
类似地有
N j xi , yi 0 , N j x j , y j 1 , N j xm , ym 0 (d) Nm xi , yi 0 , Nm x j , y j 0 , Nm xm , ym 1
式中 1、2、…6是待定常数。因三角形单元共有六个
自由度，且位移函数u、v在三个节点处的数值应该等于 这些点处的位移分量的数值。假设节点i、j、m的坐标分 别为（xi , yi ）、（xj , yj ）、（xm , ym ），代入 (b) 式， 得：
ui 1 2 xi 3 yi
有限元 2-弹性力学平面问题有限单元法(2.1三角形单元,2.2几个问题的讨论)
第2章 弹性力学平面问题有限单元法2.1 三角形单元(triangular Element)三角形单元是有限元分析中的常见单元形式之一，它的优点是:①对边界形状的适应性较好，②单刚形式及其推导比较简单,故首先介绍之。
一、结点位移和结点力列阵设图为从某一结构中取出的一典型三角形单元。
在平面应力问题中，单元的每个结点上有沿x、y两个方向的力和位移，单元的结点位移列阵规定为： 相应结点力列阵为: (式2-1-1){}⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧=⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=m j i m ed d d d m j j i v u v u v u i {}ii j j m X Y X (2-1-1)Y X Y iej m m F F F F ⎧⎫⎪⎪⎪⎪⎧⎫⎪⎪⎪⎪⎪⎪==⎨⎬⎨⎬⎪⎪⎪⎪⎩⎭⎪⎪⎪⎪⎪⎪⎩⎭二、单元位移函数和形状函数前已述及，有限单元法是一种近似方法，在单元分析中，首先要求假定（构造）一组在单元内有定义的位移函数作为近似计算的基础。
即以结点位移为已知量，假定一个能表示单元内部（包括边界）任意点位移变化规律的函数。
构造位移函数的方法是：以结点(i,j,m)为定点。
以位移(u i ,v i ,…u m v m 3)为定点上的函数值，利用普通的函数插值法构造出一个单元位移函数。
在平面应力问题中，有u,v 两个方向的位移,若假定单元位移函数是线性的,则可表示成:(,)12u u x y x yααα+46y ==+ 5(,)v v x y x ααα+==+ (2-1-2)a式中的6个待定常数α1 ,…, α6 可由已知的6个结点位移分量(3个结点的坐标)确定。
将13个结点坐标(x i,3iy y i ),(x j,y j ),(x m,y m )代入上式得如下两组线性方程: 12i i u x ααα+3=+12j j j x y u αα=+α+3m y (a)12m m u x ααα=++46i y和5i i v x αα=+α+465j j j x y v αα=+α+46m y (b)5m m v x ααα=++利用线性代数中解方程组的克来姆法则,由(a)可解出待定常数1α 、2α 、3α :211A Aα=22A 3A Aα=3Aα=式中行列式:2111i i 1i i i j m j j m m u x y A u x y u x y =j jm mu y A u y u y =3111i i j jm mx u A 2111i i j j m mAx y A x y x y x u x u ===A为△ijm 的面积,只要A不为0,则可由上式解出:112i i j j a u a u ()m m a u A α=++21(2i ij j bu b u )m m b u A α=++ （C）312i i j j c u c u ()m mc u A α=++i j a x y =−j i y x y =−m i j j i y x y 式中:m m j x y a x a x m m i =−y m y y =−m i j y ym i j b y =− b b j i =− （d）3c m i j x x =− j i c m x x =−m j i c x x =−m iy x y =−m为了书写方便，可将上式记为: a xm i j b i jy y =−(,,) i u j m uu u ruuu u r i jc m x x =−(,,)i j m uuu u r uuu u r)m m N x y u N x y u N x y u =++)m x y v 表示按顺序调换下标,即代表采用i,j,m 作轮换的方式便可得到(d)式。
第7章——有限元三角形单元程序设计
C ( 4, 1)
B(0, 2)
12:45
有限单元法
崔向阳
三角形单元程序
clear all first_time= first_time =cputime cputime; ; format long g %-------------------------------------------%input data for control parameters %-------------------------------------------lengthx=4; lengthx =4; %length of xx-axis side of problem lengthy=2; %length of yy-axis side of problem emodule=1.0; emodule=1.0; ; poisson=0.0; poisson =0.0; fload= fload =-1; %elastic modulus %Poisson's ratio % the total load
B(0, 2)
A(0, 0) D (4, 0)
F 1
C ( 4, 4 1)
12:45
有限单元法
崔向阳
三角形单元程序
lx=16; % number of element in x x-axis ly=8; ly =8; % number of element in yy-axis nel=2*lx* nel =2*lx*ly ly; ; % number of element nnel=3; nnel =3; %number of nodes per element ndof=2; ndof =2; %number of dofs per node nnode=(lx+1)*(ly+1); nnode =(lx+1)*(ly+1); (lx 1) (ly 1); %total number of nodes in system sdof= sdof =nnode* nnode*ndof; ndof; %total system dofs edof= edof =nnel* nnel*ndof; ndof; %degrees of freedom per element %------------------------------------------------------%input data for nodal coordinate values %------------------------------------------------------x0=[]; for i=1:lx+1 for j j=1:ly+1 1:ly+1 x0=[x0; (i(i-1)* 1)*lengthx lengthx/lx /lx end end 12:45
有限元教学程序及使用说明§A-1平面三角形3结点有限元程序
附录有限元教学程序及使用说明§A-1平面三角形３结点有限元程序1、程序名：，2、程序功能本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种荷载的作用，在计算自重时y轴取垂直向上为正；能处理非零已知位移，如支座沉降的作用。
主要输出的内容包括：结点位移、单元应力、主应力、第一主应力与x轴的夹角以及约束结点的支座反力。
程序采用Visual Fortran 编制而成，输入数据全部采用自由格式。
3、程序流程及框图图A-1 程序流程图图A-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 ——要计算支座约束反力的结点数目。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
. .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）最里层：区分为第 奇/偶 数个计算单元刚度的集成：[][][][][][]''''''215656665656266256561661e Z e e e Z e Z e e e e k k k K k k k k k k +⋯++=⇓=⇒==⇒==⇒=⨯⨯⨯⨯⨯⨯边界约束的处理：划0置1法 X Y P X Y P 节点编号 单元编号适用：这种方法适用于边界节点位移分量为已知(含为0)的各种约束。
做法：（1）将总刚矩阵〔K〕中相应于已知位移行主对角线元素置1，其他元素改为零；同时将载荷列阵｛R｝中相应元素用已知位移置换。
◎这样，由该方程求得的此位移值一定等于已知量。
（2）将〔K〕中已知位移相应的列的非主对角成元素也置0，以保持〔K〕的对称性。
◎当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为保证方程成立，须在方程右端减去已知位移对该方程的贡献——已知位移和相应总刚元素的乘积。
◎若约束为零位移约束时，此步则可省去。
特点：（1）经以上处理同样可以消除刚性位移(约束足够的前提下)，去掉未知约束反力。
（2）但这种方法不改变方程阶数，利于存贮。
（3）不过，若是要求出约束反力，仍要重新计算各个划去的总刚元素。
程序如下：变量说明NNODE 单元节点数NPION 总结点数NELEM 单元数NVFIX 受约束边界点数FIXED 约束信息数组NFORCE 节点力数FORCE 节点力数组COORD 结构节点坐标数组LNODS 单元定义数组YOUNG 弹性模量POISS 泊松比THICK 厚度B 单元应变矩阵(3*6)D 单元弹性矩阵(3*3)S 单元应力矩阵(3*6)A 单元面积ESTIF 单元刚度矩阵ASTIF 总体刚度矩阵ASLOD 总体荷载向量ASDISP 节点位移向量ELEDISP 单元节点位移向量STRESS 单元应力%********************************************************** %初始化clearformat short e %设定输出类型clear %清除存变量NELEM=36 %单元个数（单元编码总数）NPION=28 %结点个数（结点编码总数）NVFIX=2 %受约束边界点数NFORCE=1 %结点荷载个数YOUNG=2e11 %弹性模量POISS=0.25 %泊松比THICK=0.1 %厚度LNODS=[1 2 3;2 4 5;2 5 3;3 5 6;4 7 8;4 8 5;5 8 9;5 9 6;6 9 10;7 11 12;7 12 8;8 12 13;8 13 9;9 13 14;9 14 10;10 14 15;11 16 17;11 17 12; 12 17 18; 12 18 13;13 18 19; 13 19 14;14 19 20;14 20 15;15 20 21;16 22 23;16 23 17;17 23 24;17 24 18;18 24 25;18 25 19;25 19 26;19 26 20;20 26 27;20 27 21;21 27 28] %单元定义数组（单元结点号）%相应为单元结点号（编码）、按逆时针顺序输入COORD=[0 0;-0.75 1.5;0.75 1.5;-1.5 3;0 3;1.5 3;-2.25 4.5;-0.75 4.5;0.75 4.5;2.25 4.5;-3 6;-1.5 6;0 6;1.5 6;3 6;-3.75 7.5;-2.25 7.5; -0.75 7.5;0.75 7.5;2.25 7.5;3.75 7.5;-4.5 9;-3 9;-1.5 9;0 9;1.5 9;3 9;4.5 9] %结点坐标数组%坐标：x，y 坐标（共NPOIN 组）FORCE=[1 0 -15] %结点力数组（受力结点编号, x 方向,y 方向）FIXED=[22 1 1;28 1 1] %约束信息（约束点，x 约束，y 约束）%有约束为1，无约束为0%********************************************************** %生成单元刚度矩阵并组成总体刚度矩阵ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0%********************************************************** for i=1:NELEM%生成弹性矩阵 DD= [1 POISS 0;POISS 1 0;0 0 (1-POISS)/2]*YOUNG/(1-POISS^2)%********************************************************** %计算当前单元的面积A=-det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)])/2%********************************************************** %生成应变矩阵 Bfor j=0:2b(j+1)=COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(re m((j+2),3))+1),2);c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(r em((j+2),3))+1),1);endB=[b(1) 0 b(2) 0 b(3) 0;0 c(1) 0 c(2) 0 c(3);c(1) b(1) c(2) b(2) c(3) b(3)]/(2*A);B1( :,:,i)=B;%********************************************************** %求应力矩阵S=D*BS=D*B;ESTIF=B'*S*THICK*A; %求解单元刚度矩阵a=LNODS(i,:); %临时向量,用来记录当前单元的节点编号for j=1:3for k=1:3ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=ASTIF((a(j)*2-1) :a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2);%根据节点编号对应关系将单元刚度分块叠加到总刚%度矩阵中endendend%********************************************************** %将约束信息加入总体刚度矩阵（对角元素改一法）for i=1:NVFIXif FIXED(i,2)==1ASTIF(:,(FIXED(i,1)*2-1))=0; %一列为零ASTIF((FIXED(i,1)*2-1),:)=0; %一行为零ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1; %对角元素为 1end%********************************************************** %生成单元刚度矩阵并组成总体刚度矩阵%**********************************************************if FIXED(i,3)==1ASTIF( :,FIXED(i,1)*2)=0; %一列为零ASTIF(FIXED(i,1)*2,:)=0; %一行为零ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; %对角元素为 1endend%********************************************************** %生成荷载向量ASLOD(1:2*NPION)=0; %总体荷载向量置零for i=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end%********************************************************** %求解力ASDISP=ASTIF\ASLOD' %计算节点位移向量ELEDISP(1:6)=0; %当前单元节点位移向量for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);%取出当前单元的节点位移向量endiSTRESS=D*B1(:, :, i)*ELEDISP' %求力end（程序计算结果和有限元软件得出的结果稍有偏差，可能是程序某些地方数据输入时出了问题，还在寻找具体原因）二、有限元软件分析设置材料参数建模网格划分。