有限单元法作业非线性分析+程序Word文件下载.docx
《有限单元法作业非线性分析+程序Word文件下载.docx》由会员分享,可在线阅读,更多相关《有限单元法作业非线性分析+程序Word文件下载.docx(55页珍藏版)》请在冰点文库上搜索。
完全相同,说明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拱单元信息
约束信息(0为约束,1为放松)荷载信息(FX,FY,M)
节点信息
2.2运用ANSYS和MATLAB进行求解对比(两端铰接)
ANSYS中模型:
图2.3ANSYS模型
图2.4MATLAB和ANSYS变形图
图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)
最低点
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
从表中可以初步得出:
在一定随着矢跨比的增加,拱仍然呈现跳跃失稳的形式,拱结构的极限承载能力有大幅度的提高;
在最低处的承载力呈现出反向,相当于有一个拉力在阻止拱结构发生跳跃失稳,矢跨比越大,拱越不容易发生跳跃失稳。
当拱的矢跨比超过一定X围后,拱将发生复杂的不同于跳跃失稳的失稳形式。
2)MATLAB与ANSYS计算结果对比
图2.8ANSYS和MATLAB对比荷载位移曲线
表2各矢跨比下拱结构的极限荷载对比
F(N)MAT
F(N)ANA
误差(%)
3042.53
3.00142
0.04
2.45
7624.91
3.96303
-0.02
1.75
14974.3
5.3157
0.00
1.61
39695.7
9.59955
0.24
-1.23
从图中可以看出:
矢跨比在一定X围内,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。
实际误差相差很小,在工程允许的X围内是可以接受的。
2.5荷载位置对拱荷载位移曲线的影响
图2.9ANSYS模型简图
1)MATLAB中计算荷载位置对拱荷载位移曲线的影响
图2.10各加载点荷载位移曲线
表3改变加载点拱结构的极限荷载
加载点
26
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.11ANSYS和MATLAB对比荷载位移曲线
表4各加载点拱结构的极限荷载
3569.73
3.24865
0.01
-1.87
4728.71
2.91862
-1.34
14324.8
1.29764
-0.05
-1.17
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.12ANSYS和MATLAB固接和铰接的荷载位移曲线
拱的边界条件对其的失稳形式有很大影响。
两端固接拱的稳定性明显优于两端铰接拱,承载能力也大幅度提高。
固接拱不会发生跳跃失稳的形式,刚度在初始阶段会减小,待到达一定程度后刚度又会增加。
而两端铰接拱会发生跳跃失稳的形式。
2.7参数m对拱失稳形式的影响
文献中给出:
m是一个由拱截面在竖向平面内的回转半径r和拱的初始矢高h无确定的无量纲参数。
当m>
=1时,扁拱不会出现跳跃屈曲,当0<
m<
1时,扁拱有可能发生跳跃屈曲,而影响扁拱是否发生跳跃屈曲的主要因素是m值和荷载参数。
2.13不同m值加载点的荷载位移曲线
2.14不同m值变形后拱曲线
从MATLAB的计算结果中可以验证:
不同的m系数对应拱不同的失稳形式。
=1时,扁拱不会出现跳跃屈曲,当0<
1时,扁拱有可能发生跳跃屈曲。
但拱的最终变形图非常接近,只是此时拱的失稳形式变了。
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
F(N)ABA
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对应的值。
这可能与单元划分的个数,求解精度有关,但在工程允许的X围内,还是可以接受的。
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:
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
子程序二:
整体坐标刚度矩阵
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;
子程序三:
局部坐标刚度矩阵
functionK=beam2d_stiffness520(E,A,I,L,cs,Ele_F1);
K=(KE+KG);
End
主程序一:
Newton-Raphson法
clc
clear
Node=xlsread('
Data.xls'
'
Node'
);
Element=xlsread('
Element'
Boundary=xlsread('
Boundary'
Section=xlsread('
Section'
Force=xlsread('
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);
N_Force
ForceMatrix(Force(i,1)*3-2:
Force(i,1)*3,1)=Force(i,2:
4)'
;
%把节点荷载向量读入一个矩阵中,形成列向量=总的自由度个数
del=[];
N_Boundary;
ifBoundary(i,2)==0;
m=3*Boundary(i,1)-2;
del=[del,[m]];
%受约束节点位移为0时所对应的指标数值1
ifBoundary(i,3)==0;
m=3*Boundary(i,1)-1;
%受约束节点位移为0时所对应的指标数值2
ifBoundary(i,4)==0;
m=3*Boundary(i,1);
%受约束节点位移为0时所对应的指标数值3
%求出约束节点的标号,便于刚度、荷载矩阵归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_For