第9章有限单元法(附录)

1第9章附录9.2.2输入数据%demo_fem_2.mfunctiondemo_fem_2clc;clear;nel=15;nnel=2;ndof=1;nnode=nel+1;sdof=nnode*ndof;gcoord(1:nnode)=0.0:1.0/nel:1.0;nodes(1:nel,1)=1:nel;nodes(1:nel,2)=2:nel+1;acoef=1;bcoef=-3;ccoef=2;%΢·Ö·½³ÌϵÊýbcdof(1)=1;bcval(1)=0;%±ß½çff=zeros(sdof,1);kk=zeros(sdof,sdof);index=zeros(nnel*ndof,1);foriel=1:nel;nl=nodes(iel,1);nr=nodes(iel,2);xl=gcoord(nl);xr=gcoord(nr);eleng=xr-xl;index=feeldof1(iel,nnel,ndof);k=feode2l(acoef,bcoef,ccoef,eleng);f=fef1l(xl,xr);[kk,ff]=feasmbl2(kk,ff,k,f,index);endff(nnode)=ff(nnode)-1;[kk,ff]=feaplyc2(kk,ff,bcdof,bcval);fsol=kk\ff;c1=(1+0.5*exp(1))/(2*exp(2)-exp(1));c2=-(1+exp(2))/(2*exp(2)-exp(1));fori=1:nnode2x=gcoord(i);esol(i)=c1*exp(2*x)+c2*exp(x)+1/2;endnum=1:1:sdof;results=[num'fsolesol'];plot(num,fsol,'rs',num,esol,'k-');endfunction[kk,ff]=feaplyc2(kk,ff,bcdof,bcval)n=length(bcdof);sdof=size(kk);fori=1:nc=bcdof(i);forj=1:sdofkk(c,j)=0;endkk(c,c)=1;ff(c)=bcval(i);endendfunction[kk,ff]=feasmbl2(kk,ff,k,f,index)edof=length(index);fori=1:edofii=index(i);ff(ii)=ff(ii)+f(i);forj=1:edofjj=index(j);kk(ii,jj)=kk(ii,jj)+k(i,j);endendendfunction[index]=feeldof1(iel,nnel,ndof)edof=nnel*ndof;start=(iel-1)*(nnel-1)*ndof;fori=1:edof3index(i)=start+i;endendfunction[f]=fef1l(xl,xr)eleng=xr-xl;f=[eleng/2;eleng/2];endfunction[k]=feode2l(acoef,bcoef,ccoef,eleng)a1=-(acoef/eleng);a2=bcoef/2;a3=ccoef*eleng/6;k=[a1-a2+2*a3-a1+a2+a3;-a1-a2+a3a1+a2+2*a3];end