节点应力和单元应力的关系
一：用三结点三角形平面单元计算平面结构的应力和位移。
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 ) ;
mu = gMaterial( gElement(ie, 4), 2 ) ;
h = gMaterial( gElement(ie, 4), 3 ) ;
xi = gNode( gElement( ie, 1 ), 1 ) ;
yi = gNode( gElement( ie, 1 ), 2 ) ;
xj = gNode( gElement( ie, 2 ), 1 ) ;
yj = gNode( gElement( ie, 2 ), 2 ) ;
xm = gNode( gElement( ie, 3 ), 1 ) ;
ym = gNode( gElement( ie, 3 ), 2 ) ;
ai=xj*ym-xm*yj;
aj=xm*yi-xi*ym;
am=xi*yj-xj*yi;
bi=yj-ym;
bj=ym-yi;
bm=yi-yj;
ci=-(xj-xm);
cj=-(xm-xi);
cm=-(xi-xj);
A=(ai+aj+am)/2;
B=[bi 0 bj 0 bm 0
0 ci 0 cj 0 cm
ci bi cj bj cm bm];
B=B/2/A;
D=[1 mu 0
mu 1 0
0 0 (1-mu)/2];
D=D*E/(1-mu^2);
k=transpose(B)*D*B*h*abs(A);
return
function AssembleStiffnessMatrix( ie, k )
% 把单元刚度矩阵集成到整体刚度矩阵
% 输入参数:
% ie --- 单元号
% k --- 单元刚度矩阵
global gElement gK
for i=1:1:3
for j=1:1:3
for p=1:1:2
for q =1:1:2
m = (i-1)*2+p ;
n = (j-1)*2+q ;
M = (gElement(ie,i)-1)*2+p ;
N = (gElement(ie,j)-1)*2+q ;
gK(M,N) = gK(M,N) + k(m,n) ;
end
end
end
end
return
function es = ElementStress( ie )
% 计算单元的应力
% 输入参数
% ie ----- 节点号
% es ----- 单元应力
global gElement gNode gDelta gMaterial
es=zeros(1,6);
de=zeros(6,1);
for j=1:1:3
de(2*j-1)=gDelta(2*gElement(ie,j)-1);
de(2*j)=gDelta(2*gElement(ie,j));
end
E = gMaterial( gElement(ie, 4), 1 ) ;
mu = gMaterial( gElement(ie, 4), 2 ) ;
h = gMaterial( gElement(ie, 4), 3 ) ;
xi = gNode( gElement( ie, 1 ), 1 ) ;
yi = gNode( gElement( ie, 1 ), 2 ) ;
xj = gNode( gElement( ie, 2 ), 1 ) ;
yj = gNode( gElement( ie, 2 ), 2 ) ;
xm = gNode( gElement( ie, 3 ), 1 ) ;
ym = gNode( gElement( ie, 3 ), 2 ) ;
ai=xj*ym-xm*yj;
aj=xm*yi-xi*ym;
am=xi*yj-xj*yi;
bi=yj-ym;
bj=ym-yi;
bm=yi-yj;
ci=-(xj-xm);
cj=-(xm-xi);
cm=-(xi-xj);
A=(ai+aj+am)/2;
B=[bi 0 bj 0 bm 0
0 ci 0 cj 0 cm
ci bi cj bj cm bm];
B=B/2/A;
D=[1 mu 0
mu 1 0
0 0 (1-mu)/2];
D=D*E/(1-mu^2);
S=D*B;
es(1:3)=S*de;
es(6)=0.5*sqrt((es(1)-es(2))^2+4*es(3)^2);
es(4)=0.5*(es(1)+es(2))+es(6);
es(5)=0.5*(es(1)+es(2))-es(6);
return
3，数据文件：
输入数据：
gNode = [0.0, 2.0
0.0, 1.0
1.0, 1.0
0.0, 0.0
1.0, 0.0
2.0, 0.0] ;
gElement = [3, 1, 2, 1
5, 2, 4, 1
2, 5, 3, 1
6, 3, 5, 1];
gMaterial = [1e0, 0, 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] ;
输出数据：
节点位移
节点号x方向位移y方向位移
1 -8.79120879e-016 -3.25274725e+000
2 8.79120879e-017 -1.25274725e+000
3 -8.79120879e-002 -3.73626374e-001
4 1.17216117e-016 -8.35164835e-016
5 1.75824176e-001 -2.93040293e-016
6 1.75824176e-001 2.63736264e-016
单元应力
X-STR Y-STR
XY-STR
单元号 1 -8.79120879e-002 -2.00000000e+000 4.39560440e-001
单元号 2 1.75824176e-001 -1.25274725e+000 2.56410256e-016
单元号 3 -8.79120879e-002 -3.73626374e-001 3.07692308e-001
单元号 4 0.00000000e+000 -3.73626374e-001 -1.31868132e-001