MATLAB 有限元计算程序3

return
function FemModel(a,b,h,m,n,E,mu,load,p,bc)
% 定义有限元模型
% 输入参数：
% a --------- 板x方向宽度
% b --------- 板y方向长度
% h --------- 板z方向厚度
i*(m+1)+j+1,...
i*(m+1)+j ] ;
end
end
gElement( :, 5 ) = 1 ;
% 定义网格密度
m = 16 ; % x方向单元数目
n = 16 ; % y方向单元数目
% 定义材料性质
E = 2.1e11 ; % 弹性模量
mu = 0.3 ; % 泊松比
% 荷载参数
end
for i=1:m
gBC1( (n+1)*2+m+1+i, : ) = [n*(m+1)+i, 1, 0] ; % y=b/2的边界上挠度为零
end
elseif bc == '固支'
gBC1 = zeros( m*4+(n+1)*3+n, 3 ) ;
% step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中
[element_number,dummy] = size( gElement ) ;
for ie=1:1:element_number
disp( sprintf( '计算单元刚度矩阵并集成，当前单元: %d', ie ) ) ;
function exam7_1
% 本程序为第七章的第一个算例，采用薄板矩形单元计算矩形板的弯曲
% 输入参数：
% 无
% 定义板的尺寸
a = 1 ; % 板x方向宽度
b = 1 ; % 板y方向长度
h = 0.01 ; % 板z方向厚度
end
end
% 确定载荷
if load == '均布力'
gNF = [] ;
gDF = gElement ;
gDF(:,5) = p ;
elseif load == '集中力'
gDF = [] ;
f( node*3-1 ) = f( node*3-1 ) + edf( 2:3:11 ) ;
f( node*3 ) = f( node*3 ) + edf( 3:3:12 ) ;
end
for i=1:n+1
gBC1( m+n+(n+1)*2+i, : ) = [i*(m+1), 3, 0 ] ; % x=a/2的边界上绕y轴的转角为零
end
for i=1:m
gBC1( m+n+(n+1)*3+i, : ) = [n*(m+1)+i, 1, 0] ; % y=b/2的边界上挠度为零
% 定义单元
gElement = zeros( n*n, 5 ) ;
for i=1:n
for j=1:m
gElement( (i-1)*m+j, 1:4) = [ (i-1)*(m+1)+j, ...
(i-1)*(m+1)+j+1, ...
% m --------- x方向单元数目
% n --------- y方向单元数目
% E --------- 弹性模量
% mu --------- 泊松比
% load ------- 载荷类型
% p --------- 载荷大小
% 输出文件名称
file_out = 'exam7_1.mat' ;
% 检查输出文件是否存在
file_prompt = 0 ;
if exist( file_out ) ~= 0 & file_prompt == 1
answer = input( sprintf( '文件 %s 已经存在，是否覆盖? ( Yes / [No] ): ', file_out ), 's' ) ;
disp( sprintf( '请使用另外的文件名，或备份已有的文件' ) ) ;
disp( sprintf( '程序终止' ) );
return ;
end
end
FemModel(a,b,h,m,n,E,mu,load,p,bc) ; % 定义有限元模型
gNode = zeros( (m+1)*(n+1), 2 ) ;
for i=1:n+1
for j=1:m+1
gNode( (i-1)*(m+1)+j, : ) = [dx*(j-1),dy*(i-1)] ;
end
end
if length( answer ) == 0
answer = 'no' ;
end
answer = lower( answer ) ;
if answer ~= 'y' | answer ~= 'yes'
end
for i=1:n+1
gBC1( m+n+i, : ) = [i*(m+1), 1, 0 ] ; % x=a/2的边界上挠度为零
end
for i=1:n+1
gBC1( m+n+n+1+i, : ) = [i*(m+1), 2, 0 ] ; % x=a/2的边界上绕x轴的转角为零
for idf = 1:1:df_number
edf = EquivalentDistPressure( gDF(idf,1:4), gDF(idf,5) ) ;
node = gDF( idf, 1:4 ) ;
f( node*3-2 ) = f( node*3-2 ) + edf( 1:3:10 ) ;
load = '均布力' ; % 也可以是 集中力
%load = '集中力' ;
p = 1e6 ; % 载荷大小
% 边界条件
%bc = '简支' ; % 也可以是 固支
bc = '固支' ;
k = StiffnessMatrix( ie ) ;
AssembleStiffnessMatrix( ie, k ) ;
end
% step3. 计算分布压力的等效节点力，并集成到整体节点力向量中
[df_number, dummy] = size( gDF ) ;
end
for i=1:n+1
gBC1( (m+1)+i, : ) = [(i-1)*(m+1)+1, 3, 0] ; % y=0的边界上绕y轴的转角等于零
end
for i=1:n+1
gBC1( m+n+2+i, : ) = [i*(m+1), 1, 0 ] ; % x=a/2的边界上挠度为零
% 定义材料
gMaterial = [ E, mu, h] ;
% 确定边界条件
if bc == '简支'
gBC1 = zeros( (n+1)*2+m*2+1, 3 ) ;
for i=1:m+1
gBC1( i, : ) = [i, 2, 0] ; % x=0的边界上绕x轴的转角等于零
% step1. 定义整体刚度矩阵和节点力向量
[node_number,dummy] = size( gNode ) ;
gK = sparse( Fra bibliotekode_number * 3, node_number * 3 ) ;
f = sparse( node_number * 3, 1) ;
SolveModel ; % 求解有限元模型
SaveResults( file_out ) ; % 保存计算结果
DisplayResults(a,b,h,m,n,E,mu,load,p,bc) ; % 把计算结果显示出来
% 2. 计算单元的等效节点力，集成整体节点力向量
% 3. 处理约束条件，修改整体刚度矩阵和节点力向量
% 4. 求解方程组，得到整体节点位移向量
global gNode gElement gMaterial gBC1 gDF gNF gK gDelta gElementMoment gNodeMoment
for i=1:m
gBC1( i, : ) = [i, 2, 0] ; % x=0的边界上绕x轴的转角等于零
end
for i=1:n
gBC1( m+i, : ) = [(i-1)*(m+1)+1, 3, 0] ; % y=0的边界上绕y轴的转角等于零
% bc --------- 边界条件
% 返回值：
% 无
global gNode gElement gMaterial gBC1 gDF gNF gK
% 计算结点坐标
dx = a / m / 2 ; % 取板的1/4进行分析
dy = b / n / 2 ;