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

解：

第一步——离散

对于一维杆问题，我们先离散成单元，对每个单元作如下计算

[][

][][

][]00,,,,0

0)(22

,,,,2222222=??

????-??????-?????????=??

????-??

?

?????????-?????????=-????-???=-??=-????

?

?????ΓΓ

Γdx x N N u u K x u N N dx x N N u u dx u u N N N N u u x u

N N u u dx ux dx x

u x u x u u dx ux dx x u u dx x x u u j i x x j i j i j i x x j

i j i j i j i x x j

i j i j i x

x x x x x x x x x j i j

i

j

i

j i j i j i

j i j

i

δδδδδδδδδδδδ

其中杆被平均离散为e 个单元（有限元不一定要均分），于是有node=e+1个结点，每个单元长度len=1/e,于是第n 个单元的左端点坐标len n x i )1(-=，右端点坐标nlen x j =；

第二步——刚度矩阵

线性插值有每个单元

)

()

(j j e

j j j e

i x x e len

x x N x x e len

x x N -=-=

-=-=????

???-=e e B e

??

?

????--===?e e e e len B B dx B B K T e e T e e x

x e j i 对于整体叠加

?????

?

?

??

???--=e e e e e K 00020 （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

第三步——力矩阵

体积力

由第一步公式

[]02

=??

????-??????-???

??????Γ

dx x N N u u K x u

N N j i x x j i j i j i

其中第三项为体积力dx

x N N F j i x x j i 21??

????=?

用matlab 中int （）函数对积分进行计算并用for 函数进行循环且赋值，其中

)

()(j j e

j j j e

i x x e len

x x N x x e len x x N -=-=

-=-=，nlen

x len

n x j i =-=)1(

程序如下

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

边界力

由第一步中公式

[]02

=??

????-??????-?????????Γ

dx x N N

u u K x u

N N j i x x j i j i j i 其中第一项

i

j

x x j i x x j i j i x

u

N N x

u N N x u N N F ==Γ????????-???

?????=????????=2为边界力

因为u （0）处约束力未知为C ，u （1）处边界条件

11

=??=x x

u

各单元之间的边界力叠加的时候均抵消，所以边界力矩阵最终为????

?

???????=102 C F

第四步——解方程

由上面我们可以得到方程

02121=-????

?

???????-F u u u K F n ,代入位移边界条件01=u

先对方程进行置一处理，令F=F2-F1且的第一项置0，刚度矩阵变换成

????

?

????

???=e e K 00020001方程变换为Kd=F ， 求逆?????

???????==-n u u F K d 21

第五步——作图

用到两个作图函数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; %刚度矩阵置一

K(1,2)=0;

K(1,1)=1;

F=F2-F1; %求合力矩阵

F(1,1)=0;

X=F*inv(K); %求结点位移

X

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;

ezplot('(1/12)*x^4+2/3*x',[0,1]); %画原函数图以10个单元为例，图

数据