有限元一维杆问题解法及程序

解：第一步——离散对于一维杆问题，我们先离散成单元，对每个单元作如下计算00,,,,00)(22,,,,2222222dxxNNuuKxuNNdxxNNuudxuuNNNNuuxuNNuudxuxdxxuxuxuudxuxdxxuudxxxuujixxjijijixxjijijijixxjijijixxxxxxxxxxjijijijijijijiji其中杆被平均离散为e个单元（有限元不一定要均分），于是有node=e+1个结点，每个单元长度len=1/e,于是第n个单元的左端点坐标lennxi)1(，右端点坐标nlenxj；第二步——刚度矩阵线性插值有每个单元)()(jjejjjeixxelenxxNxxelenxxNeeBeeeeelenBBdxBBKTeeTeexxeji对于整体叠加eeeeeK00020（K为一个node×node阶矩阵）程序用for循环给K赋值K=zeros(node,node);K1=zeros(node,node);forn=1:(node-1);K1(n:n+1,n:n+1)=[e,-e;-e,e];K=K1+K;K1=zeros(node,node);End第三步——力矩阵体积力由第一步公式02dxxNNuuKxuNNjixxjijiji其中第三项为体积力dxxNNFjixxji21用matlab中int（）函数对积分进行计算并用for函数进行循环且赋值，其中)()(jjejjjeixxelenxxNxxelenxxN，nlenxlennxji)1(程序如下Symsx;F1=zeros(1,node);F11=zeros(1,node);G=zeros(1,2);forn=1:e;B=[(xj(n)*x^2-x^3)/len,(x^3-xi(n)*x^2)/len];G=int(B,x,xi(n),xj(n));G=double(G);F11(1,n:n+1)=[G(1,1),G(1,2)];F1=F11+F1;F11=zeros(1,node);End边界力由第一步中公式02dxxNNuuKxuNNjixxjijiji其中第一项ijxxjixxjijixuNNxuNNxuNNF2为边界力因为u（0）处约束力未知为C，u（1）处边界条件11xxu各单元之间的边界力叠加的时候均抵消，所以边界力矩阵最终为102CF第四步——解方程由上面我们可以得到方程02121FuuuKFn,代入位移边界条件01u先对方程进行置一处理，令F=F2-F1且的第一项置0，刚度矩阵变换成eeK00020001方程变换为Kd=F，求逆nuuFKd210第五步——作图用到两个作图函数plot、ezplot分别做原函数图和折线图折线图程序如下x=0:len:1;y=zeros(1,node);forn=1:node;y(1,n)=X(1,n);endplot(x,y,'r');holdon;（holdon可将两图画在一个坐标下）原函数图程序如下ezplot('(1/12)*x^4+2/3*x',[0,1]);附：程序clearformatlongfirst_time=cputime;e=10;%单元数node=e+1;%结点数len=1/e;%单元长度xi=0:len:(1-len);%单元左端点坐标xj=len:len:1;%单元右端点坐标K=zeros(node,node);%刚度矩阵（由于是线性插值）K1=zeros(node,node);forn=1:(node-1);K1(n:n+1,n:n+1)=[e,-e;-e,e];K=K1+K;K1=zeros(node,node);endsymsx;%F1力矩阵F1=zeros(1,node);F11=zeros(1,node);G=zeros(1,2);forn=1:e;B=[(xj(n)*x^2-x^3)/len,(x^3-xi(n)*x^2)/len];G=int(B,x,xi(n),xj(n));G=double(G);F11(1,n:n+1)=[G(1,1),G(1,2)];F1=F11+F1;F11=zeros(1,node);endF2=zeros(1,node);%F2力矩阵F2(1,node)=1;K(2,1)=0;%刚度矩阵置一K(1,2)=0;K(1,1)=1;F=F2-F1;%求合力矩阵F(1,1)=0;X=F*inv(K);%求结点位移Xx=0:len:1;%画折线图y=zeros(1,node);forn=1:node;y(1,n)=X(1,n);endplot(x,y,'r');holdon;ezplot('(1/12)*x^4+2/3*x',[0,1]);%画原函数图以10个单元为例，图数据