几何非线性大作业荷载增量法和弧长法程序设计.doc

几何非线性大作业荷载增量法
和弧长法程序设计

系（所）：建筑工程系
学 号：1432055
姓 名：焦 联 洪
培养层次：专业硕士
指导老师：吴 明 儿

2015年6月19日

一、 几何非线性大作业（ Newton-Raphson法）
用荷载增量法（Newton-Raphson法）编写几何非线性程序：
（1）用平面梁单元，可分析平面杆系
（2）算例：悬臂端作用弯矩。悬臂梁最终变形形成周长为悬臂梁长度的圆。

1.1 Newton-Raphson算法基本思想

图1.1 Newton-Raphson算法基本思想
1.2 悬臂梁参数
基本参数：L=2m, D=0.03m, A=7.069E-4m2, I=3.976E-08m4 ,E=2.0E11N/m2

图1.2 悬臂梁单元信息
将悬臂梁分成10个单元，如图1.2所示
2.1 MATLAB输入信息
材料信息 单元信息

约束信息（0为约束，1为放松） 荷载信息（FX,FY,M）

节点信息

2.2 求解过程
梁弯成圆形：理论弯矩M=EIY"=24981.944N.m ,直径为0.642m
运用ABAQUS和MATLAB进行求解对比：

图1.3 加载图

图1.4 ABAQUS变形图

图1.5 MATLAB变形曲线

ABAQUS和MATLAB变形对比，最终在理论荷载作用下都弯成了一个圆，其直径为0.64716m，与理论值相对比值为：（0.64716-0.642）/0.642=0.00804.非常接近。
2.3 加载点荷载位移曲线

图1.5 加载点Y方向的荷载位移曲线

加载点的最大竖向位移分别为1.4525m和1.45246m，相对比值（1.4525-1.45246）/1.45246=2.75395E-05。完全相同，说明MATLAB的计算结果很好。

二、 几何非线性大作业（ 弧长法）
用弧长法编写几何非线性程序，分析荷载位移全过程曲线：
1) 用平面梁单元，可分析平面杆系结构
2) 算例
(1)受集中荷载的拱：考察拱的矢跨比、荷载位置对荷载位移曲线的影响。
(2)其他有复杂平衡路径的结构
3) 将结果与相关文献进行对比

1.1 弧长法基本思想

图2.1 弧长法基本思想
1.2 拱基本参数
拱参数：L=100m, A=0.32m2 ,I=1m4 ,E=1.0e7N/m2，F=-5000N，拱曲线 y=5×sin(3.1415926*x/L)
将拱结构分成25个单元，如图2所示

图2.2 拱单元信息

2.1 MATLAB输入信息

材料信息 单元信息

约束信息（0为约束，1为放松） 荷载信息（FX,FY,M）

节点信息

2.2 运用ANSYS和MATLAB进行求解对比（两端铰接）
ANSYS中模型:

图2.3 ANSYS模型

图2.4 MATLAB和ANSYS变形图

2.3 加载点荷载位移曲线

图2.5 加载点荷载位移曲线

ANSYS求得的极限承载力3042.53,对应位移3.00142
MATLAB求得的极限承载力3043.8, 对应位移3.0768
相对误差分别为0.0417%,2.45%，模拟效果比较好。

2.4 拱的矢跨比a对拱荷载位移曲线的影响
不同矢跨比（1/20，3/40，1/10，3/20）下加载点的荷载位移曲线

1）MATLAB中计算拱的矢跨比a对拱荷载位移曲线的影响

图2.6 荷载位移曲线

图2.7 荷载位移曲线
表1 各矢跨比下拱结构的极限荷载
参数
矢高
极值点
F(N)
位移（m）
最低点
F(N)
位移（m）
5mm

3043.8
3.0768

1765.2
7.0816
7.5mm

7623.3
4.0335

-595.82
11.21
10mm

14974
5.4026

