求教非线性有限元增量迭代matlab编程问题

卢武铉

最近做非线性作业，编写了update空间梁的matlab程序，最后的求解出了问题，有哪位高手编写过增量迭代求解的，给指导指导，看看我程序哪里错了，所有程序等作业提交了一同分享给大家
%三节点梁单元主程序
clear;
clc;
%读取节点信息
%[x y z vs1 vs2 vs3 b vt1 vt2 vt3 a]
%t为高度方向，s为宽度方向，r为长度方向
point = load('node.txt');
%读取单元信息
%[1 2 3]
global element;
element=load('element.txt');
%计算节点数量
number_point=length(point(:,1));
%计算单元数量
number_element=length(element(:,1));
%材料参数
E=1.2e9;
mu=0;%泊松比
D=D(E,mu);
%初始化应力矩阵
stress_new=zeros(12,3,number_element);
%施加载荷大小
step=0.1:0.1:1;
for step_N=1:length(step)
    R=zeros(6*number_point,1);
    R(50,1)=100;
    F=zeros(6*number_point,1);
    error=R-F;
    %循环计算得到结果
    node_new=point;
    i=1;
    while i<=10
        K=zeros(6*number_point,6*number_point);
        node_old=node_new;
        stress_old=stress_new;
        %形成整体刚度矩阵
        for element_N=1:number_element  %element_N为单元号
            %形成单元刚度矩阵
            [element_K element_F] = el_K( element_N ,stress_old( :,:,element_N),D,node_old);
            for k1=1:3
                node_N1=element(element_N,k1);  %node_N为节点号
                for k2=1:3
                    node_N2=element(element_N,k2);  %node_N为节点号
                    K(6*node_N1-5:6*node_N1,6*node_N2-5:6*node_N2)=...
                    K(6*node_N1-5:6*node_N1,6*node_N2-5:6*node_N2)+element_K(6*k1-5:6*k1,6*k2-5:6*k2);
                end
                F(6*node_N1-5:6*node_N1)=element_F(6*k1-5:6*k1);
            end
        end
        %施加边界条件
        K(1:6,=0;
        K(:,1:6)=0;
        K(1:6,1:6)=eye(6,6);
        %求解节点位移矩阵
        displacement=K\(R-F);
        %更新坐标，以及更新应力矩阵
        for element_N=1:length(element(:,1))
            gauss = gauss_integral();
            for gauss_N=1:length(gauss(:,1))
                [J,detJ]= jacobian(element_N,gauss(gauss_N,1:3),node_old);
                [ BL BNL BL1 B] =element_B( element_N,gauss(gauss_N,:),stress_old(gauss_N,:,element_N),J,D,node_old);
                element_displacement=zeros(18,1);
                for k=1:3
                    node_N=element(element_N,k);
                    element_displacement(6*k-5:6*k)=displacement(6*node_N-5:6*node_N);
                end
                stress_new(gauss_N,:,element_N)=stress_old(gauss_N,:,element_N)+(D*B*element_displacement)';
            end
        end
        for node_N=1:length(node_old(:,1))
            node_new(node_N,1:3)=node_old(node_N,1:3)+(displacement(6*node_N-5:6*node_N-3))';
            node_new(node_N,4:6)=node_old(node_N,4:6)+cross(displacement(6*node_N-2:6*node_N),node_old(node_N,4:6));
            node_new(node_N,8:10)=node_old(node_N,8:10)+cross(displacement(6*node_N-2:6*node_N),node_old(node_N,8:10));
        end
        %判断循环条件，设置收敛条件
        if norm(R-F)<=0.1*norm(error)
            break;
        end
        i=i+1;
    end
end
plot(node_new(:,1),node_new(:,2));