有限元法MATLAB程序设计

whutzll

重新学习结构力学中的有限元法
------我的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=[1 0 0 x 0 0];
Hv=[0 1 x 0 x^2 x^3];
B=[diff(Hu,x,1);-y*diff(Hv,x,2)]*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=[1 0 0 0 0 0 x 0 0 0 0 0];
Hv=[0 1 0 0 0 x 0 x^2 0 0 0 x^3];
Hw=[0 0 1 0 x 0 0 0 x^2 0 x^3 0];
Hq=[0 0 0 1 0 0 0 0 0 x 0 0];
H=[diff(Hu,x,1);-y*diff(Hv,x,2);-z*diff(Hw,x,2);r*diff(Hq,x,1)];
B=H*inv(A);
D=diag([E,E,E,G]);
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

molen

你用的什么语言？

whutzll

说明：
% 节点坐标: x      y
    gNode = [0.0,   0.0          % 节点 1
             1.0,   0.0          % 节点 2
             1.0,   1.0          % 节点 3
             0.0,   1.0];        % 节点 4
[node_umber,dummy]=size(gNode)   %返回矩阵gNode的行列数
% 单元定义: 节点1  节点2  材料号
    gElement = [1,      2,      1       % 单元 1
                2,      3,      1       % 单元 2
                3,      4,      1       % 单元 3
                4,      1,      1 ];    % 单元 4
[elemnet_umber,dummy]=size(gElement) %返回矩阵gElement的行列数
% 第一类约束条件
%     节点号   自由度号    约束值
gBC1 = [ 1,        1,        0.0
              1,        2,        0.0
              1,        3,        0.0];
[bc1_umber,dummy]=size(gBC1) %返回矩阵gBC1的行列数
% 第二类约束条件
%     节点号   自由度号    约束值
gBC2 = [ 2,        1,        0.0
              2,        2,        0.0];
[bc2_umber,dummy]=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 约束条件处理
% 采用划行划列法处理第一类约束条件，以避免出现刚度运动。
    [bc1_number,dummy] = 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
% 采用划行划列法处理第二类约束条件，实现自由度耦合。
    [bc2_number,dummy] = 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 编辑 ]

whutzll

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 编辑 ]

whutzll

Question. 等参数单元刚度矩阵求法(公式稍后附上)
function k = StiffnessMatrix
    k = zeros( 16, 16 ) ;
    D = MatrixD( 1 ) ;
    % 2 x 2  高斯积分点和权系数
    x = [-0.577350269189626, 0.577350269189626] ;                    
    w = [1, 1] ;                                                     
    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* [  1  A1   0
              A1   1   0
               0   0  A2] ;
return

function B = MatrixB( xi, eta ) %  计算单元的应变矩阵B
    [N_x,N_y] = 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_x, N_y] = N_xy( xi, eta ) %  计算形函数对整体坐标的导数
    J = Jacobi( xi, eta ) ;
    [N_xi,N_eta] = N_xieta( xi, eta ) ;
    A=inv(J)*[N_xi;N_eta] ;
    N_x = A(1,: ) ;
    N_y = A(2,: ) ;
return

function [N_xi, N_eta] = 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_xi,N_eta] = 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 编辑 ]

whutzll

连续梁数值分析matlab程序[beam]
〔将结构力学学习进行到底!
上传个调试成功例子〕

[ 本帖最后由 whutzll 于 2008-2-22 16:21 编辑 ]

whutzll

平面桁架数值分析程序

[ 本帖最后由 whutzll 于 2008-2-22 11:04 编辑 ]

whutzll

平面刚架数值分析程序(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 编辑 ]

whutzll

Ans:求解三角形平面有限元模型，过程如下
%        1. 计算单元刚度矩阵，集成整体刚度矩阵
%        2. 计算单元的等效节点力，集成整体节点力向量
%        3. 处理约束条件，修改整体刚度矩阵和节点力向量
%        4. 求解方程组，得到整体节点位移向量
%        5. 计算单元应力和节点应力
Q:多贴几个弹性平面单元刚度矩阵，为以后学习打下基础
%  计算三角形平面应变单元单元刚度矩阵[输入参数: ie-单元号;返回值:k-单元刚度矩阵]
syms  E       %弹性模量         
syms  mu      %伯松比
syms  h       %厚度
syms  xi      %第i点的x坐标
syms  yi      %第i点的y坐标
syms  xj      %第j点的x坐标
syms  yj      %第j点的y坐标
syms  xm      %第m点的x坐标
syms  ym      %第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 = [bi  0 bj  0 bm  0
          0 ci  0 cj  0 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）： [sx, sy, txy, s1, s2, tmax]
    es = zeros( 1, 6 ) ;   % 单元的应力分量
    de = zeros( 6, 1 ) ;   % 单元节点位移列阵
    % 定义整体刚度矩阵和节点力向量
    [node_number,dummy] = 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 编辑 ]

whutzll

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 = [x1; x2; x3] ; y = [y1; y2; y3] ;
    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 = [1/6 1/6 1/6] ;               %权系数
    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 = [x1; x2; x3] ; y = [y1; y2; y3] ;
    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*[Nii; Njj; Nij] ;
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 编辑 ]

whutzll

Q:地震结构动力响应程序[Newmark法]

% 定义全局变量
% 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. 定义整体刚度矩阵和节点力向量
    [node_number,dummy] = size( gNode ) ;
    gK = sparse( node_number * 3, node_number * 3 ) ;
    gM = sparse( node_number * 3, node_number * 3 ) ;

    % step2. 计算单元刚度和质量矩阵，并集成到整体刚度和质量矩阵中
    [element_number,dummy] = size( gElement ) ;
    for ie=1:1:element_number
        k = StiffnessMatrix( ie ) ;
        m = MassMatrix( ie ) ;
        AssembleGlobalMatrix( ie, k, m ) ;
    end

  % step3. 处理第一类约束条件，修改刚度矩阵和质量矩阵。（采用划行划列法）
    [bc1_number,dummy] = 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阶特征值和特征向量
    [gEigVector, gEigValue] = 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. 定义整体刚度矩阵和节点力向量
    [node_number,dummy] = 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. 计算单元刚度和质量矩阵，并集成到整体刚度和质量矩阵中
    [element_number,dummy] = 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 初始计算
    [acce_number,dummy] = size( gA) ;
    gDeltaT = 0.02 ;
    gTimeEnd = acce_number*gDeltaT  ;    % 计算时间为载荷通过所需时间的两倍
    % timestep =  floor(gTimeEnd/gDeltaT)  ;
    timestep =  acce_number;
   
    gama = 0.5 ;
    beta = 0.25 ;
    [N,N] = 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进行边界条件处理
    [bc1_number,dummy] = 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
    [KL,KU] = lu(K1) ;   % 进行三角分解，节省后面的求解时间
   
    % step5.6 计算地震荷载，初始加速度
    [node_number,dummy] = size( gNode ) ;
    [acce_number,dummy] = 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进行边界条件处理
        [bc1_number,dummy] = 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 编辑 ]

mglts

太谢谢了！！

168012

牛人啊你

shikun4852

谢了,拜读一下

等待的巧

谢谢，资料还真齐全

UltraKong

多谢了，就喜欢楼主这样无私的~

feixiangqin

谢谢，佩服呀，学习了。

houqujushi

这个程序比较常用，很好。

houqujushi

我正在找类似的资料呢。怎么才能看见呢？

xtjia713

LZ大牛啊，膜拜中