-6408.1
14.886
20mm

39791
9.4831

-63049
30.513

从表中可以初步得出：在一定随着矢跨比的增加，拱仍然呈现跳跃失稳的形式，拱结构的极限承载能力有大幅度的提高；在最低处的承载力呈现出反向，相当于有一个拉力在阻止拱结构发生跳跃失稳，矢跨比越大，拱越不容易发生跳跃失稳。当拱的矢跨比超过一定范围后，拱将发生复杂的不同于跳跃失稳的失稳形式。

2）MATLAB与ANSYS计算结果对比

图2.8 ANSYS和MATLAB对比荷载位移曲线

表2 各矢跨比下拱结构的极限荷载对比
参数
矢高
F(N)MAT
位移（m）
F(N)ANA
位移（m）
误差（%）
误差（%）
5mm
3043.8
30768
3042.53
3.00142
0.04
2.45
7.5mm
7623.3
4.0335
7624.91
3.96303
-0.02
1.75
10mm
14974
5.4026
14974.3
5.3157
0.00
1.61
20mm
39791
9.4831
39695.7
9.59955
0.24
-1.23

从图中可以看出：矢跨比在一定范围内，MATLAB与ANSYS计算的荷载位移曲线非常吻合，验证了MATLAB程序的可行性。当矢跨比为0.15时，ANSYS中将跟踪不到失稳后复杂的平衡路径。
从表中可以得出：MATLAB与ANSYS计算中拱的极限荷载和极限荷载时所对应的位移非常接近，加载点均为顶点26。具体为：矢高5mm，荷载误差为0.04，位移误差为2.45；矢高7.5mm，荷载误差为-0.02，位移误差为1.75；矢高10mm，荷载误差为0，位移误差为-1.61；矢高20mm，荷载误差为0.24，位移误差为-1.23。实际误差相差很小，在工程允许的范围内是可以接受的。
2.5 荷载位置对拱荷载位移曲线的影响

图2.9 ANSYS模型简图
1）MATLAB中计算荷载位置对拱荷载位移曲线的影响

图2.10 各加载点荷载位移曲线
表3 改变加载点拱结构的极限荷载
参数
加载点
极值点
F(N)
位移（m）
最低点
F(N)
位移（m）
26

3043.8
3.0768

1765.2
7.0816
16

3570
3.1891

2258.8
6.116
11

4728
2.88

3220.5
4.7959
4

14317
1.2826

10569
1.7829
误差=100*（MATLAB-ANSYS）/ANSYS
从表中可以初步得出：随着加载点位置越靠近支座处，拱结构的极限承载能力有大幅度的提高；在最低处的承载力也大幅度提高。当加载点位置靠近支座时，拱的承载力增加幅度最大，拱的稳定性很强，不容易发生失稳。

2）MATLAB与ANSYS计算结果对比

图2.11 ANSYS和MATLAB对比荷载位移曲线

表4 各加载点拱结构的极限荷载
参数
矢高
F(N)MAT
位移（m）
F(N)ANA
位移（m）
误差（%）
误差（%）
26
3043.8
3.0768
3042.53
3.00142
0.04
2.45
16
3570
3.1891
3569.73
3.24865
0.01
-1.87
11
4728
2.88
4728.71
2.91862
-0.02
-1.34
4
14317
1.2826
14324.8
1.29764
-0.05
-1.17
误差=100*（MATLAB-ANSYS）/ANSYS
从图中可以看出：MATLAB与ANSYS计算的荷载位移曲线非常吻合，验证了MATLAB程序的可行性。从表中可以得出：MATLAB与ANSYS计算中拱的极限荷载和极限荷载时所对应的位移非常接近。具体为：26加载点，荷载误差为0.04，位移误差为2.45；16加载点，荷载误差为0.01，位移误差为-1.87；11加载点，荷载误差为-0.02，位移误差为-1.34；4加载点，荷载误差为-0.05，位移误差为-1.17。实际误差相差很小，在工程允许的范围内是可以接受的。

