平面三角形单元有限元程序的设计说明

实验三：三角形单元的形函数性质验证
一、 实验目的
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节点等参单元的有限元程序：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)的各种约束。

基于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
请注意，附件的具体文件名和数据内容将根据你的实际情况进行调整。

法律名词及注释部分需要根据实际需要进行扩充和修改。

请合理对文档结构和章节内容进行调整，以符合你的具体要求。