有限元法MATLAB程序设计
本帖最后由 whutzll 于 2010-2-7 10:38 编辑重新学习结构力学中的有限元法
------我的qq号码：25010849鲨客。欢迎志同道合的朋友！
Question.1 平面梁单元的单元刚度矩阵求法：（MATLAB语言）
syms E L x y
A=[ 1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
1 0 0 L 0 0
0 1 0 0 L^2 L^3
0 0 1 0 2*L 3*L^2];
Hu=;
Hv=;
B=*inv(A);
Ke=int(E*transpose(B)*B,x,0,L)
运行结果：
Ke =
[ E/L, 0, 0, -E/L, 0, 0]
[ 0,12*E*y^2/L^3,-6*E*y^2/L^2, 0, -12*E*y^2/L^3, 6*E*y^2/L^2]
[ 0,-6*E*y^2/L^2, 4*E*y^2/L, 0, 6*E*y^2/L^2, -4*E*y^2/L]
[ -E/L, 0, 0, E/L, 0, 0]
[ 0, -12*E*y^2/L^3, 6*E*y^2/L^2, 0,12*E*y^2/L^3,-6*E*y^2/L^2]
[ 0, 6*E*y^2/L^2, -4*E*y^2/L, 0,-6*E*y^2/L^2, 4*E*y^2/L]
请教问题：运行结果与结构力学书中公式不一致，（出入处见红色标记），百思不得其解。
************************************对比分析**********************************************************************
空间梁单元的单元刚度矩阵，（运行结果一致）
syms E G L x y z r
A=[ 1 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 0
1 0 0 0 0 0 L 0 0 0 0 0
0 1 0 0 0 L 0 L^2 0 0 0 L^3
0 0 1 0 L 0 0 0 L^2 0 L^3 0
0 0 0 1 0 0 0 0 0 L 0 0
0 0 0 0 1 0 0 0 2*L 0 3*L^2 0
0 0 0 0 0 1 0 2*L 0 0 0 3*L^2];
Hu=;
Hv=;
Hw=;
Hq=;
H=;
B=H*inv(A);
D=diag();
Ke=int(transpose(B)*D*B,x,0,L)
-----------------------------------------------------------
Question 2 计算平面梁单元的坐标转换矩阵(局部坐标 -> 整体坐标)
syms xi xj yi yj
L = sqrt( (xj-xi)^2 + (yj-yi)^2 ) ;
c = (xj-xi)/L ;
s = (yj-yi)/L ;
T=[ c-s 0 0 0 0
s c 0 0 0 0
0 0 1 0 0 0
0 0 0 c-s 0
0 0 0 s c 0
0 0 0 0 0 1] ;
Ke = T*Ke*transpose(T) ;
请教问题：在T矩阵中c与s的位置是否正确？与书中仍表示有出入
参考书目：徐荣桥《结构分析的有限元法与MATLAB程序设计》，人民交通出版社，2006.5 你用的什么语言？
有限元法MATLAB程序设计
说明：% 节点坐标: x y
gNode = [0.0, 0.0 % 节点 1
1.0, 0.0 % 节点 2
1.0, 1.0 % 节点 3
0.0, 1.0]; % 节点 4
=size(gNode) %返回矩阵gNode的行列数
% 单元定义: 节点1节点2材料号
gElement = [1, 2, 1 % 单元 1
2, 3, 1 % 单元 2
3, 4, 1 % 单元 3
4, 1, 1 ]; % 单元 4
=size(gElement) %返回矩阵gElement的行列数
% 第一类约束条件
% 节点号 自由度号 约束值
gBC1 = [ 1, 1, 0.0
1, 2, 0.0
1, 3, 0.0];
=size(gBC1) %返回矩阵gBC1的行列数
% 第二类约束条件
% 节点号 自由度号 约束值
gBC2 = [ 2, 1, 0.0
2, 2, 0.0];
=size(gBC2) %返回矩阵gBC2的行列数
% 定义整体刚度矩阵和节点力向量（稀疏矩阵表达）
gK = sparse( node_number * 3, node_number * 3 ) ; % 定义整体刚度矩阵
f = sparse( node_number * 3, 1 ) ; % 定义结点荷载列阵
%结点号为i的第p个自由度在整体刚度矩阵中的行号是m=（i-1）*3+p
-------------------------------------------------------------------------------------
Question 3 把单元刚度矩阵集成到整体刚度矩阵;输入参数:ie-单元号、k-单元刚度矩阵
（依据自由度转换公式或称“对号叠加”），怎样体会“根据结点号和自由度号来得到在整体刚度矩阵中的行号和列号”
function AssembleStiffnessMatrix( ie, k )
global gElement gK %全局变量
for i=1:1:2
for j=1:1:2
for p=1:1:3
for q =1:1:3
m = (i-1)*3+p ;
n = (j-1)*3+q ;
M = (gElement(ie,i)-1)*3+p ;
N = (gElement(ie,j)-1)*3+q ;
gK(M,N) = gK(M,N) + k(m,n) ;
end
end
end
end
return
-----------------------------------------------------------------------------------
Question 4 约束条件处理
% 采用划行划列法处理第一类约束条件，以避免出现刚度运动。
= size( gBC1 ) ;
for ibc=1:1:bc1_number
n = gBC1(ibc, 1 ) ;%提取单元ibc的结点号
d = gBC1(ibc, 2 ) ;%提取单元ibc的自由度号
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
% 采用划行划列法处理第二类约束条件，实现自由度耦合。
= 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
[ 本帖最后由 whutzll 于 2008-2-23 13:38 编辑 ] Question 5
平面梁单元质量矩阵求法(在结构动力学问题中多用)
syms L A r x
n1=1-x/L;
n2=x/L;
n3=1-3*x^2/L^2+2*x^3/L^3;
n4=x-2*x^2/L+x^3/L^2;
n5=3*x^2/L^2-2*x^3/L^3;
n6=-x^2/L+x^3/L^2;
N=[n1 0 0 n2 0 0;
0 n3 n4 0 n5 n6];
m=int(r*A*transpose(N)*N,x,0,L)
求解结果: [与书中结果一致]
m =
[ 1/3*r*A*L, 0, 0, 1/6*r*A*L, 0, 0]
[ 0, 13/35*r*A*L,11/210*r*A*L^2, 0, 9/70*r*A*L, -13/420*r*A*L^2]
[ 0,11/210*r*A*L^2, 1/105*r*A*L^3, 0,13/420*r*A*L^2,-1/140*r*A*L^3]
[ 1/6*r*A*L, 0, 0, 1/3*r*A*L, 0, 0]
[ 0, 9/70*r*A*L,13/420*r*A*L^2, 0, 13/35*r*A*L, -11/210*r*A*L^2]
[ 0, -13/420*r*A*L^2,-1/140*r*A*L^3, 0, -11/210*r*A*L^2, 1/105*r*A*L^3]
[ 本帖最后由 whutzll 于 2008-2-5 12:04 编辑 ] Question. 等参数单元刚度矩阵求法(公式稍后附上)
function k = StiffnessMatrix
k = zeros( 16, 16 ) ;
D = MatrixD( 1 ) ;
% 2 x 2高斯积分点和权系数
x = [-0.577350269189626, 0.577350269189626] ;
w = ;
for i=1:1:length(x)
for j=1:1:length(x)
B = MatrixB( x(i), x(j) ) ;
J = Jacobi( x(i), x(j) ) ;
k = k + w(i)*w(j)*transpose(B)*D*B*det(J)
end
end
return
function D = MatrixD( opt ) %计算单元的弹性矩阵D
global E mu
if opt == 1 % 平面应力的弹性常数
A1 = mu ;
A2 = (1-mu)/2 ;
A3= E/(1-mu^2) ;
else % 平面应变的弹性常数
A1 = mu/(1-mu) ;
A2 = (1-2*mu)/2/(1-mu) ;
A3 = E*(1-mu)/(1+mu)/(1-2*mu) ;
end
D = A3* [1A1 0
A1 1 0
0 0A2] ;
return
function B = MatrixB( xi, eta ) %计算单元的应变矩阵B
= N_xy( xi, eta );
B = zeros( 3, 16 ) ;
for i=1:1:8
B(1:3,(2*i-1):2*i) = [ N_x(i) 0
0 N_y(i)
N_y(i),N_x(i)];
end
return
function = N_xy( xi, eta ) %计算形函数对整体坐标的导数
J = Jacobi( xi, eta ) ;
= N_xieta( xi, eta ) ;
A=inv(J)* ;
N_x = A(1,: ) ;
N_y = A(2,: ) ;
return
function = N_xieta( xi, eta ) %计算形函数对局部坐标的导数
x = [ -1, 1, 1, -1 ] ;
e = [ -1, -1, 1, 1 ] ;
N_xi= zeros( 1, 8 ) ;
N_eta = zeros( 1, 8 ) ;
N_xi( 5 )= xi*(eta-1) ;
N_eta( 5 ) = 0.5*(xi^2-1) ;
N_xi( 6 )= 0.5*(1-eta^2) ;
N_eta( 6 ) = -eta*(xi+1) ;
N_xi( 7 )= -xi*(eta+1) ;
N_eta( 7 ) = 0.5*(1-xi^2) ;
N_xi( 8 )= 0.5*(eta^2-1) ;
N_eta( 8 ) = eta*(xi-1) ;
N_xi(1)= x(1)*(1+e(1)*eta)/4 - 0.5*( N_xi(5)+ N_xi(8) );
N_eta(1) = e(1)*(1+x(1)*xi)/4- 0.5*( N_eta(5) + N_eta(8) ) ;
N_xi(2)= x(2)*(1+e(2)*eta)/4 - 0.5*( N_xi(5)+ N_xi(6) );
N_eta(2) = e(2)*(1+x(2)*xi)/4- 0.5*( N_eta(5) + N_eta(6) ) ;
N_xi(3)= x(3)*(1+e(3)*eta)/4 - 0.5*( N_xi(6)+ N_xi(7) );
N_eta(3) = e(3)*(1+x(3)*xi)/4- 0.5*( N_eta(6) + N_eta(7) ) ;
N_xi(4)= x(4)*(1+e(4)*eta)/4 - 0.5*( N_xi(7)+ N_xi(8) );
N_eta(4) = e(4)*(1+x(4)*xi)/4- 0.5*( N_eta(7) + N_eta(8) ) ;
return
function J = Jacobi( xi, eta )%计算雅克比矩阵
global x y
= N_xieta( xi, eta ) ;
x_xi= N_xi* x ;
x_eta = N_eta * x ;
y_xi= N_xi* y ;
y_eta = N_eta * y ;
J = [ x_xi, y_xi; x_eta, y_eta ];
return
[ 本帖最后由 whutzll 于 2008-2-10 15:28 编辑 ] 连续梁数值分析matlab程序
〔将结构力学学习进行到底!
上传个调试成功例子〕
[ 本帖最后由 whutzll 于 2008-2-22 16:21 编辑 ] 平面桁架数值分析程序
[ 本帖最后由 whutzll 于 2008-2-22 11:04 编辑 ] 平面刚架数值分析程序(Matlab语言)【建议比较三楼提供的资料Exam1.rar】
内容含
1平面刚架杆件坐标定位和转轴矩阵
2.1局部坐标系下平面刚架杆端位移与杆端力
2.2整体坐标系下平面刚架杆端位移与杆端力
2.3平面刚架的转轴变换矩阵
3.1局部坐标系下平面刚架的单元刚度矩阵
3.2整体坐标系下平面刚架的单元刚度矩阵
4平面刚架的总节点位移和总支座反力
5平面刚架的总节点刚度矩阵
6平面刚架的综合节点荷载
7.1杆件上的集中荷载引起的约束杆端力
7.2杆件上的分布力荷载引起的约束杆端力
8平面刚架的最终杆端力
[ 本帖最后由 whutzll 于 2008-2-23 13:48 编辑 ] Ans:求解三角形平面有限元模型，过程如下
% 1. 计算单元刚度矩阵，集成整体刚度矩阵
% 2. 计算单元的等效节点力，集成整体节点力向量
% 3. 处理约束条件，修改整体刚度矩阵和节点力向量
% 4. 求解方程组，得到整体节点位移向量
% 5. 计算单元应力和节点应力
Q:多贴几个弹性平面单元刚度矩阵，为以后学习打下基础
%计算三角形平面应变单元单元刚度矩阵[输入参数: ie-单元号;返回值:k-单元刚度矩阵]
symsE %弹性模量
symsmu %伯松比
symsh %厚度
symsxi %第i点的x坐标
symsyi %第i点的y坐标
symsxj %第j点的x坐标
symsyj %第j点的y坐标
symsxm %第m点的x坐标
symsym %第m点的y坐标
k = zeros( 6, 6 ) ; %定义初始矩阵
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) ;
area = (ai+aj+am)/2 ; %计算面积
B = [bi0 bj0 bm0
0 ci0 cj0 cm
ci bi cj bj cm bm] ;
B = B/2/area ; %计算单元的应变矩阵B
D = [ 1-mu mu 0
mu 1-mu 0
0 0 (1-2*mu)/2] ;
D = D*E/(1-2*mu)/(1+mu) ;
k = transpose(B)*D*B*h*abs(area) %单元刚度计算
--------------------------------------------------------------------------------
Q: 计算单元应力分量ElementStress
%计算es-单元应力分量列阵（1×6）：
es = zeros( 1, 6 ) ; % 单元的应力分量
de = zeros( 6, 1 ) ; % 单元节点位移列阵
% 定义整体刚度矩阵和节点力向量
= size( gNode ) ;
gK = sparse( node_number * 2, node_number * 2 ) ;
f = sparse( node_number * 2, 1 ) ;
gDelta = gK \ f ;% 点位移向量
for j=1:1:3
de( 2*j-1 ) = gDelta( 2*gElement( ie, j )-1 ) ;
de( 2*j ) = gDelta( 2*gElement( ie, j ) ) ;
end
es(1:3) = D * B * 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) ;
[ 本帖最后由 whutzll 于 2008-2-25 14:40 编辑 ] Q：三角形薄板的刚度矩阵求解
function k = StiffnessMatrix( ie )
%计算薄板三角形单元的刚度矩阵
syms ie %ie-单元号
syms E u h %弹性摸量，泊松比，厚度
syms x1 x2 x3 y1 y2 y3 %三角形三点x,y坐标
x = sparse(3,1); y = sparse(3,1) ;%定义坐标的稀疏矩阵
x = ; y = ;
i=1; j=2; m=3;
ai = x(j)*y(m) - x(m)*y(j) ;
aj = x(m)*y(i) - x(i)*y(m) ;
am = x(i)*y(j) - x(j)*y(i) ;
delta2 = ai+aj+am ;
k = zeros( 9, 9 ) ;
L = [1/2 1/2 0;
0 1/2 1/2;
1/2 0 1/2] ;
W = ; %权系数
D = MatrixD( ie ) ; %弹性矩阵
for i=1:length(W)
B = MatrixB( ie, L( i, : ) ) ;%取应变矩阵中形函数第i个L
k = k + W(i)*transpose(B)*D*B*delta2 ;
end
return
%计算薄板三角形单元的应变矩阵
function B = MatrixB( ie, L )
syms ie %ie-单元号
syms E u h
syms x1 x2 x3 y1 y2 y3
x = sparse(3,1) ;
y = sparse(3,1) ;
x = ; y = ;
i=1; j=2; m=3;
ai = x(j)*y(m) - x(m)*y(j) ;
aj = x(m)*y(i) - x(i)*y(m) ;
am = x(i)*y(j) - x(j)*y(i) ;
bi = y(j) - y(m) ;
bj = y(m) - y(i) ;
bm = y(i) - y(j) ;
ci = -(x(j)-x(m)) ;
cj = -(x(m)-x(i)) ;
cm = -(x(i)-x(j)) ;
delta2 = ai+aj+am ;
T = [ bi^2 bj^2 2*bi*bj
ci^2 cj^2 2*ci*cj
2*bi*ci 2*bj*cj 2*(bi*cj+bj*ci) ] ;
Li = L(1); Lj = L(2); Lm = L(3) ;
Nii = [ -4*Lj-12*Li+6, ...
-6*bj*Li+2*bj-3*bj*Lj-bm*Lj, ...
2*cj-6*cj*Li-3*cj*Lj-cm*Lj, ...
-4*Lj, ...
(-bm+bi)*Lj, ...
-(cm-ci)*Lj, ...
-6+12*Li+8*Lj, ...
bi*Lj-6*bj*Li+4*bj-3*bj*Lj, ...
ci*Lj+4*cj-6*cj*Li-3*cj*Lj ] ; %形函数
Njj = [ -4*Li, ...
-(bj-bm)*Li, ...
(-cj+cm)*Li, ...
6-4*Li-12*Lj, ...
bm*Li+6*bi*Lj-2*bi+3*bi*Li, ...
cm*Li+6*ci*Lj-2*ci+3*ci*Li, ...
8*Li-6+12*Lj, ...
6*bi*Lj-4*bi+3*bi*Li-bj*Li, ...
6*ci*Lj-4*ci+3*ci*Li-cj*Li] ;
Nij = [ -4*Li-4*Lj+2, ...
-3*bj*Li-bm*Li+1/2*bj-bj*Lj-1/2*bm+bm*Lj, ...
-3*cj*Li-cm*Li+1/2*cj-cj*Lj-1/2*cm+cm*Lj, ...
-4*Li-4*Lj+2, ...
bm*Lj+3*bi*Lj+1/2*bm-bm*Li-1/2*bi+bi*Li, ...
cm*Lj+3*ci*Lj+1/2*cm-cm*Li-1/2*ci+ci*Li, ...
-4+8*Li+8*Lj, ...
3*bi*Lj-3/2*bi+bi*Li-3*bj*Li+3/2*bj-bj*Lj, ...
3*ci*Lj-3/2*ci+ci*Li-3*cj*Li+3/2*cj-cj*Lj] ;
B = -1/delta2^2*T* ;
return
%计算薄板三角形单元的弹性矩阵
function D = MatrixD( ie )
syms ie %ie-单元号
syms E u h
D=sparse(3,3);
E = E ; %弹性模量
u = u ; %泊松比
h = h ; %厚度
D = E*h^3/12/(1-u^2)*[1 , u, 0;...
u , 1, 0;...
0 , 0, (1-u)/2; ] ; %单元弹性矩阵
return
[ 本帖最后由 whutzll 于 2008-2-25 11:03 编辑 ] Q:地震结构动力响应程序
% 定义全局变量
% gNode:节点坐标; gElemen:t单元定义; gMaterial:材料性质;gBC1:第一类约束条件;
% gK:整体刚度矩阵;gDeltaT:时间步长;gTimeEnd:计算结束时刻; gDisp:位移时程响应;
% gVelo:速度时程响应;gAcce:加速度时程响应;
% 该函数求解有限元模型，过程如下
% 1. 计算单元的刚度和质量矩阵，集成整体刚度和质量矩阵
% 2. 求解特征值问题
% 3. 用Newmark法计算时程响应
function SolveModel
global gNode gElement gMaterial gBC1 gBC2 gK gM gDeltaT gTimeEnd gDisp gVelo ...
gA gAcce gLoad gEigValue gEigVector C
% step1. 定义整体刚度矩阵和节点力向量
= size( gNode ) ;
gK = sparse( node_number * 3, node_number * 3 ) ;
gM = sparse( node_number * 3, node_number * 3 ) ;
% step2. 计算单元刚度和质量矩阵，并集成到整体刚度和质量矩阵中
= size( gElement ) ;
for ie=1:1:element_number
k = StiffnessMatrix( ie ) ;
m = MassMatrix( ie ) ;
AssembleGlobalMatrix( ie, k, m ) ;
end
% step3. 处理第一类约束条件，修改刚度矩阵和质量矩阵。（采用划行划列法）
= size( gBC1 ) ;
w1max = max( diag(gK)./diag(gM) ) ;
for ibc=1:1:bc1_number
n = gBC1(ibc, 1 ) ;
d = gBC1(ibc, 2 ) ;
m = (n-1)*3 + d ;
gK(:,m ) = zeros( node_number*3, 1 ) ;
gK(m,: ) = zeros( 1, node_number*3 ) ;
gK(m,m ) = 1;
gM(:,m ) = zeros( node_number*3, 1 ) ;
gM(m,: ) = zeros( 1, node_number*3 ) ;
gM(m,m) = gK(m,m)/w1max/1e10 ;
end
% step4. 求解特征值问题
% step4.1为了使刚度矩阵和质量矩阵对称（在计算时可能引入舍入误差）
for i=1:node_number*3
for j=i:node_number*3
gK(j,i) = gK(i,j) ;
gM(j,i) = gM(i,j) ;
end
end
% step4.2 计算前6阶特征值和特征向量
= eigs(gK, gM, 3, 'SM' )
% step4.3 修改特征向量中受约束的自由度
for ibc=1:1:bc1_number
n = gBC1(ibc, 1 ) ;
d = gBC1(ibc, 2 ) ;
m = (n-1)*3 + d ;
gEigVector(m,:) = gBC1(ibc,3) ;
end
% step5. 计算时程响应(Newmark法)
% step5.1. 定义整体刚度矩阵和节点力向量
= size( gNode ) ;
gK = sparse( node_number * 3, node_number * 3 ) ;
gM = sparse( node_number * 3, node_number * 3 ) ;
C= sparse( node_number * 3, node_number * 3 ) ;
% step5.2. 计算单元刚度和质量矩阵，并集成到整体刚度和质量矩阵中
= size( gElement ) ;
for ie=1:1:element_number
k = StiffnessMatrix( ie ) ;
m = MassMatrix( ie ) ;
AssembleGlobalMatrix( ie, k, m ) ;
end
% step5.3 计算阻尼矩阵
fre_number = length(diag(gEigValue)) ;
for i=fre_number:-1:1
freq(fre_number-i+1)= sqrt(gEigValue(i,i))/2/pi;
end
alphad=2*freq(1)*freq(2)*0.05*(freq(2)-freq(1))/(freq(2)*freq(2)-freq(1)*freq(1));
betad=2*0.05*(freq(2)-freq(1))/(freq(2)*freq(2)-freq(1)*freq(1)) ;
C = alphad*gM+betad*gK;
% step5.4 初始计算
= size( gA) ;
gDeltaT = 0.02 ;
gTimeEnd = acce_number*gDeltaT; % 计算时间为载荷通过所需时间的两倍
% timestep =floor(gTimeEnd/gDeltaT);
timestep =acce_number;
gama = 0.5 ;
beta = 0.25 ;
= size( gK ) ;
alpha0 = 1/beta/gDeltaT^2 ;
alpha1 = gama/beta/gDeltaT ;
alpha2 = 1/beta/gDeltaT ;
alpha3 = 1/2/beta - 1 ;
alpha4 = gama/beta - 1 ;
alpha5 = gDeltaT/2*(gama/beta-2) ;
alpha6 = gDeltaT*(1-gama) ;
alpha7 = gama*gDeltaT ;
K1 = gK + alpha0*gM + alpha1*C;
% step5.5 对K1进行边界条件处理
= size( gBC1 ) ;
K1im = zeros(N,bc1_number) ;
for ibc=1:1:bc1_number
n = gBC1(ibc, 1 ) ;
d = gBC1(ibc, 2 ) ;
m = (n-1)*3 + d ;
K1im(:,ibc) = K1(:,m) ;
K1(:,m ) = zeros( node_number*3, 1 ) ;
K1(m,: ) = zeros( 1, node_number*3 ) ;
K1(m,m ) = 1.0 ;
end
= lu(K1) ; % 进行三角分解，节省后面的求解时间
% step5.6 计算地震荷载，初始加速度
= size( gNode ) ;
= size( gA) ;
gDisp = zeros( node_number*3, timestep ) ;
gVelo = zeros( node_number*3, timestep ) ;
gAcce = zeros( node_number*3, timestep ) ;
gA2 = zeros( node_number*3, timestep ) ;
gDisp(:,1 ) = zeros( node_number*3, 1 ) ;
gVelo(:,1 ) = zeros( node_number*3, 1 ) ;
gLoad = zeros( node_number*3, timestep ) ;
for j=1:1:acce_number
gA2(:,j) = gA(j,1) ;
end
gLoad = gM * gA2;
gAcce(:,1) = gM\(gLoad(:,1) -gK*gDisp(:,1)-C*gVelo(:,1)) ;
%step5.7 对每一个时间步计算
for i=2:1:timestep
if mod(i,100) == 0
fprintf( '当前时间步：%d\n', i ) ;
end
f1 = gLoad( :,i-1) + gM*(alpha0*gDisp(:,i-1)+alpha2*gVelo(:,i-1)+alpha3*gAcce
( :,i-1)) ...
+ C*(alpha1*gDisp(:,i-1)+alpha4*gVelo(:,i-1)+alpha5*gAcce(:,i-
1)) ;
% 对f1进行边界条件处理
= size( gBC1 ) ;
for ibc=1:1:bc1_number
n = gBC1(ibc, 1 ) ;
d = gBC1(ibc, 2 ) ;
m = (n-1)*3 + d ;
f1 = f1 - gBC1(ibc,3) * K1im(:,ibc) ;
f1(m) = gBC1(ibc,3) ;
end
y = KL\f1 ;
gDisp( :,i) = KU\y ;
gAcce( :,i) = alpha0*(gDisp(:,i)-gDisp(:,i-1)) - alpha2*gVelo(:,i-1) -
alpha3*gAcce( :,i-1) ;
gVelo( :,i) = gVelo(:,i-1) + alpha6*gAcce(:,i-1) + alpha7*gAcce(:,i) ;
end
return
[ 本帖最后由 whutzll 于 2008-2-28 15:31 编辑 ] 太谢谢了！！ 牛人啊你 谢了,拜读一下 谢谢，资料还真齐全 多谢了，就喜欢楼主这样无私的~ 谢谢，佩服呀，学习了。 这个程序比较常用，很好。 我正在找类似的资料呢。怎么才能看见呢？ LZ大牛啊，膜拜中
页:
[1]
2