2.6 两端铰接和固接对拱荷载位移曲线的影响
矢高为5mm时，拱两端为固接和铰接时的荷载位移曲线如下：

图2.12 ANSYS和MATLAB固接和铰接的荷载位移曲线

从图中可以看出：拱的边界条件对其的失稳形式有很大影响。两端固接拱的稳定性明显优于两端铰接拱，承载能力也大幅度提高。固接拱不会发生跳跃失稳的形式，刚度在初始阶段会减小，待到达一定程度后刚度又会增加。而两端铰接拱会发生跳跃失稳的形式。

2.7 参数m对拱失稳形式的影响
文献中给出:m是一个由拱截面在竖向平面内的回转半径r 和拱的初始矢高h无确定的无量纲参数。

当m>=1 时，扁拱不会出现跳跃屈曲, 当0=1 时，扁拱不会出现跳跃屈曲，当0受集中荷载。这个结构是研究分歧问题的经典题目。将半跨结构划分为8个单元, 得到图4b的基本路径和分歧路径, 并与JChreseielewski 和Rsehmiot的结果进行了比较。本文对此结构也进行了缺陷分析。
拱的基本参数：L=254cm，R=381cm，I=41.62cm4，A=6.45cm2，E=6898kN/cm2。

文献中的计算结果。

采用MATLAB和ANSYS对其进行求解，得到其荷载位移曲线：

图2.15 ABAQUS模型

图2.16 ABAQUS变形图

图2.17 ANSYS、MATLAB、ABAQUS加载点荷载位移曲线
从MATLAB、ANSYS、ABAQUS计算的荷载位移曲线可以看出：两者的荷载位移曲线基本吻合。MATLAB的计算结果可以看出在后期其承载能力会有较大提高，与文献中的荷载位移曲线趋势相同，所以验证出程序的可靠性。ABAQUS不能模拟出后续段曲线也许是单元划分过少。

图2.18 MATLAB加载点荷载位移曲线
第二次极值点会超过第一次极值点所对应的荷载，与文献一致，荷载点也接近。
加入初始缺陷：L/1000, L/2000初始缺陷后观察加载点的荷载位移曲线变化趋势。

图2.19 ANSYS加入初始缺陷后的加载点荷载位移曲线

2.20 初始缺陷为0.0001时的荷载位移曲线
加入初始缺陷后，拱的极限承载能力明显降低。且失稳形式也发生了变化，初始缺陷的大小对其荷载位移曲线有明显影响。所以在工程设计中应考虑结构或构件的初始缺陷，使结构的设计更加合理，安全。
为了研究初始缺陷对拱失稳路径的影响，应用ABAQUS和ANSYS以及MATLAB中加水平力模拟拱结构初始缺陷下的荷载位移曲线。
为了探究ABAQUS和ANSYS计算结果，取其前三阶模态进行对比分析：

2.21 一阶屈曲模态
ABAQUS和MATLAB中的一阶屈曲系数为0.55884和0.564512，对应的屈曲荷载为74325.72N 和75080.096N。

2.22 二阶屈曲模态

ABAQUS和MATLAB中的二阶屈曲系数为1.2259和1.253，对应的屈曲163044.7N 和
166649N。

2.23 三阶屈曲模态
ABAQUS和MATLAB中的三阶屈曲系数为2.166和2.255，对应的屈曲299915N 和
288078N。

从屈曲模态中可以看出，两种软件的前二阶模态趋势吻合，屈曲系数和极限荷载也是吻合的较好。第三阶模态出现不一样的变形趋势，屈曲系数和极限荷载也是也相差比较大，但计算时只引入一阶屈曲模态。

