用matlab表示e有限元,matlab网格划分程序与matlab有限元的结合

1.

2.matlab有限元可以参考徐荣桥的书

3.这里本人打算画一个园内包含一个椭圆的模型：

具体程序如下：

a.

网格划分：

>>fh=@(p) 0.05+0.3*dellipse(p,[0.5,0.2]);

>> fd=@(p) ddiff(dcircle(p,0,0,1),dellipse(p,[0.5,0.2]));

>> [p,t]=distmesh2d(fd,fh,0.05,[-1,-1;1,1],[-1,-1;-1,1;1,-1;1,1]);

b.在distmesh2d.m最后插入语句，导出数据格式(节点坐标，单元材料属性，边界约束条件)

fid = fopen('exam4_2.dat3', 'w') ;

[ilength,jlength]=size( p );

fprintf(fid,'%d\n',ilength);

for i=1:1:ilength

fprintf(fid,'%d\t%f\t%f\n',i,p(i,1),p(i,2));

end

[ilength,jlength]=size( t );

fprintf(fid,'%d\n',ilength);

for i=1:1:ilength

fprintf(fid,'%d\t%d\t%d\t%d\t%d\n',i,t(i,1),t(i,2),t(i,3),1);

end

fclose(fid);

文件'exam4_2.dat3内容如下：

5105(节点个数)

1 -0.707106 -0.707106 (每个节点坐标)

2 -0.707106 0.707107

...

9869(三角单元个数)

1 5006 4951 4934 1  (标号，i，j，k，材料属性的编号)

2 197 147 100 1

3 413 2 366 1

4 413 2 323 1

...

5(5种材料属性)

1   2.00e9    0.25  100.0  22.0e3(编号，杨氏模量，泊松比，厚度[平面应力填1]，密度)2   2.60e9    0.20  1.0  23.0e3

3   2.85e10   0.20  1.0  25.0e3

4   1.85e10   0.20  1.0  23.0e3

5   2.85e10   0.20  1.0  22.0e3

2(约束个数)

1 2  0.0 (节点，方向【2，为y方向】，0.0(0为固定))

1 1  0.0

4 2  0.0 (节点，方向【2，为y方向】，0.0(0为固定))

41  0.0

c。调用有限元程序exam_2.m(见徐荣桥书)

>>exam4_2 exam4_2.dat3

d.后处理

>> exam4_2_post

最大主应力云图如下：

最大剪应力云图如下：

最小主应力云图如下：

exam_2。m源程序如下：

function exam4_2( file_in )

% 本程序为第四章的第二个算例，采用三角形单元计算隧道在自重作用下的变形和应力

%      exam4_2( filename )

%  输入参数：

%      file_in  ---------- 有限元模型文件

% 定义全局变量

%      gNode ------------- 节点坐标

%      gElement ---------- 单元定义

%      gMaterial --------- 材料性质

%      gBC1 -------------- 第一类约束条件

%      gK ---------------- 整体刚度矩阵

%      gDelta ------------ 整体节点坐标

%      gNodeStress ------- 节点应力

%      gElementStress ---- 单元应力

global gNode gElement gMaterial gBC1 gK gDelta gNodeStress gElementStress gNF

if nargin < 1