matlab⽹格划分程序与matlab有限元的结合 1.
2.matlab有限元可以参考徐荣桥的书
3.这⾥本⼈打算画⼀个园内包含⼀个椭圆的模型：
具体程序如下：
a.
⽹格划分：
>>) 0.05+0.3*dellipse(p,[0.5,0.2]);
>> ) 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
file_in = 'exam4_2.dat' ;
end
% 检查⽂件是否存在
if exist( file_in ) == 0
disp( sprintf( '错误：⽂件 %s 不存在', file_in ) )
disp( sprintf( '程序终⽌' ) )
return ;
end
% 根据输⼊⽂件名称⽣成输出⽂件名称
[path_str,name_str,ext_str] = fileparts( file_in ) ;
ext_str_out = '.mat' ;
file_out = fullfile( path_str, [name_str, ext_str_out] ) ;
% 检查输出⽂件是否存在
if exist( file_out ) ~= 0
answer = input( sprintf( '⽂件 %s 已经存在，是否覆盖? ( Yes / [No] ):  ', file_out ), 's' ) ;        if length( answer ) == 0
answer = 'no' ;
end
answer = lower( answer ) ;
if answer ~= 'y' | answer ~= 'yes'
disp( sprintf( '请使⽤另外的⽂件名，或备份已有的⽂件' ) ) ;
disp( sprintf( '程序终⽌' ) );
return ;
end
end
% 建⽴有限元模型并求解，保存结果
FemModel( file_in ) ;          % 定义有限元模型
SolveModel ;                  % 求解有限元模型
SaveResults( file_out ) ;      % 保存计算结果
% 计算结束
disp( sprintf( '计算正常结束，结果保存在⽂件 %s 中', file_out ) ) ;
disp( sprintf( '可以使⽤后处理程序 exam4_2_post.m 显⽰计算结果' ) ) ;
return ;
function FemModel(filename)
fprintf作用
%  定义有限元模型
%  输⼊参数：
%      filename --- 有限元模型⽂件
%  返回值：
%      ⽆
%  说明：
%      该函数定义平⾯问题的有限元模型数据：
%        gNode ------- 节点定义
%        gElement ---- 单元定义
%        gMaterial --- 材料定义，包括弹性模量，梁的截⾯积和梁的抗弯惯性矩
%        gBC1 -------- 约束条件
global gNode gElement gMaterial gBC1 gNF
% 打开⽂件
fid = fopen( filename, 'r' ) ;
% 读取节点坐标
node_number = fscanf( fid, '%d', 1 ) ;
gNode = zeros( node_number, 2 ) ;
for i=1:node_number
dummy = fscanf( fid, '%d', 1 ) ;
gNode( i, : ) = fscanf( fid, '%f', [1, 2] ) ;
end
% 读取单元定义
element_number = fscanf( fid, '%d', 1 ) ;
gElement = zeros( element_number, 4 ) ;
for i=1:element_number
dummy = fscanf( fid, '%d', 1 ) ;
gElement( i, : ) = fscanf( fid, '%d', [1, 4] ) ;    end
% 读取材料信息
material_number = fscanf( fid, '%d', 1 ) ;
gMaterial = zeros( material_number, 4 ) ;
for i=1:material_number
dummy = fscanf( fid, '%d', 1 ) ;
gMaterial( i, : ) = fscanf( fid, '%f', [1,4] ) ;    end
% 读取边界条件
bc1_number = fscanf( fid, '%d', 1 ) ;
gBC1 = zeros( bc1_number, 3 ) ;
for i=1:bc1_number
gBC1( i, 1 ) = fscanf( fid, '%d', 1 ) ;
gBC1( i, 2 ) = fscanf( fid, '%d', 1 ) ;
gBC1( i, 3 ) = fscanf( fid, '%f', 1 ) ;
end
% 关闭⽂件
fclose( fid ) ;
% 集中⼒
%    节点号  ⾃由度号  集中⼒值
gNF = [  1,      1,        -800e3;
1,      2,        -800e3;
4,      1,        -800e3;
4,      2,        -800e3;
]
;
return