图2.24 ANSYS、ABAQUS、MATLAB加载点荷载位移曲线
从图中可以看出：ANSYS对缺陷的模拟效果比较好，与文献中的比较接近ABAQUS模拟的极限荷载稍低于ANSYS，而MATLAB模拟的极限荷载远低于ANSYS,但是曲线最终都在位移为300多mm时交于一点。还是有一定规律性。

图2.25 ANSYS和ABAQUS引入初始缺陷加载点荷载位移曲线
始缺陷并一定都会降低承载力，也会有对结构的承载能力有益的初始缺陷。ANSYS的计算结果可以看出，当初始缺陷为1/2000和-1/2000时，其承载能力不变。由于为对称结构，所以缺陷的正负影响不大。

图2.26 ANSYS和ABAQUS引入初始缺陷加载点荷载位移曲线

表6 对比ANNSYS和ABAQUS的极限荷载值和其对应位移值
参数
矢高
F(N)ANS
位移（m）
F(N)ABA
位移（m）
误差（%）
误差（%）
1/1000
58444.7
68.9799
55795.628
79.9318
4.53261
-15.8769
1/2000
60579.9
70.1384
57924.958
65.4542
4.38255
6.67851
误差=100*（ABAQUS-ANSYS）/ABAQUS
表中可以得出：ABAQUS求解出的机线荷载小于ANSYS，单对应的位移大于ANSYS对应的值。这可能与单元划分的个数，求解精度有关，但在工程允许的范围内，还是可以接受的。
ABAQUS中迭代步的跳跃很快，位移增加速度很快，其路径不是很准确，可能是由于其单元划分比较少。
体会：
1）注意坐标系的转化和力、位移更新时所对应的状态（C1--C2）
2) 拱是否发生跳跃失稳与矢跨比、矢高与截面旋转半径有关。矢跨比太大不会发生跳跃失稳；m大于1时不能发生跳跃失稳。
3）有些拱在ANSYS中不能得出下降段，可见ANSYS中对拱的跨度矢高、截面参数的比值有一定要求。内部计算和程序中有一些差别。

附录
子程序一：刚度组装矩阵
function K_G=Assemble(K_Element,Element,N_Node)
K_G=zeros(N_Node*3,N_Node*3);
for i=1:2
n1=Element(1,i);
for j=1:2
n2=Element(1,j);
K_G(3*n1-2:3*n1,3*n2-2:3*n2)=K_Element(3*i-2:3*i,3*j-2:3*j);
end
end
end
子程序二：整体坐标刚度矩阵
function K=beam2d_stiffness530(E,A,I,L,cs,Ele_F1);
F=Ele_F1(1,4);
M1=Ele_F1(1,3);
M2=Ele_F1(1,6);
T=[cs(1,1),cs(1,2),0,0,0,0;
-cs(1,2),cs(1,1),0,0,0,0;
0,0,1,0,0,0;
0,0,0,cs(1,1),cs(1,2),0;
0,0,0,-cs(1,2),cs(1,1),0;
0,0,0,0,0,1];
KE=[ E*A/L,0,0,-E*A/L,0,0;
0,12*E*I/L^36*E*I/L^2,0,-12*E*I/L^3,6*E*I/L^2;
0,6*E*I/L^2, 4*E*I/L, 0, -6*E*I/L^2, 2*E*I/L;
-E*A/L,0,0,E*A/L,0,0;
0,-12*E*I/L^3,-6*E*I/L^2,0,12*E*I/L^3,-6*E*I/L^2;
0,6*E*I/L^2,2*E*I/L,0,-6*E*I/L^2,4*E*I/L];

