但拱的最终变形图非常接近,只是此时拱的失稳形式变了。
2.8具有复杂失稳形式的拱
铰支圆拱
该结构及其几何参数、物理性质均示于图4a中。
中心受集中荷载。
这个结构是研究分歧问题的经典题目。
将半跨结构划分为8个单元,得到图4b的基本路径和分歧路径,并与JChreseielewski和Rsehmiot的结果进行了比较。
本文对此结构也进行了缺陷分析。
拱的基本参数:
L=254cm,R=381cm,I=41.62cm4,A=6.45cm2,E=6898kN/cm2。
文献中的计算结果。
采用MATLAB和ANSYS对其进行求解,得到其荷载位移曲线:
图2.15ABAQUS模型
图2.16ABAQUS变形图
图2.17ANSYS、MATLAB、ABAQUS加载点荷载位移曲线
从MATLAB、ANSYS、ABAQUS计算的荷载位移曲线可以看出:
两者的荷载位移曲线基本吻合。
MATLAB的计算结果可以看出在后期其承载能力会有较大提高,与文献中的荷载位移曲线趋势相同,所以验证出程序的可靠性。
ABAQUS不能模拟出后续段曲线也许是单元划分过少。
图2.18MATLAB加载点荷载位移曲线
第二次极值点会超过第一次极值点所对应的荷载,与文献一致,荷载点也接近。
加入初始缺陷:
L/1000,L/2000初始缺陷后观察加载点的荷载位移曲线变化趋势。
图2.19ANSYS加入初始缺陷后的加载点荷载位移曲线
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.24ANSYS、ABAQUS、MATLAB加载点荷载位移曲线
从图中可以看出:
ANSYS对缺陷的模拟效果比较好,与文献中的比较接近ABAQUS模拟的极限荷载稍低于ANSYS,而MATLAB模拟的极限荷载远低于ANSYS,但是曲线最终都在位移为300多mm时交于一点。
还是有一定规律性。
图2.25ANSYS和ABAQUS引入初始缺陷加载点荷载位移曲线
始缺陷并一定都会降低承载力,也会有对结构的承载能力有益的初始缺陷。
ANSYS的计算结果可以看出,当初始缺陷为1/2000和-1/2000时,其承载能力不变。
由于为对称结构,所以缺陷的正负影响不大。
图2.26ANSYS和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中对拱的跨度矢高、截面参数的比值有一定要求。
内部计算和程序中有一些差别。
附录
子程序一:
刚度组装矩阵
functionK_G=Assemble(K_Element,Element,N_Node)
K_G=zeros(N_Node*3,N_Node*3);
fori=1:
2
n1=Element(1,i);
forj=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
子程序二:
整体坐标刚度矩阵
functionK=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^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,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^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=T'*(KE+KG)*T;
end
子程序三:
局部坐标刚度矩阵
functionK=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,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^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);
fori=1:
N_Force
ForceMatrix(Force(i,1)*3-2:
Force(i,1)*3,1)=Force(i,2:
4)';%把节点荷载向量读入一个矩阵中,形成列向量=总的自由度个数
end
del=[];
fori=1:
N_Boundary;
ifBoundary(i,2)==0;
m=3*Boundary(i,1)-2;
del=[del,[m]];%受约束节点位移为0时所对应的指标数值1
end
ifBoundary(i,3)==0;
m=3*Boundary(i,1)-1;
del=[del,[m]];%受约束节点位移为0时所对应的指标数值2
end
ifBoundary(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;
forn=0:
Allstep-1
n=n+1
K_Global=zeros(N_Node*3,N_Node*3);%总刚矩阵
fori=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);
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);%形成总刚K0
end%整体刚度矩阵
Delta_Force=ForceMatrix/Allstep;%初始荷载