一用三结点三角形平面单元计算平面结构的应力和位移讲解

一：用三结点三角形平面单元计算平面结构的应力和位移。

1，设计说明书

计算简图，网格划分，单元及结点的编号如下图所示。由于结构对称，去四分之一结构分析。其中

E=2e10pa,mu=0.167,h=1m.

变量注释：

Node ------- 节点定义

gElement ---- 单元定义

gMaterial --- 材料定义，包括弹性模量，泊松比和厚度

gBC1 -------- 约束条件

gNF --------- 集中力

gk------------总刚

gDelta-------结点位移

子程序注释：

PlaneStructualModel ———定义有限元模型

SolveModel ———————求解有限元模型

DisplayResults ——————显示计算结果

k = StiffnessMatrix( ie )———计算单元刚度

AssembleStiffnessMatrix( ie, k )—形成总刚

es = ElementStress( ie )————计算单元应力

function exam1

% 输入参数：无

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

PlaneStructualModel ; % 定义有限元模型

SolveModel ; % 求解有限元模型

DisplayResults ; % 显示计算结果

return ;

function PlaneStructualModel

% 定义平面结构的有限元模型

% 输入参数：无

% 说明：

% 该函数定义平面结构的有限元模型数据：

% gNode ------- 节点定义

% gElement ---- 单元定义

% gMaterial --- 材料定义，包括弹性模量，泊松比和厚度% gBC1 -------- 约束条件

% gNF --------- 集中力

global gNode gElement gMaterial gBC1 gNF

% 节点坐标

% x y

gNode = [0.0, 2.0 % 节点1

0.0, 1.0 % 节点2

1.0, 1.0 % 节点3

0.0, 0.0 % 节点4

1.0, 0.0 % 节点5

2.0, 0.0] ; % 节点6

% 单元定义

% 节点1 节点2 节点3 材料号

gElement = [3, 1, 2, 1 % 单元1

5, 2, 4, 1 % 单元2

2, 5, 3, 1 % 单元3

6, 3, 5, 1]; % 单元4 % 材料性质

% 弹性模量泊松比厚度

gMaterial = [1e0, 0, 1] ; % 材料1

% 第一类约束条件

% 节点号自由度号约束值

gBC1 = [ 1, 1, 0.0

2, 1, 0.0

4, 1, 0.0

4, 2, 0.0

5, 2, 0.0

6, 2, 0.0] ;

% 集中力

% 节点号自由度号集中力值

gNF = [ 1, 2, -1] ;

return

function SolveModel

% 求解有限元模型

% 输入参数：无

% 说明：

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

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

% 2. 计算单元的等效节点力，集成整体节点力向量% 3. 处理约束条件，修改整体刚度矩阵和节点力向量% 4. 求解方程组，得到整体节点位移向量

global gNode gElement gMaterial gBC1 gNF gK gDelta

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

[node_number,dummy] = size( gNode ) ;

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

f = sparse( node_number * 2, 1 ) ;

% step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中[element_number,dummy] = size( gElement ) ;

for ie=1:1:element_number

k = StiffnessMatrix( ie ) ;

AssembleStiffnessMatrix( ie, k ) ;

end

% step3. 把集中力直接集成到整体节点力向量中

[nf_number, dummy] = size( gNF ) ;

for inf=1:1:nf_number

n = gNF( inf, 1 ) ;

d = gNF( inf, 2 ) ;

f( (n-1)*2 + d ) = gNF( inf, 3 ) ;

end

% step4. 处理约束条件，修改刚度矩阵和节点力向量。采用乘大数法[bc_number,dummy] = size( gBC1 ) ;

for ibc=1:1:bc_number

n = gBC1(ibc, 1 ) ;

d = gBC1(ibc, 2 ) ;

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

f(m) = gBC1(ibc, 3)* gK(m,m) * 1e15 ;

gK(m,m) = gK(m,m) * 1e15 ;

end

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

gDelta = gK \ f ;

return

function DisplayResults

% 显示计算结果

% 输入参数：无

global gNode gElement gMaterial gBC1 gNF gK gDelta

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

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

[node_number,dummy] = size( gNode ) ;

for i=1:node_number

fprintf( '%6d %16.8e %16.8e\n',...

i, gDelta((i-1)*2+1), gDelta((i-1)*2+2)) ;

end

fprintf( '\n\n单元应力\n' ) ;

fprintf( ' X-STR Y-STR XY-STR\n' ) ;

[element_number, dummy] = size( gElement ) ;

for ie = 1:element_number

es = ElementStress( ie ) ;

fprintf( '单元号%6d %16.8e %16.8e %16.8e\n', ...

ie, es(1), es(2), es(3));

end

return

function k = StiffnessMatrix( ie )

% 计算单元刚度矩阵

% 输入参数:

% ie ------- 单元号

global gNode gElement gMaterial

k = zeros( 6, 6 ) ;

E = gMaterial( gElement(ie, 4), 1 ) ;