KG=[ F/L,0,-M1/L,-F/L,0,-M2/L;
0,12*F*I/A/L^3+6*F/5/L, 6*F*I/A/L^2+F/10, 0,-(12*F*I/A/L^36*F/5/L), 6*F*I/A/L^2+F/10;
-M1/L, 6*F*I/A/L^2+F/10, 4*F*I/A/L+2*F*L/15, M1/L, -(6*F*I/A/L^2+F/10), 2*F*I/A/L-F*L/30;
-F/L,0, M1/L,F/L,0,M2/L;
0,-12*F*I/A/L^3-6*F/5/L,-6*F*I/A/L^2-F/10,0,12*F*I/A/L^3+6*F/5/L,-6*F*I/A/L^2-F/10;
-M2/L, 6*F*I/A/L^2+F/10, 2*F*I/A/L-F*L/30, M2/L, -(6*F*I/A/L^2+F/10) 4*F*I/A/L+2*F*L/15];
K=T'*(KE+KG)*T;
end
子程序三：局部坐标刚度矩阵
function K=beam2d_stiffness520(E,A,I,L,cs,Ele_F1);
F=Ele_F1(1,4);
M1=Ele_F1(1,3);
M2=Ele_F1(1,6);
KE=[ E*A/L,0,0,-E*A/L,0,0;
0,12*E*I/L^3,6*E*I/L^2,0,-12*E*I/L^3,6*E*I/L^2;
0,6*E*I/L^2, 4*E*I/L, 0, -6*E*I/L^2, 2*E*I/L;
-E*A/L,0,0,EA/L,0,0;
0,-12*E*I/L^3,-6*E*I/L^2,0,12*E*I/L^3,-6*E*I/L^2;
0,6*E*I/L^2,2*E*I/L,0,-6*E*I/L^2,4*E*I/L];

KG=[ F/L,0,-M1/L,-F/L,0,-M2/L;
0,12*F*I/A/L^3+6*F/5/L, 6*F*I/A/L^2+F/10, 0,-(12*F*I/A/L^3+6*F/5/L), 6*F*I/A/L^2+F/10;
-M1/L, 6*F*I/A/L^2+F/10, 4*F*I/A/L+2*F*L/15, M1/L, -(6*F*I/A/L^2+F/10), 2*F*I/A/L-F*L/30;
-F/L,0, M1/L,F/L,0,M2/L;
0,-12*F*I/A/L^3-6*F/5/L,-6*F*I/A/L^2-F/10,0,12*F*I/A/L^3+6*F/5/L,-6*F*I/A/L^2-F/10;
-M2/L, 6*F*I/A/L^2+F/10, 2*F*I/A/L-F*L/30, M2/L, -(6*F*I/A/L^2+F/10), 4*F*I/A/L+2*F*L/15];
K=(KE+KG);
End
主程序一：Newton-Raphson法
clc
clear
Node=xlsread('Data.xls','Node');
Element=xlsread('Data.xls','Element');
Boundary=xlsread('Data.xls','Boundary');
Section=xlsread('Data.xls','Section');
Force=xlsread('Data.xls','Force');
%读取相关数据
Allstep=1000; %迭代步数
N_Node=size(Node,1); %节点个数
N_Element=size(Element,1); %单元个数
N_Force=size(Force,1); %节点力个数
N_Boundary=size(Boundary,1); %约束节点个数
Displacement=zeros(N_Node,3); %位移向量（a0）
%初始位移及转角为0
All_Disp=zeros(N_Node,3); %初始位移和转角为零（迭代后的节点位移）
All_F=zeros(N_Node*3,1); %初始荷载向量为零 （存放节点荷载向量）
Internal_F=zeros(N_Node*3,1); %节点内力向量

