MATLAB有限元计算程序2

MATLAB有限元计算程序2

function exam3_2

% 本程序为第三章的第二个算例，采用平面梁单元计算斜腿刚架桥的变形和内力

% 输入参数： 无

% 输出结果： 节点位移和单元节点力

% 检查数据文件是否存在

file_in = 'exam3_2.dat' ;

if exist( file_in ) == 0

disp( sprintf( '错误：文件 %s 不存在', file_in ) )

disp( sprintf( '程序终止' ) )

return ;

end

% 读入模型，求解并显示结果

PlaneFrameModel(file_in) ; % 定义有限元模型

SolveModel ; % 求解有限元模型

DisplayResults ; % 显示计算结果

return ;

function PlaneFrameModel( file_in )

% 定义平面杆系的有限元模型

% 输入参数：

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

% 返回值：

% 无

% 说明：

% 该函数定义平面杆系的有限元模型数据： % gNode ------- 节点定义

% gElement ---- 单元定义

% gMaterial --- 材料定义，包括弹性模量，梁的截面积和梁的抗弯惯性矩

% gBC1 -------- 第一类边界条件

% gBC2 -------- 第二类约束条件, 处理铰接的结点

global gNode gElement gMaterial gBC1 gBC2

% 打开数据文件

fid_in = fopen( file_in, 'r' ) ;

% 读取节点坐标和节点温度

node_number = fscanf( fid_in, '%d', 1 ) ;

gNode = zeros( node_number, 3 ) ;

for i=1:node_number

dummy = fscanf( fid_in, '%d', 1 ) ;

gNode(i,:) = fscanf( fid_in, '%f', [1 3] ) ;

end

% 读取单元定义

element_number = fscanf( fid_in, '%d', 1 ) ;

gElement = zeros( element_number, 3 ) ;

for i=1:element_number

dummy = fscanf( fid_in, '%d', 1 ) ;

gElement(i,:) = fscanf( fid_in, '%d', [1 3] ) ;

end

% 读取材料性质

material_number = fscanf( fid_in, '%d', 1 ) ; gMaterial = zeros( material_number, 5 ) ;

for i=1:material_number

dummy = fscanf( fid_in, '%d', 1 ) ;

gMaterial(i,:) = fscanf( fid_in, '%f', [1 5] ) ;

end

% 读取第一类约束条件

bc1_number = fscanf( fid_in, '%d', 1 ) ;

gBC1 = zeros( bc1_number, 3 ) ;

for i=1:bc1_number

gBC1(i,1) = fscanf( fid_in, '%d', 1 ) ;

gBC1(i,2) = fscanf( fid_in, '%d', 1 ) ;

gBC1(i,3) = fscanf( fid_in, '%f', 1 ) ;

end

% 读取第二类约束条件

bc2_number = fscanf( fid_in, '%d', 1 ) ;

gBC2 = zeros( bc2_number, 3 ) ;

for i=1:bc2_number

gBC2(i,:) = fscanf( fid_in, '%d', [1 3] ) ;

end

% 关闭数据文件

fclose( fid_in ) ;

return

function SolveModel

% 求解有限元模型

% 输入参数：

% 无 % 返回值：

% 无

% 说明：

% 该函数求解有限元模型，过程如下

% 1. 计算单元刚度矩阵，集成整体刚度矩阵

% 2. 计算单元的等效节点力，

集成整体节点力向量

% 3. 处理约束条件，修改整体刚度矩阵和节点力向量

% 4. 求解方程组，得到整体节点位移向量

global gNode gElement gMaterial gBC1 gBC2 gK gDelta

% step1. 定义整体刚度矩阵和节点力向量

[node_number,dummy] = size( gNode ) ;

gK = sparse( node_number * 3, node_number * 3 ) ;

f = sparse( node_number * 3, 1 ) ;

% step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中

[element_number,dummy] = size( gElement ) ;

for ie=1:1:element_number

k = StiffnessMatrix( ie, 1 ) ;

AssembleStiffnessMatrix( ie, k ) ;

end

% step3. 计算自重温度变化产生的等效节点力

for ie=1:1:element_number

egf = EquivalentGravityForce( ie ) ;

etf = EquivalentThermalForce( ie ) ;

i = gElement( ie, 1 ) ;

j = gElement( ie, 2 ) ; f( (i-1)*3+1 : (i-1)*3+3 ) = f( (i-1)*3+1 : (i-1)*3+3 ) + egf( 1:3 )

+ etf( 1:3 );

f( (j-1)*3+1 : (j-1)*3+3 ) = f( (j-1)*3+1 : (j-1)*3+3 ) + egf( 4:6 )

+ etf( 4:6 );

end

% step4. 处理第二类约束条件，修改刚度矩阵和节点力向量。（采用划行划列法）

[bc2_number,dummy] = size( gBC2 ) ;

for ibc=1:1:bc2_number

n1 = gBC2(ibc, 1 ) ;

n2 = gBC2(ibc, 2 ) ;

d = gBC2(ibc, 3 ) ;

m1 = (n1-1)*3 + d ;

m2 = (n2-1)*3 + d ;

f(m2) = f(m2)+f(m1) ;

f(m1) = 0 ;

gK(m2,:) = gK(m2,:)+gK(m1,:) ;

gK(:,m2) = gK(:,m2)+gK(:,m1) ;

gK(m1,:) = zeros(1, node_number*3) ;

gK(:,m1) = zeros(node_number*3, 1) ;

gK(m1,m1) = 1.0 ;

gK(m1,m2) = -1 ;

gK(m2,m1) = -1 ;

gK(m2,m2) = gK(m2,m2)+1 ;

end

% step5. 处理第一类约束条件，修改刚度矩阵和节点力向量。（采用划行划列法）

[bc1_number,dummy] = size( gBC1 ) ;

for ibc=1:1:bc1_number

n = gBC1(ibc, 1 ) ;

d = gBC1(ibc, 2 ) ;

m = (n-1)*3 + d ;

f = f - gBC1(ibc,3) * gK(:,m) ;

f(m) = gBC1(ibc,3) ;

gK(:,m) = zeros( node_number*3, 1 ) ;

gK(m,:) = zeros( 1, node_number*3 ) ;

gK(m,m) = 1.0 ;

end

% step 7. 求解方程组，得到节点位移向量

gDelta = gK \ f ;

return

function DisplayResults

% 显示计算结果

% 输入参数：

% 无

% 返回值：

% 无

global gNode gElement gMaterial gBC1 gNF gDF gK gDelta

fprintf( '节点位移\n' ) ;

fprintf( ' 节点号 x方向位移 y方向位移 转角\n' ) ;

[node_number,dummy] = size( gNode ) ;