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

解： 第一步——离散 对于一维杆问题，我们先离散成单元，对每个单元作如下计算 其中杆被平均离散为e个单元（有限元不一定要均分），于是有node=e+1个结点，每个单元长度len=1/e,于是第n个单元的左端点坐标，右端点坐标； 第二步——刚度矩阵 线性插值有每个单元 对于整体叠加 （K为一个node×node阶矩阵） 程序用for循环给K赋值 K=zeros(node,node); K1=zeros(node,node); for n=1:(node-1); K1(n:n+1,n:n+1)=[e,-e;-e,e]; K=K1+K; K1=zeros(node,node); End 第三步——力矩阵 体积力 由第一步公式 其中第三项为体积力 用matlab中int（）函数对积分进行计算并用for函数进行循环且赋值，其中 ， 程序如下 Syms x; F1=zeros(1,node); F11=zeros(1,node); G=zeros(1,2); for n=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 边界力 由第一步中公式 其中第一项为边界力 因为u（0）处约束力未知为C，u（1）处边界条件 各单元之间的边界力叠加的时候均抵消，所以边界力矩阵最终为 第四步——解方程 由上面我们可以得到方程 ,代入位移边界条件 先对方程进行置一处理，令F=F2-F1且的第一项置0，刚度矩阵变换成 方程变换为Kd=F， 求逆 第五步——作图 用到两个作图函数plot、ezplot 分别做原函数图和折线图 折线图 程序如下 x=0:len:1; y=zeros(1,node); for n=1:node; y(1,n)=X(1,n); end plot(x,y,r); hold on; （hold on可将两图画在一个坐标下） 原函数图 程序如下 ezplot((1/12)*x^4+2/3*x,[0,1]); 附： 程序 clear format long first_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); for n=1:(node-1); K1(n:n+1,n:n+1)=[e,-e;-e,e]; K=K1+K; K1=zeros(node,node); end syms x; %F1力矩阵 F1=zeros(1,node); F11=zeros(1,node); G=zeros(1,2); for n=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 F2=zeros(1,node); %F2力矩阵 F2(1,node)=1; K(2,1)=0;