ForceMatrix=zeros(N_Node*3,1); %总荷载向量
C=zeros(N_Element,2);
L=zeros(N_Element,1);
for i=1:N_Force
ForceMatrix(Force(i,1)*3-2:Force(i,1)*3,1)=Force(i,2:4)'; %把节点荷载向量读入一个矩阵中,形成列向量=总的自由度个数
end
del=[];
for i=1:N_Boundary;
if Boundary(i,2)==0;
m=3*Boundary(i,1)-2;
del=[del,[m]]; %受约束节点位移为0时所对应的指标数值1
end
if Boundary(i,3)==0;
m=3*Boundary(i,1)-1;
del=[del,[m]]; %受约束节点位移为0时所对应的指标数值2
end
if Boundary(i,4)==0;
m=3*Boundary(i,1);
del=[del,[m]]; %受约束节点位移为0时所对应的指标数值3
end
end
%求出约束节点的标号，便于刚度、荷载矩阵归0
Update_Node=Node+Displacement(:,1:2); %更新后的节点坐标向量（x，y坐标）
Ele_F=zeros(N_Element,6); %单元节点荷载向量
ELEDISP=zeros(N_Element,6); %单元节点位移向量
Ele_F1=zeros(N_Element,6); %单元节点荷载刚度矩阵中向量
Ele1_F=zeros(1,6);
ELEDISP1=zeros(1,6);
qq(1,1)=0;
pp(1,1)=0;
for n=0:Allstep-1
n=n+1
K_Global=zeros(N_Node*3,N_Node*3); %总刚矩阵
for i=1:N_Element
N1=Element(i,1); %i节点编号
N2=Element(i,2); %j节点编号
N_Section=Element(i,3); %截面的形状控制参数
C(i,:)=Update_Node(N2,:)-Update_Node(N1,:); %a0下坐标向量增量
L(i)=norm(C(i,:)); %变形后的长度
cs=C(i,:)/L(i); %杆件的cos和sin值
E=Section(N_Section,1);
ASection(N_Section,2);
I=Section(N_Section,3);
K_Element=beam2d_stiffness530(E,A,I,L(i),cs,Ele_F1(i,:));
K_Global=K_Global+Assemble(K_Element,Element(i,:),N_Node); %形成总刚K0
end%整体刚度矩阵
Delta_Force=ForceMatrix/Allstep;%初始荷载向量（Q0）
Equation=[K_Global,Delta_Force];
%方程组
Disp_Transefer=ones(N_Node*3,1); %建立调节位移矩阵的位移向量
Disp_Transefer(del,:)=0; %将约束节点的位移值定为0,其他的定位1
Equation(del,:)=[];%把方程中约束节点所对应的行和列去掉
Equation(:,del)=[];
%引入约束条件修改方程组
n1=size(Equation,2); % 方程组中列数
dis1=(Equation(:,1:n1-1))\Equation(:,n1); % 刚度矩阵求逆后乘以荷载向量, Equation(:,n1)荷载向量，得到节点的位移（da0）
%求解方程组
zz=1;
%识别约束
for i=1:N_Node*3;
if Disp_Transefer(i,1)==1;
Disp_Transefer(i,1)=dis1(zz,1); %总的位移向量
zz=zz+1;
end
end
for i=1:N_Node
Displacement(i,:)=Disp_Transefer(i*3-2:i*3,1); % 位移增量，形成n*3的位移向量（da0）
end
All_Disp=All_Disp+Displacement %位移向量更新得到a1（总的位移增量）
All_F=All_F+Delta_Force; %外荷载向量更新Q1
Internal_F1=zeros(N_Node*3,1); %节点内力向量
Update_Node1=Update_Node; %C1状态
Update_Node=Node+All_Disp(:,1:2); %C2状态坐标位置更新（改）（迭代后的位置
for i=1:N_Element
F_Global=zeros(N_Node*3,1); %全局的荷载向量
for j=1:2
ELEDISP1(j*3-2:j*3)=Displacement(Element(i,j),:); %整体 取出当前单元节点位移向量
end
N1=Element(i,1); %i节点编号
N2=Element(i,2); %j节点编号
C1(i,:)=Update_Node1(N2,:)-Update_Node1(N1,:); %a0下坐标向量增量
L1(i)=norm(C1(i,:)); %变形后的长度
cs1=C1(i,:)/L1(i); %杆件的cos和sin值
T1=[cs1(1,1),cs1(1,2),0,0,0,0;
-cs1(1,2),cs1(1,1),0,0,0,0;
0,0,1,0,0,0;
0,0,0,cs1(1,1),cs1(1,2),0;
0,0,0,-cs1(1,2),cs1(1,1),0;
0,0,0,0,0,1];
ELEDISP(i,:)=T1* ELEDISP1(:);
X1=L1(i)+ ELEDISP(i,4)- ELEDISP(i,1);
Z1= ELEDISP(i,5)-ELEDISP(i,2);
LN=(X1^2+Z1^2)^0.5;
sin1=Z1/LN;
cos1=X1/LN;
Citaloca=atan2(sin1,cos1);
Ub=LN-L1(i);
THRA=ELEDISP(i,3)-Citaloca;
THRB=ELEDISP(i,6)-Citaloca;
ELEDISP(i,:)=[0,0,THRA,Ub,0,THRB];
K_Element=beam2d_stiffness520(E,A,I,L1(i),cs1,EleF1(i,:)); %L（i）为a0对应下的
Ele_F(i,:)=K_Element*ELEDISP(i,:)'; %局部坐标系下荷载（Q0）作用下的节点力
Ele_F1(i,:)= Ele_F1(i,:)+ Ele_F(i,:);
C2(i,:)=Update_Node(N2,:)-Update_Node(N1,:); %a0下坐标向量增量
L2(i)=norm(C2(i,:)); %变形后的长度
cs2=C2(i,:)/L2(i); 杆件的cos和sin值
T2=[cs2(1,1),cs2(1,2),0,0,0,0;
-cs2(1,2),cs2(1,1),0,0,0,0;
0,0,1,0,0,0;
0,0,0,cs2(1,1),cs2(1,2),0;
0,0,0,-cs2(1,2),cs2(1,1),0;
0,0,0,0,0,1];
Ele1_F(:)=T2'*Ele_F1(i,:)';行向量
N1=Element(i,1);
N2=Element(i,2);
F_Global(3*N1-2:3*N1,1)=Ele1_F(1:3); %i节点荷载
F_Global(3*N2-2:3*N2,1)=Ele1_F(4:6); %j节点荷载
Internal_F1=Internal_F1+F_Global; %a1对应下的P1
end
Val=Internal_F1-All_F %求出上次迭代完成时的残余应力Q1-P1
Correc_dis=zeros(N_Node,3); %构造新的节点位移向量，每次更新
Correc_dis1=zeros(N_Node,3); %构造新的节点位移向量，叠加位移增量
N_dis=size(dis1,1); %未受约束的节点位移数目,不为零的节点位移数目
dis3=zeros(N_dis,1); %构造一个向量
k=0;
%修改位移矩阵形式
while norm(Val)>1e-7 & k<500 %在一个荷载增量步下进行的对此迭代
k=k+1;
K_Global=zeros(N_Node*3,N_Node*3);
for i=1:N_Element
N1=Element(i,1);
N2=Element(i,2);
N_Section=Element(i,3);
C(i,:)=Update_Node(N2,:)-Update_Node(N1,:); %a1下坐标向量增量
L(i)=normC(i,:)); %变形后的长度
cs=C(i,:)/L(i);
E=Section(N_Section,1);
A=Section(N_Section,2);
I=Section(N_Section,3);
K_Element=beam2d_stiffness530(E,A,I,L(i),cs,Ele_F1(i,:));
K_Global=K_Global+Assemble(K_Element,Element(i,:),N_Node); %形成总刚k1
end
Equation=[K_Global,Val]; %在残余应力下的位移求解
Disp_Transefer=ones(N_Node*3,1);
Disp_Transefer(del,:)=0;
Equation(del,:)=[];
Equation(:,del)=[];
n1=size(Equation,2);
Dis2=-(Equation(:,1:n1-1))\Equation(:,n1) %荷位移增量da1