有限单元法作业非线性分析 程序Word格式文档下载.docx
《有限单元法作业非线性分析 程序Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《有限单元法作业非线性分析 程序Word格式文档下载.docx（47页珍藏版）》请在冰豆网上搜索。
（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.00142MATLAB求得的极限承载力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
从表中可以初步得出：
在一定随着矢跨比的增加，拱仍然呈现跳跃失稳的形式，拱结构的极限承载能力有大幅度的提高；
在最低处的承载力呈现出反向，相当于有一个拉力在阻止拱结构发生跳跃失稳，矢跨比越大，拱越不容易发生跳跃失稳。
当拱的矢跨比超过一定范围后，拱将发生复杂的不同于跳跃失稳的失稳形式。
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
从图中可以看出：
矢跨比在一定范围内，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.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
误差=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.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对应的值。
这可能与单元划分的个数，求解精度有关，但在工程允许的范围内，还是可以接受的。
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_Force=ForceMatrix/Allstep;
%初始荷载
升级会员 