fem.zip_FEM_fem matlab_matlab 有限元_三角形单元_有限元编程

平面应力问题

1.算例

一简支梁如图所示：跨度 10m,梁高 1m，梁宽 1m。受均布荷载作用，荷载大小为 10kN/m。

材料为钢材，弹性模量取 2.0E11，泊松比取 0.3。不考虑重力作用。分别利用有限元编程、

手算以及商业软件 Abaqus 进行计算。

2.有限元编程

利用 Matlab 编写平面问题的有限元程序，该程序可以通过读入不同的数据文件对不同的算

例进行求解。该程序现存的问题是：绘制云图时，云图形状为矩形。也就是说，该程序输出

的云图只能应用于矩形平面算例。

3 3 4 15

4 4 5 16

5 5 6 17

6 6 7 18

7 7 8 19

8 8 9 20

9 9 10 21

10 10 11 22

11 12 1 13

12 13 2 14

19 20 9 21

20 21 10 22

结点信息：

NodeInfo=

1 0 0

2 1 0

3 2 0

4 3 0

5 4 0

6 5 0

13 1 1

14 2 1

15 3 1

16 4 1

17 5 1

18 6 1

19 7 1

20 8 1

21 9 1

22 10 1

16 0 -10000

17 0 -10000

18 0 -10000

19 0 -10000

20 0 -10000

21 0 -10000

22 0 -5000

X方向的约束信息

ResX=

1 0

%求解前提：计算者手动划分网格，定节点和单元，手动将荷载信息简化为点荷载

%------------------------------------------------------------------------------

-----

%P 荷载列向量； K 刚度矩阵； U 位移向量； Ua 已知位移向量（边界条件）；Ub 未知

位移; Pb 已知力；Pa 未知力（反力）

%ResX: NodeNum RestrainX;

%ResY: NodeNum RestrainY;

%ResN 被限制的位移总量; ResX X方向边界信息矩阵； ResY Y方向边界信息矩阵；

%BIG 大数

%------------------------------------------------------------------------------

E0=2.0E11;

pRatio0=0.3;

D=E0/(1-pRatio0^2)*[1 pRatio0 0;pRatio0 1 0;0 0 (1-pRatio0)/2];

BIG=1e10;

%-----------------------------------------------------------------------------

K=zeros(2*NodeN);

for Ele_n=1:EleN

%节点信息 i(xi,ji) (i,j,m)

i=EleInfo(Ele_n,2);

j=EleInfo(Ele_n,3);

m=EleInfo(Ele_n,4);

xi=NodeInfo(i,2);

yi=NodeInfo(i,3);

xj=NodeInfo(j,2);

bi=yj-ym;

ci=xm-xj;

aj=xm*yi-xi*ym;

bj=ym-yi;

cj=xi-xm;

am=xi*yj-xj*yi;

bm=yi-yj;

cm=xj-xi;

%应变矩阵

A=(bi*cj-bj*ci)/2;

iter11=K_serial(iter1);

for iter2=1:6

iter22=K_serial(iter2);

K(iter11,iter22)=K(iter11,iter22)+Ke(iter1,iter2);

end

end

end

PNum=Ps(iter,1);

P(2*PNum-1)=P(2*PNum-1)+Ps(iter,2);

P(2*PNum)=P(2*PNum)+Ps(iter,3);

end

Judge1=input('Consider gravity: 1 OR 0\n');

Judge2=input('Please input the value of the volume force (N/m3):');

if Judge1==0

g=0;

end

m=EleInfo(Ele_n,4);

xi=NodeInfo(i,2);

yi=NodeInfo(i,3);

xj=NodeInfo(j,2);

yj=NodeInfo(j,3);

xm=NodeInfo(m,2);

ym=NodeInfo(m,3);

ai=xj*ym-xm*yj;

bi=yj-ym;

am=xi*yj-xj*yi;

bm=yi-yj;