90 lines
2.7 KiB
Matlab
90 lines
2.7 KiB
Matlab
function SolveModel
|
|
% 求解有限元模型
|
|
% 该函数求解有限元模型，过程如下
|
|
% 1. 计算单元刚度矩阵，集成整体刚度矩阵
|
|
% 2. 计算单元的等效节点力，集成整体节点力向量
|
|
% 3. 处理约束条件，修改整体刚度矩阵和节点力向量
|
|
% 4. 求解方程组，得到整体节点位移向量
|
|
% 5. 计算单元应力和节点应力
|
|
|
|
global gNode gElement gMaterial gBC1 gK gDelta gNodeStress gElementStress gDF
|
|
|
|
% step 1. 定义整体刚度矩阵和节点力向量
|
|
[node_number,dummy] = size( gNode ) ;
|
|
gK = sparse( node_number * 2, node_number * 2 ) ;
|
|
f = sparse( node_number * 2, 1 ) ;
|
|
|
|
% step 2. 计算单元刚度矩阵，并集成到整体刚度矩阵中
|
|
[element_number, dunmmy] = size( gElement ) ;
|
|
for ie=1:1:element_number
|
|
k = StiffnessMatrix( ie ) ;
|
|
AssembleStiffnessMatrix( ie, k ) ;
|
|
end
|
|
|
|
% step 3. 计算地面超载产生的等效节点力
|
|
[df_number,dummy] = size( gDF ) ;
|
|
for idf = 1:1:df_number
|
|
enf = EquivalentNodeForce( gDF(idf,1), gDF(idf,2), gDF(idf,3), gDF(idf,4) ) ;
|
|
i = gElement( gDF(idf,1), 1 ) ;
|
|
j = gElement( gDF(idf,1), 2 ) ;
|
|
m = gElement( gDF(idf,1), 3 ) ;
|
|
f( (i-1)*2+1 : (i-1)*2+2 ) = f( (i-1)*2+1 : (i-1)*2+2 ) + enf( 1:2 ) ;
|
|
f( (j-1)*2+1 : (j-1)*2+2 ) = f( (j-1)*2+1 : (j-1)*2+2 ) + enf( 3:4 ) ;
|
|
f( (m-1)*2+1 : (m-1)*2+2 ) = f( (m-1)*2+1 : (m-1)*2+2 ) + enf( 5:6 ) ;
|
|
end
|
|
|
|
% step 4. 处理约束条件，修改刚度矩阵和节点力向量。采用乘大数法
|
|
[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) * 1e20 ;
|
|
gK(m,m) = gK(m,m) * 1e20 ;
|
|
end
|
|
|
|
% step 5. 求解方程组，得到节点位移向量
|
|
gDelta = gK \ f ;
|
|
|
|
% step 6. 计算单元应力
|
|
gElementStress = zeros( element_number, 6) ;
|
|
for ie=1:element_number
|
|
es = ElementStress( ie ) ;
|
|
gElementStress( ie, : ) = es ;
|
|
end
|
|
|
|
% % step 7. 计算节点应力(采用绕节点加权平均)
|
|
gNodeStress = zeros( node_number, 6 ) ;
|
|
for i=1:node_number
|
|
S = zeros( 1, 3 ) ;
|
|
A = 0 ;
|
|
for ie=1:1:element_number
|
|
for k=1:1:3
|
|
if i == gElement( ie, k )
|
|
area= ElementArea( ie ) ;
|
|
S = S + gElementStress(ie,1:3 ) * area ;
|
|
A = A + area ;
|
|
break ;
|
|
end
|
|
end
|
|
end
|
|
gNodeStress(i,1:3) = S / A ;
|
|
gNodeStress(i,6) = 0.5*sqrt( (gNodeStress(i,1)-gNodeStress(i,2))^2 + 4*gNodeStress(i,3)^2 ) ;
|
|
gNodeStress(i,4) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)) + gNodeStress(i,6) ;
|
|
gNodeStress(i,5) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)) - gNodeStress(i,6) ;
|
|
end
|
|
return
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|